\section{Introduction}

One of the most significant discoveries in modern cosmology is the observational evidence of the accelerated expansion of the present Universe, which has been independently inferred from the observations of distant type Ia supernovae \cite{SupernovaSearchTeam:1998fmf,SupernovaCosmologyProject:1998vns}. Such an accelerated expansion is thought to be caused by an unknown yet dominant component of the Universe called dark energy. Within the regime of general relativity (GR), the Standard Model of Cosmology, also known as the $\Lambda$-Cold Dark Matter ($\Lambda$CDM) model, can explain this phenomenon with the help of the cosmological constant $\Lambda$. Despite being the simplest model describing the Universe, the $\Lambda$CDM framework is plagued by fundamental issues like the cosmological constant problem and lack of physical understanding regarding late-time acceleration of the Universe, which might be arising due to some unknown dynamics at large scales. There are mainly two approaches that have been followed in order to address these issues -- either by considering GR to be the accepted theory of gravity and modifying the dark energy components \cite{Copeland:2006wr} or by modifying the gravitational theory itself in a way that ensures the accelerated expansion of the Universe, while reducing to GR in certain limits \cite{Clifton:2011jh}.

In the former case, one disregards the cosmological constant $\Lambda$ and instead introduces a scalar field that mimics the dynamical dark energy, thereby leading to the late-time accelerated expansion.  On the other hand, a modification of the underlying gravitational theory is usually done by introducing some extra degrees of freedom (DoF) in addition to Einstein's gravity or GR. One of the ways to manifest this is via the scalar-tensor theories of gravity \cite{Fujii:2003pa} wherein the extra DoF are mixed with curvature, and on subjecting the action to a conformal transformation, they get non-minimally coupled with the matter. The idea is that one should be able to make use of such modified gravity theories to generate extra DoF that mimic the late time acceleration at large scales and hence resolve the problems that the $\Lambda$CDM model faces. However, with respect to the local measurements, these DoF should be suppressed so that Einstein's gravity, which agrees very well with observations, becomes the underlying gravitational theory. The required local suppression of the degrees of freedom is done via various screening mechanisms such as the chameleon \cite{Khoury:2003rn} and the Vainshtein mechanisms \cite{Vainshtein:1972sx}.

Over the years, various cosmological tensions have also been observed in different datasets, one of the most interesting being the discrepancy in the measurement of the Hubble constant $H_0$ at early and late times.  The Hubble tension thus raises a direct question about our understanding of the expansion of the Universe since the prediction of the $\Lambda$CDM model inferred from the Planck
CMB data leads to $H_0=(67.4 \pm 0.6)$ km/s/Mpc \cite{Planck:2018vyg}, while  its determination using type Ia supernovae by the SH0ES collaboration  \cite{Riess:2021jrx} implies $H_0=(73.3\pm 1.04)$ km/s/Mpc which indicates a tension of $\sim 5\sigma$ with respect to $\Lambda$CDM. The physical implication of this would be a faster expansion of the Universe than the prediction from $\Lambda$CDM. Besides, a lot of progress has been made in measuring $H_0$ with alternative methods, but all of them seem to indicate some tension with the $\Lambda$CDM inferred value. 
Apart from the $H_0$ tension,  the $\sigma_8$ (the cosmological parameter related to the clustering of matter in an $8h^{-1}$Mpc radius sphere) tension  refers to a discrepancy between the estimation from the CMB and the measurement from surveys like BOSS \cite{BOSS:2016wmc, Philcox:2021kcw}.

Recently, a large number of attempts have been made to alleviate the Hubble tension using different methods. A summary of these approaches can be found in \cite{Verde:2019ivm, Freedman:2021ahq, Knox:2019rjx, DiValentino:2021izs,Perivolaropoulos:2021jda,Zumalacarregui:2020cjh,Abdalla:2022yfr,Shah:2021onj,Schoneberg:2021qvd}. Some of these approaches include combining the Baryon Acoustic Oscillations (BAO) and the Big Bang Nucleosynthesis (BBN) data \cite{Schoneberg:2019wmt, Schoneberg:2022ggi}, imposing constraints on early dark energy \cite{Karwal:2016vyq, Agrawal:2019lmo, Smith:2019ihp, Lin:2019qug, Poulin:2018cxd, Niedermann:2019olb, Niedermann:2020dwg, Herold:2022iib}, hot new early dark energy \cite{Niedermann:2021ijp, Niedermann:2021vgd}, creating an early dark energy scenario using thermal friction \cite{Berghaus:2019cls, Berghaus:2022cwf}, using new gravitational scalar-tensor theories that are free of ghosts by definition \cite{Banerjee:2022ynv}, using massive neutrinos \cite{Poulin:2018zxs, Kreisch:2019yzn, Sakstein:2019fmf},  employing ``quintessence" potentials to solve the $H_0$ and $\sigma_8$ tensions simultaneously \cite{Adil:2022hkj}, considering a quintessence field that tracks the equation of state during recombination and then boosts the late-time expansion \cite{DiValentino:2019exe}, tuning the late-time acceleration of the Universe to cause a monotonous increase in the expansion rate \cite{Adil:2021zxp}, employing a phenomenological two-parameter family of dark energy models \cite{Lee:2022cyh},  interacting dark energy models \cite{DiValentino:2017iww,Gariazzo:2021qtg,Yao:2023jau} and many more \cite{Alexander:2019rsc, Krishnan:2020obg, Adi:2020qqf, Braglia:2020iik}.

In this work, we focus our attention on the alleviation of the $H_0$ tension using an approach that involves the modification of the background evolution at late times with the help of a dynamical scalar field in the framework of Horndeski gravity \cite{Horndeski:1974wa}, which is one of the most general scalar-tensor theories involving second-order field equations. Within the Horndeski framework, we construct two different models: Model I with a non-trivial $G_5$ term, and Model II with both the $G_4$ and $G_5$ terms being non-trivial,  which enable the dynamics allowing for the late-time accelerated expansion of the Universe and also simultaneously help alleviate the $H_0$ tension. Such models resort to late dark energy \cite{Cai:2021weh,Gangopadhyay:2022bsh,Efstathiou:2021ocp}, contrary to the usual early dark energy approach \cite{Poulin:2018cxd,Herold:2022iib,Niedermann:2021ijp,Niedermann:2021vgd}. For these models to be consistent and well-behaved, we also examine various issues, such as gradient and ghost instabilities and superluminal propagation of the scalar and tensor perturbations. 
 
For each model, we consider two different classes of scalar field potentials, namely freezing and thawing potentials \cite{Caldwell:2005tm,Amendola:2015ksp}. In the first case, even before the acceleration sets in, the scalar field has already been rolling down towards the minimum of the potential and comes to a stop as it starts dominating the Universe, thus exhibiting a tracking behaviour. On the other hand, thawing potentials are the ones in which the field has been slightly displaced from the minimum and got frozen and starts rolling down to the minimum only at late times.
We find that Model I successfully alleviates the Hubble tension but suffers from the superluminal propagation of tensor modes and running of the effective Planck mass.
On the other hand, in addition to  alleviating the Hubble tension, Model II also remains free from superluminal propagation as well as the running of the effective Planck mass. We thus establish that Model II is a better  scenario for resolving the Hubble tension in this framework.  

Our paper is organised as follows: in Section \ref{sec:horndeski_intro}, we briefly introduce the Horndeski theory and its dynamics and state the conditions for the avoidance of any gradient and ghost instabilities. In Section \ref{sec:model_I}, we describe and analyze the results for our first model, pertaining to a generalized non-minimal derivative coupling (GNMDC) involving only the $G_5$ term of the Horndeski Lagrangian. In Section \ref{sec:model_II}, we first argue the reasons for going beyond Model I and introduce Model II, which includes, in addition to the $G_5$ term, the $G_4$ term as well and show that this model is indeed far better than Model I. Finally, in Section \ref{sec:conclusion}, we summarize our results and conclude with some discussions and future directions.

