\section{Introduction}

Deep learning is expected to be an important enabler for many wireless communication challenges in 6G. Deep neural networks (DNNs) are being proposed to handle a wide range of wireless communication tasks including encoding/decoding operations, spectrum sensing and RF signal classification. Power allocation for D-MIMO networks is one of the challenging tasks that we see the utilization of DNNs \cite{hexax}.

D-MIMO is a new network type considered for 6G communication systems where many radio units (RUs) are geographically distributed in a region to increase the coverage and reliability. The results obtained in \cite{demir2021foundations} show the benefits of D-MIMO compared to traditional uncoordinated small-cells which are used in all previous generations and collocated massive MIMO which is currently used in 5G systems.

In D-MIMO, a precoding is applied at each RU to spatially focus the jointly transmitted signal on the desired user equipment (UE). There are different precoding techniques studied in the literature. Maximal ratio transmission (MRT) is one of the most attractive techniques where each RU applies a local precoding. To mitigate the multi-user interference, some power control techniques are applied together with precoding. According to the slowly changing channel statistics, RUs allocate their power among users to optimize the system's performance. Max-min fairness (MMF) power control aims to maximize the minimum spectral efficiency (SE) of users, and it is shown that the result satisfies uniform quality-of-service to all users. Future generation communication systems are required to satisfy uniformly good services to all users in an environment with ubiquitous demand, making MMF highly attractive for power control. On the other hand, the analytical solution of MMF requires high complexity operations.

To overcome the complexity problem, Bashar et al. \cite{bashar2020exploit} propose deep learning methods to approximate the analytical solution using a low-complexity artificial intelligence (AI) model. Although AI models can accurately approximate the analytical result, it is known that they have some intrinsic properties which make them vulnerable to adversarial attacks~\cite{szegedy, Tuna2022}. Adversarial machine learning has lately attracted significant attention in wireless security domain as DNNs have become more commonly deployed.

In this study, we investigate the potential effects of adversarial attacks targeting AI-driven power control systems in D-MIMO. We explain the main constraints of the adversary resulting from the distributed nature of wireless domain and focus only on the possible practical scenarios to observe the severity of adversarial attack threats. We work on attacks based on universal adversarial perturbation (UAP) which are not based on the exact input samples of the DNN but instead generalize the characteristics of a randomly selected subset of input samples. We propose a novel modified UAP (m-UAP) technique that crafts a specific perturbation for each input where there is only a partial knowledge about some of input entries. By performing several numerical simulations, we show that m-UAP with a surrogate AI model has much larger disruptive effect than standard Gaussian noise perturbation. To the best of our knowledge, this is the most practical work on the security of AI-driven power allocation in D-MIMO.

The organization of the paper is as follows. Section II describes the system model. In Section III, we present possible adversarial attacks. Section IV includes detailed numerical simulations, and Section V concludes the paper.

\section{System Model}
We consider a D-MIMO network with $M$ single antenna RUs which are connected to a central processor (CP) via wired fronthaul links. We assume that $K$ single antenna UEs are jointly served by all RUs in a given time/frequency resource block. In Fig. 1, we present an example D-MIMO network with 16 RUs, 4 UEs, and a CP.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.35\textwidth]{d_mimo_fig1.png}
    \caption{An example D-MIMO network.}
    \label{fig1}
\end{figure}

RUs apply MRT precoding where conjugate of the channel vectors are used to precode the user data. In this case, the transmitted signal from $m$-th RU can be written as
\begin{equation}
    x_m = \displaystyle\sum_{k=1}^K w_{m,k}s_k
\end{equation}
where $w_{m,k}=\sqrt{P_t\eta_{m,k}}h_{m,k}^{*}$ is the precoding vector of the $m$-th RU for the $k$-th user, $s_k$ is the information symbol of the $k$-th user satisfying $\mathbb{E}[|s_k|^2]=1$ for all $k$, $P_t$ is the maximum total transmit power of each RU, $\eta_{m,k}$ is the power control coefficient, and $h_{m,k}$ is the channel coefficient between RU $m$ and UE $k$. To satisfy transmit power constraints, we need
\begin{equation}
    \mathbb{E}[|x_m|^2] = \displaystyle\sum_{k=1}^K \mathbb{E}[|w_{m,k}|^2] = \displaystyle\sum_{k=1}^K \eta_{m,k}\beta_{m,k} \leq P_t, \quad \forall m
