Scalar modes in extended hybrid metric-Palatini gravity: weak field phenomenology

TL;DR

This paper explores the weak field phenomenology of extended hybrid metric-Palatini gravity, revealing scalar modes that influence gravitational potential corrections and wave polarizations, with implications for dark matter and gravitational wave observations.

Contribution

It introduces and analyzes two coupled scalar modes in extended hybrid metric-Palatini gravity, detailing their effects on Newtonian potential and gravitational wave polarizations.

Findings

Yukawa corrections to Newtonian potential from scalar fields

Scalar modes can mimic dark matter effects

Scalar gravitational waves exhibit breathing and longitudinal polarizations

Abstract

We investigate the nature of additional scalar degrees of freedom contained in extended hybrid metric-Palatini gravity, outlining the emergence of two coupled dynamical scalar modes. In particular, we discuss the weak field limit of the theory, both in the static case and from a gravitational waves perspective. In the first case, performing an analysis at the lowest order of the post parameterized Newtonian (PPN) structure of the model, we stress the settling of Yukawa corrections to the Newtonian potential. In this respect, we show that one scalar field can have long range interactions and used in the principle for mimicking dark matter effects. Concerning the gravitational waves propagation, instead, we demonstrate that is possible to have well-defined physical degrees of freedom, provided by suitable constraints on model parameters. Moreover, the study of the geodesic deviation points out the presence of breathing and longitudinal polarizations due to these novel scalar waves, which on peculiar assumptions can give rise to beating phenomena during their propagation.

Equations

$$S=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}f(R,\mathcal{R})+S_M$$

$$R=g^{\mu \nu }{R}_{\mu \nu }(\Gamma (g)),$$

$$\Gamma \text{indices}{}_{\mu \nu }^{\rho }=\frac{g^{\rho \tau }}{2}\left({\partial }_{\mu }{g}_{\nu \tau }+{\partial }_{\nu }{g}_{\mu \tau }-{\partial }_{\tau }{g}_{\mu \nu }\right),$$

$$\mathcal{R}=g^{\mu \nu }{\mathcal{R}}_{\mu \nu }(\tilde{\Gamma }).$$

$${\tilde{\nabla }}_{\rho }\left(\sqrt{-g}{f}_{\mathcal{R}}{g}^{\mu \nu }\right)=0,$$

$$\tilde{\Gamma }\text{indices}{}_{\mu \nu }^{\rho }=\frac{{\tilde{g}}^{\rho \tau }}{2}\left({\partial }_{\mu }{\tilde{g}}_{\nu \tau }+{\partial }_{\nu }{\tilde{g}}_{\mu \tau }-{\partial }_{\tau }{\tilde{g}}_{\mu \nu }\right),$$

$$\mathcal{R}=R+\frac{3}{2f_{\mathcal{R}}^2}{\partial }_{\mu }{f}_{\mathcal{R}}{\partial }^{\mu }{f}_{\mathcal{R}}-\frac{3\square f_{\mathcal{R}}}{f_{\mathcal{R}}}.$$

$$f_R{R}_{\mu \nu }+f_{\mathcal{R}}{\mathcal{R}}_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }f-({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square )f_R=\kappa T_{\mu \nu },$$

$$f_{R}R+f_{\mathcal{R}}\mathcal{R}-2f+3\square f_R=\kappa T.$$

$$\{\begin{array}{ll}& 3\square f_{\mathcal{R}}-f_{\mathcal{R}}R+f_{\mathcal{R}}\mathcal{R}-\frac{3}{2f_{\mathcal{R}}}{\partial }_{\mu }{f}_{\mathcal{R}}{\partial }^{\mu }{f}_{\mathcal{R}}=0\\ & 3\square f_R+f_{R}R+f_{\mathcal{R}}\mathcal{R}-2f=\kappa T.\end{array}$$

$$f_{\mathcal{R}}\mathcal{R}-2f=\kappa T+R,$$

$$f_{2\mathcal{R}}\mathcal{R}-2f_2=\kappa T+2f_1-3\square f_{1R}-f_{1R}R,$$

$$f_{RR}{f}_{\mathcal{RR}}-f_{R\mathcal{R}}^2\ne 0,$$

$$S=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}\left(\psi R+\xi \mathcal{R}-V(\psi ,\xi )\right)+S_M,$$

$$V(\psi ,\xi )\equiv \psi R(\psi ,\xi )+\xi \mathcal{R}(\psi ,\xi )-f(R(\psi ,\xi ),\mathcal{R}(\psi ,\xi )).$$

$$S_g=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}\left((\psi +\xi )R+\frac{3}{2\xi }{\partial }_{\mu }\xi {\partial }^{\mu }\xi -V(\psi ,\xi )\right),$$

$$S_g=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}\left(\varphi R+\frac{3}{2\xi }{\partial }_{\mu }\xi {\partial }^{\mu }\xi -W(\varphi ,\xi )\right),$$

$$\begin{array}{rl}{R}_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R& -\frac{1}{\varphi }({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square )\varphi +\frac{3}{2\xi \varphi }{\partial }_{\mu }\xi {\partial }_{\nu }\xi \\ & -\frac{1}{2\varphi }{g}_{\mu \nu }\left(\frac{3}{2\xi }{\partial }_{\rho }\xi {\partial }^{\rho }\xi -W(\varphi ,\xi )\right)=\frac{\kappa }{\varphi }{T}_{\mu \nu },\end{array}$$

$$R=\frac{\partial W(\varphi ,\xi )}{\partial \varphi }$$

$$\frac{3}{2{\xi }^2}{\partial }_{\mu }\xi {\partial }^{\mu }\xi -\frac{3\square \xi }{\xi }-\frac{\partial W(\varphi ,\xi )}{\partial \xi }=0.$$

$$R=\frac{3\square \varphi }{\varphi }-\frac{3{\partial }_{\mu }\xi {\partial }^{\mu }\xi }{2\xi \varphi }+\frac{2W}{\varphi }-\frac{\kappa T}{\varphi },$$

$$\square \varphi -\frac{1}{2\xi }{\partial }_{\mu }\xi {\partial }^{\mu }\xi +\frac{2W-\varphi W_{\varphi }}{3}=\frac{\kappa }{3}T$$

$$\square \xi -\frac{1}{2\xi }{\partial }_{\mu }\xi {\partial }^{\mu }\xi +\frac{\xi W_{\xi }}{3}=0,$$

$$\square \psi +\frac{2V-\psi V_{\psi }-\xi V_{\xi }}{3}=\frac{\kappa }{3}T$$

$$\square \xi -\frac{1}{2\xi }{\partial }_{\mu }\xi {\partial }^{\mu }\xi +\frac{\xi (V_{\xi }-V_{\psi })}{3}=0.$$

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

$$\varphi ={\varphi }_0+\delta \varphi \xi ={\xi }_0+\delta \xi .$$

$$\begin{array}{rl}W(\varphi ,\xi )\simeq & W_0+W_{0,\varphi }\delta \varphi +W_{0,\xi }\delta \xi +\\ & +\frac{1}{2}\left(W_{0,\varphi \varphi }\delta {\varphi }^2+W_{0,\xi \xi }\delta {\xi }^2+2W_{0,\varphi \xi }\delta \varphi \delta \xi \right),\end{array}$$

$$\begin{array}{rl}& \det \text{H}(W_0)\equiv W_{0,\varphi \varphi }{W}_{0,\xi \xi }-W_{0,\varphi \xi }^2>0,\\ & W_{0,\varphi \varphi },W_{0,\xi \xi }>0.\end{array}$$

$$R^{(0)}=\frac{2W_0}{{\varphi }_0},$$

$$\begin{array}{rl}{R}_{\mu \nu }=& \frac{\kappa }{\varphi }\left(T_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }T\right)+\frac{1}{2\varphi }{g}_{\mu \nu }\left(\square \varphi +W\right)\\ & +\frac{1}{\varphi }{\nabla }_{\mu }{\nabla }_{\nu }\varphi -\frac{3}{2\varphi \xi }{\partial }_{\mu }\xi {\partial }_{\nu }\xi ,\end{array}$$

