1 Introduction↩︎

The \(\Lambda\)CDM concordance cosmological model, based on general relativity, associated with the cold dark matter sector and the cosmological constant, is a well-defined model that has been highly effective in describing the evolution of the Universe, accurately reproducing observations such as the cosmic microwave background and the distribution of large structures [1]. However, the \(\Lambda\)CDM model presents theoretical and phenomenological limitations that challenge its consistency and motivate new investigations [2]. From a theoretical point of view, fundamental questions such as the origin and nature of the cosmological constant and the lack of an explanation for dark matter remain open. The model also presents difficulties in explaining certain observational discrepancies, such as the tensions in the determination of the Hubble constant, the behavior of structures on small scales, and several other less discussed tensions [2], [3]. These limitations reinforce the need to study beyond the standard model, driving the development and exploration of alternative theories and extensions to the \(\Lambda\)CDM paradigm.

One of the main tensions related to \(\Lambda\)CDM concerns the current value of the Hubble parameter \(H_{0}\). According to the most recent global measurements by the Planck Collaboration in 2020 (P20), \(H_0 = (67.40 \pm 0.50)~\mathrm{km} \cdot \mathrm{s}^{-1} \cdot \mathrm{Mpc}^{-1}\) was estimated at the confidence level \(68\%\) [4]. In contrast, the latest local measurements performed by the SH\(0\)ES collaboration in 2021 (R21), based on Cepheid-calibrated supernovae, provided a significantly larger value, estimated at \(H_0 = (73.04 \pm 1.04)~\mathrm{km} \cdot \mathrm{s}^{-1} \cdot \mathrm{Mpc}^{-1}\) at confidence level \(68\%\), representing a discrepancy of about \(5.0 \sigma\) between both \(H_{0}\) values [5]. Although extensive discussions have sought to determine whether this tension can be attributed to unidentified systematic errors, several observations with other alternative methods have also indicated a tension in \(H_{0}\) inferred from \(\Lambda\)CDM [6]–[10]. Thus, there is growing evidence that the discrepancy in the \(H_{0}\) values may indeed be an indication of new physics beyond the Standard Model. Recently, several studies have tried to solve this question, using different methods [2], [3], [11]–[16]. Since the tension in \(H_{0}\) suggests a faster expansion of the universe than predicted by \(\Lambda\)CDM, a promising approach to mitigate this discrepancy involves the use of a modified gravity theory. Such a theory should qualitatively induce a reduction in gravitational intensity during the middle and late epochs of the cosmic expansion [17]–[19].

In recent years, several studies have investigated specific subclasses of Horndeski theories in an effort to reconcile conflicting measurements of the Hubble constant \(H_{0}\) [17]–[20]. Horndeski gravity is regarded as one of the most general scalar-tensor theories that preserve second-order equations of motion, thereby avoiding dynamical instabilities. This framework allows for the inclusion of a scalar field that interacts directly with gravity, providing the flexibility needed to modify the Universe expansion dynamics across different cosmological epochs. In this work, we explore the application of Horndeski gravity [21], focusing on a specific subclass known as the John sector, which is a non-minimal derivative coupling between the scalar field and the Einstein tensor, by adopting an approach that aims to mitigate the observational tension around \(H_{0}\). Two scenarios are considered. In the first, dark matter is described by a scalar field, as discussed in [22] and in the second, the scalar field plays the role of dark energy. Through this, estimates are possible for the values of some important cosmological quantities, such as the value of \(H_0\) itself and the current matter energy density parameter \(\Omega_m\).

This sector of Horndeski gravity allows a specific interaction between the scalar field and the spacetime curvature, enabling variations in the strength of gravity over time. Such behavior can contribute to the accelerated expansion of the Universe without relying exclusively on the cosmological constant. In recent years, this sector of theory has been widely investigated in some cosmological contexts, covering several areas of interest, such as the construction of domain wall solutions [23], simulations involving dark energy [24], and dark matter [22], [25]. Moreover, numerous studies have been carried out in the context of black holes, including locally asymptotically AdS and planar solutions [26], [27], rotating black holes with probe strings [28], and thermodynamic analyzes of static solutions [29]. In the realm of compact objects, neutron star configurations have been developed, even accounting for slow rotation [30], [31]. Finally, this sector has also been explored in the search for braneworld solutions [32]–[34]. Thus, this approach has proven to be both promising and relevant in the literature, being a possible viable path toward alleviating the Hubble tension, by adjusting the late-time expansion rate of the universe and bringing it into better agreement with both local and global measurements of \(H_{0}\).

The paper is organized as follows. In Sec. 2 we briefly review Horndeski gravity, its field equations, and some conditions for its cosmological viability. In Sec. 3 we address the non-minimal coupling sector between the cosmological solutions of this class of Horndeski theories. To do so, we start from the flat FRW metric and find the differential equations to be solved in this dynamic. Next, in Sec. 4 we present the numerical and statistical analysis performed on the investigated model, obtaining statistical results in relation to some important cosmological quantities, in addition to the free parameters of the model, addressing issues such as a possible alternative for the relief of the Hubble tension. In Sec. 5 we address the issues of Laplacian instabilities and ghosts in the model. Finally, in Sec. 6 we have the conclusions.