We work with natural units such that $\hbar=c=1$, and set the reduced Planck mass $M_{Pl}=(8\pi G)^{-1/2}=1$.
Our metric signature is ~$(-,+,+,+)$.

\section{Horndeski theory and its dynamics}
\label{sec:horndeski_intro}

We work in the framework of the Horndeski theory which is a generalised scalar-tensor theory in four dimensions whose Lagrangian is constructed out of the metric tensor and a scalar field and leads to second order equations of motion, thereby free of Ostrogradsky instability. The complete Lagrangian for the Horndeski theory (or equivalently, for the generalized Galileons) is given by 
\begin{equation}
\mathcal{L}= \sum_{i=2}^5 \mathcal{L}_{i}\,,
\label{e: horndeski_lag}
\end{equation}

where
\begin{eqnarray}
     && \mathcal{L}_2\ =\  G_2(\phi, X),\nonumber\\
     && \mathcal{L}_3\ =\  -G_3(\phi, X)\Box{\phi},\nonumber\\
     && \mathcal{L}_4\ =\  G_4(\phi, X)R + G_{4,X}(\phi, X)\Big[(\Box{\phi})^2 - (\nabla_{\mu}\nabla_{\nu}\phi)^2\Big],\nonumber\\
     && \mathcal{L}_5\ =\  G_5(\phi,X)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi - \frac{1}{6}  G_{5,X}(\phi, X)\Big[(\Box{\phi})^3-3\Box{\phi}(\nabla_{\mu}\nabla_{\nu}\phi)^2+2(\nabla_{\mu}\nabla_{\nu}\phi)^3\Big], \nonumber