$$h\text{indices}{}_{\nu ,\mu }^{\mu }-\frac{1}{2}h\text{indices}{}_{\mu ,\nu }^{\mu }=\frac{{\partial }_{\nu }\delta \varphi }{{\varphi }_0},$$

Full text

I Introduction

During last years, modified theories of gravity have been intensively studied in order to address problems of modern cosmology. Indeed, current evidences of a phase of accelerated expansion for the Universe Riess:1998cb ; Perlmutter:1998np ; Knop:2003iy ; Amanullah:2010vv ; Weinberg:2012es , along with structure dynamics in astrophysical scenarios, e.g. galaxy rotation curves or clusters properties Persic:1995ru ; Wu:1998ju ; Firmani:2000ce , represent inescapable issues for any reliable attempt of providing a unitary theoretical picture of gravitational interaction on different scales. In fact, our present description of Universe evolution, based on the so-called $\Lambda$CDM model, requires the uncomfortable introduction of two unspecified dark components into the matter-energy budget of the Universe. Dark energy, responsible for a de Sitter phase of accelerated expansion and comparable with a cosmological constant term in Einstein equations ($\Lambda$), and cold dark matter (CDM), thought as non relativistic particles interacting with ordinary matter mostly gravitationally Overduin:2004sz ; Baer:2014eja ; Bernal:2017kxu . This model however, even though phenomenologically well-grounded, is not capable of a satisfactory theoretical justification for its additional dark elements. Especially, it is still object of debate the process originating the effective cosmological constant Martel:1997vi ; Carroll:2000fy ; Peebles:2002gy ; Padmanabhan:2002ji , whose observed value is in contrast with predictions of quantum field theory, or the proper nature of dark matter particles Navarro:1995iw ; Jungman:1995df ; Bertone:2004pz ; ArkaniHamed:2008qn . With this regard, a different perspective is therefore offered by the possibility of modifying the nature of the gravitational interaction as predicted by General Relativity, with the aim of accounting for these exotic phenomena as purely dynamical effects, e.g. introducing additional degrees of freedom Nojiri:2005jg ; DeLaurentis:2015fea ; Joyce:2016vqv or modified stress-energy couplings to geometry Nicolis:2008in ; Deffayet:2009wt ; Harko:2011kv ; Odintsov:2013iba ; Wu:2018idg ; Barrientos:2018cnx . Of course, a large number of choices for an extended Einstein-Hilbert action is actually feasible, involving different contributions in metric derivatives Bergmann:1968ve ; Lovelock:1971yv ; Horndeski:1974wa ; Whitt:1984pd ; Schmidt:2006jt ; Bahamonde:2015zma , as well as gauge theory approaches for the gravitational field Hehl:1994ue ; Jackiw:2003pm ; Blagojevic:2012bc ; Hehl:2013qga ; Ashtekar:1986yd ; Immirzi:1996di . Among these available models, $f(R)$ theories stand for their relevance and simplicity Sotiriou:2008rp , where a new degree of freedom is introduced by replacing the Ricci scalar $R$ of the standard General Relativity action with a generic function of it, leading to fourth order equations of motion for the metric field. Cosmological scenarios stemming from such revisited theoretical framework have been deeply investigated, with dark energy-like solutions widely discussed Capozziello:2005ku ; Cognola:2007zu ; Nojiri:2010wj ; Nojiri:2017ncd , and dark matter issue addressed by means of the additional scalar mode featuring this reformulation Boehmer:2007kx ; Stabile:2013jon , made manifest in its scalar tensor restatement Flanagan:2003iw ; Olmo:2005hc ; Capone:2009xk ; ST . In this respect, however, the requirement of preserving solar system local dynamics Will:2005va ; Zakharov:2006uq ; Chiba:2006jp ; Schmidt:2008qi ; Berry:2011pb , consisting in very short range scalar interaction, turned out to be inconsistent with demands of late time expansion, involving instead astrophysical range deformations of gravitational force, and led to introduce peculiar screening mechanisms Khoury:2003rn ; Brax:2008hh ; Capozziello:2007eu .

Another source of ambiguity for the gravitational field dynamics is offered by the nature of the metric field and the affine connection, which could be considered in principle as independent variables. Such an approach, corresponding to the so-called Palatini (or first order) formulation, appears very promising especially for its implication in the quantization of gravity as a gauge field theory Ashtekar:1986yd ; Immirzi:1996di . However, even if the Hilbert-Palatini action is in vacuum at all equivalent to the metric Einsten-Hilbert analogous Holst:1995pc , it outlines significant differences for instance when fermions are included in the dynamics Mercuri:2007ki . In fact, spinor fields couple to connection and induce a non-vanishing torsion in spacetime structure, so that the equivalence with the second order approach is intrinsically lost, forcing us to deal instead with Einstein-Cartan geometry Hehl:1976kj ; Shapiro:2001rz . Similar issues arise when Palatini scheme is applied to $f(\mathcal{R})$ models Olmo:2011uz , and several discrepancies emerge with respect to the corresponding metric (or second order) analysis. Especially, the connection turns out to be an auxiliary field devoid of proper dynamics, whose expression depends on the form of the function $f(\cdot)$, and Palatini $f(\mathcal{R})$ gravity can be conveniently restated as a metric theory endowed with torsion Olmo:2011uz ; Bombacigno:2018tbo ; Bombacigno:2018tih . Therefore, the additional scalar degree is not dynamical, but it affects the way matter sources and spacetime curvature interact, and also in vacuum the two reformulations are not equivalent, being Palatini case featured by an effective cosmological constant, inherently related to the form of the $f(\cdot)$ function.