\end{equation}
where $\beta_{m,k}$ is the large-scale fading coefficient of the channel between RU $m$ and UE $k$, i.e., $\mathbb{E}[|h_{m,k}|^2]=\beta_{m,k}$.

The received signal for the $k$-th user can be written as
\begin{equation}
    r_k = \displaystyle\sum_{m=1}^M h_{m,k} \displaystyle\sum_{\ell=1}^K w_{m,\ell}s_{\ell} + z_k
\end{equation}
where $r_k$ is the received signal and $z_k\sim \mathcal{C}\mathcal{N}(0, \sigma_k^2)$ is the receiver noise for the $k$-th user. Here, we assume that the channel coefficients $h_{m,k}$ are perfectly known at RU side.  The channel estimation can be performed using uplink pilots in time-division duplex (TDD) mode and when the channel coherence time is large enough, the channel coefficients can be estimated accurately.

We assume that the users only have statistical channel knowledge, and they know the mean value of the effective channels. This assumption is widely used in the literature \cite{demir2021foundations},\cite{ngo2017cell}, \cite{nayebi2017precoding}. When the number of RUs is large, by means of channel hardening, the effective channels become nearly deterministic and hence, statistical knowledge is sufficient. Under this assumption, the received signal of the $k$-th user is given by
\begin{equation}
    r_k = r_{k, \text{desired}}+r_{k, \text{mismatch}}+r_{k, \text{interference}}+r_{k, \text{noise}}
\end{equation}
where
\begin{align}
    r_{k, \text{desired}} &= \displaystyle\sum_{m=1}^M\sqrt{P_t\eta_{m,k}}\mathbb{E}[|h_{m,k}|^2] \\
    r_{k, \text{mismatch}} &= \displaystyle\sum_{m=1}^M\sqrt{P_t\eta_{m,k}}|h_{m,k}|^2 - \sqrt{P_t\eta_{m,k}}\mathbb{E}[|h_{m,k}|^2] \notag \\
    r_{k, \text{interference}} &= \displaystyle\sum_{\ell \neq k}\displaystyle\sum_{m=1}^M\sqrt{P_t\eta_{m,\ell}}h_{m,k}h_{m,\ell}^{*}, \: \: r_{k, \text{noise}} = z_k. \notag
\end{align}
Here $r_{k, \text{desired}}$ is the desired signal part, $r_{k, \text{mismatch}}$ is the part involving the uncertainty of the effective channel, $r_{k, \text{interference}}$ includes multi-user interference, and $r_{k, \text{noise}}$ is the noise part.

Using the approach given in \cite{demir2021foundations} and \cite{medard2000effect}, achievable user SEs can be given as \cite{nayebi2017precoding}
\begin{equation}
\begin{aligned}
\text{SINR}_k &= \dfrac{P_t\left(\displaystyle\sum_{m=1}^M \sqrt{\eta_{m,k}}\beta_{m,k}\right)^2}{P_t\displaystyle\sum_{\ell=1}^K\displaystyle\sum_{m=1}^M \eta_{m,\ell}\beta_{m,\ell}\beta_{m,k}+\sigma_k^2} \\
\text{SE}_k &= \log_2(1+\text{SINR}_k), \quad \forall k.
\end{aligned}
\end{equation}
Energy efficiency (EE) is another important parameter of networks which is aimed to be as large as possible. It is generally defined as \cite{ngo2017cell}
\begin{equation}
\text{EE} = \dfrac{\displaystyle\sum_{k=1}^K \text{BW}_k\cdot \text{SE}_k}{P_t\displaystyle\sum_{m=1}^M \displaystyle\sum_{k=1}^K\eta_{m,k}\beta_{m,k}}
\end{equation}
where EE is the average energy efficiency, $\text{BW}_k$ is the frequency bandwidth allocated for the $k$-th user. Here we ignore the energy consumption related to the fronthaul links and the other parts of the network, and focus on the energy efficiency of the access link only. EE shows the average number of bits transmitted per consumed energy unit.