2 Horndeski Gravity↩︎

In this section, we briefly discuss the Horndeski gravitational theory, introduced into the literature in 1975 by Gregory Horndeski [21]. In recent years, this theory has been widely applied in several studies involving modified gravity, particularly in topics related to cosmology. This theory was recently rediscovered in the context of generalizations of the Galileon models and represents the most general action for a scalar-tensor theory in a four-dimensional spacetime, with second-order field equations.

2.1 The Lagrangian↩︎

The Horndeski theory is characterized by the action \[\label{acaohorndeski} S_{H}[g_{\mu \nu}, \phi] = \int d^{4}x \sqrt{-g} \, \mathcal{L}_{H}[g_{\mu \nu}, \phi],\tag{1}\] where \(g\) is the determinant of the metric \(g_{\mu \nu}\) and \(\mathcal{L}_{H}\) is the Horndeski Lagrangian [21], given by \[\label{eqlagHorn} \mathcal{L}_{H} = \sum_{i=2}^{5} \mathcal{L}_{i},\tag{2}\] where \[\begin{align} \mathcal{L}_2 &= G_{2}(\phi, X), \tag{3} \\ \mathcal{L}_3 &= -G_{3}(\phi, X) \Box \phi, \tag{4} \\ \mathcal{L}_4 &= G_{4}(\phi, X) R + G_{4,X} \left[ (\Box \phi)^2 - (\nabla_{\mu} \nabla_{\nu} \phi) (\nabla^{\mu} \nabla^{\nu} \phi) \right], \tag{5} \\ \mathcal{L}_5 &= G_{5}(\phi, X) G_{\mu \nu} (\nabla^{\mu} \nabla^{\nu} \phi) - \frac{1}{6} G_{5,X} \big[ (\Box \phi)^3 - 3 (\Box \phi) (\nabla_{\mu} \nabla_{\nu} \phi) (\nabla^{\mu} \nabla^{\nu} \phi) \nonumber \\ &\quad + 2 (\nabla^{\mu} \nabla_{\alpha} \phi) (\nabla^{\alpha} \nabla_{\beta} \phi) (\nabla^{\beta} \nabla_{\mu} \phi) \big]. \tag{6} \end{align}\] Here, \(G_{i} \, (i=2,3,4,5)\) are arbitrary functions of the scalar field \(\phi\) and their canonical kinetic term \(X\equiv-\frac{1}{2}\nabla^{\mu}\phi \nabla_{\mu} \phi\), with \(\Box \phi=\nabla_{\mu} \nabla^{\mu} \phi\) and partial derivatives \(G_{j,X}(\phi, X)=\partial G_{j}(\phi, X)/\partial X\) with \(j=4,5\), \(R\) is the Ricci scalar and \(G_{\mu \nu}\) is the Einstein tensor. Thus, the complete action of the Horndeski theory is written as \[\label{acaohorndeskigeral} S[g_{\mu \nu}, \phi]=\int d^{4}x\sqrt{-g}(\mathcal{L}_{H}+\mathcal{L}_{M}),\tag{7}\] with \(\mathcal{L}_{M}\) representing the matter and radiation content of the universe, which corresponds to a perfect fluid with power density \(\rho_M\) and pressure \(p_M\).

2.2 Background Equations of Motion↩︎