As originally proposed in Harko:2011nh , an intriguing perspective is constituted by the possibility of combining both the approaches, considering actions which contain Palatini $f(\mathcal{R})$ modifications to ordinary Einstein-Hilbert metric Lagrangian. Particularly, these theories successfully accomplish the result of providing long range scalar mode, able to reproduce dark matter effects Capozziello:2013uya ; Capozziello:2012qt ; Capozziello:2013yha , without violating solar system observational constraints and invoking the so-called “chameleon mechanism”Khoury:2003rn ; Brax:2008hh ; Capozziello:2007eu . Furthermore, cosmological solutions have been extensively investigated, obtaining accelerated expansion scenarios Carloni:2015bua ; Leanizbarrutia:2017xyd , and studies about compact objects and spherically symmetric static configuration have been performed Danila:2016lqx ; Danila:2018xya .

Here, we deal with a further generalization of these mixed models, and consider a scalar action as in Tamanini:2013ltp ; Rosa:2017jld , where the function $f$ is assumed to depend on both the Ricci scalars, metric and Palatini ones. Especially, we analyze in detail the features of the theory in the weak field limit. It is easy to recognize that in such a type of theory, the scalar-tensor representation is still possible, but now two distinct scalar degrees of freedom come out. These non-minimally coupled scalar fields are dynamically characterized by the form of the potential term they obey as a result of the form of the original function $f({R,\mathcal{R}})$. Critical points of the potential (minima, maxima and saddle points) are relevant for the local gravitational field dynamics, as it is concerned in the PPN limit or when the propagation of gravitational waves is taken into account.

We analyze situation in which the two emerging massive scalar modes have, in the diagonal representation, well-defined masses, ruling out of the theory the non-physical situations in which tachyon modes are present (see Koivisto:2013kwa for a comparison).

On the level of PPN analysis, we show how General Relativity can be still recovered with high degree of precision in the Solar System, as far as the theory parameters are suitably constrained, also in the presence of long range scalar interaction, which can be adopted in principle to reproduce dark matter effects. Then, we analyze the propagation of the gravitational waves in the presence of the two additional massive scalar modes. The deformation of the standard wave polarizations is investigated in some detail for a rather general spanning of the parameter space. In particular, we discuss the intriguing case of nearly degenerate massive modes, and we study the very peculiar phenomenon of wave beating. Such a beating mode is a very striking track of the considered modified theory of gravity and it suggests that upper limits on the existence of mixed $f(R,\mathcal{R})$ model can be experimentally put via present and incoming interferometer devices TheLIGOScientific:2016src ; Abbott:2017tlp ; Abbott:2018utx ; Capozziello:2008rq ; Chatziioannou:2012rf ; Isi:2015cva ; Maselli:2016ekw ; Zhang:2017sym ; Blaut:2019fxb .

The paper is organized as follows. In Sec. II extended hybrid metric-Palatini models are briefly discussed and their scalar tensor representation introduced; in Sec. III we analyze the first PPN corrections in the static weak field limit, pointing out the appearing of Yukawa corrections to gravitational potential given by both the additional scalar degrees; in Sec. IV we address the propagation of gravitational waves in vacuum, investigating to some extent the theory structure in order to have well-defined physical modes; in Sec. V we study the effects on geodesic deviation equation of scalar fields, tracing analogies with metric $f(R)$ theories, and discussing the settling of beating phenomena; in Sec. VI we refine the analysis of Secs. III and IV putting several constraints on the form of the function $f(R,\mathcal{R})$; in Sec. VII conclusions are drawn.

II Extended hybrid metric-Palatini theories

Formerly introduced in Flanagan:2003iw ; Tamanini:2013ltp , extended hybrid metric-Palatini theories are described by the action111We set $\kappa=8\pi G$ and $c=1$.

[TABLE]

where $S_{M}$ stands for the generic matter contribution and the function $f(R,\mathcal{R})$ is assumed to depend on both the metric and Palatini Ricci scalars, denoted by $R$ and $\mathcal{R}$, respectively. Accordingly, we deal with two different kind of affine connections, i.e. the standard Levi-Civita connection related to the metric curvature scalar222We adopted the mostly plus spacetime signature $(-1,+1,+1,+1)$ and the following convention for the Riemann tensor: $R\indices{{}^{\mu}_{\nu\rho\sigma}}=\partial_{\rho}\Gamma\indices{{}^{\mu}_{\nu\sigma}}-\partial_{\sigma}\Gamma\indices{{}^{\mu}_{\nu\rho}}+\Gamma\indices{{}^{\mu}_{\tau\rho}}\Gamma\indices{{}^{\tau}_{\nu\sigma}}-\Gamma\indices{{}^{\mu}_{\tau\sigma}}\Gamma\indices{{}^{\tau}_{\nu\rho}}$, with $R_{\mu\nu}=R\indices{{}^{\rho}_{\mu\rho\nu}}$. $R$:

[TABLE]

with

[TABLE]

and the independent connection $\tilde{\Gamma}\indices{{}^{\rho}_{\mu\nu}}$ defining the Palatini Ricci scalar

[TABLE]

The form of $\tilde{\Gamma}\indices{{}^{\rho}_{\mu\nu}}$ can be dynamically determined by evaluating its equation of motion from (1). Indeed, under the assumption that the matter fields only minimally couple with the metric, it results in

[TABLE]

where $f_{\mathcal{R}}\equiv\partial_{\mathcal{R}}f$ (similarly for $R$ and higher order derivatives), and $\tilde{\nabla}_{\rho}$ denotes the covariant derivative from $\tilde{\Gamma}\indices{{}^{\rho}_{\mu\nu}}$.

Especially, if we neglect the issue concerning the role of the torsion in Palatini $f(\mathcal{R})$ models (see Olmo:2011uz for a review and Bombacigno:2018tbo ; Bombacigno:2018tih for specific applications), the solution of (5) is given by

[TABLE]

which represents the Levi-Civita analougous for the conformal metric $\tilde{g}_{\mu\nu}\equiv f_{\mathcal{R}}\,g_{\mu\nu}$, and allows us to recast the Palatini Ricci scalar as

[TABLE]

Of course, since the function $f_{\mathcal{R}}$ can contain in principle both the curvature scalars, in (7) is actually embedded a differential equation relating $R$ and $\mathcal{R}$. Therefore, in order to make definition (6) well-grounded, we have to provide a further relation between the metric Ricci scalar $R$ and the Palatini curvature $\mathcal{R}$. That can be accomplished by tracing the equation of motion for the metric field $g_{\mu\nu}$, i.e.