To optimize the power allocation, we consider MMF problem (P0) which can be defined as
\begin{equation}
(\text{P}0) \: \: \underset{\eta_{m,k}}{\text{max}} \: \min_k \text{SE}_k \quad \text{s.t.} \quad \displaystyle\sum_{k=1}^K \eta_{m,k}\beta_{m,k} \leq 1, \: \: \forall m.
\end{equation}
The optimal solution of the problem (P0) requires iterative quadratic programming which is of high complexity. To find an approximate and low-complexity solution deep learning method is proposed in the literature \cite{bashar2020exploit}. In this study, we assume that a DNN model is implemented at CP. The model is trained using an input data consisting of $\beta_{m,k}$'s and the output $\eta_{m,k}$'s used for training are calculated by solving (P0) analytically and in an offline fashion. The model details and the channel models to produce $\beta_{m,k}$'s are described in section \ref{simulation_result}.

\section{Adversarial Attacks in D-MIMO}

\subsection{Adversarial Attacks}

The goal of the adversarial attacks is to craft a perturbation $\delta$ under given constraints ($\delta<\epsilon$) which yields to an incorrect prediction as $f(x+\delta)=y^{adv}$ which differs from a prediction on a clean sample $f(x)=y$. The success criteria of the attack might change depending on the type of task. For a classification task, the attack can be considered successful if the model predicts a class other than the actual class. However, for regression tasks, many different objectives can be considered. In this study, our objective is to minimize user SEs.

Fast Gradient-Sign Method (FGSM) \cite{goodfellow2015explaining} is one of the most popular adversarial attacks in literature which utilizes the derivative of the model's loss function with respect to the input to determine in which direction the feature values of the input vector should be changed to minimize the loss function of the model. Once this direction is extracted, it changes all features simultaneously in the opposite direction to maximize the loss. Later, Kurakin et al. \cite{kurakin2017adversarial} proposed a small but effective improvement to the FGSM, known as Basic Iterative Method (BIM). In this approach, rather than taking only one step of size $\epsilon$ in the gradient sign's direction, the attacker takes several but smaller steps $\alpha$, and use the given $\epsilon$ value to clip the result. Crafting adversarial samples under $L_{\infty}$ norm for BIM attack is given by (\ref{eq:bim}).
\begin{equation}
\begin{aligned}
\mathbf{x}_1^* & = \mathbf{x} \quad \text{and for all} \: \: i=1, 2, \ldots, i_{\text{max}} \\
\mathbf{x}_{i+1}^* & = \text{Clip}_{\epsilon} \{ \mathbf{x}_i^* + \alpha \cdot \text{sign} \left( \nabla_{\mathbf{x}} \left(J(\mathbf{x}_i^*, \bm{\theta}) \right)\right) \}
\end{aligned}
\label{eq:bim}
\end{equation}
where $\mathbf{x}$ is the input sample, $\mathbf{x}_i^*$ is the crafted adversarial sample (to be added to the input of the victim) at $i$\textsuperscript{th} iteration, $\bm{\theta}$ is the vector containing AI model weights, $J$ is the objective function, $\epsilon$ is a tunable parameter, limiting maximum level of perturbation for $L_{\infty}$ norm, $\alpha$ is the step size, and $\text{Clip}_{\epsilon}\{\cdot\}$ is the clipping operator that clips entries of the argument larger than $\epsilon$ to $\epsilon$ and smaller than -$\epsilon$ to -$\epsilon$.

\subsection{Sources of Adversarial Attack Threats in D-MIMO}

