Testing f(R) Gravity Using Gravitational-Wave Signals from Binary Mergers

1. Introduction

Nowadays, General Relativity (GR) remains one of the foundational frameworks of our description of nature. The theory has withstood numerous experimental and observational tests, from Solar-System measurements to the recent detections of gravitational waves (GW) by ground-based interferometers. Nonetheless, theoretical motivations arising from attempts to construct a complete quantum theory of gravity, together with cosmological observations and the problem of dark energy, continue to motivate tests of GR in regimes where its predictions may be affected by new physics [1,2]. In particular, the strong-field, highly dynamical regime of compact-binary coalescences constitutes a unique laboratory for probing the limits of GR and constraining alternative gravitational theories [3].

The advent of GW astronomy [4,5], particularly in the multimessenger context [6,7], has opened a novel observational window onto these extreme environments [8]. The signals emitted during the inspiral, merger, and ringdown of binary black holes and neutron stars encode detailed information about the underlying theory of gravity. Potential deviations from GR would appear as modifications of the waveform phasing, amplitude, polarization content, or propagation properties. Consequently, GW observations furnish a direct and powerful probe to test GR and to constrain or exclude specific theoretical extensions.

Among the proposed modifications of GR, $f(R)$ theories constitute one of the simplest and most extensively investigated classes [9,10]. In these models the Einstein–Hilbert Lagrangian is replaced by a general function $f(R)$ of the Ricci scalar R, thereby introducing an additional scalar degree of freedom that can produce observable deviations from GR in both cosmological and astrophysical contexts. Although many $f(R)$ models are constructed to reproduce the observed cosmic acceleration without a cosmological constant, they are required to satisfy stringent constraints from Solar-System experiments, binary pulsar timing, and GW observations.

In contrast to $f(R)$ models, other gravity extension theories, such as $f(R,T)$ [12] and $f(R,Q)$ models [13,14], introduce explicit couplings between matter and curvature, leading to much stronger observational constraints. In $f(R,T)$ gravity, the dependence on the trace of the stress-energy tensor generally implies modification of the standard conservation rules of $T_{\mu \nu }$ producing force terms that violate the Weak Equivalence Principle. These effects are tightly constrained by various Solar System tests. Besides, in $f(R,Q)$ models, with $Q=R_{\mu \nu }{T}^{\mu \nu }$, matter-curvature couplings can modify the propagation of GW and generate ghost-like instabilities. The multimessenger observation of GW170817 [8] imposed extremely stringent limits on deviations of the GW speed from the speed of light, excluding large regions of parameter space for these models. Additional bounds arise from binary pulsars, cosmological structure formation, and GW observations from the LIGO-Virgo-KAGRA collaboration [15,16,17,18].

For these reasons, in this work we will focus on $f(R)$ gravity, since it is considerably less constrained by current experimental observations. This scenario still provides a viable framework in which deviations from GR may be experimentally tested with next-generation GW detectors such as Einstein Telescope (ET) [19]. In this work, we investigate the predictions of $f(R)$ gravity for the emission of GW from coalescing compact binaries. We focus on the implications of the modified field equations for the GW signal and its dependence on the functional form of $f(R)$.

We then assess the prospects for testing this class of models with future GW observatories. By forecasting the sensitivity of next-generation interferometers to the deviations predicted by $f(R)$ gravity, we evaluate their ability to constrain model parameters and to discriminate these theories from GR. Our results underscore the role of advanced GW facilities in probing the strong-field regime and in testing GR as the low-energy limit of more general gravitational theories. In this context, we place particular emphasis on the ET sensitivity, comparing the capabilities of its different planned configurations.

The paper is organized as follows. In Section 2, we introduce $f(R)$ theories starting from the modified Einstein–Hilbert action, deriving the corresponding field equations and the massive scalar mode that is absent in GR. In Section 3, we describe the propagation of GW in the $f(R)$ framework, highlighting the differences with respect to the standard GR scenario. In Section 4, we discuss the main effects introduced by $f(R)$ gravity that can be tested experimentally through GW observations. In Section 5, we consider future interferometric GW detectors, focusing mainly on ET and illustrating which detector configurations are best suited to probe the various effects predicted by $f(R)$ models. Finally, in Section 6, we summarize the main results presented in this paper.