[TABLE]

yielding to

[TABLE]

Then, since $f_{R}$ is function of both the metric and the Palatini scalars as well, relation (9) constitutes a second differential equation for $R$ and $\mathcal{R}$, which along with (7) forms a set of highly coupled differential equations for the two different curvatures, here rewritten for the sake of clarity:

[TABLE]

This suggests that in extended hybrid models we actually deal with two independent additional degrees of freedom, somehow connected to the two type of scalar curvatures the theory is equipped with. It is worth noting that this property depends crucially on the form of (9), with special focus on $f_{R}$ contributes. Indeed, in hybrid models originally discussed in Harko:2011nh we simply have $f_{R}=1$, and (9) reduces to an algebraic constraint relating the Palatini curvature both to the metric Ricci scalar and the trace of the stress energy tensor, i.e.

[TABLE]

which can be solved in principle for $\mathcal{R}=\mathcal{R}(R,T)$. This in turn implies that the first of (10) boils down to a differential equation for the metric scalar R, in the presence of non trivial stress energy source terms, and we just retain an additional degrees of freedom. We point out that such an outcome is to some extent preserved also in extended hybrid theories, when metric and Palatini terms in the function $f$ are actually separable, i.e. $f(R,\mathcal{R})=f_{1}(R)+f_{2}(\mathcal{R})$. In this case, in fact, constraint (11) generalizes to

[TABLE]

and the Palatini scalar $\mathcal{R}$ turns out to depend also on the derivatives of $R$, i.e. $\mathcal{R}=\mathcal{R}(R,\nabla R,\Box R,T)$. Again, the only additional dynamical degree is still $R$, even though its evolution is now encoded by a higher order differential equation.

In the following we will restrict our attention to the general case $f_{R\mathcal{R}}\neq 0$, where $R$ and $\mathcal{R}$ represent truly independent degrees of freedom, whose dynamics is described by fourth-order differential equations as displayed in (10).

II.1 Scalar-tensor formulation

As discussed in Tamanini:2013ltp , if the determinant of the Hessian matrix for $f(R,\mathcal{R})$ is not vanishing, i.e.

[TABLE]

action (1) can be rearranged in the scalar-tensor form

[TABLE]

where we introduced the scalar fields $\psi\equiv f_{R},\;\xi\equiv f_{\mathcal{R}}$, together with the potential term

[TABLE]

Hence, taking into account (7) and the definition of $\xi$, (14) can be rearranged in its gravitational part as

[TABLE]

Then, defining a new scalar field $\phi=\psi+\xi$, we can finally rewrite (16) in the form

[TABLE]

where $W(\phi,\xi)\equiv V(\phi-\xi,\xi)$ and only the scalar $\phi$ is coupled to the metric Ricci scalar.

Then, varying (16) with respect to the metric field we get

[TABLE]

while equations for $\phi$ and $\xi$ are given by, respectively

[TABLE]

Now, we can evaluate the trace of (18), resulting in

[TABLE]

which plugged into (19)-(20) leads to the following set of coupled differential equations

[TABLE]

where $W_{\phi}$ and $W_{\xi}$ are defined by analogy with $f_{R,\mathcal{R}}$. This set of equations is the scalar-tensor equivalent of (10), with the additional scalar degrees of freedom now embodied in the independent fields $\phi,\,\xi$. With this respect, even if the transformation $\phi=\psi+\xi$ seems to artificially relate the degrees $\phi$ and $\xi$, it actually preserves the dynamical content of the theory, a fact that can be further appreciated by evaluating the equations of motion for the original fields $\psi$ and $\xi$ directly from (14), i.e.

[TABLE]

This in turn guarantees that when interested in perturbation theory, the departure of $\phi$ and $\xi$ from background values could be considered truly independent, as long as they represent proper dynamical variables, as previously discussed when we look at (10).

III Post parameterized Newtonian corrections

It is a well-established result (see Olmo:2005hc ; Harko:2011nh ) that additional scalar degrees can remarkably affect the dynamics of gravitating system in weak field and slow motion case. In particular, Yukawa corrections are usually obtained for the Newtonian potential, and the requirement of reproducing local experiment results allows in general to put several constraints on theory parameters Zakharov:2006uq ; Chiba:2006jp ; Berry:2011pb . In this regard, the easiest way of determining the effects of the fields $\phi$ and $\xi$ in a slightly curved spacetime is to consider a quasi-Minkowskian system of local coordinates where the metric can be put into the form

[TABLE]

with $\left|h_{\mu\nu}\right|\ll 1$, and the scalar fields $\phi,\,\xi$ are given by

[TABLE]

Here $\phi_{0}$ and $\xi_{0}$ represent background values fixed by cosmological boundary conditions, which evolve adiabatically in time according the cosmological background curvature. Local fluctuations are denoted by $\delta\phi,\,\delta\xi\sim\mathcal{O}(h)$, which we assume to vanish outside the region described by (26). Now, by virtue of (27) we can expand the potential $W(\phi,\xi)$ as:

[TABLE]

where the subscript [math] denotes evaluation at the point $(\phi_{0},\,\xi_{0})$, which we require to be located in the neighbourhood of a stable minimum for $W(\phi,\xi)$, so that the smallness of corrections $\delta\phi$ and $\delta\xi$ be preserved by the dynamics. We assume therefore that

[TABLE]

where we introduced $\text{H}(W_{0})$, the Hessian matrix for the potential $W(\phi,\xi)$ evaluated at the point $(\phi_{0},\xi_{0})$. Moreover, given that the value $W_{0}$ is related to the background curvature by

[TABLE]

as it can be inferred by taking the lowest order in (21), quasi-Minkowskian conditions imply that locally we can set $R^{(0)}=\epsilon$, with $\epsilon$ a small parameter quantifying departure from spacetime flatness. In particular, since it turns out to be actually responsible for a divergent term in the PPN corrections to the gravitational potential (see later), we constraint it to be small enough so that its contribution is negligible on the considered scales, i.e. wherever (26) is valid.

Now, let us recast (18) into the form

[TABLE]

which, once we fixed the Nutku gauge Nutku

[TABLE]

gives at the lowest order in perturbation the following equations for the metric $h_{\mu\nu}$ components:

[TABLE]

where we neglected time derivatives, we set $\bigtriangleup=\nabla^{2}$ for the Laplacian operator and $T_{00}=-T\approx\rho$, $T_{ij}\approx 0$. Analogously, we can rearrange (22) and (23) as

[TABLE]

with

[TABLE]

and zero order terms satisfying

[TABLE]

Now, in order to solve (35)-(36) is useful to find a suitable change of variables with the aim of decoupling the equations of motion for $\delta\phi,\,\delta\xi$. That can be accomplished by introducing the matrix