In Fig. \ref{fig2}, we demonstrate a D-MIMO network with potential attack sources. The first possible attack scenario is related to malicious RUs or fronthaul links. A man-in-the-middle (MITM) attack can be applied to change the channel related data transmitted from RUs to CP by an adversary. This type of attack enables to perturb the channel data of RUs on which the attack is applied. The second scenario is related to the UEs. There may be some malicious UEs in the network with modified RF and baseband components so that the corresponding channel information of those UEs are perturbed via applying some perturbation to pilot signals. In TDD systems, the channel knowledge is obtained at RU side using uplink pilots received from UEs. Therefore, a change in pilot signalling can produce false channel knowledge at RU side. We will analyze the effects of both these attack types in section \ref{simulation_result}.
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.35\textwidth]{d_mimo_fig2.png}
    \caption{D-MIMO network with potential attacks.}
    \label{fig2}
    \vspace{-5mm}
\end{figure}
\subsection{Practical Limitations of the Adversary}
From the adversary's perspective, there are three important constraints which limits the success of the adversarial attack in a D-MIMO network. Firstly, the adversary mostly does not have access to the details (architecture and weights) of the original AI model, therefore cannot use it in a white-box setting for crafting adversarial samples. Secondly, the adversary may not have complete knowledge of the input features of the AI model. Because, it is almost impossible for the adversary to know the channel information of each UE. Lastly, in a practical scenario, the adversary does not have the capability to introduce perturbations to all parts of the input vector, even if the channel information is known beforehand. However, despite all these limitations, there are proven ways in literature which increase the success of the attacker. Regarding the first limitation, it has been shown that a surrogate AI model might be sufficient to launch an effective attack \cite{transferability} due to the transferability nature of the adversarial samples. Regarding the other limitation, the universal adversarial perturbation (UAP) method \cite{uap} is proposed for cases where the complete input knowledge is not available.

\subsection{Proposed modified UAP (m-UAP) Method}

To craft a perturbation when there is partial input knowledge, we propose a modified UAP (m-UAP) method which is based on the version suggested by Santos et al. \cite{santos} as in Algorithm 1 where principal component analysis (PCA) is applied to a set of output vectors obtained by the attack method.

\begin{algorithm}[!ht]
\DontPrintSemicolon
  \KwInput{$\textbf{x}, \epsilon$}
  \KwData{$\{\textbf{x}_1, \textbf{x}_2, \ldots, \textbf{x}_N\}, \widetilde{f}(\cdot,\bm{\theta})$}
  \KwOutput{$\bm{\delta}$}
  Define a matrix $\widehat{\textbf{X}}^{N \times MK}$ using the known entries of the input $\textbf{x}$ and the vectors $\{\textbf{x}_1, \textbf{x}_2, \ldots, \textbf{x}_N\}$ taken from the test set. Firstly, initialize the matrix $\widehat{\textbf{X}}$ as $\widehat{\textbf{X}}=[\textbf{x}_1 \: \textbf{x}_2 \: \cdots \textbf{x}_N]^T$. Then update $\widehat{\textbf{X}}$ using the known entries of $\textbf{x}$, i.e., $[\widehat{\textbf{X}}]_{i,j}=[\textbf{x}]_j$ for all known term indices $j$ and for all $i=1, 2, \ldots, N$.  \\
  For each row $\widehat{\textbf{x}}_i^T$ of $\widehat{\textbf{X}}$, apply the BIM to generate the matrix $\bm{P}^{N \times MK}=[\bm{\rho}_{\widehat{\textbf{x}}_1}, \bm{\rho}_{\widehat{\textbf{x}}_2}, \ldots, \bm{\rho}_{\widehat{\textbf{x}}_N}]^T=[\nabla_{\widehat{\textbf{x}}_1}J (\widetilde{f}(\widehat{\textbf{x}}_1, \bm{\theta})), \nabla_{\widehat{\textbf{x}}_2}J (\widetilde{f}(\widehat{\textbf{x}}_2, \bm{\theta})), \ldots, \nabla_{\widehat{\textbf{x}}_N}J (\widetilde{f}(\widehat{\textbf{x}}_N, \bm{\theta}))]^T$ \\
  Compute the principal right singular vector $\textbf{v}_1$ of $\textbf{X}$ as $\bm{P}=\textbf{U}\bm{\Sigma}\textbf{V}^H$, and $\textbf{v}_1$ is the first column of $\textbf{V}$. \\
  Compute two perturbations $\bm{\delta}_1=\epsilon\cdot \text{sign}(\textbf{v}_1), \: \bm{\delta}_2=-\bm{\delta}_1$ and calculate the sum objectives corresponding to the perturbed inputs, i.e.,  $J_u=\displaystyle\sum_{i=1}^N J(\widehat{\textbf{x}}_i+\bm{\delta}_u, \bm{\theta})$ for $u=1, 2$. Find the index $u_0 \in \{1, 2\}$ such that $u_0 = \underset{u \in \{1, 2\}}{\text{argmax}} \: J_u$. \\
  $\bm{\delta}=\bm{\delta}_{u_0}$. \\
  \Return $\bm{\delta}$