2. $f(R)$ Theory Formulation

The motivation underlying the formulation of $f(R)$ theories lies in the fact that the Einstein–Hilbert action may be not unique, but rather the simplest covariant choice for describing gravity. In this framework, one generalizes the Lagrangian from a linear dependence on the Ricci scalar R to an arbitrary function $f(R)$, while still preserving general covariance and locality. Such extensions are well motivated from both theoretical and phenomenological perspectives. Indeed, various types of higher-order curvature corrections naturally appear in effective actions of QG and in semiclassical expansions, suggesting the possibility of terms beyond the linear R contribution. From the phenomenological side, suitable choices of $f(R)$ can account for late-time cosmic acceleration or early-time inflation without introducing additional matter fields. In the weak-curvature limit, consistency with GR is recovered by requiring $f(R)\to R$, ensuring that standard gravitational physics is preserved.

As stated before, the formulation of $f(R)$ extensions of GR is based on a modification of the Einstein–Hilbert action [9,10,11]:

$$S={\displaystyle \frac{1}{16\pi G}}\int \sqrt{|g|}f(R)d^{4}x,$$

(1)

where the Ricci scalar contribution R is replaced by a generic function $f(R)$. Considering the variation of the function $f(R)$,

$$\delta f(R)=f^{\prime }(R)\delta R,$$

(2)

where $f^{\prime }(R)=df(R)/dR$, together with the variation of the Ricci scalar R,

$$\delta R=R_{\mu \nu }\delta g^{\mu \nu }+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\delta g^{\mu \nu },$$

(3)

the equations of motion can be obtained by varying the action with respect to the metric:

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

(4)

Computing the trace of Equation (4) one can obtain the following relation:

$$f^{\prime }(R)R-2f(R)+3\square f^{\prime }(R)=8\pi GT.$$

(5)

which is fundamental because it explicitly reveals the presence of an additional propagating scalar degree of freedom that is not foreseen in the standard formulation of GR.

On a Minkowski background we consider the perturbative expansion of the metric:

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

(6)

where the perturbative condition $|h_{\mu \nu }|\ll 1$ must hold. As in standard GR, the weak-field regime is investigated by introducing the harmonic gauge. Posing:

$${\overline{h}}_{\mu \nu }=h_{\mu \nu }-{\displaystyle \frac{1}{2}}{\eta }_{\mu \nu }h,$$

(7)

the gauge condition is imposed requiring:

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

(8)

Expanding $f(R)$ perturbatively,

$$f(R)=f(0)+f^{\prime }(0)R+{\displaystyle \frac{1}{2}}{f}^{\prime \prime }(0)R^2+\dots ,$$

(9)

and assuming a flat background with $R_0=0$, one must impose $f^{\prime }(R_0)=f^{\prime }(0)=1$ in order to recover standard GR in the linearized limit. Thus one obtains:

$$f(R)\simeq R+{\displaystyle \frac{1}{2}}{f}^{\prime \prime }(R_0)R^2.$$

(10)

The fundamental linearized geometric quantities are obtained as in the context of the linearization of standard GR. The connection (or Christoffel symbol) can be written in the form:

$${\Gamma }_{\mu \nu }^{\alpha }={\displaystyle \frac{1}{2}}{\eta }^{\alpha \beta }\left({\partial }_{\mu }{h}_{\nu \beta }+{\partial }_{\nu }{h}_{\beta \mu }-{\partial }_{\beta }{h}_{\mu \nu }\right),$$

(11)

the linearized Ricci tensor as:

$$R_{\mu \nu }^{(1)}={\displaystyle \frac{1}{2}}\left({\partial }_{\alpha }{\partial }_{\mu }{h}_{\nu }^{\alpha }+{\partial }_{\alpha }{\partial }_{\nu }{h}_{\mu }^{\alpha }-\square h_{\mu \nu }-{\partial }_{\mu }{\partial }_{\nu }h\right),$$

(12)