\end{eqnarray}
where $R$ is the Ricci scalar, $G_i$ are four independent arbitrary functions of $\phi$ and $X$, and $G_{i,Y} =\partial{G}_i/\partial{Y}$ with $Y=\{\phi,X\}$.
Thus, in the Horndeski gravity, the complete action can be given as,
\begin{equation}
\mathcal{S} = \int \mathrm{d}^4x \sqrt{-g}\, \left(\sum_{i=2}^5 \mathcal{L}_{i}+\mathcal{L}_M\right)\,,
\label{e:horndeski_action}
\end{equation}
where,
$\mathcal{L}_M$ accounts for the matter (and the radiation) component. Such a theory can well explain the dark energy phase of Universe modelled by a scalar field $\phi$ coupled to gravity (minimal or non-minimal coupling) depending on the choice of $G_i(\phi,X)$. In a flat FRW background with $ds^2=-dt^2+a(t)^2 d\Bar{x}^2$, the equations of motion are given by \cite{DeFelice:2011bh},
\begin{eqnarray}
    2X G_{2,X}-G_2+6X\dot\phi H G_{3,X}-2 X G_{3,\phi}-6H^2 G_4+24H^2 X(G_{4,X}+X G_{4,XX}-12H X \dot\phi G_{4,\phi X} \nonumber\\
    -6H \dot\phi G_{4,\phi}+2H^3 X \dot\phi(5G_{5,X}+2X G_{5,XX})-6H^2 X(3G_{5,\phi}+2X G_{5,\phi X})=-\rho_M
\label{e:friedmann_first}
\end{eqnarray}
    
\begin{eqnarray}
    & G_{2}-2X(G_{3,\phi}+\ddot\phi G_{3,X})+2(3H^2+2\dot H)G_4-12H^2 X G_{4,X}-4H \dot X G_{4,X}-8\dot H X G_{4,X}\nonumber\\& -8H X\dot X G_{4,XX} +2(\ddot \phi+2H \dot \phi)G_{4,\phi}+4X G_{4,\phi\phi}+4X(\ddot \phi-2H\dot \phi)G_{4,\phi X}-2X(2H^3 \dot \phi+ 2H \dot H \dot \phi +3H^2 \ddot\phi)G_{5,X}\nonumber\\ & -4H^2 X^2 \ddot\phi G_{5,XX}+ 4H X(\dot X-H X)G_{5,\phi X}+2[2(\dot H X+H \dot X)+3H^2 X]G_{5,\phi}+4H X \dot\phi G_{5,\phi \phi}=-p_M
\label{e:friedmann_second}
\end{eqnarray}
    
\begin{equation}
    \frac{1}{a^3}\frac{d}{dt}(a^3
    \mathcal{J})=\mathcal{P}_\phi
\label{e:phi_evolution}
\end{equation}
where,
\begin{eqnarray}
    \mathcal{J}=\dot\phi G_{2,X}+6H X G_{3,X}-2 \dot\phi G_{3,\phi}+6H^2 \dot\phi(G_{4,X}+2X G_{4,XX})-12H X G_{4,\phi X} \nonumber\\
    +2H^3 X(3G_{5,X}+2X G_{5,XX})-6H^2 \dot\phi(G_{5,\phi}+X G_{5,\phi X}),
\label{e:J}
\end{eqnarray}

\begin{eqnarray}
    \mathcal{P}_\phi=G_{2,\phi}-2X(G_{3,\phi\phi}+\ddot\phi G_{3,\phi X})+6(2H^2+\dot H)G_{4,\phi}+6H(\dot X+2H X)G_{4,\phi X}\nonumber\\
    -6H^2 X G_{5,\phi \phi}+2H^3 X \dot \phi G_{5,\phi X}.
\label{e:P_phi}
\end{eqnarray}
Here, Eqs. \ref{e:friedmann_first} and \ref{e:friedmann_second} are the two Friedmann equations while Eq. \ref{e:phi_evolution} gives the evolution of the dark energy field $\phi$. Also, $\rho_M$ and $p_M$ are the energy density and pressure of the matter field (and radiation), which satisfy the continuity equation,
\begin{equation}
    \dot \rho_M+3 H \rho_M (1+w_M)=0,
\label{e:continuity}
\end{equation}
where $w_M=p_M/\rho_M$ is the equation of state for the fluid. The complete background evolution of the Universe can be obtained by solving the above set of equations. However, due to the non-trivial and non-canonical structure of the Horndeski Lagrangian, sometimes these models can have instability issues in the evolution of the perturbations, rendering the background evolution inappropriate. It is therefore necessary to keep a check on the parameters of the theory, even if one is not concerned with the evolution of perturbations for a given dark energy model. The Laplacian or gradient instability is related to the propagation speed of the scalar and tensor perturbations which arise in the regime where the square of sound speed of perturbations becomes negative. This leads to an unstable growth of the perturbation modes on small scales. Another crucial one is the ghost instability which arises when the sign of the kinetic term of the perturbations goes negative. We shall work in the parameter space wherein we avoid all these instabilities so as to have a well-behaved theory. Following the standard technique of linear cosmological perturbation theory, one can obtain the second-order action for the scalar and tensor perturbations in Horndeski theory as \cite{Bellini:2014fua},
\begin{equation}
    \mathcal{S}_2=\int \mathrm{d}t \mathrm{d}^3 x a^3 \left[\mathcal{Q}_s\left(\mathcal{\dot R}^2-\frac{c_s^2}{a^2}(\partial_i \mathcal{R})^2\right)+\mathcal{Q}_T\left(\dot h_{ij}^2-\frac{c_T^2}{a^2}(\partial_k h_{ij})^2\right)\right],
\label{e:action_scalartensor}
\end{equation}
where $\mathcal{R}$ is the scalar curvature perturbation while $h_{ij}$ are the tensor modes or gravitational waves (GWs). Similarly $c_s$ and $c_T$ are the propagation speeds of the scalar and tensor modes, respectively and are given as,
\begin{eqnarray}
    c_s^2=-\frac{(2-\alpha_B)\big[\dot H-\frac{1}{2}H^2 \alpha_B(1+\alpha_T)-H^2(\alpha_M-\alpha_T)\big]-H\dot\alpha_B+\rho_M+p_M}{H^2 D},\hspace{.5cm}c_T^2=1+\alpha_T
\label{e:cs_cT}
\end{eqnarray}
where,
\begin{equation}
\mathcal{Q}_s=\frac{2 M_{*}^2 D}{(2-\alpha_B)^2},\hspace{.5cm}
D=\alpha_K+\frac{3}{2}\alpha_B^2, \hspace{.5 cm}
\mathcal{Q}_T=\frac{M_{*}^2}{8}
\label{e:Qs_D_QT}
\end{equation}
and $M_{*}$ is the effective Planck mass and the Horndeski functions $\alpha_B,\alpha_K,\alpha_M,\alpha_T$ are usually   defined as follows \cite{Bellini:2014fua}

\begin{eqnarray}
    M_{*}^2=2\big(G_4-2X G_{4,X}+X G_{5,\phi}-\dot\phi H X G_{5,X}\big)
\label{e:effective planck mass}    
\end{eqnarray}

\begin{eqnarray}
    H M_{*}^2 \alpha_M=\frac{d}{dt}M_{*}^2
\label{e:alphaM}    
\end{eqnarray}

\begin{eqnarray}
    H^2 M_{*}^2 \alpha_K=2X\big(G_{2,X}+2X G_{2,XX}-2G_{3,\phi}-2X G_{3,\phi X}\big)+12\dot\phi X H\big(G_{3,X}+X G_{3,XX}-3G_{4,\phi X}\nonumber\\ -2X G_{4,\phi XX}\big)
    +12X H^2\big(G_{4,X}+8X G_{4,XX}+4X^2 G_{4,XXX}\big)-12X H^2 \big(G_{5,\phi}\nonumber\\+5X G_{5,\phi X}+2X^2 G_{5,\phi X X}\big)
    +4\dot \phi X H^3\big(3G_{5,X}+7X G_{5,XX}+2X^2 G_{5,XXX}\big)
\label{e:alphaK}    
\end{eqnarray}

\begin{eqnarray}
    H M_{*}^2 \alpha_B=2\dot\phi\big(X G_{3,X}-G_{4,\phi}-2X G_{4,\phi X}\big)+8H X\big(G_{4,X}+2X G_{4,XX}-G_{5,\phi}-X G_{5,\phi X}\big)\nonumber\\+2\dot \phi X H^2\big(3G_{5,X}+2X G_{5,XX}\big)
\label{e: alphaB}    
\end{eqnarray}

\begin{eqnarray}
    M_{*}^2 \alpha_T=2X\big(2G_{4,X}-2G_{5,\phi}-(\ddot\phi-\dot\phi H)G_{5,X}\big)
\label{e: alphaT}    
\end{eqnarray}
In the above relations, $\alpha_M$ denotes the running of the effective Planck mass $M_*$, $\alpha_K$ is the kineticity contributing to the kinetic energy of scalar perturbations, $\alpha_B$ gives the braiding or the amount of kinetic mixing between the scalar and tensor perturbations, and $\alpha_T$ is related to the speed of the tensor modes as given in Eq. \ref{e:cs_cT} \cite{Noller:2018wyv,Matos:2022uew}.
For a consistent dynamics, free of any gradient and ghost instabilities for both scalar and tensor modes, one must satisfy

\begin{eqnarray}
    c_s^2>0, \hspace{.5cm} \mathcal{Q}_s>0, \hspace{.5cm} c_T^2>0, \hspace{.5cm} \mathcal{Q}_T>0
\end{eqnarray}

Another non-trivial characteristic of Horndeski theories is that the speed of GWs can evolve with time for models containing $G_4(\phi,X)$ and $G_5(\phi,X)$, as can be seen from Eqs. \ref{e:cs_cT} and \ref{e: alphaT}. The recent observations of GWs from LIGO-VIRGO collaboration and their electromagnetic counterparts put a very stringent bound on speed of GWs, such that
\cite{Gong:2017kim, LIGOScientific:2017zic}
\begin{align}
-3\times 10^{-15}<c_T-1<7\times10^{-16} \,\label{e:cTbound}
\end{align}
which implies that a consistent dark energy model should not violate Eq. \ref{e:cTbound} in any regime. For Horndeski theories, this bound can be stated in terms of $\alpha_T$ as \cite{Kreisch:2017uet},
\begin{align}
-6\times 10^{-15}<\alpha_{T}<1.4\times10^{-15} \,\label{e:alphaTbound}
\end{align}
thereby putting strong constraints on the same.


\section{Model I}
\label{sec:model_I}
In order to address the $H_0$ tension, we first construct a model within the framework of Horndeski Lagrangian, which involves a non-zero $G_5(\phi,X)$. This acts as an extension to 
the quintessence dark energy scenario \cite{Tsujikawa:2013fta}, such that
\begin{equation}
    G_2=X-V(\phi);\quad G_3=0;\quad G_4=\frac{1}{2};\quad 
    G_5=G_5(\phi)
\label{e:G5_setup}
\end{equation}
where $\phi$ is the scalar field which drives the dark energy phase, $X=\dot\phi^2/2$ is the kinetic term, $V=V(\phi)$ is the potential of the field and we choose $G_5$ to be a function of $\phi$ only.
Clearly, $G_5=0$ corresponds to quintessence, while $\Lambda$CDM can be obtained by choosing $G_2=-2\Lambda$, $G_4=1/2$ and $G_3=G_5=0$. The choice of a non-zero $G_5$, as in Eq. \ref{e:G5_setup} is often referred to as a generalized non-minimal derivative coupling (GNMDC), which involves a coupling between the Einstein tensor and the derivative of the scalar field. Usually, such an interaction acts as an additional friction to the scalar field evolution and, depending on the choice of the coupling function, can affect the dynamics of the scalar field in a controlled way leading to distinguishable signatures. In the context of early Universe physics, such type of couplings are well explored, leading to interesting implications like the generation of non-trivial features in primordial power spectrum to explain large-scale CMB anomalies \cite{Tiwari:2022zzz} and formation of primordial black holes \cite{Fu:2019ttf,Dalianis:2019vit}. Also, due to the enhanced friction effect, such couplings make certain inflationary models with steep potentials like Higgs inflation \cite{Germani:2010gm}, compatible with observational data, which are otherwise ruled out by Planck \cite{Planck:2018jri}. Apart from the early Universe, GNMDC theories are also explored in the context of quintessence models of dark energy as well as phantom cosmology \cite{Saridakis:2010mf,DiValentino:2019jae,Quiros:2017gsu}. The freedom to choose $G_5$, makes this class of models phenomenologically attractive. Recently in the context of the late Universe, such an interaction has been explored to resolve the $H_0$ tension \cite{Petronikolou:2021shp}, by choosing $G_5$ as a function of the kinetic energy of the scalar field. Here, as our first model, we choose $G_5$ to be a function of only the dark energy scalar field, such that the strength of this interaction comes to play a role only in the very late Universe while leaving negligible effects at high redshifts. Using the GNMDC setup outlined in Eq. \ref{e:G5_setup}, the Eqs. \ref{e:friedmann_first}, \ref{e:friedmann_second} and \ref{e:phi_evolution} for Model I become
\begin{equation}
    3H^2=\frac{1}{2} \dot\phi^2 \left(9 \theta(\phi)  H^2+1\right)+\text{$\rho _M$}+V(\phi)
\label{e:friedmann_first_modelI}
\end{equation}
\begin{equation}
   -\dot\phi^2 \left(3 \theta(\phi)  H^2-\theta(\phi) \dot H+1\right)+2 \theta(\phi)  H
  \ddot\phi\dot\phi+H \theta^\prime(\phi) \dot\phi^3-2
   \dot H-\rho_M-w_M\rho_M=0
   \label{e:friedmann_second_modelI}
\end{equation}
\begin{equation}
    \ddot\phi \left(3 \theta(\phi)  H^2+1\right)+3 H \dot\phi \left(\theta(\phi) 
   \left(3 H^2+2 \dot H\right)+1\right)+\frac{3}{2} H^2 \theta^\prime(\phi)
   \dot\phi^2+V^\prime(\phi)=0
  \label{e:kg_modelI}
\end{equation}
 where $\theta(\phi)=-2G_{5,\phi}$ and $V'(\phi)=dV/d\phi$. We refer to $\theta(\phi)$ as the GNMDC function and we choose it to be a monotonically increasing function of dark energy scalar field $\phi$. To be specific, we choose $\theta(\phi)=c\phi^p$, $c$ being a quantity that numerically establishes the magnitude of the coupling, and $p$ being any positive number that gives a polynomial form to $\theta(\phi)$.
In this work, we fix $p$, the choice of which will be discussed later, and use $c$ as a tunable parameter. We will show further that by increasing $c$, and thereby increasing the strength of the coupling in a controlled manner, one can have strong deviations from the $\Lambda$CDM value of $H_0$, while being indistinguishable from the former at high redshifts, thereby helping to alleviate the $H_0$ tension. As can be obtained from Eq. \ref{e:friedmann_first}, the energy density of the dark energy field for this model can be given as,
\begin{equation}
    \rho_\phi=V(\phi)+\frac{1}{2} \dot\phi^2 +\frac{9}{2} \dot\phi^2 \theta(\phi)  H^2
    \label{e:rho_phi}
\end{equation}

\begin{equation}
    p_\phi=-\frac{1}{2} \dot\phi \left(\dot\phi \left(3 \theta(\phi)
    H^2+2 \theta(\phi) \dot H-1\right)+4 \theta(\phi)  H \ddot\phi+2 H \theta^\prime(\phi) \dot\phi^2\right)-V(\phi).
\label{e:p_phi}
\end{equation}
With these expressions, we can now define the equation of state parameter $w_\phi$ for the dark energy scalar field as 
\begin{equation}
w_{\rm DE}=w_\phi = \frac{p_\phi}{\rho_\phi}
\end{equation}
Evidently, for $\Lambda$CDM cosmology, $w_{\rm DE} =-1$ 
while for the quintessence scalar field models, $w_{\phi}$ does strongly depend on the choice of the potential. It is also evident that the 
conservation equation (Eq. \ref{e:phi_evolution}) will simply reduce to the following equation
\begin{equation}
    \dot \rho_{\phi} + 3H (1+w_{\phi}) \rho_{\phi} =0.
\end{equation}
In order to get a significant effect from the GNMDC term on the scalar field dynamics, we notice from Eqs. \ref{e:friedmann_first_modelI} and \ref{e:rho_phi}, that the function $\theta(\phi)$ should be chosen such that 
\begin{equation}
    \theta(\phi)H^2\sim \mathcal{O}(1)
    \label{e:x4_condition}
\end{equation} 
at late times. This condition ensures that the contribution to the energy density arising from the GNMDC term becomes significant in comparison to the canonical kinetic term which eventually leads to a larger value of the Hubble constant today, thereby alleviating the tension with $\Lambda$CDM.

\subsection{Methodology and analysis}
\label{sec:model_I_method}


In order to proceed with solving these differential equations for the background dynamics, we introduce a set of dynamical variables following the approach of \cite{Amendola:2015ksp}. Similar approaches have also been discussed in \cite{Bhattacharya:2015wlz,Hrycyna:2015vvs,Matsumoto:2017gnx,Humieja:2019ywy} in different contexts.
\begin{eqnarray}
    x_1=\frac{\dot\phi}{\sqrt{6}H};\quad x_2=\frac{\sqrt{V(\phi)}}{\sqrt{3}H}; \quad x_3=\frac{\sqrt{\rho_r}}{\sqrt{3}H};\quad \lambda=\frac{-V^\prime(\phi)}{V(\phi)};\nonumber\\
    x_4=\theta(\phi) H^2;\quad x_5=\theta^\prime(\phi) H^2;\quad \alpha_I=\frac{\theta^{\prime\prime}(\phi)}{\theta^\prime(\phi)}
\label{e:dynamical_variables_G5}
\end{eqnarray}
where the subscript $I$ corresponds to Model I.
For the purpose of computational convenience, we define a function $N$ of the redshift $z$ such that $N=-\ln(1+z)$. Upon  differentiating the above dynamical variables with respect to $N$, we get the following first-order differential equations, which are much easier to solve than the second-order differential equations that we encountered prior to introducing these dynamical variables.

\begin{eqnarray}
    \frac{\text{d}x_1}{\text{d}N}&=&\frac{\ddot\phi}{\sqrt{6}H^2}-x_1\frac{\dot H}{H^2}
    \label{e:eq_x1_dynam_var_G5}\\
    \frac{\text{d}x_2}{\text{d}N}&=& -\sqrt{\frac{3}{2}}\lambda x_1 x_2 - x_2\frac{\dot H}{H^2}\\
    \frac{\text{d}x_3}{\text{d}N}&=& -2 x_3-x_3 \frac{\dot H}{H^2}\\
    \frac{\text{d}x_4}{\text{d}N}&=& \sqrt{6}x_1 x_5+2 x_4 \frac{\dot H}{H^2}\\
    \frac{\text{d}x_5}{\text{d}N}&=& \sqrt{6} \alpha_I x_5 x_1 + 2 x_5 \frac{\dot H}{H^2}\\
     \frac{\text{d}\lambda}{\text{d}N}&=& \sqrt{6}x_1 \lambda^2(1-\Gamma),\quad \text{where}\quad \Gamma=\frac{V^{\prime\prime}(\phi) V(\phi)}{V^\prime(\phi)^2}\\
    \frac{\text{d}\alpha_I}{\text{d}N}&=& \sqrt{6} x_1\alpha_I^2(\beta_I-1),\quad \text{where}\quad \beta_I=\frac{\theta^{\prime\prime\prime}(\phi)\theta^\prime(\phi)}{(\theta^{\prime\prime}(\phi))^2}
\label{e:eqns_dynam_var_G5}
\end{eqnarray}
and where,
\begin{equation}
    \frac{\dot H}{H^2}=\frac{3 \left(w_M (3 x_4+1) \left(x_1^2 (9
   x_4+1)+x_2^2-1\right)+2 \sqrt{6} x_1^3 x_5-x_1^2 (3
   x_4+1) (9 x_4+1)+2 \sqrt{6} \lambda x_1 x_2^2
   x_4+\left(x_2^2-1\right) (3 x_4+1)\right)}{6 \left(x_1^2
   x_4 (9 x_4-1)+x_4\right)+2}
\label{e:HdotbyH2_G5}
\end{equation}

\begin{equation}
    \frac{\ddot \phi}{H^2}=-\frac{9 \sqrt{6} x_1^3 x_4 (9 w_M
   x_4+w_M-2)+3 \sqrt{6} x_1 \left(3 x_4 \left(w_M
   \left(x_2^2-1\right)+x_2^2\right)+1\right)+81
   x_1^4 x_4 x_5+9 x_5^2 \left(\lambda
    x_2^2 x_4+x_5\right)-3 \lambda 
   x_2^2}{3 \left(x_1^2 x_4(9
   x_4-1)+x_4\right)+1}
\label{e:phiddotbyH2_G5}
\end{equation}
In order to avoid encountering any undefined expressions, as in for $\beta_I$ in Eq. \ref{e:eqns_dynam_var_G5}, we restrict ourselves to $p\geq 2$, and choose $p=4$ although the model works fine for any $p\geq 2$.
We consider two models of dark energy potential belonging to two classes of quintessence, namely thawing and freezing, and analyse how they fare when alleviating the $H_0$ tension under the domain of Model I. We work with potentials where $\phi$ increases with time while rolling towards its minima, so that the strength of the coupling also increases with time, as the coupling function is a monotonic increasing function of $\phi$.
Hence we consider two potentials for the dark energy scalar field
\begin{equation}
    V_1(\phi) = \frac{V_0}{\phi}
    \label{e:V1}
\end{equation}
which is a  simple example of a tracker potential for the dark energy scalar field \cite{Zlatev:1998tr}, and belongs to the class of freezing models of quintessence. 
This potential does not have a minimum and therefore, the field keeps rolling down the potential towards infinity.
The other potential that we choose is
\begin{equation}
    V_2(\phi) = V_0(1-\phi)
    \label{e:V2}
\end{equation}
which belongs to the class of thawing models \cite{Kallosh:2003bq}. 
For this case, we restrict our analysis to the field values such that the potential always remains positive.

 In order to understand the dynamics of the system, we solve  Eqs. \ref{e:eq_x1_dynam_var_G5}$-$\ref{e:eqns_dynam_var_G5}, using appropriate initial conditions in the early radiation domination epoch. The kinetic energy of the scalar field is very small at early times, which sets the initial value of $x_1\simeq 0$. Since the potential of the dark energy field governs the dynamics of expansion in the late Universe, the initial condition for $x_2$ is fixed accordingly. The contribution of radiation to the total energy density is quantified by $x_3$, hence  $x_3\simeq 0.99$ (or $\rho_m, \rho_\phi \ll \rho_r$) in the very early Universe. Using the arguments pertaining to Eq. \ref{e:x4_condition}, we tune the parameter $c$ to set the initial conditions on $x_4$ and $x_5$ such that $x_4\sim \mathcal{O}(1)$ at late times. Finally the initial conditions for $\lambda$ depend on the shape of the potential, given that the expansion phase occurs around the minimum of the potential.
 
 \begin{figure}[t]
    \centering
    \includegraphics[width=.48\linewidth]{images/trackerG5_Hz.pdf}
    \hspace{6pt}\includegraphics[width=.48\linewidth]{images/trackerG5_wphi.pdf}
    \caption{Evolution of the normalized Hubble parameter (left) and the equation of state of the dark energy scalar field (right) for the tracker potential $V_1(\phi)$ of Model I, and for different values of $\tilde c$, where $\tilde c$ is given by $c=\tilde c\times 10^{120}$ such that $\theta(\phi)H^2\sim \mathcal{O}(1)$ at late times in Planck units, where $c$ has a dimension of $1/M_{Pl}^{8}$. Evidently, at high redshifts, the Hubble parameter is the same in all the different cases. However, it rises towards larger values at the present epoch in the presence of the GNMDC coupling, thereby alleviating the tension with $\Lambda$CDM. The equation of state also evolves to the present epoch and becomes more negative  for the  GNMDC coupling, as compared to the quintessence case. At high redshifts, the GNMDC term is subdominant and thus, the model remains indistinguishable from the quintessence.}

   \label{fig:trackerG5}
\end{figure}

 In Fig. \ref{fig:trackerG5}, we show the evolution of the normalized Hubble parameter (left panel) and the evolution of the dark energy equation of state (right panel), with redshift for the potential $V_1(\phi)$. As we increase the strength of the GNMDC term by tuning $c$, the value of $H_0$ increases at present, thereby relieving the tension with the $\Lambda$CDM value of $H_0$. We also show the evolution of $\Lambda$CDM and the quintessence model ($c=0$) for reference. The quintessence model, denoted by a vanishing coupling parameter agrees well with $\Lambda$CDM in the early Universe as expected, and also at late times for suitably chosen initial conditions. However, the dynamics of the scalar field mimicking accelerated expansion in the case of quintessence can slightly deviate from $\Lambda$CDM at some intermediate epochs given that the quintessence has an evolving equation of state, in contrast to $\Lambda$CDM. From the plot in the right panel, it is clear that as the strength of the coupling increases, $w_\phi$ approaches smaller values today ($z=0$), which shows that introducing this coupling can also help us get better values for $w_\phi$ \cite{Planck:2018vyg}. 
 More importantly, in this case, the equation of state never becomes phantom-like in order to alleviate the Hubble tension, as is usually the case with the late dark energy models.

Similar to Fig. \ref{fig:trackerG5}, we plot the normalized Hubble parameter and the  equation of state as a function of redshift for $V_2(\phi)$ in Fig. \ref{fig:linearG5}. From the plot in the left panel, we see that in case of linear potential, the setup is more sensitive to the effect of GNMDC, for the alleviation of $H_0$ tension  as compared to the case of $V_1(\phi)$ [c.f. Fig. \ref{fig:trackerG5}]. Apart from this, $w_\phi$ (right panel) evolves to values much closer to $-1$ at $z=0$ compared to the case of $V_1(\phi)$, which is in better agreement with observations \cite{Planck:2018vyg,Riess:2019cxk}.  In general, for tracker potentials of the form of Eq. \ref{e:V1} one obtains the present value of  $w_\phi\sim -0.8$ \cite{Caldwell:2005tm,Clemson:2008ua,Durrive:2018quo}. However, various attempts have been made to formulate generalised tracker models that give rise to more realistic values of $w_\phi$ closer to $-1$ \cite{Urena-Lopez:2020npg}.

\begin{figure}[t]
      \includegraphics[width=.48\linewidth]{images/linearG5_Hz.pdf}
    \hspace{6pt}
\includegraphics[width=.48\linewidth]{images/linearG5_wphi.pdf}
 \centering   
 \caption{Evolution of the normalized Hubble parameter (left) and the equation of state of the dark energy scalar field (right) for the linear potential in case of Model I, and for different values of $\tilde c$, where $\tilde c$ is given by $c=\tilde c\times 10^{124}$ such that $\theta(\phi)H^2\sim \mathcal{O}(1)$ at late times in Planck units, where $c$ has a dimension of $1/M_{Pl}^{8}$. At high redshifts, the Hubble parameter is the same in all the different cases. However, it rises towards larger values at the present epoch in the presence of the GNMDC coupling, thereby alleviating the tension with $\Lambda$CDM. The equation of state also evolves to the present epoch and becomes much closer to $-1$ for the  GNMDC coupling, as compared to the quintessence case. At high redshifts, as one would expect, the GNMDC model remains indistinguishable from $\Lambda$CDM or the quintessence.}
    \label{fig:linearG5}
\end{figure}
As discussed in Section \ref{sec:horndeski_intro}, for consistent model building, it is necessary to ensure that the sound speed of perturbation modes is free of gradient instabilities as well as superluminal propagation. As we show in Fig \ref{fig:stabilityG5}, for Model I with the potential $V_2(\phi)$, though the scalar modes remain subluminal, tensor modes become superluminal and violate the observational bounds given by Eq.  \ref{e:cTbound} at late times. We also show in Fig \ref{fig:stabilityG5} the evolution of the Horndeski functions for this setup, where we can see that such type of models give  non-zero values of $\alpha_T$ at low redshifts and can, in principle,  violate the bound in Eq. \ref{e:alphaTbound}. 
As evident from Fig. \ref{fig:stabilityG5}, all the Horndeski functions evolve to non-zero values in the late Universe. Such models with $\alpha_M \neq 0$ exhibit the running of the effective Planck mass due to the evolution of $\alpha_M$ in the late Universe, and therefore, could be possibly constrained by the future gravitational wave distance measurements \cite{Matos:2022uew,Amendola:2017ovw,Maggiore:2019uih}. 
\begin{figure}[t]
    \centering
    \includegraphics[width=.48\linewidth]{images/cscTplot_G5.pdf}
    \hspace{6pt}\includegraphics[width=.48\linewidth]{images/alphasG5.pdf}
    \caption{The sound speed of scalar and tensor perturbations (left) and the Horndeski functions (right) for Model I with the potential $V_2(\phi)$. At high redshift, both the sound speeds are equal to unity but the tensor perturbations become superluminal at late times while the scalar perturbations remain subluminal. Also, all the Horndeski functions vanish at early times but evolve to different non-zero values at late times. In particular, $\alpha_T > 0$ leads to superluminal tensor modes while $\alpha_M \neq 0$ induces the running of the effective Planck mass with redshift.}
    \label{fig:stabilityG5}
\end{figure}


\section{Model II}
\label{sec:model_II}
As discussed in the previous section, although Model I succeeds in the alleviation of the $H_0$ tension, it encounters the issue of superluminality of tensor modes since $\alpha_T$ [c.f. Eq. \ref{e: alphaT}] is found to be quite large ($\alpha_T \sim 0.04$), which is in strong violation of the observational constraints given in Eq. \ref{e:alphaTbound}. In order to solve this problem that renders an otherwise successful model inconsistent, we propose a new model that is constructed by setting $\alpha_T=0$, to begin with. This is done by choosing $G_{4,X}=G_{5,\phi}$ as well as $G_5=G_5(\phi)$ as already considered in Model I, thus making $\alpha_T=0$ following Eq. \ref{e: alphaT}. As is evident, $G_{4,X}=G_{5,\phi}$ gives $G_4= \frac{1}{2}+X G_{5,\phi}+A(\phi)$, where $A(\phi)$ is an integration constant that we set to zero for simplicity.
Hence the setup for Model II is as follows,

\begin{equation}
    G_2=X-V(\phi);\quad G_3=0; \quad G_{4,X}=G_{5,\phi}=-\frac{\theta(\phi)}{2};\quad G_4=\frac{1}{2}+XG_{5,\phi};\quad G_{5,X}=0
\label{e:G4+G5_setup}
\end{equation}

\subsection{Methodology and analysis}
\label{sec:model_II_method}
In a manner similar to that adopted in Section \ref{sec:model_I_method}, we can obtain the Friedmann and Klein-Gordon equations in the case of Model II as 
\begin{equation}
    3H^2=\frac{1}{2}\dot\phi^2\left(1+9H\dot\phi\theta^\prime(\phi)\right)+\rho_M+V(\phi)
    \label{e:friedmann_first_modelII}
\end{equation}
\begin{equation}
    \dot\phi^2+\frac{9}{2}\dot\phi^3 H \theta^\prime(\phi)-\frac{3}{2}\dot\phi^2 \ddot\phi \theta^\prime(\phi)-\frac{1}{2}\theta^{\prime\prime}\dot\phi^4+\rho_M+w_M\rho_M+2\dot H=0
    \label{e:friedmann_second_modelII}
\end{equation}
\begin{equation}
    -18 H^3 \theta \dot\phi-6 H^2 \theta  \ddot\phi+\frac{15}{2} H^2 \theta^\prime(\phi) \dot\phi^2-12 H \dot H \theta  \dot\phi+9 H
   \theta^\prime(\phi) \ddot\phi\dot\phi+3 H
   \theta^{\prime\prime}(\phi)\dot\phi^3+3 H \dot\phi+\frac{9}{2} \dot H \theta^\prime(\phi) \dot\phi^2+V^\prime(\phi)+\ddot\phi=0
   \label{e:kg_ModelII}
\end{equation}
The energy density and pressure of the dark energy scalar field in this setup are given by,
\begin{eqnarray}
    \rho_\phi&=&V(\phi)+\frac{1}{2}\dot\phi^2\left(1+9H\dot\phi\theta^\prime(\phi)\right)
    \label{e: rho_phi_G4G5}\\
    p_\phi&=&-\frac{1}{2} \dot\phi^2 \left(3 \theta^\prime(\phi) \ddot\phi+\theta^{\prime\prime}(\phi) \dot\phi^2-1\right)-V(\phi)
    \label{e:p_G4G5}
\end{eqnarray}
As we see from Eqs. \ref{e:friedmann_first_modelII} and \ref{e: rho_phi_G4G5}, for an appropriate choice of $\theta(\phi)$ such that 

\begin{equation}
    \theta^\prime(\phi)H\dot\phi\sim \mathcal{O}(1)
    \label{e:x4_condition_G4G5}
\end{equation} 
at late times, the contribution in the energy density arising from the combination of $G_4$ and $G_5$ terms, as given by Eq. \ref{e:G4+G5_setup}, becomes significant in comparison to the canonical kinetic term. This additional contribution modifies the dynamics of the background evolution at late times as compared to the quintessence model and leads to a large $H_0$ thereby alleviating the Hubble tension.

\begin{figure}[t]
    \centering
    \includegraphics[width=.48\linewidth]{images/trackerG4G5_Hz.pdf}
    \hspace{6pt}
    \includegraphics[width=.48\linewidth]{images/trackerG4G5_wphi.pdf}
    \caption{Evolution of the normalized Hubble parameter (left) and the dark energy equation of state (right) for the tracker potential $V_1(\phi)$ in case of Model II, and for different values of $\tilde c$, where $\tilde c$ is given by $c=\tilde c\times 10^{119}$ such that $\theta^\prime(\phi)\dot\phi H \sim \mathcal{O}(1)$ at late times in Planck units, where $c$ has a dimension of $1/M_{Pl}^{7}$. The left figure is broadly the same as the one in Fig. \ref{fig:trackerG5}. However, the effect of increasing $\tilde c$ on $w_\phi$ is opposite to the corresponding behaviour of Model I.}
    \label{fig:trackerG4G5}
\end{figure}
While the dynamical variables $x_1$ through $x_5$ remain the same as Eq. \ref{e:dynamical_variables_G5} of Model I, we further introduce $x_6$ and $\alpha_{II}$, as the setup of Model II involves a higher derivative of $\theta(\phi)$ as compared to Model I. Hence, in case of Model II, the dynamical variables are defined as,
\begin{eqnarray}
    x_1=\frac{\dot\phi}{\sqrt{6}H};\quad x_2=\frac{\sqrt{V(\phi)}}{\sqrt{3}H}; \quad x_3=\frac{\sqrt{\rho_r}}{\sqrt{3}H};\quad \lambda=\frac{-V^\prime(\phi)}{V(\phi)};\nonumber\\
    x_4=\theta(\phi) H^2;\quad x_5=\theta^\prime(\phi) H^2;\quad x_6=\theta^{\prime\prime}(\phi)H^2; \quad \alpha_{II}=\frac{\theta^{\prime\prime\prime}(\phi)}{\theta^{\prime\prime}(\phi)}
\label{e:dynamical_variables_G5}
\end{eqnarray}
where the subscript II corresponds to Model II. Thus the differential equations in Eqs. \ref{e:eq_x1_dynam_var_G5} through \ref{e:eqns_dynam_var_G5} stay the same for $x_1$ till $x_4$, as well as for $\lambda$ and the others get modified as follows
\begin{eqnarray}
    \frac{\text{d}x_5}{\text{d}N}&=& \sqrt{6} x_6 x_1 + 2 x_5 \frac{\dot H}{H^2}\\
    \frac{\text{d}x_6}{\text{d}N}&=& \sqrt{6} \alpha_{II} x_6 x_1 + 2 x_6 \frac{\dot H}{H^2}\\
    \frac{\text{d}\alpha_{II}}{\text{d}N}&=& \sqrt{6} x_1\alpha_{II}^2(\beta_{II}-1),\quad \text{where}\quad \beta_{II}=\frac{\theta^{\prime\prime\prime\prime}(\phi)\theta^{\prime\prime}(\phi)}{(\theta^{\prime\prime\prime}(\phi))^2}
\label{e:diffeq_var_G4G5}
\end{eqnarray}
whereas $\dot H/H^2$ and $\ddot\phi/H^2$ are given as

\begin{equation}
    \frac{\dot{H}}{H^2} = \frac{\splitfrac{3(w_M (-54 \sqrt{6} x_1^3 x_4^2+3 x_4 (x_1^2 (3 x_1(54 x_1 x_5+\sqrt{6})-2)-2 x_2^2+2)+(x_1^2+x_2^2-1) (9 \sqrt{6} x_1 x_5+1))-3 x_1^4}{ \left(2 (6 x_4-1) x_6+45 x_5^2\right)+18
   \sqrt{6} x_1^3 (3 x_4-1) x_5+x_1^2 \left(9 \lambda  x_2^2 x_5+6
   x_4-1\right)+9 \sqrt{6} x_1 \left(x_2^2-1\right)
   x_5-\left(x_2^2-1\right) (6 x_4-1))}}{-243 x_1^4 x_5^2+18
   \sqrt{6} x_1 x_5 \left(6 x_1^2 x_4-1\right)+12 x_4-2}
\end{equation}
\begin{equation}
    \frac{\ddot\phi}{H^2}= \frac{\splitfrac{27 \sqrt{6} x_1^5 (27 w_M x_5^2-8 x_4 x_6)-81
   x_1^4 x_5 (w_M (24 x_4-1)+1)+36 \sqrt{6} x_1^3 (-w_M
  x_4+x_4+x_6)+9 x_1^2 x_5 (9 (w_M+1) x_2^2}{-9
   w_M+1)+6 \sqrt{6} x_1 \left(1-6 x_4 \left(w_M
   \left(x_2^2-1\right)+x_2^2\right)\right)+486 x_1^6 x_5 x_6-6 \lambda 
   x_2^2}}{-243 x_1^4 x_5^2+18
   \sqrt{6} x_1 x_5 \left(6 x_1^2 x_4-1\right)+12 x_4-2}
\end{equation}

Similar to Model I, we evolve the above equations by providing suitable initial conditions in the early radiation domination epoch. The motivation for choosing initial $x_1$, $x_2$, $x_3$ and $\lambda$ is the same as discussed earlier in Section \ref{sec:model_I_method}. The initial conditions on $x_4$, $x_5$ and $x_6$ are chosen appropriately by tuning the parameter $c$, such that $x_5 \sim \mathcal{O}(1)$ at late times. Here the expression for $\beta_{II}$ in Eq. \ref{e:diffeq_var_G4G5} restricts us to $p\geq 3$ to avoid any inconsistencies or unphysical quantities. Thereby we choose $p=4$ once again, as done in the case of Model I.


We now proceed to discuss our results for the two classes of potentials exemplified by Eqs. \ref{e:V1} and \ref{e:V2} similar to Model I.
In the case of $V_1(\phi)$, the left panel of Fig. \ref{fig:trackerG4G5} shows the evolution of the normalized Hubble parameter for different values of $c$, and as is evident, increasing $c$ leads to more acceptable values of the Hubble parameter at present, which is similar to the behaviour of Fig. \ref{fig:trackerG5} of Model I. The right panel shows the evolution of the dark energy equation of state. Unlike Model I, in this case, increasing the strength of the coupling leads to larger values of $w_\phi$ in comparison to the minimal quintessence model, corresponding to $c=0$. In other words, the hierarchy of $w_\phi$ for different values of $c$ is inverted, as compared to the behaviour of Fig. \ref{fig:trackerG5}.

\begin{figure}[t]
    \centering
    \includegraphics[width=.48\linewidth]{images/linearG4G5_Hz.pdf}
    \hspace{6pt}\includegraphics[width=.48\linewidth]{images/linearG4G5_wphi.pdf}
    \caption{Evolution of the normalized Hubble parameter (left) and the dark energy equation of state (right) for the linear potential in case of Model II, and for different values of $\tilde c$, where $\tilde c$ is given by $c=\tilde c\times 10^{122}$ such that $\theta^\prime(\phi)\dot\phi H \sim \mathcal{O}(1)$ at late times in Planck units, where $c$ has a dimension of $1/M_{Pl}^{7}$. These figures are also similar to Fig. \ref{fig:linearG5} of Model I, except for the small changes in the behaviour of $w_\phi$.}
    \label{fig:linearG4G5}
\end{figure}

In the case of $V_2(\phi)$, the evolution of the normalized Hubble parameter and the dark energy equation of state are shown in Fig. \ref{fig:linearG4G5}. The behaviour is nearly the same as in the case of Model I. Increasing the strength of the coupling leads to an increased value of $H_0$, as well as smaller values of $w_\phi$ going close to $-1$ at present. 

Finally, we show the evolution of the sound speed of perturbations and the four Horndeski functions in the case of Model II and for the potential $V_2$ in Fig. \ref{fig:stabilityG4G5}. It is clear that by choosing the suitable combination of $G_4$ and $G_5$ terms, as discussed earlier in this section, the sound speed of tensor modes $c_T$ turns out to be unity (left panel) and $\alpha_T$ is zero throughout the evolution (right panel). Hence, such a setup resolves the $H_0$ tension and is also consistent with the observational bounds (Eqs. \ref{e:cTbound} and \ref{e:alphaTbound}). The scalar perturbations remain subluminal at late times, similar to the case of Model I, which could be attributed to the non-zero evolution of $\alpha_K$ and $\alpha_B$. Although for  consistency, one requires $0 \le c_s^2  \le 1$, the present observational  data does not constrain very well the propagation speed of scalar perturbations. Previous studies have considered the effect of $c_s < 1$, which leads to the clustering of dark energy and thereby can leave observational imprints in the matter and the CMB power spectrum \cite{Linton:2017ged, Xia:2007km}. At the same time, various studies have also shown that these observational signals due to subluminal $c_s$ are important whenever $c_s^2 \sim \mathcal{O}(10^{-3})$ \cite{Basse:2012wd, dePutter:2010vy}. It is to be noted that the suitable combination of the $G_4$ and $G_5$ terms given by Eq. \ref{e:G4+G5_setup}, also renders the function $\alpha_M$ zero throughout the evolution. This is straightforward from Eqs. \ref{e:effective planck mass} and \ref{e:alphaM} as $M_*^2=1$ in such a setup and thus, no running of the effective Planck mass. Hence, even though Model II is very different from the minimal quintessence setup, it gives the same observational results for gravitational wave distance measurements as the latter.
The clustering of dark energy can thus be used as a novel probe to differentiate this scenario from the other models \cite{Batista:2021uhb}.


\begin{figure}[t]
    \centering
    \includegraphics[width=.48\linewidth]{images/cscTplot_G4G5.pdf}
    \hspace{6pt}\includegraphics[width=.48\linewidth]{images/alphasG4+G5.pdf}
    \caption{The sound speed of scalar and tensor perturbations (left) and the Horndeski functions (right) for Model II with the potential $V_2(\phi)$. At high redshift, both the sound speeds are equal to unity and the tensors remain luminal even later along with the scalars being subluminal as earlier. Also, $\alpha_M$ and $\alpha_T$ vanish at all times for this setup but the other two functions evolve to non-zero values at late times. }
    \label{fig:stabilityG4G5}
\end{figure}

\section{Conclusion and Discussions}
\label{sec:conclusion}

The Hubble tension is certainly one of the most intriguing problems in cosmology today and has recently received a lot of attention to address it. Numerous attempts have been made as well as ongoing in order to bridge the gap between the $\Lambda$CDM predicted value of $H_0$ and the values obtained through the late-time observations. While the mismatch can certainly be attributed to the systematic errors associated with different measurements, the tension may also be a true reflection of reality, requiring exotic new physics and, perhaps, a dramatic revision to our current understanding of cosmic evolution.

In this work, we propose two dark energy models within the framework of Horndeski gravity which are found to be quite successful in the alleviation of the $H_0$ tension. However, in order for a model to be consistent, one needs to keep a check on the propagation speed of the scalar and tensor perturbation modes in order to avoid issues such as gradient instability or superluminal propagation, even if the background evolution is fairly stable. Using these arguments, we establish the preference of one of our proposed models over the other. Our first model (Model I) incorporates an evolving $G_5$ term of the Horndeski Lagrangian as a function of the scalar field $\phi$, which involves a coupling between the Einstein tensor and the derivative of the scalar field. The strength of this coupling can be tuned in a controlled manner using a dimensionful parameter $c$ such that we can obtain a larger value of $H_0$ at late times. In order to solve the Friedmann and Klein-Gordon equations, we employ a set of dynamical variables that reduce the second-order differential equations to first-order and enable us to put suitable initial conditions on these variables. We analyze our results for two different classes of scalar field potentials and examine the alleviation of the $H_0$ tension for different values of $c$. Our results show that by increasing the coupling strength $c$, the value of $H_0$ increases, while the dark energy equation of state parameter decreases at present. Although this model works well for resolving the $H_0$ tension, it is plagued by the superluminal propagation of tensor modes (or GWs), and is thereby in strong disagreement with the recent GW observations from LIGO/VIRGO and their electromagnetic counterparts. For these reasons, Model I is not a compatible model, which motivates us to propose another model (Model II) that avoids the superluminal propagation of tensor modes by construction. For Model II, we choose a particular combination of the $G_4$ and $G_5$ terms of the Horndeski Lagrangian in a way such that the function $\alpha_T$ that quantifies the tensor speed excess is zero. Thereafter, we use an approach similar to that used for Model I and proceed to obtain $H_0$ values that are in better agreement with observations. Hence, Model II not only successfully alleviates the $H_0$ tension but is also free from any disagreements or inconsistencies. 

\begin{table}[t]
\begin{center}
	\begin{adjustbox}{width=0.6\textwidth}
	\renewcommand{\arraystretch}{2.2}
		\begin{tabular}{ |c|c|c| } 	
			\hline
		    \textbf{Horndeski functions} & \textbf{Model I} & \textbf{Model II}\\
			\hline
               
		    $\alpha_M$ & \large{$-\frac{\dot\phi\ddot\phi\theta(\phi)+\frac{1}{2}\dot\phi^3\theta^\prime(\phi)}{HM_*^2}$} & 0 \\
		    $\alpha_K$ & \large{$\frac{\dot\phi^2+3\dot\phi^2 H^2\theta(\phi)}{H^2 M_*^2}$} & \large{$\frac{\dot\phi^2+9\dot\phi^3 H\theta^\prime(\phi)}{H^2M_*^2}$}\\
		    $\alpha_B$ & \large{$\frac{2\dot\phi^2\theta(\phi)}{M_*^2}$} & \large{$\frac{3\dot\phi^3\theta^\prime(\phi)}{2H M_*^2}$}\\
		    $\alpha_T$ & \large{$\frac{\dot\phi^2\theta(\phi)}{M_*^2}$} & 0\\
		    \hline
		    \textbf{Sound speeds} &  & \\
		    $c_s$ & $0 \leq c_s^2 \leq 1$ & $0 \leq c_s^2 \leq 1$ \\
		    $c_T$ & $c_T^2 > 0,\quad c_T\geq 1$ & $c_T^2>0, \quad c_T=1$\\
			\hline
		 \end{tabular}
    \end{adjustbox}
\end{center}
\caption{We tabulate all four Horndeski functions, and the scalar and tensor sound speeds for the two models that we have studied in this work. While Model I is plagued by the superluminal propagation of tensor modes and the running of the effective Planck mass, Model II leads to luminal tensor modes and vanishing running of the effective Planck mass. Therefore, Model II is a better working scenario than Model I.}
\label{t:alphas}
\end{table}

The expressions for various  Horndeski functions, as well as the behaviour of the scalar and tensor speeds for the two models, have been summarised in Table \ref{t:alphas}. For Model I, all four Horndeski functions are non-vanishing. In particular, the non-zero evolution of $\alpha_M$ induces a running of the effective Planck mass at late times. Such a scenario also leads to additional friction in the propagation of GWs and, thereby, can constrain $\alpha_M$ by future GW distance measurements. On the other hand, the $\alpha_M$ in Model II is always zero, which can act as a null test for modified gravity models which constrain $\alpha_M$. The functions $\alpha_K$ and $\alpha_B$ are non-zero for both models and affect the speed of propagation of scalar modes, leading them to subluminal values at late times. Moreover, $\alpha_T=0$ in Model II also ensures that the tensor modes propagate at the speed of light, contrary to being superluminal in Model I.

In terms of the dynamics of such dark energy models, the presence of non-canonical coupling terms (i.e. evolving $G_4$ and $G_5$) at late times, especially the $G_5$ term, acts as friction which hinders the rolling of the scalar field down the potential, thereby suppressing the field velocity. Therefore, the evolution of the dark energy equation of state is modified at late times, leading to an increase in the energy density of the dark sector in a controlled way. This is reflected in a large value of the Hubble constant today, without disturbing the early Universe physics.  A very interesting aspect of such models is that they don't run into a phantom regime to solve the Hubble tension, which is considered a necessity in the case of late dark energy models trying to address the same \cite{Alestas:2020mvb,DiValentino:2020naf,Vagnozzi:2019ezj}. However, due to the constraints from other large-scale structure and late Universe observations, one may find that such models can not possibly alleviate the tension completely. Nevertheless, these can still be considered an interesting and viable class of models that help relieve the tension and hence assure a step in the right direction. Needless to mention, an ideal dark energy model would be one which does not affect the pre-recombination physics, thereby governing the dynamics of accelerated expansion at late times and must also be in agreement with various cosmological observations. 

In fact, the next step in this direction would be to study the evolution of linear density perturbations in such Horndeski models to see the effects on CMB and matter power spectrum. Following that, it would be ideal to find the limitations and constraints on the parameter space by employing a detailed Monte Carlo analysis taking into account the observational data  from  CMB, BAO, Redshift Space Distortions,
Cosmic Chronometers and SNIa combined together. We consider this as a future prospect to perform a complete data analysis and explore the parameter space of such late Universe modifications arising within the Horndeski theory, along with examining the possibility of alleviating the $\sigma_8$ tension simultaneously. Finally, it would be very interesting to construct consistent models within this framework, which can simultaneously alleviate the $H_0$ and $\sigma_8$ tensions. We leave this for future work. 

\acknowledgments
BG acknowledges support from the DST-INSPIRE Faculty Fellowship grant No. DST/INSPIRE/04/2020/001534. RKJ acknowledges support from the Science and Engineering 
Research Board, Department of Science and Technology, Government of India, 
through the Core Research Grant No.~CRG/2018/002200.
RKJ also acknowledges financial support from the new faculty seed start-up 
grant of the Indian Institute of Science, Bengaluru, India and the Infosys Foundation, Bengaluru, India through the Infosys Young Investigator Award. BG and YT thank Ruchika Kaushik and Farbod Hassani for clarifying some doubts, and Ujjwal Kumar Upadhyay for useful comments.

\bibliographystyle{JHEP}