\caption{PCA-based modified UAP (m-UAP) method under $L_{\infty}$ norm}
\end{algorithm}

In our algorithm, $\widetilde{f}(\cdot, \bm{\theta})$ is the surrogate AI model used by the adversary, $\epsilon$ is the maximum allowed perturbation amount per each entry, $J$ is the objective function which is defined as the sum SE of users, i.e., $J(\textbf{x}, \bm{\theta}) = \displaystyle\sum_{k=1}^K \text{SE}_k$, and $\{\textbf{x}_1, \textbf{x}_2, \ldots, \textbf{x}_N\}$ are the $N$ randomly selected input samples from the test data set. We generate a matrix $\widehat{\textbf{X}}$ using the known entries of the input vector $\textbf{x}$ and other unknown elements are selected from $\textbf{x}_i$'s. Because, in a practical scenario, the adversary will most likely has access to the AI input data of some UEs or RUs only. Therefore, some part of the AI model's input is already known and modifiable by the adversary. After the generation of the matrix $\widehat{\textbf{X}}$, we apply BIM to each row to obtain the accumulated gradient ($\bm{\rho}_{\widehat{\textbf{x}}_i}$) of the objective function for all $i=1, 2, \ldots, N$. The BIM results are collected in a matrix $\textbf{P}$ and PCA is applied to find the principal right singular vector $\textbf{v}_1$. By this, we are able to obtain a generalized perturbation direction that reflects $\bm{\rho}_{\widehat{\textbf{x}}_i}$'s common characteristics. As $-\textbf{v}_1$ is also a valid right singular vector corresponding to the same singular value, we check which one to choose by evaluating the corresponding sum SEs of users. Here we evaluate the SE values using the input vectors obtained by the rows of $\widehat{\textbf{X}}$. Perturbations are calculated using the direction of the corresponding singular vectors and the perturbation limit $\epsilon$ to limit $L_{\infty}$ norm of the perturbation vector.

Lastly, it is important to choose the right objective function in a practical manner. In our simulations, our objective was to minimize the sum SE of UEs. If we instead chose to maximize sum of allocated powers above a predetermined limit as in \cite{santos}, then the attack would easily be mitigated by the network by normalizing the estimated power values.

\section{Simulation Results}\label{simulation_result}

\subsection{Experimental Setup}

We setup a simulation environment in MATLAB to generate our experimental training and test data sets to be used in model training and attack scenarios. We assume that there are 16 RUs distributed on a uniform grid in a $500 \text{m} \times 500 \text{m}$ area. There are 4 UEs and their locations are randomly chosen in the same region at each trial. All RU-UE channels have some large and small-scale fading parts, and they are assumed to be independent and constant within a coherence block. Large-scale part consists of path-loss with a 3-slope model and log-normal shadowing with a standard deviation $8$ dB, and small-scale part is assumed to be Rayleigh. We take $P_t=0.2$ W, $\sigma_k^2=-92$ dBm and $\text{BW}_k=20$ MHz for all $k$. These parameters are taken from \cite{ngo2017cell} and \cite{nayebi2017precoding}. We generate 400K pair of data consisting of large-scale fading coefficients ($\beta_{m,k}$'s) and their associated optimum power control coefficients ($\eta_{m,k}$'s). We use 390K pairs for model training and 10K pair of data for our tests.