The next step is to consider an expanding Universe with a spatial flat, homogeneous, and isotropic geometry, described by the Friedmann-Robertson-Walker (FRW) metric system, written in the form \[\label{eqmetfrw} ds^2=-dt+a^2(t)\delta_{ij}dx^{i}dx^{j},\tag{8}\] from which the Friedmann equations assume the form [35], \[\label{eqFRI1} \begin{align} & 2 X G_{2, X}-G_2+6 X \dot{\phi} H G_{3, X}-2 X G_{3, \phi}-6 H^2 G_4+24 H^2 X\Big(G_{4, X}+X G_{4, X X}\Big) -12 H X \dot{\phi} G_{4, \phi X}\\ &-6 H \dot{\phi} G_{4, \phi}+2 H^3 X \dot{\phi}\Big(5 G_{5, X}+2 X G_{5, X X}\Big) -6 H^2 X\Big(3 G_{5, \phi}+2 X G_{5, \phi X}\Big)=-\big(\rho_A+\rho_B\big), \end{align}\tag{9}\] and \[\label{eqFRI2} \begin{align} G_2- & 2 X\Big(G_{3, \phi}+\ddot{\phi} G_{3, X}\Big)+2\Big(3 H^2+2 \dot{H}\Big) G_4 -12 H^2 X G_{4, X}-4 H \dot{X} G_{4, X}-8 \dot{H} X G_{4, X}\\ & -8 H X \dot{X} G_{4, X X}+2\Big(\ddot{\phi}+2 H \dot{\phi}\Big) G_{4, \phi}+4 X G_{4, \phi \phi}+4 X\Big(\ddot{\phi}-2 H \dot{\phi}\Big) G_{4, \phi X} \\ & -2 X\Big(2 H^3 \dot{\phi}+2 H \dot{H} \dot{\phi}+3 H^2 \ddot{\phi}\Big) G_{5, X}-4 H^2 X^2 \ddot{\phi} G_{5, X X}+4 H X\Big(\dot{X}-H X\Big) G_{5, \phi X}\\ & +2\left[2\Big(\dot{H} X+H \dot{X}\Big) +3 H^2 X\right] G_{5, \phi}+4 H X \dot{\phi} G_{5, \phi \phi}=-\big(p_A+p_B\big). \end{align}\tag{10}\] The subscripts \(A\) and \(B\) in densities and pressures represent two perfect fluids, which are generally attributed to amounts of matter and radiation, and both equations 9 and 10 can be written in their usual forms, as follows: \[\label{eqfri1} H^2=\frac{8 \pi G}{3}\big(\rho_M+\rho_R +\rho_{\phi}\big),\tag{11}\] and \[\label{eqfri2} 2\dot{H}+3H^2=-8\pi G\big(p_M+p_R+p_{\phi} \big).\tag{12}\] The modification terms of the theory can all be compactly expressed in the energy density and pressure quantities associated with the scalar field in Eqs. 11 and 12 . We vary the action 7 with respect to \(\phi\), thus obtaining the evolution of the scalar field given in the form [35], \[\label{eqfri3} \frac{1}{a^3}\frac{d}{dt}\big(a^3 J\big)=P_{\phi},\tag{13}\] where \[\begin{align} J= & \dot{\phi} G_{2, X}+6 H X G_{3, X}-2 \dot{\phi} G_{3, \phi}+6 H^2 \dot{\phi}\Big(G_{4, X}+2 X G_{4, X X}\Big)-12 H X G_{4, \phi X} \\ & +2 H^3 X\Big(3 G_{5, X}+2 X G_{5, X X}\Big)-6 H^2 \dot{\phi}\Big(G_{5, \phi}+X G_{5, \phi X}\Big) \\ \end{align}\] and \[\begin{align} P_{\phi}= & G_{2, \phi}-2 X\left(G_{3, \phi \phi}+\ddot{\phi} G_{3, \phi X}\right)+6\left(2 H^2+\dot{H}\right) G_{4, \phi}\\ & +6 H\left(\dot{X}+2 H X\right) G_{4, \phi X}-6 H^2 X G_{5, \phi \phi}+2 H^3 X \dot{\phi} G_{5, \phi X}, \end{align}\] being \(J\) a current and \(P_{\phi}\) the scalar source. In the equations presented above, the two perfect fluids satisfy the following continuity equations for matter and radiation, respectively, \[\dot{\rho}_M+3H\rho_M \big(1+\omega_M \big)=0,\] and \[\dot{\rho}_R+3H\rho_R \big(1+\omega_R \big)=0.\] Through this, with the set of solutions provided by these equations listed so far, it is possible to obtain the full background evolution of the universe.

2.3 Perturbations, Instabilities and Gravitational Waves Constraints↩︎

In the context of perturbations [35], [36], the interest is focused on scalar and tensor perturbations, verifying the conditions of absence of ghost and Laplacian instabilities, which, if satisfied, will guarantee the cosmological viability of the proposed model. In particular, for Horndeski theory to be free of Laplacian instabilities associated with propagation speed of the scalar and tensor (gravitational waves) fields, we must, respectively, have the following expressions \[\label{eqcs2} c_S^2 \equiv \frac{3\big(2 w_1^2 w_2 H-w_2^2 w_4+4 w_1 w_2 \dot{w}_1-2 w_1^2 \dot{w}_2\big)-6w_1^2\big[(1+\omega_A)\rho_A +(1+\omega_B)\rho_B\big]}{w_1\big(4 w_1 w_3+9 w_2^2\big)} \geq 0,\tag{14}\] and \[\label{eqcT} c_T^2 \equiv \frac{ w_4}{w_1} \geq 0,\tag{15}\] while for the absence of ghost instabilities associated with the kinetic energy of the scalar and tensor perturbations, we must have, respectively, \[\begin{align} \label{eqqs} Q_S \equiv \frac{w_1\big(4 w_1 w_3+9 w_2^2\big)}{3 w_2^2}>0 \quad \text{and} \quad Q_T \equiv \frac{w_1}{4}>0 . \end{align}\tag{16}\] In the Horndeski gravity formalism applied to a FRW cosmological background, these coefficients \(w_{i}\,(i=1,2,3\) and \(4)\) are obtained from the perturbed Lagrangian, and their general expression can be found in [35].