[TABLE]

and the vectors

[TABLE]

which allow us to rewrite the set (35)-(36) as

[TABLE]

where $\text{I}_{2\times 2}$ denotes the identity matrix of dimension 2.

Then, the system (41) can be decoupled, that is A turned into diagonal form, by simply evaluating the matrix P of its eigenvectors. Therefore, let us rearrange (41) like

[TABLE]

with $\text{A}_{\text{D}}\equiv\text{P}^{-1}\text{A}\text{P}$ diagonal and $\mathbf{\Phi_{D}}\equiv\text{P}^{-1}\mathbf{\Phi}$. We observe that the stress energy contributions to (35)-(36) are shuffled, so that we expect that matter sources could enter now both the equations for the decoupled scalar fields.

Let us denote the elements of P and $\text{P}^{-1}$ with

[TABLE]

and $\mathbf{\Phi}_{D}=(\delta\phi_{D},\delta\xi_{D})$. It is then possible to rewrite (42) as

[TABLE]

where $M_{\phi,\xi}$ are the masses of the decoupled scalar fields which we require to be positive and that still have to be determined explicitly. Solutions for (44)-(45) can be easily obtained, i.e.

[TABLE]

where the integration is performed over the matter source. Hence, the solution for (33)-(34) can be written down, noting that the field $\delta\phi$ is actually a linear combination by means of P of the decoupled modes $\delta\phi_{D},\,\delta\xi_{D}$, i.e.

[TABLE]

which leads to

[TABLE]

where $x_{S}$ is an integration constant related to the source. Now, in spherically symmetric case and far away from the source, $\delta\phi(x)$ and the metric perturbations take the simpler form

[TABLE]

with $M_{\odot}$ the Newtonian mass of the central body and the modified gravitational constant defined as

[TABLE]

We also introduced the PPN $\gamma$, given by

[TABLE]

It is now evident from (53)-(54) that the parameter $\epsilon$ must satisfy

[TABLE]

as long as the point $r$ lays in the region described by (26). Analogously, solar system measurements Will:2005va constraint $\gamma\simeq 1$ and in contrast with standard metric $f(R)$ predictions Olmo:2005hc but by close analogy with Harko:2011nh , we see that in principle such a requirement can be fulfilled also in the presence of long range scalar interactions. Indeed, Yukawa contributes to (56) are tuned by coefficients $p_{11}\bar{p}_{11},\,p_{12}\bar{p}_{21}$, related to the potential expansion (28), which can make the corrections due to the scalar field negligible also for very low masses. However, by virtue of $\text{P}\text{P}^{-1}=\text{I}$, they are not truly independent but are compelled to satisfy the condition

[TABLE]

so that we can a priori clearly distinguish two different scenarios. In the first case, taking $p_{11}\bar{p}_{11}$ and $p_{12}\bar{p}_{21}$ of the same magnitude we are forced to consider large masses for both the scalar fields in order to recover $\gamma\simeq 1$. In the second case, on the contrary, setting one of the coefficient nearly vanishing, we can retain a low mass mode which can affect astrophysical and cosmological scales. Lastly, we note that to preserve the attractive behaviour of gravity at the leading order implies, by virtue of (55), that condition $\phi_{0}>0$ be satisfied.

III.1 Case $W_{0,\xi}=0$

In order to see if a configuration characterized by a low mass mode is actually attainable, we have to go back to conditions (38), which besides fixing the value of $W_{0,\phi}=\epsilon$, also implies either $\xi_{0}=0$ or $W_{0,\xi}=0$333We do not take into account the very special case $\xi_{0}=W_{0,\xi}=0$, when the matrix A has a vanishing row and the scalar fields are already decoupled, with $\delta\xi$ massless.. Especially, when $W_{0,\xi}=0$ the matrices P, $\text{P}^{-1}$ turn out to be, respectively

[TABLE]

and

[TABLE]

where $U$ is defined as

[TABLE]

Then, combining (39) and (59)-(60), the masses of the decoupled scalar fields can be obtained, i.e

[TABLE]

and the coefficients ruling the Yukawa corrections evaluated

[TABLE]