We use DNN-based regression models to learn the mapping between the large-scale fading coefficient vector $\textbf{x}=[\beta_{1,1} \: \beta_{1,2} \: \ldots \: \beta_{M,K}]^T \in\mathbb{R}^{MK}$ and the power control coefficient vector $\bm{\eta}=[\eta_{1,1} \: \eta_{1,2} \: \ldots \: \eta_{M,K}]^T \in\mathbb{R}^{MK}$. We train two different AI models which we denote as $f$ and $\widetilde{f}$ for the original AI model used in D-MIMO system and the surrogate model that the adversary uses to craft adversarial samples, respectively. Detailed model architectures are given in Fig. \ref{fig3} and \ref{fig4}. In Fig. \ref{fig5}, we present the per-user CDFs for the original and surrogate AI models, together with the result of analytical solution. We observe that surrogate model performance is slightly worse than the original one. Furthermore, both models can accurately find an approximate MMF solution as there is a small gap between the CDF curves of AI models and the CDF curve of the analytical result.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.9\textwidth]{Figures/orig_new.png}
\caption{Original AI Model}
\label{fig3}
\vspace{-5mm}
\end{figure*}
\begin{figure*}[t!]
\centering
\includegraphics[width=0.9\textwidth]{Figures/surrogate_new.png}
\caption{Surrogate AI Model}
\label{fig4}
\vspace{-5mm}
\end{figure*}

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.38\textwidth]{paper_fig0.eps}
    \caption{Performance of original and surrogate models.}
    \label{fig5}
    \vspace{-2mm}
\end{figure}

\subsection{Experimental Results}
In this section, we perform several simulations to see the effects of adversarial attacks under several constraints.

We assume that the original model $f$ is used in CP to determine the power allocation. However, we assume that the adversary mainly uses the surrogate model $\widetilde{f}$ to craft adversarial samples. To see the performance gap from adversary's point of view, we also present the results in the case where the adversary has access to the original model. In our simulations, we present the effects of the adversarial attacks on per-user SEs and average EE with different amounts of perturbations and different levels of input information.

\subsubsection{The ultimate performance of adversarial attacks}
We begin our simulations by firstly showing the extreme scenario with the most devastating consequence where the adversary has access (read/modify) to all the parts of input vector fed to the AI model in CP. We set $\epsilon = 8$ dB which is equal to the standard deviation of the shadowing in the system. In Fig. \ref{fig6}, we see the CDFs of per-user SEs under different attack types.
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.38\textwidth]{paper_fig1.eps}
    \caption{Comparison of different attack types ($\epsilon = 8$ dB).}
    \label{fig6}
    \vspace{-5mm}
\end{figure}
It is clear that m-UAP attack results in much more devastating consequences than standard Gaussian perturbation. Under the attack, the performance of the analytical solution also degrades showing that adversarial training is not suitable for this regression task. Because applying adversarial training will degrade the natural (clean) performance of the system. This shows the necessity of developing new defense solutions for such kind of smart attack threats. Finally, we observe that surrogate model (indicated as blackbox) has slightly worse disruptive performance than the original model (indicated as whitebox). This result proves that the adversary does not need to have access to the original AI model for crafting effective adversarial samples.

\subsubsection{The performance with partial input knowledge}

To cover more practical cases which might be seen in a real-world scenario, we assume that the adversary might have a partial channel knowledge and partial perturbation capability on the input. In Fig. 7-8, we present the median (the SE value corresponding to $\text{CDF} = 0.5$) and the 5th percentile (the SE value corresponding to $\text{CDF} = 0.05$) per-user SEs for various levels of input information and perturbation capability. Fig. 7 shows the results of the first type of attack scenario which might be originated by malicious RUs. Fig. 8 shows the results of the second type of attack which might be launched by malicious UEs. In each of these cases, we assume that some portion of the UEs or RUs are malicious and only their channel information can be known and perturbed by the adversary. In both scenarios, the adversary uses the surrogate model with $\epsilon = 8$ dB.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.37\textwidth]{paper_fig2.eps}
    \caption{The effect of attack on RUs ($\epsilon = 8$ dB).}
    \label{fig7}
    \vspace{-2mm}