The propagation of tensor perturbations, Eq. 15 , is of fundamental interest, as it imposes serious constraints on Horndeski gravity, particularly on the functions \(G_4(X, \phi)\) and \(G_5(X, \phi)\), with respect to the speed of gravitational waves \((c_T)\). According to the GW170817 observation, conducted by the LIGO/Virgo collaboration [37], together with its electromagnetic counterpart GRB 170817A [38]–[40] — in which a gravitational wave signal originating from the merger of binary neutron stars was detected — a strict bound was imposed on \(c_T\) [41]–[43], indicating that the speed of gravitational waves in the late universe must satisfy \[|c_{T}^2-1|\lesssim 10^{-15}\label{ct95cond}.\tag{17}\] This constraint implies that the speed of gravitational waves must be practically identical to that of electromagnetic waves (\(c_T = c\)). This becomes particularly relevant because the arbitrariness in the functions \(G_i\) is drastically reduced, especially in the terms \(G_4(X, \phi)\) and \(G_5(X, \phi)\). In many Horndeski models, \(G_5(X, \phi)\) is simply neglected in the description of the current accelerated expansion of the Universe, as well as the terms proportional to \(G_{4,X}\), \(G_{5,X}\), and \(G_{5,\phi}\) [41]–[48]. The vanishing of these terms becomes necessary to ensure that \(c_T^2 = 1\) in Eq. 15 . However, there are well-established limits where it is still possible to consider \(G_{4,X} \neq 0\) and \(G_{5,\phi} \neq 0\) without violating the constraints imposed on \(c_T\), as done in [49], for a model with non-minimal derivative coupling between the scalar field and the Einstein tensor, showing that a wide range of values for the coupling parameter \(\eta\) is compatible with the stringent bounds on \(c_T\). Indeed, there are also other possibilities, such as in ‘extended Horndeski gravities’, where this constraint can be easily respected. For a recent investigation on this issue, see, e.g., Ref.  [50].

3 The Non-minimal Derivative Coupling Sector↩︎

In this model of Horndeski theory, our main interest lies in the non-minimal coupling between the derivative term of the scalar field and the Einstein tensor \(G_{\mu \nu}\). These models are widely applied in cosmological scenarios, especially within the framework of Horndeski theories, which employ this and other more general couplings between the scalar field and curvature terms [51]–[54]. They are also explored in other extended gravity theories such as, for example, those with additional coupling to the four-dimensional Gauss–Bonnet invariant, giving rise to interesting cosmological phases [55], L.N.?, Granda_2012?. Moreover, this coupling has been investigated in the context of the curvature model [56], in discussions of inflation during the rapid oscillation of a scalar field [57], and finally, within the framework of four-dimensional \(N=1\) minimal supergravity [58].

3.1 The Model↩︎

In our model, we consider the sector of Horndeski gravity that can be obtained by defining the functions \(G_{i}\) as follows \[\label{eqmodel} G_{2}(\phi, X) = \alpha X -2\kappa \Lambda - V(\phi), \qquad G_{3}(\phi, X) = 0, \qquad G_{4}(\phi, X) = \kappa \qquad \text{and} \qquad G_{5}(\phi, X) = -\frac{1}{2} \eta \phi.\tag{18}\] The model examined in this work belongs to a subclass of Horndeski gravity, commonly known as the John sector [59], [60], as framed within the Fab Four (F4) formulation [61], [62]. The corresponding action for this sector is written in the form [22]–[33], \[\label{eqacao} S[g_{\mu \nu}, \phi] = \int d^4x \sqrt{-g} \left[ \kappa( R-2\Lambda) - \frac{1}{2} (\alpha g_{\mu \nu} - \eta G_{\mu \nu}) \nabla^{\mu} \phi \nabla^{\nu} \phi - V(\phi) \right] + S_m[g_{\mu \nu}],\tag{19}\] with \(\kappa = (16 \pi G)^{-1}\). The parameters that control the intensity of the couplings are \(\alpha\) and \(\eta\) the first is the dimensionless parameter and the second having dimensions of \((mass)^{-2}\). It is observed that, by defining \(\alpha = 1\) and \(\eta = 0\) in the action above, we recovered the usual Einstein theory for scalar fields, with minimally coupled gravity to the scalar field \(\phi\) with potential \(V(\phi)\).

3.2 Cosmological Dynamics of the Model↩︎