Now, in order to assure the existence of the function $U$ together with the positivity of (62)444We are disregarding configurations where $M^{2}_{\phi},\,M^{2}_{\xi}<0$, which would lead to an oscillatory behaviour for $\gamma$ (see Olmo:2005hc ) and to tachyonic instabilities in gravitational waves propagation (see Sec. IV., the following set of inequalities must hold

[TABLE]

which combined yields to

[TABLE]

In particular, the reality of $U$ implies

[TABLE]

which is always satisfied for $\xi_{0}<0$, while for a positive $\xi_{0}$ can be rewritten as

[TABLE]

Instead, the inequality $U(\phi_{0},\xi_{0})\leq 3\left(m_{\phi}^{2}+m_{\xi}^{2}\right)$ states that

[TABLE]

whereas the fact that the sum of the masses of the coupled modes is compelled to be positive implies

[TABLE]

In the case $\xi_{0}>0$ condition (69) yields to

[TABLE]

which, combined with (70), implies the redundancy of (68), hence the latter shall not be considered in the following. Now, according the value of $U$ the mass spectrum can exhibit quite different behaviour. Indeed, when $U\neq 0$ the masses of the decoupled modes are distinct, i.e.

[TABLE]

and the spectrum is not degenerate. Then, since we are interested in peculiar scenarios where at least one scalar field is long range, we can take $U\simeq U_{max}\equiv 3(m_{\phi}^{2}+3m_{\xi}^{2})$, for which the scalar field $\delta\xi_{D}$ is very tiny while the mode $\delta\phi_{D}$ is endowed with the mass $M_{\phi}^{2}\simeq m_{\phi}^{2}+m_{\xi}^{2}$. We also note that for $U\simeq U_{max}$, relation (69) implies

[TABLE]

which taken into account (29) leads to $W_{0}>0$.

Now, for this arrangement of the masses the coefficients (64) boils down to

[TABLE]

and choosing $m_{\phi}^{2}\gg m_{\xi}^{2}$, i.e.

[TABLE]

we can reproduce the conditions $p_{11}\bar{p}_{11}\simeq 1$ and $p_{12}\bar{p}_{21}\simeq 0$, which are compatible with the requirement of preserve the observed dynamics at local scales, provided we properly set $m_{\phi}^{2},\,m_{\xi}^{2}$.

III.2 Case $\xi_{0}=0$

If $\xi_{0}=0$ the function $f(R,\mathcal{R})$ has no linear contribution in the Palatini scalar $\mathcal{R}$ and the matrix A turns in triangular form. In this special case, the matrices P, $\text{P}^{-1}$ read as

[TABLE]

where $p\equiv\frac{2W_{0,\xi}-\phi_{0}W_{0,\phi\xi}}{\epsilon-W_{0,\xi}-\phi_{0}W_{0,\phi\phi}}$. Since now $\bar{p}_{11}=1,\;\bar{p}_{21}=0$, the matter source is not shuffled and (44)-(45) take the simpler form

[TABLE]

with

[TABLE]

The condition for the masses of the decoupled fields to be real leads in this case to

[TABLE]

Solutions can be written as

[TABLE]

where we normalized conveniently the expression for $\delta\xi_{D}$. Then, far away from the central source solutions are still given by (53)-(54), where now

[TABLE]

and

[TABLE]

In this situation, since the Yukawa correction due to the scalar mode $\delta\phi_{D}$ cannot be tuned, likewise ordinary metric $f(R)$ case, we are compelled to consider configurations in which it is very massive and its contribute appreciable only at short scales. Conversely, properly setting the parameter $p$ we can still have a light $\delta\xi_{D}$ mode, provided $2W_{0,\xi}\simeq\phi_{0}W_{0,\phi\xi}$, where $p\simeq 0$.

IV Gravitational waves propagation

Now, we consider (26) exact and we study the propagation of gravitational degrees of freedom on a globally flat spacetime. In this case, the consistency at the lowest order for equations (18), (22) and (23), requires that the configuration $(\phi_{0},\xi_{0})$ be also a zero, beside a stable minimum, for the potential $W$, i.e.

[TABLE]

Hence, restricting our attention to the vacuum case ($T_{\mu\nu}=0$), the equation of motion for the metric field is given by

[TABLE]

where we did not fix a priori any gauge conditions and $R_{\mu\nu}^{(1)}$ and $R^{(1)}$ are the Ricci tensor and the Ricci scalar expressed at first order in $h_{\mu\nu}$. The linearized equations for $\phi$ and $\xi$ instead turn out to be, respectively

[TABLE]

where $\Box$ is the D’Alambert operator $\Box\equiv\partial_{\mu}\partial^{\mu}$.

We point out that (86) features the same form of metric $f(R)$ theories Moretti:2019yhs , even if the dynamical degree $\phi$ actually satisfies the remarkably different equation (87), which still at the linear order is coupled with the corresponding equation (88) for $\xi$. They represent a pair of coupled wave equations for massive fields, so that extended hybrid metric-Palatini gravity seems to be characterized in vacuum by two further propagating degrees of freedom in addition to the ordinary tensorial ones. Of course, the theory could be in principle affected by instabilities concerning possible tachyonic modes, and in this respect we will see that there exists a suitable region into the parameter space of the theory, where both the modes are allowed to propagate.

IV.1 Decoupling of the wave equations

Following the analysis made in Sec. III, the set (87)-(88) can be rearranged into the form

[TABLE]

where now B is given by

[TABLE]

For $W_{0,\phi\xi}\neq 0$ the set must be turned to diagonal form, and that can be accomplished provided $\phi_{0},\,\xi_{0}\neq 0$. Again, it is possible to rearrange (89) into the form

[TABLE]

with $\text{B}_{\text{D}}\equiv\text{P}^{-1}\text{A}\text{P}$, where the matrix P is still given by (59), provided we replace

[TABLE]

Then, the following set of decoupled equations for $\mathbf{\Phi}_{D}$ can be written down

[TABLE]

with $M_{\phi}^{2},\,M_{\xi}^{2}$ as in (62), taken into account (92). Now, in order to assure that (93)-(94) actually describe propagating physical fields, the set of inequalities (67)-(65c) can be restated, by writing explicitly $U$ and squaring (66), as

[TABLE]

Since we are interested in stable minimum configurations (29), from (95b) it follows that $\phi_{0}$ and $\xi_{0}$ have to exhibit opposite sign. Thus, considering (70) for $W_{0,\phi}=0$ (i.e. $\epsilon=0$), this in turn implies that the only possible case satisfying all the criteria is

[TABLE]

In fact, when (96) holds, relations (70) and (95b) are strictly satisfied and also the squared masses of the non diagonal modes (92) turn out to be positive. There exist, then, suitable configurations of the theory, corresponding to peculiar minima for the potential $W(\phi,\xi)$, characterized by two additional scalar degrees of freedom which propagate like linear waves on a Minkowski background. Of course, since the potential $W$ is ultimately related to the functional form $f(\cdot)$ by means of (15), this selects specific classes of $f(R,\mathcal{R})$ models and the existence of one or more propagating scalar degrees could be not in general guaranteed (we remind the reader to Sec. VI for details). Thus, if we restrict our attention to $f$ functions able to produce these scalar waves, we see that for $U\neq 0$ the masses corresponding to the scalar modes are distinguished, with $M_{\phi}^{2}>M_{\xi}^{2}$ for every value of $U$ between [math] and $U_{max}$. The specific configuration $U=U_{max}$, where the mode $\xi_{D}$ is predicted to become massless, has to be instead disregarded. Indeed, in that condition (95b) would imply $\det\text{H}(W_{0})=0$, where as discussed in Sec. VI the scalar-tensor representation is not valid, being $\det\text{H}(f)=\infty$. We have to restrict therefore the study to the case $\det\text{H}(W_{0})=\Delta U$, with $\Delta U$ a positive small parameter, where with a bit a manipulation can be shown that

[TABLE]

and the decoupled scalar modes are endowed with the squared masses

[TABLE]

When $U=0$, instead, by virtue of (61) and (96), the following constraints have to be separately satisfied:

[TABLE]

and it follows from (39) that we actually deal with a system already decoupled. The procedure involving the definition of $U$, therefore, is not well grounded, and we cannot simply perform the limit of $U\rightarrow 0$ in (62), that would result in the degenerate spectra

[TABLE]

Rather, if $W_{0,\phi\xi}=0$, we just retain (92), where the masses could be in principle different: The mass spectrum is affected by a discontinuity for $U=0$, where the masses of the actual physical modes do not coincide with the values predicted by (62).

V Geodesic deviation

In order to analyze the phenomenology of gravitational waves in extended hybrid theories, we can evaluate, via the geodesic deviation equation, the perturbations induced by the scalar modes on a sphere of test masses. These are displayed along with the tensorial degrees in

[TABLE]

where we introduce the vector

[TABLE]

denoting the separation between two nearby geodesics, with $x_{0}$ and $\delta_{x}$ indicating the rest position and the displacement of order $\mathcal{O}(h)$ induced by waves, respectively555Analogously for $y,z$.. Then, following Flanagan:2005yc , we introduce for the metric perturbation $h_{\mu\nu}$ the generic decomposition for a symmetric tensor of rank two, i.e.

[TABLE]

with $\delta_{ij}$ the Kronecker delta, $\bigtriangleup=\partial_{i}\partial^{i}$ the Laplacian operator and symmetrization given by $A_{(ij)}\equiv\frac{1}{2}(A_{ij}+A_{ji})$. The irreducible parts introduced in (103) are accompanied by the conditions

[TABLE]

which, as stressed in Flanagan:2005yc (see also Weinberg:2008zzc for the curved background case), are required in order to preserve the uniqueness and the consistency of the procedure. By means of these quantities we can then introduce the set of variables

[TABLE]

which turns out to be invariant, together with $h_{ij}^{TT}$, under a linear gauge transformation.

Now, as it was outlined in Moretti:2019yhs for metric $f(R)$ theories, it is possible to rearrange the linearized equation for the metric (86) into the form

[TABLE]

where we introduced the modified static degrees

[TABLE]

From (106), it is evident that beyond the scalar degrees discussed in Sec. IV, we retain the standard tensorial modes for the metric $h_{\mu\nu}$. Moreover, it is worth noting that with respect to the discussion in Moretti:2019yhs , the scalar field involved into the definition (107) does not represent a proper degree of freedom. Indeed, the quantity $\delta\phi$ is actually related by means of (59) to the diagonal scalar modes $\left(\delta\phi_{D},\,\delta\xi_{D}\right)$, i.e.

[TABLE]

Therefore, when we look at the components of the linearized Riemann tensor entering (101), these can be expressed in gauge invariant variables as

[TABLE]

which can be rewritten, neglecting the static contributions and taking into account (107) and (108), like:

[TABLE]

with $\mathcal{U}^{(+)},\,\mathcal{U}^{(-)}$ given by, respectively

[TABLE]

Now, leaving aside the tensorial degrees contained in $h_{ij}^{TT}$ and choosing the $z$ axis coincident with the direction of propagation of the waves, the scalar degrees can be described by

[TABLE]

with frequencies

[TABLE]

wave vectors fixed in

[TABLE]

and $A_{\phi},\,A_{\xi}$ the amplitudes of the waves. Thus, the geodesic deviation equation takes the form:

[TABLE]

where we disregarded terms of order $\mathcal{O}(h^{2})$.

By close analogy with the discussion in Moretti:2019yhs ; Liang:2017ahj , we see that both the scalar degrees are able to induce two type of polarizations. In fact, they are separately responsible for a breathing mode on the transverse plane $xy$, as well as for a longitudinal excitation along the direction of propagation of the wave. Moreover, the corresponding polarizations are modulated for $U\neq 0$ by factors of distinct magnitude, and $\delta\phi_{D},\,\delta\xi_{D}$ propagate with different speed. In particular, when $U\rightarrow U_{max}$, by virtue of (97) and (98) the longitudinal contribute of $\delta\xi_{D}$ turns out to be or order $\varepsilon$, i.e.

[TABLE]

In this case, therefore, it mostly affects the geodesic deviation as a breathing on the plane $xy$ and the longitudinal polarization is almost entirely due to the massive mode $\delta\phi_{D}$.

Conversely, when $U\rightarrow 0$ the angular frequencies of the scalar modes are very close to each other, and we expect that typical interference patterns between waves, i.e. beatings, could take place. That can be considered a very distinctive marker of gravitational wave propagation in generalized hybrid metric-Palatini theories, absent in ordinary metric $f(R)$ gravity, with specific phenomenological implications.

Thus, let us write for $U\simeq 0$ the solution of (115) as

[TABLE]

where we set $z=0$ and effective amplitudes $A_{\phi,\xi}^{B,L}$ given by

[TABLE]

with $A_{\phi}\sim A_{\xi}=A$ and $K_{\phi}\sim K_{\xi}=k$. After a bit of manipulation, it is possible to recast666We just report the result for $\delta_{x}$. Similar considerations hold for $\delta_{y},\,\delta_{z}$. (117) like

[TABLE]

where we defined $\Delta\Omega\equiv\frac{\Omega_{\phi}-\Omega_{\xi}}{2}$ and $\bar{\Omega}\equiv\frac{\Omega_{\phi}+\Omega_{\xi}}{2}$. Then, the perturbation described by (119) represents a superposition of two waves of frequencies $\bar{\Omega}$ with a phase shift of $\pi/2$, both modulated by the beating frequency $\Delta\Omega$ (Fig. 1), and relative amplitude $A_{\phi}^{B}-A_{\xi}^{B}>A_{\phi}^{B}+A_{\xi}^{B}$.

Finally, when $W_{0,\phi\xi}=0$ the set of equations (87)-(88) is naturally decoupled, and the transformation (108) is no longer necessary. In this case the relevant components of the Riemann are given by

[TABLE]

and we see that only $\delta\phi$ enters the geodesic deviation. Therefore, the phenomenology described is identical to that descending from the scalar-tensor formulation of metric $f(R)$ theories, i.e.

[TABLE]

Nevertheless, even if $\delta\xi$ does not appear explicitly in (121), we cannot infer that the functional dependence of $f(\cdot)$ on $\mathcal{R}$ have no phenomenological implications. Indeed, since $\phi$ is actually the combination of $\psi$ and $\xi$, we clearly see that both contributions of $f$ from $R$ and $\mathcal{R}$ concur in determining the effects of (121).

VI Constraints on the form of $f(R,\mathcal{R})$

In this section we analyze in detail the implications onto the form of the function $f(R,\mathcal{R})$ of conditions discussed in Sec. IV-V. In particular, we are interested in establishing clear relations between derivatives of the potential $W$ and corresponding derivatives of the function $f$ with respect to the curvatures $R$ and $\mathcal{R}$. Then, in order to do that, it is useful to express the derivatives of $W(\phi,\xi)$ in terms of derivatives of $V(\psi,\xi)$, i.e.

[TABLE]

when we considered $W(\phi,\xi)=V(\psi(\phi,\xi),\xi)$. It follows that second order derivatives are given by

[TABLE]

which can be further rewritten, taking into account definitions of the potential $V$, as

[TABLE]

Now, since $\psi$ and $\xi$ are function of $R$ and $\mathcal{R}$ by means of the first derivatives of $f$, in evaluating (127) we can apply the inverse function theorem for the two dimensional case, leading to

[TABLE]

where we use the fact that the Jacobian matrix of the transformation relating $(\psi,\,\xi)$ to $(R,\,\mathcal{R})$ coincides with the Hessian matrix for $f$. We can then express the determinant of the Hessian matrix of $W$ in terms of derivatives of $f$, that is

[TABLE]

In the continuing we will evaluate these quantities at background values $\phi_{0},\xi_{0}$ and we see that condition $\det\text{H}(f)\neq 0$, required for scalar-tensor representation to exist, guarantees that $(R,\mathcal{R})$ could be solved for $(\phi,\xi)$ and therefore computed in $(\phi_{0},\xi_{0})$. Furthermore, from (129) it is evident that we have to disregard from the analysis configurations with $\det\text{H}(W_{0})=0$, corresponding to $\det\text{H}(f)=\infty$, where a propagating massless mode is theoretically predicted for gravitational waves (see Sec. V).

VI.1 Post parameterized Newtonian corrections

In studying the PPN corrections we made the following general assumptions

[TABLE]

which can be translated by means of (127) and (128) in conditions on the derivatives of the function $f$. This results in

[TABLE]

where we used the definitions of $\phi$ and introduced, by analogy with $W$, the subscript 0 also for $f$. We see that second order derivatives constitute an independent subsystem of inequalities, whose solution is given by

[TABLE]

while first order derivatives do not require further manipulations. In the continuing, we will investigate in detail other conditions required in each case discussed in Sec. III.

VI.1.1 Case $W_{0,\xi}=0$, $\xi_{0}<0$

The condition of reality for the function $U$ (67) always holds, while inequalities (69) and (70) can be reformulated as

[TABLE]

where we have outlined the negative sign of $f_{0,\mathcal{R}}=\xi_{0}$. The requirement $m^{2}_{\phi}\gg m^{2}_{\xi}$ that we impose in order to correctly reproduce the Newtonian limit can be restated as

[TABLE]

VI.1.2 Case $W_{0,\xi}=0$, $\xi_{0}>0$

Conditions (70), (71) and (75) can be reformulated as follows

[TABLE]

where in the first inequality we have omitted $\epsilon^{2}$ terms. It is immediate to translate this set of inequalities in terms of constraints on the derivatives of the function $f$, yielding to

[TABLE]

VI.1.3 Case $\xi_{0}=0$

Conditions (80) reads, in the special setting $p\simeq 0$, as

[TABLE]

In order to reproduce PPN corrections which are compatible with local measurements we have to impose $M_{\phi}\bar{r}\gg 1$, with $\bar{r}$ much smaller than the size of the source. By combining this request with (57) we get

[TABLE]

VI.2 Gravitational wave modes

As we saw in Sec. V, for stable minima of $W(\phi,\xi)$ two additional massive scalar modes are expected to propagate, provided the set (130) of inequalities hold, along with the additional condition $\xi_{0}<0$. In terms of $f$ this leads to

[TABLE]

and as in the PPN case the constraints on the second derivatives are satisfied if

[TABLE]

Finally, for the special setting $W_{0,\phi\xi}=0$, where the scalar modes are already decoupled and equipped with massed as in (92), conditions (139) are endowed with the further requirement $f_{0,R\mathcal{R}}+f_{0,\mathcal{R}\mathcal{R}}=0$. That, plugged back into (139), leads then to $f_{0,RR}>f_{0,\mathcal{R}\mathcal{R}}=-f_{0,R\mathcal{R}}>0$.

VII Concluding remarks

Extended hybrid metric-Palatini theories are promising generalizations of the two main approaches to the study of $f(R)$ gravity, and the most intriguing feature of these models is certainly the presence of two dynamical scalar fields non minimally coupled to gravity. In fact, this enrichment of the dynamical structure can be used, in principle, to remove technical and conceptual problems, common to both the metric and the Palatini approach, that arise when one tries to mimic dark matter or dark energy effects without spoiling Solar System tests. In this work we investigated the weak field limit of the theory in its scalar-tensor formulation, by analyzing the first PPN order and the gravitational wave propagation. In both cases we found the two scalar fields solve coupled equations, and the decoupled scalar fields are shown to be massive, with masses that vary in a range determined by the potential $W$. Particularly, the mass spectrum spans continuously the interval that goes from two modes having nearly the same mass to the case in which one field has the maximum mass (again determined by $W$) while the other field is massless. In this respect, we clarified how this peculiar setting depends crucially on the parameter $\epsilon$ quantifying the departure at the background level from Minkowski spacetime. We showed, in fact, that when the propagation of gravitational scalar waves is addressed, corresponding to $\epsilon=0$, the configuration where $\delta\xi_{D}$ is massless is not actually feasible, in that scalar tensor representation is not attainable. With regard to the PPN expansion we found that the presence of the massive fields implies that the parameters $G_{eff}$ and $\gamma$ acquire Yukawa-like corrections. The intensity of these modifications is governed by coefficients that can be tuned through specific constraints on the potential function $W$. We claim that it is possible to make corrections in the expressions of the PPN parameters small enough to stay within the constraints of current Solar System tests, still having the presence of a scalar massive field light enough to act as dark matter on galactic scales. This can be accomplished by choosing $W$ such that the masses of the scalar fields are widely separated: with this expedient the mass of the heavier can be set to a value that implies the suppression of the relative exponential factor over a convenient scale, while the lighter can be forced to have a decay length comparable with galactic scales. The correction relative to the light scalar in the expressions of $G_{eff}$ and $\gamma$ can be made small enough through a precise choice on the corresponding coefficient. For what concerns the gravitational wave study, we performed a linear metric approximation around a Minkowski background, hence we restricted the dynamics of the scalar fields to small oscillations around a local minimum of $W$. We showed that the decoupled fields solve two independent Klein-Gordon equations with masses varying in the above mentioned range. The analysis of the phenomenology associated to the scalar fields, performed via the geodesic deviation equation for a sphere of test particles, demonstrated that each field is detectable as the superposition of two independent polarizations, namely a breathing plus a longitudinal mode. It must be stressed that such a finding cannot be claimed to be a real marker for this specific model: indeed in Liang:2017ahj ; Moretti:2019yhs is shown that in the metric formalism the only additional scalar field is responsible for the same mixture of polarizations, whereas in Montani:2018iqd is demonstrated that in General Relativity gravitational waves travelling in molecular media, like galaxies, are expected to contain additional mode characterized by the same feature. The striking peculiarity of this model is instead the fact that the two scalar fields can mutually interact and produce beatings. We studied this phenomenon in the special case of nearly degenerate masses, in which the beating frequency is much smaller than the signal frequency. However, the same feature should be detectable for any value of the masses in the allowed range, at least until the lighter scalar can be properly distinguished from the massless tensorial degrees, in which case we expect their mutual interaction as in Bombacigno:2018tih . Finally, we established precise relations between the constraints on the potential $W(\phi,\xi)$, obtained in analyzing PPN and gravitational wave settings, and the form of the function $f(R,\mathcal{R})$. These resulted in a set on inequalities connecting second and first derivatives of $f$ with respect both the curvatures. In particular, they could be used in principle for the construction of a definite model meeting the requirements that are necessary to mimic dark matter effects and pass the Solar System tests, paying special attention to already known potentials (as discussed in Rosa:2017jld ) characterized by accelerating cosmological solutions as well.