\end{figure}
\vspace{-3mm}
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.37\textwidth]{paper_fig3.eps}
    \caption{The effect of attack on UEs ($\epsilon = 8$ dB).}
    \label{fig8}
    \vspace{-2mm}
\end{figure}

According to the results obtained in Fig. 7-8, we conclude that m-UAP attack outperforms conventional approach of applying standard Gaussian noise in all cases. The gap gets larger when the ratio of involved malicious actors in the network increases. When the level of information and perturbation capability about the input become larger, the adversary can degrade the performance more, as expected. In all these adversarial attack cases, we observe a significant decrease in the user SE performance of the D-MIMO network.

\subsubsection{The effect of amount of perturbation}

The results obtained so far assumes $\epsilon=8$ dB. Considering the shadowing standard deviation (which is equal to $8$ dB), it is theoretically very hard for the system to detect an attack with $\epsilon = 8$ dB as the probability of observing that much variation in large-scale fading coefficients is roughly $32\%$, which is not negligible. On the other hand, for $\epsilon>16$ dB, the same probability decreases down to $5\%$ making the detection possible.

In Fig. 9, we observe the effect of $\epsilon$ for values up to $16$ dB. We consider the scenario where half of the RUs are malicious and employ m-UAP method using a surrogate model.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.36\textwidth]{paper_fig4.eps}
    \caption{The effect of different $\epsilon$ values.}
    \label{fig9}
    \vspace{-1mm}
\end{figure}

The results show that adversarial attacks are more effective than standard Gaussian noise, and the attack becomes more disruptive as $\epsilon$ increases.
\subsubsection{The effect on energy efficiency}

Lastly, we consider the impacts of the adversarial attacks on the EE in the network as sustainability is one of the important concerns in 6G. The power consumption may differ under attack, so we need to analyze the EE separately. The average EE for various $\epsilon$ are analyzed for the case where half of the RUs are malicious and we employ m-UAP method using a surrogate model.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.37\textwidth]{paper_fig5.eps}
    \caption{The effect on energy efficiency for various $\epsilon$ values.}
    \label{fig10}
    \vspace{-1mm}
\end{figure}

In Fig. 10, we see the comparison of no attack, adversarial attack and standard Gaussian noise attack cases. The adversary can decrease EE more than $8$ percent in $\epsilon=8$ dB case and the effect becomes more severe with increasing $\epsilon$ values. We conclude that adversarial attacks can also degrade the energy efficiency of the system.

\section{Conclusion}
In this study, we investigate the potential effects of adversarial attacks on AI-driven power control in D-MIMO. We propose m-UAP method to craft input specific yet PCA-based perturbations by making use of the partial knowledge about the channel. We consider different practical attack types and analyze the disruptive impacts of adversarial attacks under various levels of knowledge and capabilities. The results show that adversarial attacks with optimized perturbations might have a potential to degrade the performance of the network in terms of both spectral and energy efficiency. Smart defense techniques should not be neglected when deploying D-MIMO networks.

\section*{Acknowledgment}
This work was supported by The Scientific and Technological
Research Council of Turkey (TUBITAK) through the 1515
Frontier Research and Development Laboratories Support Program
under Project 5169902, and has been partly funded by
the European Commission through the H2020 project Hexa-X
(Grant Agreement no. 101015956).

\bibliographystyle{IEEEtran}


\section{Introduction}

Recently, the Large Hadron Collider (LHC) at CERN has reached a new, unprecedented center-of-mass energy of $\sqrt{s} = 13.6 \, \mathrm{TeV}$, starting LHC Run 3. This presents the opportunity to re-measure relevant physical quantities at the new energy frontier, thereby checking the predictions of the standard model (SM).