and the Ricci scalar at the linear perturbative order becomes:

$$R^{(1)}={\partial }_{\mu }{\partial }_{\nu }{h}^{\mu \nu }-\square h,$$

(13)

where $h={\eta }^{\mu \nu }{h}_{\mu \nu }$. Therefore, at first perturbative order, one finds:

$$f^{\prime }(R)={\displaystyle \frac{df(R)}{dR}}=1+f^{\prime \prime }(R_0)R^{(1)}.$$

(14)

3. Gravitational Waves in the $f(R)$ Scenario

To illustrate how Equation (4) introduces an additional scalar degree of freedom, we consider the linearization of the Ricci scalar [20]. The linearized approximation can be obtained posing:

$$\begin{array}{cc}\hfill & R=R_0+\delta R\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & f^{\prime }(R)=f^{\prime }(R_0)+f^{″}(R_0)\delta R\hfill \end{array}$$

(16)

In vacuum with $T=0$, Equation (5) becomes:

$$3f^{″}(R_0)\square \delta R+(f^{″}(R_0)R_0-f^{\prime }(R_0))\delta R=0$$

(17)

since the background contribution vanish because the background satisfies the zeroth-order field equations, yielding an overall zero contribution. The background curvature $R_0$ is fixed by the zeroth-order trace equation, which requires:

$$f^{\prime }(R_0)R_0-2f(R_0)=0.$$

(18)

This condition ensures that the background spacetime is a constant curvature solution of the field equations, allowing a clean separation between background and perturbations and leading to a well-defined Klein–Gordon equation for the scalar degree of freedom [21].

The (17) can be recast in the form of a Klein-Gordon Equation:

$$\square \delta R-m_s^2\delta R=0,$$

(19)

which at first order gives the relation:

$$\square R^{(1)}-m_s^2{R}^{(1)}=0.$$

(20)

$m_s$ represents the effective mass of the propagating massive scalar mode of GW:

$$m_s^2={\displaystyle \frac{f^{\prime }(R_0)-R_0{f}^{″}(R_0)}{3f^{″}(R_0)}}.$$

(21)

In the Minkowski case, where $R_0=0$ and $f^{\prime }(R_0)=1$ the effective mass becomes:

$$m_s^2={\displaystyle \frac{1}{3f^{″}(R_0)}}.$$

(22)

The different GW emission modes can be decomposed separating the tensorial perturbation from the scalar one:

$$h_{\mu \nu }=h_{\mu \nu }^{TT}+h_{\mu \nu }^{(S)},$$

(23)

introducing the tensorial traceless perturbation of the metric $h_{\mu \nu }^{TT}$ and the scalar perturbation $h_{\mu \nu }^{(S)}$. The tensorial mode satisfies the relation:

$$\square h_{\mu \nu }^{TT}=0.$$

(24)

This relation admits the usual solution that depends on the polarization ${ϵ}_{\mu \nu }$:

$$h_{\mu \nu }^{TT}={ϵ}_{\mu \nu }{e}^{i{k}_{\alpha }{x}^{\alpha }}.$$

(25)

The resulting tensor mode propagates exactly as in GR, with dispersion relation given by:

$$k^{\alpha }{k}_{\alpha }=0.$$

(26)

In the harmonic traceless gauge the tensorial mode propagates with two different degrees of freedom, related to the usual GW polarizations $h_{+}$ and $h_{\times }$. The propagation group velocity of the tensorial perturbation is the light speed $c=1$, as expected for a non-massive propagating solution.

On the other hand, the scalar mode behaves as a massive mode and the solution of Klein-Gordon Equation (20) gives the relation:

$$R^{(1)}=A{e}^{i{k}_{\alpha }{x}^{\alpha }},$$

(27)

and the dispersion relation:

$$k^{\alpha }{k}_{\alpha }-m_s^2=0.$$

(28)

As a result, the group velocity associated with the massive scalar mode results:

$$v_g^{(S)}={\displaystyle \frac{\partial \omega }{\partial k}}={\displaystyle \frac{k}{\sqrt{k^2+m_s^2}}}.$$

(29)