A possible approach to this topic consists in varying the action in Eq. 19 with respect to the metric and the scalar field, thus obtaining the field equations of the Horndeski theory, which can be solved using the FRW metric presented in Eq. 8 . However, this algebraic procedure can be bypassed by directly applying the model equations 18 to the expressions 9 and 10 . In doing so, we obtain the two Friedmann equations in their standard forms \[H^2 = \frac{8 \pi G}{3} \big(\rho_m + \rho_r + \rho_{\phi} + \rho_{\Lambda}\big),\] and \[2\dot{H} + 3H^2 = -8\pi G \big(p_m + p_r + p_{\phi} + p_{\Lambda}\big).\] The equations corresponding to the energy density and pressure of the scalar field, associated with the model defined in 18 , now take the following forms \[\label{eqdensidphi} \rho_{\phi}=\dfrac{\dot{\phi}^{2}}{2}\left(\alpha+9 \eta H^2 \right)+V(\phi)\tag{20}\] and \[\label{eqpphi} p_{\phi}=\frac{\alpha \dot{\phi}^2}{2}-\frac{\eta\dot{\phi}^2}{2} \big(2\dot{H}+3H^2\big) -2\eta H \dot{\phi} \ddot{\phi}+V(\phi).\tag{21}\] We take the derivative of \(\rho_{\phi}\) with respect to time, thus obtaining its continuity equation for the scalar field, written in the form \[\label{eqdceponto} \dot{\rho_{\phi}}-\phi \ddot{\phi} \big(\alpha + 9\eta H^2\big) - 9 \eta \dot{\phi}^{2} H \dot{H}-\dot{\phi} V_{\phi}(\phi)=0,\tag{22}\] while the continuity equations for the other components of matter and radiation, which in turn satisfy their usual forms, are written in the following way \[\label{eqdm} \dot{\rho}_m+3H \rho_m=0 \qquad \text{and} \qquad \dot{\rho}_r+4H \rho_r=0.\tag{23}\] Another important equation that drives the cosmological evolution is the equation of motion for \(\phi\), which can be obtained directly by applying 18 to 13 , obtaining the following expression \[\label{eqceponto} \Ddot{\phi}+3H \dot{\phi}+\frac{6 \eta \dot{\phi}H \dot{H}}{\alpha +3\eta H^2}+\frac{V_{\phi}(\phi)}{\alpha +3\eta H^2}=0.\tag{24}\] For \(V_{\phi}=0\), we return to the form found in [22]. Finally, we will compute equations for the energy density parameters whose components satisfy \[\Omega_m+\Omega_r+\Omega_{\Lambda}+\Omega_{\phi}=1.\] The cosmic time \(t\), is related to the redshift \(z\) by the relation \(1+z=a(t_0)/a(t)\), with \(t_0\) being the current time such that \(a(t_0)=1\). The differential equations can be written in terms of \(z\) by using the following relation \[\label{eqtransf} \frac{d}{dt}=-(1+z)H(z)\frac{d}{dz}.\tag{25}\] With this, we rewrite the differential equations in terms of \(z\). The Hubble function is expressed in the form \[\label{funhubble} H(z)= \sqrt{H_{0}^{2}\big[\Omega_{m0}(1+z)^{3} +\Omega_{r0}(1+z)^{4} +\Omega_{\Lambda 0}\big]+\dfrac{8\pi G\rho_{\phi}}{3}},\tag{26}\] where \(H_0\) is the Hubble parameter today, while \(\Omega_{i0}\) stands for the present values of the energy density parameters of each component, and \(\Omega_{\phi}\) represents the quantity associated with the relative energy density of the scalar field. In the following sections, we will consider two specific cases for the field \(\phi\).