Here, we present the first ever measurement of the top quark pair ($\mathrm{t\bar{t}}$) production cross section in proton-proton collisions \cite{CMS:TOP-22-012}. Compared to the previous center-of-mass energy of $\sqrt{s} = 13 \, \mathrm{TeV}$, the cross section is expected to rise by about 10\% to $921\,^{+18}_{-16}\,(\mathrm{scale})\pm4\,(\mathrm{PDF+}\alpha_s) \, \mathrm{pb}$. The measurement was performed using the CMS detector \cite{CMS:Detector-2008}, which underwent several upgrades inbetween the end of Run 2 and the beginning of Run 3. Since the detection of top quarks covers a wide range of physics objects reconstructed in CMS, this measurement can also be seen as a first test of the CMS calibration and performance on the new data.

\section{Measurement strategy}

Data corresponding to an integrated luminosity of $1.20 \pm 0.07 \, \mathrm{fb}^{-1}$ are analyzed. Events are selected with two charged leptons (electrons and muons) of opposite charge (dilepton channel), and with a single lepton (lepton+jets channel). In order to improve the separation between signal and background events, the presence of hadronic jets, as well as jets identified as originating from b quarks (b-tagged jets), is required. The inclusion of both decay channels in the analysis allows for constraining lepton selection efficiencies, making the measurement independent from external input on calibrations. Events are categorized by the number and flavours of the leptons, number of jets and b-tagged jets.

The $\mathrm{t\bar{t}}$ signal as well as background contributions are estimated using Monte Carlo simulation, with the exception of background originating from non-prompt leptons in QCD events, which is estimated using a data-driven method. The hadronic decay $\mathrm{t} \rightarrow \mathrm{b W} \rightarrow \mathrm{b q q'}$ in the lepton+jets channels is utilized to derive a coarse calibration of the jet energy by constructing a dijet mass distribution and using it to correct the simulation to data. Figure \ref{fig:control_plots} shows, as an example, the agreement between data and simulation in the dimuon and muon+jets channels.

\begin{figure}[!tb]
      \centering
      \includegraphics[width=0.49\textwidth]{njet_mm_all.pdf}
      \hfill
      \includegraphics[width=0.49\textwidth]{njet_mj_all.pdf}
      \caption{%
      The number of jets for data and simulation in the dimuon (left) and the muon+jets channel (right) \cite{CMS:TOP-22-012}.
      }
      \label{fig:control_plots}
\end{figure}

To carry out the measurement of the $\mathrm{t\bar{t}}$ cross section, events are further binned by the number of jets. A maximum likelihood fit over all event categories is performed, and systematic uncertainties are included in the fit as nuisance parameters. This way, both lepton selection and b-tagging efficiencies are determined directly in-situ together with the cross section. The uncertainty on the luminosity is not directly considered in the fit, and instead treated externally. Figure \ref{fig:postfit} shows the post-fit agreement of all categories in the final fit binning.

\begin{figure}[!tbp]
      \centering
      \includegraphics[width=0.8\textwidth]{postfithist.pdf}
      \caption{
            Resulting post-fit agreement and uncertainties in the fit binning after performing the fit \cite{CMS:TOP-22-012}.
      }
      \label{fig:postfit}
\end{figure}

\section{Result}

From the likelihood fit, a cross section of $887^{+43}_{-41}(\mathrm{stat+syst}) \pm 53 (\mathrm{lumi}) \, \mathrm{pb}$ is measured. A comparison to previous measurements from CMS at lower center-of-mass energies, as well as the SM prediction, can be seen in Fig. \ref{fig:xsec}. The statistical and systematic error encompass the full uncertainty of the likelihood fit, while the error on the luminosity is estimated from emittance scans \cite{CMS:DP-2018-011}. The result is in good agreement with the SM, and shows excellent performance from the first CMS data in Run 3.

\begin{figure}[!ht]
      \centering
      \includegraphics[width=0.8\textwidth]{tt_curve_cms_Sep2022.pdf}
      \caption{%
          The $\mathrm{t\bar{t}}$ cross section as a function of $\sqrt{s}$, as obtained in previous measurements by the CMS experiment in pp collisions (blue markers) and at the Tevatron in $\mathrm{p\bar{p}}$ collisions (empty black inverted triangle).
          The red bullet shows the result from this analysis at $\sqrt{s} = 13.6 \, \mathrm{TeV}$ \cite{CMS:TOP-22-012}.
      }

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


\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}