The propagating scalar mode has only one degree of freedom, so it can be written as

$$h_{ij}^{(S)}={\delta }_{ij}\Phi ,$$

(30)

with $\Phi \propto R^{(1)}$, after introducing the polarization structure in the transverse plane, which takes the form of the Euclidean metric ${\delta }_{ij}$ in the appropriate gauge. In the massive case, the scalar mode can induce a longitudinal component due to the non-null dispersion relation.

4. Physical Effects Caused by the Massive Scalar Mode Emission

The difference caused by this dispersive propagation can be potentially observed in the GW emission of coalescing objects. In the following we list the different effects that can be potentially detected in the GW sector.

4.1. Different Propagation Velocity

The first effect concerns the time delay accumulated by the scalar mode during propagation [20].

The emission of the scalar mode from a system of coalescing astrophysical objects can be obtained by introducing a source term in Equation (20), represented by the trace of the stress-energy tensor:

$$(\square -m_s^2)R^{(1)}={\displaystyle \frac{8\pi G}{3}}T.$$

(31)

The solution of this equation can be formally expressed through the retarded Green function of the Klein–Gordon operator:

$$R^{(1)}(x)={\displaystyle \frac{8\pi G}{3}}\int d^4{x}^{\prime }{G}_{\mathrm{ret}}^{(m_s)}(x-x^{\prime })T(x^{\prime }),$$

(32)

where $G_{\mathrm{ret}}^{(m_s)}$ satisfies:

$$(\square -m_s^2)G_{\mathrm{ret}}^{(m_s)}(x-x^{\prime })={\delta }_{(4)}(x-x^{\prime }).$$

(33)

In four-dimensional spacetime, the retarded Green function takes the explicit form

$$G_{\mathrm{ret}}^{(m_s)}(t,\mathbf{r})=\Theta (t)\left[{\displaystyle \frac{\delta (t-r)}{4\pi r}}-\Theta (t-r){\displaystyle \frac{m_s}{4\pi }}{\displaystyle \frac{J_1\left(m_s\sqrt{t^2-r^2}\right)}{\sqrt{t^2-r^2}}}\right],$$

(34)

where $r=|\mathbf{r}|$, $\Theta$ denotes the Heaviside step function, and $J_1$ is the Bessel function of the first kind.

The structure of the Green function highlights the fundamental difference between the massive and massless cases. For a massless field ($m_s=0$), one recovers:

$$G_{\mathrm{ret}}^{(0)}={\displaystyle \frac{\delta (t-r)}{4\pi r}},$$

(35)

consequently, the perturbation propagates exactly at the speed of light and all the information reaches the observer at the retarded time $t=r$. On the other hand, for a massive scalar field ($m_s\ne 0$), an additional contribution appears:

$$-\Theta (t-r){\displaystyle \frac{m_s}{4\pi }}{\displaystyle \frac{J_1\left(m_s\sqrt{t^2-r^2}\right)}{\sqrt{t^2-r^2}}}.$$

(36)

This term is non-vanishing for $t>r$. As a result, the signal is no longer confined to the light cone, and part of the perturbation is delayed. This behaviour is the mathematical origin of the dispersive propagation of the scalar mode and ultimately leads to the modified dispersion relation, which is responsible for the time delays and phase shifts.

Considering the limit of negligible mass of the scalar mode compared to the propagation energy $m_s\ll k_0$, the group velocity (Equation (29)) of the scalar mode can be written as:

$$v_g^{(S)}\simeq 1-{\displaystyle \frac{m_s^2}{2k_0^2}}.$$

(37)

The propagation time of the tensorial non-massive mode:

$$t_T=D$$

(38)

can be compared with the time related to the scalar mode:

$$t_S={\displaystyle \frac{D}{v_g^{(S)}}}\simeq D\left(1+{\displaystyle \frac{m_s^2}{2k_0^2}}\right),$$

(39)

obtaining in the ultrarelativistic approximation ($k_0\gg m_s$) the time delay in the form:

$$\Delta t=D{\displaystyle \frac{m_s^2}{2k_0^2}}.$$

(40)