We rewrite the equation 20 in terms of \(z\) as follows \[\label{eq1461fried} \rho_{\phi}=\dfrac{(1+z)^{2}\phi'^{2}H^{2}\left(\alpha+9 \eta H^2 \right)}{2}+V(\phi).\tag{27}\] In the same way, we have that equations 22 and 24 are given respectively by \[\label{eqdif5} \rho'_{\phi} + \phi (1 + z)(H'\phi' + H\phi'')(\alpha - 9\eta H^2) - 9\eta (1 + z)^2 H^2 H' \phi'^2 - \phi' V(\phi) = 0,\tag{28}\] and \[\label{eqdif6} (1 + z)^2 H(H'\phi' + H\phi'') - 3(1 + z) H^2 \phi' + \frac{6\eta (1 + z)^2 H^3 H'\phi'}{\alpha + 3\eta H^2} + \frac{V_{\phi}(\phi)}{\alpha + 3\eta H^2} = 0.\tag{29}\] Combining equations (28 ) and (29 ), we obtain a first-order differential equation in \(\phi\), given by \[\label{eqGG} (\alpha + 3\eta H^2)\rho'_{\phi} - 3(1 + z) H^2 \phi'^2 \big[ (\alpha + 3\eta H^2)(\alpha + 9\eta H^2) + \eta (1 + z) H H'(\alpha - 9\eta H^2) \big] + 6\eta H^2 \phi' V'(\phi) = 0.\tag{30}\]

4 Numerical and statistical analysis↩︎

In this section, we describe the procedures performed for the numerical and statistical analysis of the cosmological background given by the equations 18 . The equations of motion of the model are solved numerically using the residual method implemented in the Mathematica software. In this way, we obtain the solutions of the differential equations involved. For statistical analysis, we employ the Markov Chain Monte Carlo (MCMC) method, which determines the parameter space for the free parameters in the model. In this way, we obtain some relevant quantities, such as the Hubble parameter \(H(z)\) and the energy density parameters \(\Omega_i\).

We apply Bayesian sampling of the posterior probability distribution of these parameters, performing an MCMC analysis implemented in the Wolfram Mathematica software, obtaining the best fitting constraints for the following free parameters: \(\alpha\), \(\eta\), \(H_0\) and \(\Omega_{i 0}\). This is done by confronting the Hubble function with experimental data obtained through the Cosmic Chronometer (CC), Baryon Acoustic Oscillations (BAO), and the SH\(0\)ES Collaboration methods for measurements of \(H(z)\) at low redshifts, as we specify below. The agreement between the results of \(H(z)\) from numerical integration and the observational data is evaluated with the following chi-square function \[\begin{align} \chi^{2} = \sum_{i=1}^{N}\dfrac{\big[H(x,z_{i}) - H_{obs}(z_{i})\big]^{2}}{\sigma^{2}(z_{i})}, \end{align}\] where \(N\) is the maximum number of observational data and \(x\) in the \(H\) function are the free parameters of the model. In our analyses, we employ an exponential-type scalar potential, expressed in the form \[V(\phi)=V_0 e^{-\lambda\phi}.\] The exponential potential has been long well-justified in the literature. For recent motivations, see, for example, in dilatonic [63], [64], quintessence [65] and deformed Starobinsky [66] models.

We establish initial conditions for the scalar field and potential amplitude \(V_{0}\), and thus numerically solve the set of differential equations given by the equations 27 and 30 . In this paper, we consider two cases for this model: In the first case, the scalar field plays the role of dark matter, while in the second case, we attribute to the scalar field the description of the amount of dark energy.

4.1 About \(H(z)\) data↩︎

Here we discuss the observational datasets of \(H(z)\) used to fit the model. We used the expansion rate value at redshift \(z=0\), one of the latest results obtained for the Hubble constant, estimated at \(H_0 = 73.04 \pm 1.04~(\mathrm{km/s/Mpc})\). This value comes from the Hubble Space Telescope (HST) measurement together with the Supernova Collaboration \(H_0\) for the Equation of State (SH\(0\)ES) [5]. Regarding the baryon acoustic oscillation (BAO) measurements of \(H(z)\) from the Sloan Digital Sky Survey Collaboration (SDSS), we use the dataset specified and cited in [4], where we can find them in references [67]–[70]. For cosmic chronometer (CC) data, we used a list of measurements with 33 results of \(H(z)\), which can be found in Table III of [71], together with their respective references for each point. In the results shown below, we made some combinations between the aforementioned datasets.

4.2 Case I: scalar field as dark matter↩︎

For this case, we have the contribution of the cosmological constant \(\Lambda\) playing the role of dark energy, while the scalar field plays the role of dark matter, and thus the quantity \(\Omega_{m}\) represents only the baryonic matter. In addition, we have the non-minimal coupling parameters of the model \(\alpha\) and \(\eta\). For statistical analysis using the following values for the potential parameters \(V_{0}=10.0 \times 10^{-124}\,M_{\mathrm{Pl}}^{4}\) and \(\lambda = 0.5\), we assume a small initial value for \(\phi'(z=0) = 0.5\times10^{-6}\, M^{2}_{\mathrm{Pl}}\) such that it guaranties condition 17 , taking into account equations 15 , 31 , and 34 . Furthermore, the initial condition \(\phi(z=0)\) can be obtained so that it satisfies 27 for \(z=0\) using the aforementioned values for the potential parameters and \(\phi'(z=0)\). We present in Table 1 the best fit (with \(1\sigma\) constraints) for the free parameters of the model, while in Table 2 we have the mean values (with \(2\sigma\) constraints) considered for this case of the non-minimal derivative coupling model of Horndeski gravity. These results were obtained by MCMC sampling, where we used different \(H(z)\) data at low and intermediate redshifts.

We performed an individual analysis for the parameter \(H(z)\), obtaining the best fit \(H_{0} = 69.9\pm 0.9~(\mathrm{km/s/Mpc})\) within the confidence level \(68\%\) of the CC measurements, which is in tension with R22 at \(2.4\sigma\). This result is midway between the R21 value and the Planck value. For the analysis of the combined measurements of CC\(+\)BAO\(+\)SH\(0\)ES, we obtain the best fit \(H_{0}=74.0\pm 0.6~(\mathrm{km/s/Mpc})\) within the \(68\%\) confidence level, presenting a value very close to R22 with a tension \(0.8\sigma\), which is particularly interesting due to the precision level and proximity to the measured value locally. However, it raises the tension with P20.

Thus, we present in Fig. 1 the graphs of \(H(z)\) and its normalization \(H(z)/(1+z)\) with their respective confidence bands. In Fig. 2, we have the results of the analysis as posterior probability distribution of the model parameters \(\alpha\) and \(\eta\), in addition to the other parameters of the cosmological background, with their respective contour plots referring to the confidence regions \(1\sigma\) and \(2\sigma\) of the MCMC sampling, obtaining results for the CC data set and for the combination CC\(+\)BAO\(+\)SH\(0\)ES.

4.3 Case II: Scalar field as dark energy↩︎

In this case, we consider the scalar field \(\phi\) playing the role of dark energy. Now, we have that \(\Omega_{m}\) represents the contribution of dark matter and baryonic matter, while \(\Omega_{\phi}\) corresponds to the amount of dark energy. The procedure is similar to that carried out for Case I, with one less free parameter. We assume the following initial conditions \(\phi'(z=0) =0.5\times 10^{-6}\, M^{2}_{\mathrm{Pl}}\), \(V_0=10 \times 10^{-124}\, M_{\mathrm{Pl}}^{4}\) and \(\lambda = 0.5\). For this analysis, we have the addition of statistics using the individual BAO data, shown in Tables 3 and 4

We perform the same individual analysis for the parameter \(H(z)\), where we obtain the best fit \(H_{0} = 70.8 \pm 0.8~(\mathrm{km/s/Mpc})\) within the confidence level of \(68\%\), with points of \(H(z)\) from the combined CC\(+\)BAO\(+\)SH\(0\)ES measurements, which is in tension with R21 at \(1.7\sigma\), is particularly interesting since it is much smaller than the tension existing between R22 and P20 of \(5\sigma\), in addition to several other dark energy models, which provides considerable relief in the Hubble tension with the model based on Horndeski gravity. Regarding the average values presented in Table 4, all remain close to this restriction.

Thus, in Fig. 3 we present the graphs of \(H(z)\) and its normalization \(H(z)/(1+z)\) with their respective confidence bands. In Fig. 4, the results of the analysis as posterior probability distribution of the parameters \(\alpha\) and \(\eta\) of the model, in addition to the other cosmological quantities, together with their respective contour plots with confidence regions 1\(\sigma\) and 2\(\sigma\) of the MCMC sampling. For this case, we analyzed the individual data set of BAO and CC, and finally the combination of the data CC\(+\)BAO\(+\)SH\(0\)ES.

Other approaches based on Horndeski gravity aim to treat, for example, the Hubble tension, in addition to cosmological and astrophysical approaches [17], [19], [20]. The model discussed here presents flexibility in adjusting the fundamental cosmological quantities and the parameters \(\alpha\) and \(\eta\) of the model, producing estimates with good agreement with the standard model, allowing to better accommodate the discrepancies between local and global measurements of \(H_0\). Thus, the numerical and statistical results obtained indicate that the non-minimal derivative coupling model of Horndeski gravity has promising potential in the description of cosmological dynamics, in addition to treating the issue of the tension in \(H_0\), reconciling the different observational regimes of cosmic expansion.

5 Stability Conditions for the Model↩︎

We apply the relations of the functions \(G_{i}\) expressed in 18 to obtain the equations of \(w_{i}\) for the model, written in the form

\[\begin{align} w_1 &= \frac{4\kappa-\eta (1+z)^2 H^2 \phi'^2}{2}, \tag{31} \\ w_2 &= 4\kappa H-3\eta (1+z)^2 H^3 \phi'^2, \tag{32} \\ w_3 &= \frac{3}{2}\alpha (1+z)^2H^2\phi'^2-18 \kappa H^2 +27\eta (1+z)^2H^4\phi'^2, \tag{33} \\ w_4 &= \frac{4\kappa +\eta (1+z)^2 H^2 \phi'^2}{2}. \tag{34} \end{align}\]

Next we write the equations of the parameters of scalar and tensor perturbations expressed in 14 , 15 and 16 in terms of the new \(w_i\) of the model. In doing this, we examine the stability of the solutions obtained for the model, investigating the square of the velocities of the scalar and tensor perturbations, in addition to the kinetic energy associated with the scalar \(Q_S\) and tensor \(Q_T\) perturbations. Thus, for consistent dynamics, free of Laplacian and ghost instabilities for the scalar and tensor modes, the following conditions must be satisfied, \[\label{eqcond} c_S^2\geq0, \qquad Q_S>0, \qquad c_T^2 \geq 0\qquad \text{and} \qquad Q_T>0.\tag{35}\] In this way, we analyze the evolutions of 14 , 16 and 15 as a function of redshift for the background solution given by the model of the equations 18 . In Fig. 5, we evolve the propagation velocity (left panel) and kinetic energy (right panel) associated with the scalar perturbations. In the figures in this section, we use the best fit results for the combination CC\(+\)BAO\(+\)SH\(0\)ES data. For \(Q_{s}\), the stability conditions are satisfied throughout the evaluated interval, showing that the solutions obtained for the model are free of ghost instabilities. When we investigate the behavior of \(c_{S}\), the behavior of the curves admits Laplacian stability within the analyzed redshift range \(0<z<5000\) (This upper bound redshift is far beyond typical astronomical observations). The curves remain above zero for both case: the scalar field playing the role of dark energy (dashed black line) and dark matter (dotted gray line). In particular, in both cases we note that at a moment in late time Universe the curve approaches zero, to \(z \approx 2\) for the scalar field as dark matter and \(z \approx 12\) for dark energy, but any case does not violate stability. This is precisely the regime where the Horndeski gravity is completely safe, as shown in the extension of the theory presented in Ref. [72].

One of the recurring discussions within the framework of Horndeski theories concerns the possibility of a time variation in the speed of gravitational waves in models that involve the terms \(G_4(\phi, X)\) and \(G_5(\phi, X)\). In our study, we consider \(G_5(\phi, X) = -\eta \phi/2\), which, in principle, does not guarantee that the speed of gravitational waves is exactly luminal (\(c_T^2 = 1\)). For this reason, contributions from \(G_5(\phi, X)\) are commonly neglected in many approaches based on this class of theories, as previously discussed. However, in the model that we are analyzing, by properly adjusting the free parameters and making appropriate initial conditions choices, it is possible to obtain a velocity \(c_T^2 \approx 1\) at \(z \approx 0\), while remaining positive throughout the entire evolution interval considered. This behavior can be verified in the graph on the left of Fig. 6. Regarding the kinetic energy quantity \(Q_T\), shown on the right of Fig. 6, it also remains positive throughout the entire evolution. These results ensure that the model satisfies the stability conditions and is, therefore, free from ghost and Laplacian instabilities associated with tensor modes. We can also study an important cosmological quantity, the parameter associated with the effective equation of state for the scalar field given by \(\omega_{\phi} = p_{\phi}/\rho_{\phi}\), which can be obtained by using the equations 20 and 21 .

We can see in Fig. 7 the behavior of the equation of state for both cases using the best-fit results for the CC\(+\)BAO\(+\)SH\(0\)ES data. For the scalar field as dark energy (dashed black line), we have a slight peak where \(\omega_{\phi}>1\) (‘super stiff matter’ regime), reflecting the fact that \(\omega_\phi\) is indeed an effective equation of state.

6 Conclusions↩︎

In this work, we explore the effects of Horndeski gravity integrated with the theoretical and observational foundations of the standard cosmological model. We analyze the statistics, using the observational data of \(H(z)\) from the SH\(0\)ES, BAO, and CC datasets, applied to a specific model featuring a non-minimal derivative coupling between the scalar field and Einstein tensor. Two scenarios were investigated: in the first (Case I), the scalar field acts as dark matter; in the second (Case II), the scalar field is responsible for the dark energy component of the Universe, replacing the cosmological constant in driving cosmic expansion.

This approach allowed us to obtain estimates for the values of the main cosmological parameter, as can be seen in the posterior distributions and in the tables. In this way, limits were set for the coupling parameters \(\alpha\) and \(\eta\) of the model, yielding values that fit well within the current cosmological framework. Particular attention was paid to the \(H_0\) values obtained in both cases. We found that using different \(H(z)\) datasets, as well as their combinations, resulted in Hubble constant values within the ranges reported by P20 and R22. This result is particularly relevant, as it contributes to alleviating the tension between local and global determinations of \(H_0\). Additionally, we observed that the behavior of \(H(z)\) remains within the expected range when evolved over \(0 \leq z \leq 2.5\).

In Case I the best mitigation of the tension was achieved using combined CC\(+\)BAO\(+\)SH\(0\)ES data, reducing the tension with R22 to only \(0.8\sigma\), while in Case II, the reduction reached \(1.7\sigma\) using the combined data, which provides relief of the tension with R22. Overall, the model studied, which has been extensively investigated in recent literature in various contexts of gravity and cosmology proved to be both promising and consistent with the results of our analysis.

Finally, concerning the effective equation of state, in neither the aforementioned case, was the scalar field able to cross the phantom divide at large redshifts as recently pointed and explored by the DESI collaboration [73], [74]. This is a point to be explored in our setup in the realm of the DESI dataset in upcoming investigations, along the lines of [65], [75], [76] in the context of “beyond Horndeski” physics.

We thank CNPq and CAPES for partial financial support. FAB acknowledges support from CNPq (Grant No. \(309092/2022-1\)). JAVC thanks the Paraíba State Research Support Foundation (FAPESQ) (Grant No. \(22/2025\)) for financial support. The authors thank Amilcar R. Queiroz for his helpful comments.

In memory of our late collaborator and friend Prof. Raimundo Silva.

References↩︎