The observation of a retarded signal can be therefore associated with the emission of the massive scalar mode predicted in the context of $f(R)$ theories.

The modified propagation velocity of the massive scalar mode can lead to several observable effects in GW signals.

4.2. Difference in the Observed GW Phase

The difference between the group velocities of the tensorial non-massive mode and the scalar massive mode can induce a difference in the observed GW phase [22].

The phase of the detected GW is:

$$\varphi =k_{0}t.$$

(41)

Thus, introducing the propagation time associated with the tensorial and the scalar modes, one can obtain the GW phases related to the different propagating modes:

$$\begin{array}{cc}\hfill & {\varphi }_T=k_{0}D,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & {\varphi }_S=k_{0}D+{\displaystyle \frac{m_s^{2}D}{2k_0}}.\hfill \end{array}$$

(43)

It therefore appears evident that the difference in propagation speeds introduces a dispersive phase difference given by:

$$\Delta \varphi ={\displaystyle \frac{m_s^{2}D}{2k_0}}.$$

(44)

The presence of a dispersive difference in the GW phase can manifest itself modifying the waveform of the propagating perturbation:

$$\tilde{h}(f)=A(f)e^{i\varphi (f)}=A(f)exp\left[i\left({\varphi }_T+{\displaystyle \frac{m_s^{2}D}{4\pi f}}\right)\right],$$

(45)

where the phase is expressed as a function of the frequency f, that is proportional to the energy $k_0$. This result implies that different frequencies arrive with different time delays, introducing a dephasing effect. Moreover, the two different propagating modes can interfere:

$$h(f)=h_T(f)+h_S(f)=A_T(f)e^{i{\varphi }_T(f)}+A_S(f)e^{i{\varphi }_S(f)},$$

(46)

producing a small modulation of the observed waveform, phase drift and slow modulation.

The observational limits that place stringent constraints on the mass of the scalar mode $m_s$ strongly restrict the possibility of observing echoes in the signal due to the time delay. Indeed, the accumulated delay is severely limited and the scalar mode is not sufficiently intense. What remains observable is a dephasing effect that can impact the waveform.

4.3. GW Emission

The emission of GWs in the context of (f(R)) theories is modified. The introduction of the scalar mode can indeed lead to dipole emission, which is not predicted in standard GR [22].

It is useful the introduction of the expansion of the source term in multipoles to deal with the emission solution. In this way one obtains monopole, dipole and quadrupole contributions. The monopole term does not radiate for isolated systems because of energy conservation, while the dipole term can become non-vanishing, differently from what predicted in the context of GR. Introducing the useful definition of scalar charges $q_A={\alpha }_A{m}_A$, the scalar dipole moment of a binary system becomes:

$$\overrightarrow{D}(t)=\sum _A{q}_A{\overrightarrow{x}}_A(t),.$$

(47)

In the case of a binary system made of two compact objects possessing different scalar charges the emission from dipole is allowed.

For a binary system of coalescing objects in quasi-circular orbit, using the center-of-mass condition, the dipole moment becomes:

$$\overrightarrow{D}(t)=\mu ({\alpha }_1-{\alpha }_2)\overrightarrow{r}(t),$$

(48)

where $M=m_1+m_2$ is the total mass and $\mu =m_1{m}_2/M$ is the reduced mass.

The dipolar emission power radiated is proportional to the square of the second time derivative of the dipole moment:

$$P_{\mathrm{dip}}\propto 〈{\ddot{\overrightarrow{D}}}^2〉.$$

(49)

For circular orbits with orbital frequency $\Omega$, one obtains:

$$P_{\mathrm{dip}}\sim G{\mu }^2{({\alpha }_1-{\alpha }_2)}^2{r}^2{\Omega }^4.$$

(50)

This contribution must be added to the standard quadrupole emission predicted by GR:

$$P_{\mathrm{GR}}={\displaystyle \frac{32}{5}}G{\mu }^2{r}^4{\Omega }^6.$$

(51)

The dipolar emission represents an additional energy dissipation mechanism, that can influence the inspiral phase of the coalescence. Indeed, from the Kepler law one can obtain that the dipole contribution scales as $v^8$, whereas the quadrupole term scales as $v^{10}$. Therefore:

$${\displaystyle \frac{P_{\mathrm{dip}}}{P_{\mathrm{GR}}}}\sim v^{-2},$$

(52)

showing that dipole radiation becomes particularly important during the early inspiral phase, where the orbital velocity is still relatively small.

The additional energy loss accelerates the inspiral evolution and modifies the GW phase. In the frequency domain, the phase receives an additional contribution due to scalar dipole emission:

$$\Psi (f)={\Psi }_{\mathrm{GR}}(f)-\beta {(\pi \mathcal{M}f)}^{-7/3}.$$

(53)

This characteristic negative post-Newtonian correction represents one of the main observational signatures of $f(R)$ gravity in GW interferometers, particularly those of the future generation such as ET.

4.4. GW Polarizations

In $f(R)$ gravity, the appearance of additional GW polarizations follows directly from the linearization of the modified field equations [23,24]. Starting from the trace equation in vacuum (Equation (17)):

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

(54)

and expanding around a constant-curvature background, one obtains at first order a dynamical Klein-Gordon Equation (20) for the Ricci scalar perturbation. Thus, $R^{(1)}$ behaves as a propagating scalar field. In parallel, the metric perturbation is decomposed as Equation (23):

$$h_{\mu \nu }=h_{\mu \nu }^{TT}+h_{\mu \nu }^{(S)},$$

(55)

where $h_{\mu \nu }^{TT}$ satisfies the standard GW propagation Equation (24), while the scalar sector is sourced by $R^{(1)}$.

At linear order, the scalar contribution to the spatial metric perturbation takes the form of Equation (30). Substituting $h_{ij}^{(S)}=\Phi {\delta }_{ij}$ with $\Phi \propto R^{(1)}$ into the geodesic deviation equation:

$${\ddot{\xi }}^i=-{\displaystyle \frac{1}{2}}{\ddot{h}}_{ij}{\xi }^j,$$

(56)

yields

$${\ddot{\xi }}^i=-{\displaystyle \frac{1}{2}}\ddot{\Phi }{\delta }_{ij}{\xi }^j,$$

(57)

which produces an isotropic expansion and contraction of a ring of test particles in the transverse plane, which is named breathing polarization mode.

In addition, when $m_s\ne 0$, the dispersion relation

$${\omega }^2=k^2+m_s^2$$

(58)

implies a non-null longitudinal response, leading to a small but non-vanishing deformation along the propagation direction. Therefore, $f(R)$ gravity predicts, in addition to the standard tensor polarizations $h_{+}$ and $h_{\times }$, a scalar breathing mode (and a possible longitudinal component in the massive case), which directly originates from the dynamical nature of the Ricci scalar.

5. Experimental Perspectives with Next-Generation Interferometers

In the framework of $f(R)$ gravity, compact-binary GW signals provide a particularly promising probe for next-generation ground-based interferometers, in particular ET. The additional massive scalar degree of freedom modifies propagation via the dispersion relation (Equation (58)). This yields a frequency-dependent group velocity (Equations (29) and (37)) and therefore a relative propagation delay between tensor and scalar components. Current empirical bounds constrain the scalar mass to very small values [25]. For concreteness, effects in ground-based detectors are relevant for scalar masses approximately $m_s\ll {10}^{-13}-{10}^{-15}\mathrm{eV}/c^2$, which correspond to wavelengths and frequency scales in the few-hertz regime. From a detection point of view, this does not generally produce clearly separated echo-like signals (which are strongly suppressed by current bounds on $m_s$), but instead manifests as a coherent deformation of the waveform, encoded as a frequency-dependent phase shift and a mild reshaping of the observed spectrum [26]. Propagation effects accumulate over astrophysical distances and can be parametrized to give simple order-of-magnitude estimates. The modifications introduced in the group velocity Equation (37) lead to the propagation delay Equation (40) that depends on the propagation length. Thus the induced phase shift grows with the propagation path and is larger at low frequencies, which makes low-frequency sensitivity crucial, because group-velocity differences scale inversely with the square of the energy, that is the inverse of the square of the frequency $f^{-2}$ (Equation (37)). Moreover, the detectability of such dephasing benefits from the large number of inspiral cycles in band: for fixed signal-to-noise ratio (SNR), the phase precision improves with the effective number of cycles, so early inspiral–dominated signals give the strongest constraints on dispersive propagation and dipole-like emission.

Propagation effects such as time delays (Equations (39) and (40)), dispersion-induced dephasing, and possible spectral distortions are most effectively constrained when the detector configuration maximizes sensitivity to phase coherence across spatially separated baselines. In this respect, a geographically separated double-L configuration (two geographically separated L-shaped interferometers [27]) may offer advantages for measurements relying on timing triangulation and inter-site consistency checks, potentially helping to constrain propagation-induced effects. However, the overall sensitivity to dispersive propagation is determined primarily by the low-frequency performance and signal-to-noise ratio of the detector network. Since propagation effects accumulate over cosmological distances, they manifest themselves as frequency-dependent phase shifts and timing offsets in the observed waveform. Networks with multiple widely separated detectors may help reduce parameter degeneracies and improve consistency tests of such effects. As a consequence, a double-L network can be a useful observational handle on dispersive propagation effects and residual time-delay signatures, particularly through improved timing triangulation and consistency checks between sites.

By contrast, detection of GW polarizations is optimally achieved with a triangular configuration. In $f(R)$ gravity, the scalar mode generates a breathing polarization (Equation (30)) that cannot be uniquely disentangled by a single L-shaped interferometer because of degeneracies in the detector antenna response. The triangular geometry, consisting of three co-located interferometers rotated by ${60}^{\circ }$ [27], provides independent projections of the same signal at a single site. This enables an almost complete reconstruction of the polarization content, greatly improving the ability to disentangle tensorial + and × modes and the scalar breathing (and, in the massive case, possible longitudinal) components. Therefore, while the double-L ET configuration can offer advantages for propagation-related observables such as time delay (Equations (39) and (40)), dispersion (Equations (29) and (37)), and spectral deformation, the triangular ET design is fundamentally superior for polarization reconstruction and model-independent tests of additional gravitational degrees of freedom. Practical searches must also account for potential degeneracies and analysis choices. Effects that can partially mimic dispersive or scalar-induced waveform deformations include orbital eccentricity, spin precession, higher-mode content, environmental interactions, and instrumental calibration errors.

6. Conclusions

Since $f(R)$ theories still exhibit a region of parameter space not yet ruled out by current experiments, GW detected with the precision expected from future detectors such as ET can provide an ideal framework for placing stringent constraints on this class of theories. We have shown that compact-binary GW signals offer a sensitive probe of additional scalar degrees of freedom through both propagation and emission effects. In particular, we have considered the signatures associated with the scalar degree of freedom predicted by $f(R)$ models, namely propagation effects such as a frequency-dependent group velocity, relative time delays between tensor and scalar components, and dispersion-induced phase shifts, as well as scalar emission effects such as dipole radiation. These effects accumulate over astrophysical distances and become more pronounced at low frequencies. Moreover, scalar emission is also more significant during the early inspiral phase.

The detectability of $f(R)$-induced modifications depends jointly on detector bandwidth, network geometry, and overall sensitivity, as well as on the ability to control astrophysical degeneracies. Consequently, the optimal experimental strategy combines complementary design and analysis features, such as high sensitivity at low frequencies (in the few-hertz regime) to maximize phase accumulation during the inspiral. Long-baseline networks, such as a double-L configuration, can be useful to probe propagation delays and dispersive signatures, offering the opportunity to strongly constrain this class of theories. On the other hand, co-located triangular configurations can enable a full reconstruction of the polarization content. In both cases we have illustrated the importance of having very high sensitivity at a few Hz, since this is essential for probing the region of parameter space most sensitive to the modifications induced by $f(R)$ models.

Implementing this strategy requires extended low-frequency performance in future detectors, coordinated multi-site observations, and, where possible, multi-band observations. In this context, the construction of next-generation GW detectors such as ET will be fundamental for testing this class of GR-modifying theories.