\section{Introduction}
Human motion prediction involves forecasting future human poses given past motion. For enabling efficient Human-Robot Collaboration, a crucial aspect of robot perception is real-time anticipatory modeling of human motion \cite{rudenko2020human,hoffman2019evaluating,Iqbal2017HRI-Bookchapter, Iqbal2016T-RO}. Fluid tasks such as collaborative assembly, handovers, and navigating through moving crowds require combining aspects of perception, representation, and motion analysis to accurately and timely predict probable human motion \cite{biswas2019fast, fosapt, mit_ucsd, islam2020hamlet, islam2020RAL, Iqbal2017RA-L, Iqbal2015TAC}. This would enable the robot to anticipate the human pose and intent and plan accordingly around the human partner without disturbing the natural flow of the human’s motion. However, accurate and timely prediction of human motion remains a non-trivial problem due to the complex and interpersonal nature of human behavior \cite{yasar2021RAL, Iqbal2014Roman}.

To address the aperiodic and stochastic nature of human motion, prior work has framed the problem of predicting future poses like that of sequence learning, modeling the spatio-temporal aspect of human motion using Recurrent Neural Networks  \cite{fragkiadaki2015recurrent,martinez2017human, aksan2019structured,yasar2021RAL}. These approaches aim to learn a unified representation from training samples that are expected to generalize for test data. However, generalization comes at the cost of learning individual subtleties of motion, which is crucial for human-robot collaboration. When training these networks, the core assumption is that the given data points are realizations of independent and identically distributed (i.i.d) random variables. However, this assumption is often violated, e.g., when training and test data come from different distributions (dataset bias or domain shift) or the data points are highly interdependent (e.g., when the data exhibits temporal or spatial correlations) \cite{darrell2015machine}. Both these cases are observed in human motion prediction, making it challenging to deploy models trained on benchmark models to the real world.

While generalization at the cost of learning individual preferences is sub-optimal, there is also a need to learn a robust representation over a diverse range of training samples. As such, training and generalizing over a benchmark dataset cannot be discarded and is, in fact, necessary as the first step to accurate motion prediction. Prior work on language modeling has demonstrated the benefit of learning a rich representation on a large training data followed by fine-tuning on a target task \cite{devlin2018bert,yang2019xlnet,brown2020language}. For human motion prediction, this can be posed as a continual learning problem whereby a motion prediction model acquires prior knowledge by observing a large range of human activities. This is followed by fine-tuning its parameters to accurately capture the subtleties of motion prediction for a particular individual. Such a learning setup, however, brings additional challenges to an already non-trivial problem, with prior work on continual learning demonstrating the risk of \emph{catastrophic forgetting} \cite{lesort2020continual, parisi2019continual}.

To address the challenges mentioned above, we propose a continual learning scheme that can improve human motion prediction accuracy while reducing the risk of catastrophic forgetting. Our framework is modular and is developed to acquire new knowledge and refine existing knowledge based on the new input. In line with prior work on computational neuroscience, which states that the brain must carry out two complementary tasks: generalize across experiences, and retain specific episodic-like events \cite{kumaran2016learning,parisi2019continual}; we utilize a two-phase learning scheme. Our framework aims to learn a robust representation of past observations by training on a benchmark dataset in the first phase. This is achieved by using a modular encoder-decoder architecture with adversarial regularization \cite{yasar2021RAL}, that has state-of-the-art performance on benchmark datasets. In the second phase, we use the representation learning aspect of the framework to condition future poses and fine-tune only the decoder module on new samples in a curriculum learning setup \cite{bengio2009curriculum}. This mitigates the problem of training from scratch while also providing performance gains, both quantitatively over short, mid, and long-term horizons and qualitatively in terms of generating motion that is perceptibly similar to the ground-truth.

\section{Problem Formulation}
Formally defined, human motion prediction is the problem of predicting the future human pose over a horizon, given their past pose and any additional contextual information. In this paper, we assume that there is only one agent in the scene. For any particular scenario, the input to our model is the past or \emph{observed} trajectory frames, spanning time $t=1$ to $\tau$,: $\textbf{X} = \{x_1,\ldots,x_\tau \}$. Each frame $x_t \in \mathbb{R}^N$ denotes the $N$-dimensional body pose. $N$ depends on the number of joints in the skeleton, $J$ and the dimension of the joints $D$, where \textcolor{black}{$N =J \times D$}. The expected output of the model is the future trajectory frames over horizon $H$, i.e. the ground truth pose over the horizon $t=\tau+1$ to $\tau + H$:  $\textbf{Y} = \{ y_{\tau+1},\ldots,y_{\tau+H} \}$.

Our first objective is to learn the underlying representation which would allow the model to generate feasible and accurate human poses $\textbf{\^{Y}} = \{ \hat{y}_{\tau+1},\ldots,\hat{y}_{\tau+H} \}$.
We assume that future human pose is conditioned on the past observed or generated poses and predict each frame in an auto-regressive manner as formulated below:
\vspace{-0.5em}
\begin{equation}
\label{equation: motion prediction}
p_{\theta} (\textbf{\^{Y}}) = \prod^{\tau+H}_{\delta=\tau+1}  p_{\theta}(\hat{y}_\delta | \hat{y}_{\tau:\delta-1}, x_{1:\tau})
\end{equation}
where the joint distribution is parameterized by $\theta$.

Next, we use these learned parameters to fine-tune for a specific agent who was not observed during the training phase, using a continual learning setup. Instead of updating all the model parameters, we update a specific module, say the decoder module, with corresponding parameters $\theta^*$. We formulate this as follows, similar to prior work in continual learning \cite{kirkpatrick2017overcoming}:
\begin{equation}
 \log p\;(\theta^*|D ) = \log p\;(D_B|\theta) + \log p\;(\theta|D_A) - \log p\;(D_B)
\end{equation}
where $D_A$ represents the first phase's training data, which involves learning a representation from the large data distribution. $D_B$ represents the second phase's training data, whereby we aim to learn the parameters for a specific human. $logp(\theta|D_A)$ embeds all the prior information learned during the training phase.

\section{Continual Learning for Human Motion Prediction}
The collective goal of our approach is to accurately predict human motion while being flexible to parameter or architectural updates, given new data. Our overall framework is comprised of an encoder and decoder, trained end-to-end with adversarial regularization on the latent variables, building on top of our prior work \cite{yasar2021RAL}. The encoder aims to learn a rich representation over past trajectories, which the decoder can use to condition its prediction. To improve model stability and robustness of the latent space, we use adversarial regularization through discriminators. This acts as a regularizer during training and can improve the network's stability during parameter updates over new data. We will first describe the overall model for motion prediction and then discuss its flexibility for fine-tuning on a particular agent.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{workshop_arch1.pdf}
    \caption{Motion prediction architecture}
    \label{fig:motion prediction}
\end{figure}

\subsection{Overall model for motion prediction}
\textbf{Motion Encoder:}
The encoder learns a representation over the high-dimensional observed trajectory, projecting the input to a low-dimensional latent space. To obtain a rich and more robust representation over the past trajectories, we extract the past velocity and acceleration features along with the provided positional values, in line with prior work on motion prediction \cite{yasar2021RAL}. The velocity and acceleration features are first and second-order derivatives of the position values for each skeleton joint.

For encoding spatio-temporal representation from the position, velocity, and acceleration data, we employ Recurrent Neural Networks, in particular Gated Recurrent Units (GRU). We use unidirectional GRUs, as we wish to predict human motion in real-time. For each stream, the stream-specific GRU aims to extract the spatio-temporal representation that summarizes the input sequence, with the operation formulated as follows:
\begin{equation}
    h_{s,t} = GRU(h_{s,t-1}, x_{s,t}, \phi_s)
\end{equation}
where $s$ represents the specific stream: position, velocity or acceleration, $x_{s,t}$ denotes the input to the GRU at $t$, $h_{s,t-1}$ corresponds to the past hidden state and $\phi_s$ represent the parameters of the GRU. The output from each GRU represents disparate information corresponding to the past trajectory and needs to be fused adaptively. As such, we use a multi-head self-attention mechanism \cite{vaswani2017attention} which is tasked to disentangle and extract relevant stream-specific representation.
\begin{equation}
    \label{equation: attention mechanism}
    \begin{aligned}
        h_t &= Concat(h_{pos,t}; h_{vel,t}; h_{acc,t}); \\
        h_{att,t} &= Attention (h_t; \phi_{att})
    \end{aligned}
\end{equation}
where $h_{att,t}$ is the output of the attention mechanism, and $\phi_{att}$ represents the parameters. The output $h_{att,t}$ is used to obtain the latent representation, which is tasked to characterize the observed trajectory.

\noindent \textbf{Latent Representation:}
The latent representation aims to capture relevant spatial and temporal semantics from the observed data, which can then be used to condition motion prediction. The latent representation is comprised of a continuous random variable and a categorical random variable. The motivation behind using both continuous and categorical variables is to jointly model the continuous aspect of human motion, such as the spatial semantics of a particular activity and the discrete characteristics of human motion such as the class activity or segment.

To obtain the continuous latent variable $z_t$, the output from the self-attention module is passed through a linear layer. In the case of the categorical latent variable $c_t$, the output from the self-attention module is passed through a linear layer followed by a softmax layer.
\begin{equation}
    \label{equation: latent variables}
    \begin{split}
     z_t &= Linear(h_{att,t}) \\
     h_{c,t} &= Linear(h_{att,t}) \\
     c_t &= softmax(h_{c,t}) \\
\end{split}
\end{equation}

\noindent \textbf{Adversarial Regularization:}
To enforce a prior on the latent space, we use adversarial learning, similar to the Adversarial Autoencoders (AAE) \cite{makhzani2015adversarial} framework. This serves the purpose of a regularizer as there is a modification to the overall objective function: the objective function now consists of a reconstruction loss and an adversarial loss.  We reason that this helps improve the stability of the overall framework for continual learning, as the parameters are updated based on two competing objectives: the reconstruction loss and the discriminator loss.

Similar to the GAN \cite{goodfellow2014generative} and AAE \cite{makhzani2015adversarial} setups, the encoder aims to confuse the discriminators by trying to ensure that its output is similar to the aggregated prior. The discriminators are trained to distinguish the true samples generated using a given prior, from the latent space output of the encoder, thus establishing a min-max adversarial game between the networks \cite{goodfellow2014generative, makhzani2015adversarial}.

We use two discriminators, one for the continuous latent variable and the other for the categorical latent variable. The discriminators compute the probability that a point $z_t$ or $c_t$ is a sample from the prior distribution that we are trying to model (positive samples), or from the latent space (negative sample).
We use a Gaussian prior for continuous latent variables and a uniform distribution prior for categorical latent variables.

\noindent \textbf{Decoder:}
The decoder uses the latent representation and the past generated pose to predict the future pose for each time step. It is auto-regressive, i.e., it uses the output of the previous timestep to predict the current pose and has only one stream: position as the expected output is future joint positions of the human.

The input to the decoder is the latent representation: $z_t$ and $c_t$ and the past generated pose, or the seed pose at time $t$ if it is predicting the first time-step, $t+1$. This is then passed to an attention mechanism that allows the decoder to adaptively condition its output on the latent variables that provide long-term information over the observed frames and the immediate generated frame. The output from this attention mechanism is next passed to a GRU cell, similar to the one at the encoder. This is followed by a Structured Prediction Layer (SPL) \cite{aksan2019structured}, which predicts each joint hierarchically following a skeleton tree, thus allowing the decoder to enforce structural prior on its final output. The operations at the decoder are formulated as follows:
\begin{equation}
    \begin{split}
    p_t &= Concat(z_t, c_t, h_{dec, t-1}) \\
    p_{att,t} &= Attention(p_{t}, \phi_{att}) \\
    h_{dec,t} &= GRU(S_{t-1}, p_{att,t}, \phi_{pos})    \\
    S_{t+1} &= \gamma(h_{dec,t})
    \end{split}
    \label{equation: decoder}
\end{equation}
\subsection{Curriculum learning for the decoder:}
The encoder-decoder architecture with adversarial regularization is trained to convergence on the training set. This training is followed by providing the overall architecture with unseen but small samples of motion data. This aims to relax the i.i.d assumption of the training procedure as our framework now has access to limited motion samples of the agent that it is trying to model.

Our choice of continual learning scheme is the curriculum learning setup \cite{bengio2009curriculum}, whereby we first train the network on a comparatively simpler task of representation learning, followed by a relatively difficult task of fine-tuning its parameters for a specific human subject. Our implementation is based on findings in connectionist models \cite{richardson2008critical, thomas2006computational}, in particular self-organizing maps which reduce the levels of functional plasticity (i.e., ability to acquire knowledge in neural networks) through a two-phase training of the topographic neural map \cite{kohonen1997exploration, kohonen1982self}. The first phase is the organization phase, where the neural network is trained with a high learning rate and large spatial neighborhood size, allowing the network to reach an initial rough topological organization. The second phase is referred to as the tuning phase, where the learning rate and the neighborhood size are iteratively reduced for fine-tuning \cite{parisi2019continual}. We aim to adopt these findings to a sequence learning framework.

Following prior work on developmental and curriculum learning \cite{bengio2009curriculum, graves2016hybrid}, we fine-tuned the architecture on the new data. We adopt techniques that will allow us to retain previous knowledge and avoid catastrophic forgetting during fine-tuning. In particular, we rely on \emph{discriminative fine-tuning} \cite{howard2018universal}, whereby we fine-tune only the decoder network at a different learning rate while freezing the encoder and the discriminator networks.

\begin{table*}[t]
    \centering
    \caption{MSE (in cm\textsuperscript{2}) comparison of fune-tuning vs no fine-tuning on UTD-MHAD for different test subjects) (Lower is better)}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{c| c c c c c c| c c c c c c| c c c c c c}
        \hline
         & \multicolumn{6}{c|}{Subject 2} & \multicolumn{6}{c|}{Subject 4} & \multicolumn{6}{c}{Subject 6} \\
        \hline
         Frames  & 2 & 4 & 8 & 10 & 13 & 15 & 2 & 4 & 8 & 10 & 13 & 15 & 2 & 4 & 8 & 10 & 13 & 15 \\
        \hline
        Zero-Velocity & 11.20 & 27.37 & 66.85 & 86.23 & 112.39 & 127.64
        & 13.11 & 32.33 & 77.22 & 97.80 & 123.99 & 138.53
        & 10.37 & 25.56 & 64.62 & 85.29 & 115.33 & 135.17
        \\
         Our method without Curriculum-Learning & \textbf{5.41} & \textbf{14.75} & 33.87 & 41.89 & 51.33 & 56.12
        & 7.78 & 18.26 & 41.66 & 51.46 & 62.67 & 68.1
        & 6.99 & 16.18 & 38.51 & 49.40 & 63.7 & 71.68
         \\
         \textbf{Our method with Curriculum-Learning} & 7.62 & 16.17 & \textbf{32.35} & \textbf{38.63} & \textbf{45.72} & \textbf{49.13}  &
         \textbf{6.98} & \textbf{14.79} & \textbf{30.29} & \textbf{36.51} & \textbf{43.81} & \textbf{47.61} &
         \textbf{6.51} & \textbf{13.44} & \textbf{27.49} & \textbf{33.66} & \textbf{41.74} & \textbf{46.41}\\
         \hline
    \end{tabular}
    }
    \label{tab:Finetuning all methods}
\end{table*}

\noindent \textbf{Fine-tuning the decoder:}
In line with equation~\ref{equation: latent variables}, the input to the model is the sequence of observed poses: $\textbf{X} = \{x_1,\ldots,x_\tau \}$, with the output being of the encoder being $z_\tau$ and $c_\tau$. However, instead of imposing a prior on the latent space and training the encoder-decoder end-to-end, we only update the decoder's parameters. We also use a lower learning rate and rely on a small number of training samples to improve model stability and reduce the likelihood of catastrophic forgetting.

We leverage the representation learning capability of the framework that it attained when training on a large and diverse dataset. The encoder network is tasked to provide a representation summarizing the past observation that is used by the decoder to condition its prediction, similar to equation \ref{equation: decoder}. The pre-trained weights from the training set are used to initialize the overall architecture and act as prior knowledge. The decoder weights are updated based on the reconstruction loss on the new data.

\section{Experimental Setup}
\subsection{Dataset}
We evaluated the performance of our approach on the widely used human-activity dataset: UTD-MHAD \cite{chen2015utd}. The dataset comprises 27 action classes covering activities from hand gestures to training exercises and daily activities, thus providing relevant activities for human-robot collaboration. Each activity was performed by 8 different subjects, with each subject repeating the activity 4 times. In our experiments, we use only Skeleton data for predicting human motion, following previous work in this domain \cite{jain2016structural,fragkiadaki2015recurrent, martinez2017human,butepage2018anticipating,aksan2019structured}, and considered each of the 20 provided joints. For all experiments, the model predicted output for the next 15 frames, using observation over the past 15 frames.

\subsection{Generalized Representation Learning}
We used the cross-subject evaluation scheme, training and validating on odd-numbered subjects for the first phase, thus providing the framework with a large training sample and maximizing the likelihood of encountering diverse demonstrations. To evaluate the performance, we hold out a section of the data for the validation set and early stopping. This reduces the likelihood of overfitting on the training data while also provide the mechanism for stopping the training procedure.

\subsection{Curriculum Learning for a specific subject}
Having learned a generalized representation, the second phase involved training the framework in a curriculum learning setup. Here, the experiments are conducted on a particular held-out even-numbered subject. We fine-tuned only the decoder using a reduced learning rate, with the encoder weights initialized from the first phase. As each subject has 4 trials, we trained on one trial and tested on the other 3 trials.

\subsection{State-of-the-art method and baseline}
For evaluating the efficacy of our curriculum learning setup, we compared against a non-curriculum learning framework, and the zero-velocity baseline \cite{martinez2017human}. The first benchmark \cite{yasar2021RAL} is comprised of an encoder-decoder framework, with adversarial regularization, but with no provision for curriculum learning. The zero-velocity baseline assumes that all the future predictions are identical to the last observed pose and is challenging to outperform for short-term prediction \cite{martinez2017human, aksan2019structured}. It also allows us to gauge the movement dynamics, with a lower MSE for zero-velocity suggesting less movement and vice-versa for higher MSE.

\subsection{Evaluation Metric}
We evaluated the performance of all models using the Mean Squared Error (MSE), which is the $l_2$ distance between the ground-truth and the predicted poses at each timestep, averaged over the number of joints and sequence length, in line with prior work \cite{butepage2017deep, adeli2020socially, butepage2019imitating, yasar2021RAL}. The MSE is calculated as:

\begin{equation}
    \mathcal{L (X, \hat{X})} =  \frac{1}{T.K}\sum^T_{t=1} \sum^K_{i=1} (x_{t}^{i} - \hat{x}_{t}^{i})^2 \\
\end{equation}

where, T and K are the total number of frame and joints respectively.

\section{Results and Discussion}
\textbf{Results:} We present the results of all approaches on the UTD-MHAD on table \ref{tab:Finetuning all methods}.
We report the performance of all approaches at distinct frame intervals to circumvent the problem of frame drops during data collection and subsequent evaluation \cite{yasar2021RAL}. Our frame intervals aim to evaluate all models on short (2 \& 4), mid (8 \& 10), and long-term motion prediction (13 \& 15). Table I depicts the performance of all approaches with respect to the test subjects. The results in Table I suggest that fine-turning the framework allows it to outperforms all other methods and the zero-velocity baseline for short, mid, and long-term prediction.

\noindent \textbf{Discussion:} Our proposed approach outperformed the prior state-of-the-art approach and baseline both quantitatively by having lower MSE and qualitatively in terms of generating motion closer to the ground-truth pose (see Fig. \ref{fig:qualitative comparison}). This shows the benefit of the curriculum learning approach while also suggesting that our overall framework is robust to catastrophic forgetting. The performance gain is especially significant over the mid and long-term, as the decoder is trained to learn the spatial-temporal movement pattern of a specific subject and can generate the future pose with higher accuracy.

Using a curriculum-learning setup, albeit on a small training sample, allows the framework to capture individual human motion subtleties, as seen by the lower MSE, particularly over the mid and long-term horizons. The approach is particularly useful when there is significant movement over the given horizon, as seen for Subjects 4 and 6 (table \ref{tab:Finetuning all methods}), who have higher MSE loss on the zero-velocity baseline. For Subject 2, there is overall less movement as seen by the zero-velocity MSE loss, and hence the performance gain is not significant over the mid and long-term and even worse over the short-term. Overall, the results are particularly promising as we did not fine-tune the encoder, instead only focusing on the decoder. Further improvement can be attained by fine-tuning the encoder.

\section{Conclusion}
In this work, we present a curriculum learning approach that opens the possibility of continual learning for human motion prediction, especially if the model is especially \emph{deployed in the wild}. Our framework first learns a general representation over diverse training samples before fine-tuning on a target human subject. Our experiments suggest the feasibility of curriculum learning with performance gains over non-curriculum learning approaches. Future work will focus on fine-tuning to activities and domains that were not observed in training in a zero-shot learning setup.

\section{Acknowledgement}
This work is supported by the CCAM Innovation Award. The authors thank Tim Bakker, Tomonari Furukawa, and Qing Chang for their support.

\bibliographystyle{ACM-Reference-Format}


\section{Introduction}\label{sec:intro}
Measurements of cosmic microwave background (CMB) anisotropies have disclosed that the cosmos can be described remarkably by a six-parameter model, with the current energy content dominated by cold dark matter (CDM) and dark energy consistent with a cosmological constant as the standard model of cosmology \cite{wmap, plank}. CMB photons along their path to the earth are subject to various interactions, namely weak gravitational lensing, integrated Sachs-Wolfe and Rees-Sciama effect, and Shapiro effect, which depend directly on the gravitational field.
If Photons cross a time-evolving gravitational potential, the energy gained by photons falling into a potential well cannot cancel out the energy loss as photons climb out of the well \cite{Sachs:1967er, 0801.0642}. For photons passing through a growing potential well, the light is redshifted, and it seems blueshifted if the potential well is decaying. This effect is known as the integrated Sachs-Wolfe (ISW) effect \cite{2007.04968}. This effect can give an important observable that constrains dark energy models. In other words, the ISW effect leads to a new direction in observational cosmology taken for constraining the cosmological models \cite{Cai:2013toa, Nishizawa:2014vga, Manzotti:2014kta, Mostaghel:2018pia, Beck:2018owr, Dong:2020fqt}.
It is known that the simplest way to study the so-called dark energy, responsible for the observed cosmic acceleration, is to postulate the existence of a cosmological constant. Nevertheless, this explanation has some serious problems; above all, the dark energy scale is incredibly smaller than the typical particle physics vacuum energy \cite{Weinberg:1988cp}. This is one of the reasons why a lot of alternative models have been proposed, providing different models of dark energy \cite{Peebles:2002gy, Padmanabhan:2002ji, Copeland:2006wr, Lombriser:2019jia}. These discussions (similar problems for dark matter) lead to a crucial idea of modifying the general relativity (GR), which cannot be well explained in the framework of GR plus $\Lambda$-cold dark matter ($\Lambda$CDM), where $\Lambda$ is the cosmological constant.
There is an important way to add dark energy to the GR, i.e., promoting the cosmological constant to a dynamical field. This gives that the basic theory governing the gravitational dynamics of the Universe may not be GR and could be an alternative gravitational scenario, which can help to realize the dark side better. One way to generalize the Ricci scalar in the Einstein-Hilbert (EH) action is to add some higher-order terms made from Ricci and Riemann tensors. Lovelock's terms will give a well-known model, which cancels the ghost degrees of freedom \cite{Lovelock1, Lovelock2, Lovelock3} if the following four assumptions can be satisfied: ({\it i}) metricity, ({\it ii}) diffeomorphism invariance, ({\it iii}) second-order equations of motion, ({\it iv}) $(3+1)$-dimensional spacetime. One of the non-trivial Lovelock terms is the Gauss-Bonnet (GB) invariant which is a modification of GR, which leads to the GB gravity \cite{Gauss-Bonnet1, Gauss-Bonnet2} and satisfies the former assumptions in higher dimensions, i.e., in $D>4$. Since GB invariant is a total derivative in 4D spacetime, it can not contribute to gravitational dynamics. Interest in the dimensional regularization of gravity with the GB term was motivated by the development of 4D EGB gravity \cite{Glavan:2019inb}. Since due to the presence of an overall factor $(D-4)$, the GB invariant vanishes in 4D, authors in Ref.~\cite{Glavan:2019inb} solved this issue by rescaling the coupling constant $ \alpha\rightarrow \alpha/(D-4) $ in order to produce non-trivial contributions to gravitational dynamics. Consequently, many aspects of this theory have been investigated in black holes and compact objects physics~\cite{Kumar:2020owy, Zhang:2020qew, Aragon:2020qdc, Yang:2020czk, Lin:2020kqe, Wei:2020poh, Konoplya:2020qqh, Konoplya:2020juj, Yang:2020jno, Heydari-Fard:2021ljh, HosseiniMansoori:2020yfj, Wei:2020ght}, and in cosmology~\cite{Li:2020tlo, Aoki:2020iwm, MohseniSadjadi:2020qnm, Narain:2020qhh, Aoki:2020ila, Shahidi:2021rnu}.
In this work, we calculate the ISW effect in the 4D EGB gravity. In this respect, the outline of the work is as follows. In Sec.~\ref{sec:ii}, we review the theoretical framework of the 4D EGB gravity. Then, in Sec.~\ref{sec:iii}, we discuss the cosmological consequences of the 4D EGB gravity concerning the homogeneous and inhomogeneous Universe. Moreover, in Sec.~\ref{sec:iv}, we calculate the ISW effect in the 4D EGB model, and compare it with the corresponding result obtained from the $\Lambda$CDM model. Finally, we scrutinize the results and summarize the findings in Sec.~\ref{sec:v}.
\section{Theoretical Framework}\label{sec:ii}
The Lovelock's theorem ensures that for a gravitational theory with four assumptions mentioned in Sec.~\ref{sec:intro}, the EH term together with the cosmological constant term lead to the unique theory for any generic form of symmetric rank-2 tensor $g_{\mu \nu}$. In other words, this theorem guarantees the uniqueness of Einstein's field equations in 4D spacetime. The well-known EH action with cosmological constant is
\begin{equation}
S_{\rm EH} = \int d^Dx \sqrt{-g}\,\kappa^{2}\,(R - 2\Lambda),
\end{equation}
where $D$ is the dimension of spacetime, $\kappa$ is Einstein's gravitational constant, $\Lambda$ is the cosmological constant, and $R$ is the scalar curvature. Also, $g={\rm det}(g_{\mu \nu})$ is the determinant of the metric tensor matrix. Nevertheless, in higher-dimensional spacetime, i.e., $D>4$, the EH action with second-order field equations is no longer unique. In this regard, a unique action, i.e., the GB action, has been proposed in higher dimensions that include quadratic corrections of the curvature tensors~\cite{Lovelock1, Lovelock2, Lovelock3, Zwiebach1985}. The GB action in $D$-dimensional spacetime is
\begin{equation}
S_{\rm GB} = \int d^Dx \sqrt{-g}\,\alpha\,\mathcal{G},
\end{equation}
where $\alpha$ is the GB coupling parameter. In this action, $\mathcal{G}$ is the GB invariant and has the following form
\begin{equation}
\mathcal{G} = R^{\mu \nu \rho \sigma} R_{\rm \mu \nu \rho \sigma} - 4R^{\mu \nu} R_{\rm \mu \nu} + R^{2},
\end{equation}
where $R_{\mu \nu \rho \sigma}$ is the Riemann curvature tensor, $R_{\mu \nu}$ is the Ricci curvature tensor and $R^{2}$ is the squared of the scalar curvature. In 4-dimensional spacetime, the GB invariant is a total derivative and does not alter the dynamics of the gravitational system~\cite{Felice}. In other words, variation of the GB action with respect to the metric yields
\begin{equation}
\frac{g_{\rm \mu \nu}}{\sqrt{-g}} \frac{\delta S_{\rm GB}}{\delta g_{\rm \mu \nu}} = \frac{(D-4)}{2}\alpha \mathcal{G},
\end{equation}
which is identically zero in $D=4$. While in higher dimensions, i.e., $D>4$, it does not vanish and appears as a dynamical term in the field equations. Recently, 4D EGB has been proposed to rescale GB coupling parameter as $\alpha \rightarrow \alpha / (D-4)$ in a way that it cancels the factor of $(D-4)$ in the GB action~\cite{Glavan:2019inb}. Under this assumption in the limit of $D \rightarrow 4$, one can expect to have a non-negligible dynamical contribution of the GB action in the field equations. Therefore, by combining the EH action, the rescaled GB one, and the action of matter field, $S_{\rm m}$, one can achieve the total action that can intelligently circumvent the Lovelock's theorem and its consequences. With these considerations, the total action can be specified as
\begin{multline} \label{tot_action}
S = \int d^Dx \sqrt{-g} \{ \kappa^2 (R - 2 \Lambda) + \frac{\alpha}{D-4} \mathcal{G} + \mathcal{L}_{\rm m} \},
\end{multline}
where $\mathcal{L}_{\rm m}$ is the Lagrangian of matter field. As can be seen, by rescaling the GB coupling parameter, the factor of $(D-4)$ in the GB term will be eliminated. Hence, variation of the total action with respect to the metric leads to the modified Einstein's field equations, namely
\begin{multline}\label{VarLE$GB$}
2 \kappa^2 (G_{\rm \mu \nu} + \Lambda g_{\rm \mu \nu}) + \frac{\alpha}{d-4} (2 R R_{\rm \mu \nu} - 4 R_{\rm \mu}^{\alpha} R_{\rm \nu \alpha} \\
- 4R^{\alpha \beta} R_{\rm \mu \alpha \nu \beta} - g_{\rm \mu \nu} \mathcal{G}) = T_{\rm \mu \nu},
\end{multline}
where $T_{\rm \mu \nu}$ is the energy-momentum tensor of the matter field and can be calculated by the following defination
\begin{equation}
T_{\rm \mu \nu} =-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{\rm m})}{\delta g^{\rm \mu \nu}}.
\end{equation}

The second parenthesis in Eq.~(\ref{VarLE$GB$}) is known as the GB tensor, which indicates the contribution of the GB action in the equations of motion while considering the rescaled GB coupling parameter. Specifically, the total action introduced in Eq.~(\ref{tot_action}) represents a classical modified theory of gravity that will reduce to the general relativity without any Ostrogradski instability if $\alpha\rightarrow 0$. Moreover, the number of degrees of freedom remains constant within this limit. Furthermore, by applying the Bianchi identity to the GB tensor, one can quickly investigate the conservation of energy-momentum tensor by the following relation
\begin{equation} \label{EM-full}
\nabla_{\mu} T^{\mu \nu} = 0,
\end{equation}
where $\nabla_{\mu}$ represents the covariant derivative. Note that the conservation of energy-momentum tensor is trivial in the absence of non-minimal coupling between the matter and geometry. In the following, we will discuss the cosmological implications of the 4D EGB model.
\section{Cosmology}\label{sec:iii}
\subsection{Homogeneous Universe}
It is known that the Universe is assumed to be isotropic and homogeneous on sufficiently large scales. In this regard, we intend to study the cosmological consequences of the 4D EGB theory by applying the homogeneous and flat spacetime to the geometry of the Universe. For this purpose, one can adopt the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric as
\begin{equation}\label{frw}
ds^2 = a^2(\tau) \eta_{\mu \nu} dx^{\mu} dx^{\nu},
\end{equation}
where $\eta_{\mu \nu}=\rm diag(-1, +1, +1,\cdots)$ is the Minkowski metric and $a(\tau)$ is the scale factor that is as a function of conformal time $\tau$. In D-dimensional spacetime, one can expand relation~(\ref{frw}) as
\begin{equation}\label{exp_frw}
ds^2=-a^{2}(\tau)d\tau^{2} + a^{2}(\tau) ( dx_1^{2} + dx_2^{2} + ... + dx_{D-1}^{2} ).
\end{equation}
Therefore, by inserting relation~(\ref{exp_frw}) into Eq.~(\ref{VarLE$GB$}) the Riemann curvature tensor, the Ricci curvature tensor, and the scalar curvature can be calculated to be
\begin{equation}
\begin{array}{l}
R_{i0i0} = -a \ddot{a}+ \dot{a}^2,	\\ \\
R_{ijij} = \dot{a}^2	, \\ \\
R_{00} = - (D-1)\dfrac{a \ddot{a} - \dot{a}^2}{a^2}, \\ \\
R_{ii} = \dfrac{(D-3) \dot{a}^2 + a \ddot{a}}{a^2}, \\ \\
R = \dfrac{(D-1) [2 a \ddot{a} + (D-4) \dot{a}^2]}{a^4},
\end{array}
\end{equation}
where dot indicates derivative with respect to the conformal time $\tau$. On the other hand, the perfect fluid is the best approximation proposed for the energy-momentum tensor in the homogeneous and isotropic Universe, which has the following form
\begin{equation}
T_{\mu \nu} = (\rho + p) u_{\mu} u_{\nu} - p g_{\mu \nu}.
\end{equation}
In this relation, $\rho$, $p$ and $u^{\mu}$ are the energy density, the pressure and the four-velocity vector field of the cosmological fluid, respectively. Additionally, we have assumed the barotropic equation of state (EoS) for the matter field, i.e. $p=\omega\rho$, where $\omega$ is a dimensionless parameter.
Under the above assumptions, the Friedmann and Raychaudhuri equations for the 4D EGB gravity can be specified as
\begin{equation}
6 \mathcal{H}^2 (\alpha \mathcal{H}^2 + a^2 \kappa^2) = a^4 (2 \kappa^2 \Lambda + \rho) ,
\label{Friedmann-eq}
\end{equation}
and
\begin{equation}
4 \mathcal{\dot{H}} (2 \alpha \mathcal{H}^2 + a^2 \kappa^2) - 2 \mathcal{H}^2 (\alpha \mathcal{H}^2 - a^2 \kappa^2) = a^4 (2 \kappa^2 \Lambda - p),
\end{equation}
where $\mathcal{H} = \dot{a}/a$ is the Hubble parameter. In the Minkowski spacetime, one can easily show that the energy-momentum tensor is conserved, namely
\begin{equation} \label{EMC}
\nabla_{\mu} T^{\mu}_{\nu} = \partial_{\mu} T^{\mu}_{\nu} + \Gamma^{\mu}_{\mu \lambda} T^{\lambda}_{\nu} - \Gamma^{\lambda}_{\mu \nu} T^{\mu}_{\lambda} = 0.
\end{equation}
Hence, by taking $\nu = 0$ in Eq.~(\ref{EMC}) the continuity equation can be obtained as follow
\begin{equation}
\dot{\rho} + 3 \mathcal{H} (\rho + p) = 0.
\end{equation}
On the other hand, the dimensionless parameter of the EoS for the ordinary non-relativistic matter with zero pressure is $\omega_{\rm m} = 0$, whereas it takes the form of $\omega_{\rm r} = 1/3$ for the case of radiation. Therefore, for a Universe that includes non-relativistic matter and radiation, density and pressure can be defined as follows
\begin{equation}
\rho = \rho_{\rm m} + \rho_{\rm r}=6 \kappa^2 H_0^2 (\dfrac{\Omega_{\rm r0}}{a^4} + \dfrac{\Omega_{\rm m0}}{a^3}), \hspace*{0.5cm} p = p_{\rm r} = \frac{\rho_{\rm r}}{3},
\label{rho-p}
\end{equation}
where $H_0$ is the present-time Hubble parameter, and $\rho_{\rm m}$ and $\rho_{\rm r}$ denote the matter and radiation contribution of the energy density, respectively. Also, $\Omega_{\rm r0}$ is the density parameter for radiation, and $ \Omega_{\rm m0}$ the density parameter for matter, and $ \Omega_{\Lambda 0}$ is the density parameter for dark energy at the present-time Universe. By inserting Eq.~(\ref{rho-p}) together with the following dimesionless parameters
\begin{equation}
\begin{array}{l}
E(a)=\dfrac{\mathcal{H}(a)}{H_0},	\kern 2pc	\beta = \dfrac{H_0^2 \alpha}{\kappa^2}, \kern 2pc \Omega_{\Lambda 0} = \dfrac{\Lambda}{3 H_0^2},
\label{parameters}
\end{array}
\end{equation}
one can rewrite Eq.~(\ref{Friedmann-eq}) to be
\begin{equation}
E^2(a) + \dfrac{\beta}{a^2} E^4(a)= a^2 (\dfrac{\Omega_{\rm r0}}{a^4} + \dfrac{\Omega_{\rm m0}}{a^3} + \Omega_{\Lambda 0}).
\label{friedmann-E-eq}
\end{equation}
This relation explains the evolution of the homogeneous Universe while considering the 4D EGB model. Moreover, at the present-time Universe (i.e., at $a=1$), a relation between those density parameters can be extracted from Eq.~(\ref{friedmann-E-eq}), namely
\begin{equation}
\Omega_{\rm r0} + \Omega_{\rm m0} + \Omega_{\rm \Lambda 0} - \beta = 1 .
\end{equation}
We have strong evidence that the Universe is Euclidean, and total density parameter is $\Omega \equiv \Omega_{\rm r}+\Omega_{\rm m}+\Omega_{\Lambda} = 1$ \cite{Dodelson2020}. Correspondingly, the value of the dimensionless coupling parameter, $\beta$, should be very small. Various approaches have been performed to find the observational constraint on $\beta$. These constraints come from gravitational waves, cosmic microwave background, baryon acoustic oscillation, type Ia supernovae, and orbit of Mercury~\cite{2006.15017, 2006.16751, 2103.12358}. By analyzing the data obtained from these phenomena, the strongest constraint on the 4D EGB coupling parameter has been obtained as $\alpha = (2.69 \pm 11.67) \times 10^{48}~ \rm eV^{-2}$, which corresponds to $\beta = (1.2 \pm 5.2) \times 10^{-17}$ \cite{2103.12358}. In this work, we have considered these strong constraints placed on the coupling parameter to calculate the ISW effect in 4D EGB model.

Up to here, we have discussed the evolution of the homogeneous Universe within the context of the 4D EGB model, while a clearer picture of how structures could have been formed and evolved is achieved by considering a Universe that is inhomogeneous. In the next section, we will study the evolution of the inhomogeneous Universe for the 4D EGB model within the context of the perturbation theory and will discuss its cosmological consequences.

\subsection{Inhomogeneous Universe}
There are different kinds of structures in the Universe like galaxies, galaxy clusters, and filaments. Although the Universe seems homogeneous on large scales, one may need some inhomogeneities in small scales to explain the formation conditions of these structures. It should be noted that these inhomogeneities in the cosmic fluid have to be negligible. For this purpose, we have employed a perturbed metric as
\begin{equation}
g_{\rm \mu \nu} = \bar{g}_{\rm \mu \nu} + \delta g_{\rm \mu \nu}, \hspace*{0.5cm} |\delta g_{\rm \mu \nu}|\ll1,
\end{equation}
where $\bar{g}_{\mu \nu}$ is the background metric and $\delta g_{\mu \nu}$ is the spacetime perturbations. Furthermore, in order to maintain the general covariance, a convenient gauge must be considered. To achieve this task, we have chosen the Newtonian gauge that leads to the following line element of spacetime
\begin{equation}
ds^2 = a^2(\tau) [-(1+2\phi) d\tau^2 + (1-2\psi) d \mathbf{x} ^2],
\end{equation}
which is known as perturbed conformal FLRW metric. In this relation, $\phi(\tau, \mathbf{x})$ and $\psi(\tau, \mathbf{x})$ are special combinations of metric perturbations and do not change under a coordinate transformation. Therefore, one can easily calculate the corresponding perturbed connection coefficients as
\begin{equation}
\begin{array}{l}
\Gamma^0_{00} = \mathcal{H} + \dot{\phi},	\kern 2.3pc	\Gamma^0_{0i} = \partial_i \phi	, \\	\\
\Gamma^i_{00} = \partial_j \phi \delta^{ij},	\kern 2pc	\Gamma^i_{j0} = (\mathcal{H} - \dot{\psi}) \delta^i_j ,	\\	\\
\Gamma^0_{ij} = \{-2\mathcal{H} (\phi + \psi) + \mathcal{H} - \dot{\psi}\}\delta_{ij}, \\ \\
\Gamma^i_{jk} = -2 \delta^i_{(j} \partial^{\ }_{k)} \psi + \delta_{jk} \delta^{il} \partial_l \psi ,
\label{Pert-Chr}
\end{array}
\end{equation}
where $\delta_{ij}$ denotes the Kronecker delta. Additionally, one can consider the small perturbations of the stress-energy tensor as
\begin{equation}\label{energy-momentum-full}
T^{\mu}_{\nu} = \bar{T}^{\mu}_{\nu} + \delta T^{\mu}_{\nu},
\end{equation}
where $\bar{T}_{\mu \nu}$ represents the energy-momentum tensor in a homogeneous and isotropic Universe and takes the form of a perfect fluid. Also, $\delta T^{\mu}_{\nu}$ is the perturbation of the energy-momentum tensor that is
\begin{equation}
\delta T^{\mu}_{\nu} = (\delta \rho + \delta p) \bar{u}^{\mu} \bar{u}_{\nu} + (\bar{\rho} + \bar{p})(\delta u^{\mu} \bar{u}_{\nu} + \bar{u}^{\mu} \delta u_{\nu}) - \delta p \delta^{\mu}_{\nu}
\end{equation}
where $\delta \rho$, $\delta p$ and $\delta u^{\mu}$ represent the perturbation of density, pressure, and four-vector velocity, and bar has been employed for the background quantities. Note that the anisotropic part of the above relation has not been considered, because its spatial part can be traceless. As a result, by choosing the anisotropic part of the perturbation of the energy-momentum tensor to be zero, the generality of the problem does not lose.
Now, assuming a non-relativistic and pressureless matter, i.e., $\omega_{\rm m} = p = 0$, as the dominant factor managing the dynamics of the Universe, non-zero components of Eq.~(\ref{energy-momentum-full}) are~\cite{Haghani}
\begin{equation} \label{pert-EM}
\bar{T}^0_0 = \rho,	\kern 2pc	\delta T^0_0 = - \rho \delta_{\rm m} ,	\kern 2pc \delta T^0_i = \rho \partial_i v,
\end{equation}
where $\delta_m=\delta \rho / \rho$ is the density contrast of matter and $v$ is the scalar mode of the velocity perturbation. Therefore, by inserting Eq.~(\ref{pert-EM}) into Eq.~(\ref{EM-full}), one can calculate the perturbed temporal and spatial components of the conserved energy-momentum tensor as
\begin{equation}
\xi = 3 \dot{\psi} - \dot{\delta}_{\rm m},
\label{nu0}
\end{equation}
and
\begin{equation}
\dot{\xi} + \mathcal{H} \xi - k^{2} \phi = 0,
\label{nui}
\end{equation}
where $\xi = \nabla_i \nabla^i v$ is the divergence of the velocity and $k$ is the wave number. Nevertheless, the definition of perturbed metric and energy-momentum tensor lead to the perturbed Einstein's field equations. In this regard, the (00) component of Eq.~(\ref{VarLE$GB$}) can be specified as
\begin{equation}
a^4 \rho \delta_{\rm m} + 4 \mathcal{P}(a) (k^2 \psi + 3 \mathcal{H}^2 \phi +3 \mathcal{H} \dot{\psi}) = 0.
\label{00component}
\end{equation}
Moreover, the spatial off-diagonal components of Eq.~(\ref{VarLE$GB$}) take the following form
\begin{equation}
\mathcal{P}(a) \phi + (\mathcal{Q}(a) - 4 \alpha \dot{\mathcal{H}}) \psi = 0.
\label{off-diagonal}
\end{equation}
and, the remaining components yield
\begin{equation}
\begin{array}{l}
2 \mathcal{H} \mathcal{P}(a) \dot{\phi} + \{ 4 \dot{\mathcal{H}} (2 \mathcal{P}(a) - a^2 \kappa^2) - k^2 \mathcal{P}(a) - 2 \mathcal{H}^2 \mathcal{Q}(a) \}\phi \\
+ 2 \mathcal{P}(a) \ddot{\psi} + 4 \mathcal{H} (a^2 \kappa^2 + 2 \alpha \dot{\mathcal{H}}) \dot{\psi} + k^2 ( 4 \alpha \dot{\mathcal{H}} - \mathcal{Q}(a) ) \psi = 0.
\end{array}
\end{equation}
where we have defined $\mathcal{P}(a)$ and $\mathcal{Q}(a)$ to be
\begin{equation}
\begin{array}{l}
\mathcal{P} (a)= 2 \alpha \mathcal{H}^2 + a^2 \kappa^2 , \\ \\
\mathcal{Q} (a)= 2 \alpha \mathcal{H}^2 - a^2 \kappa^2.
\label{PQ}
\end{array}
\end{equation}
It is worth noting that while we ignore the anisotropic stress tensor, one can expect that by taking the limit of $\alpha \rightarrow 0$, Eq.~(\ref{off-diagonal}) reduces to $\phi = \psi$, similar to the $\Lambda$CDM model.

It is believed that the Fourier modes can be divided into two parts~\cite{Uzan}. Super-Hubble modes, i.e., $k \ll H$, and sub-Hubble modes, i.e., $k \gg H$. Since the growing modes are of interest, we have considered the sub-Hubble limit. With this consideration, by substituting Eq.~(\ref{00component}) into Eq.~(\ref{off-diagonal}) the Poisson equation of the 4D EGB model takes the following form
\begin{equation}
4 \mathcal{P}(a)^2 k^2 \phi = a^4 (\mathcal{Q}(a) - 4 \alpha \dot{\mathcal{H}}) \rho \delta_m.
\label{Poisson-eq}
\end{equation}
Actually, by using this relation one can specify the evolution of the gravitational potential, $\phi$, for the 4D EGB model as a function of redshift. In Fig.~\ref{GP-fig}, we have depicted the redshift-evolution of the gravitational potential for the 4D EGB model, and have compared it with the corresponding result obtained from the $\Lambda$CDM model. As it is clear, the gravitational potential of the 4D EGB model is lower than the corresponding one obtained from the $\Lambda$CDM model for all values of redshifts. Since the ISW effect corresponds to the integrated change of the gravitational potentials, it is expected that this difference leads to a change in the related cosmological parameters to the 4D EGB model.
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-1.pdf}
\caption{The redshift-evolution of dimensionless gravitational potential for the 4D EGB and $\Lambda$CDM models. The dot-dashed (black) line indicates the calculations for the 4D EGB, while the solid (red) line shows the corresponding results for the $\Lambda$CDM model.}
\label{GP-fig}
\end{figure}
In addition, one can express the matter density contrast in terms of decaying modes and growing modes. On the large and small scales, we are only interested in growing modes, because the constructive and destructive interference of these modes lead to clusters or voids, respectively.
Additionally, all modes of interest have entered the Hubble horizon in the late-time Universe. In models with $\Omega_{\rm m0} \neq 1$, and with $\Omega_{\Lambda0}$, all modes experience the same growth factor in a way that the result leads to the same evolution for those modes~\cite{Dodelson}. Moreover, the linear growth factor can be defined as
\begin{equation}
D_{+} (z) = \frac{\delta_{\rm m} (z)}{\delta_{\rm m} (z=0)}.
\label{Growth-factor}
\end{equation}
By substituing Eq.~(\ref{off-diagonal}) into Eq.~(\ref{nu0}) and Eq.~(\ref{nui}), one can achieve a second-order differential equation for the linear growth factor as
\begin{equation}
\begin{array}{l}
D'' + \dfrac{E'}{E} D' \\ \\
\kern 1.4pc + \dfrac{3 \Omega_{m0} \{ 2 \beta E^2 (1+z)^2 + 4 \beta (1+z)^3 E E' -1 \}}{2 (1+z) \{2 \beta E^3 (1+z)^2 + E \}^2}D = 0,
\end{array}
\end{equation}
where prime denotes derivative with respect to the redshift $z$. To solve such an equation in general form, a numerical solution with proper initial conditions is required. For this purpose, we have considered the following initial condition
\begin{equation}
\frac{d D}{d \ln a} |_{z_*} = \gamma D |_{z_*},
\end{equation}
where $\gamma$ and $z_*$ are arbitrary constant and redshift. Note that the observational data must set $ \gamma$ and $\beta$. It should also be mentioned that we have used the best-fit value of parameters $ \gamma$ and $z_{*}$ obtained in Ref.~\cite{Haghani}, and have presented those along with other cosmological parameters in Table~\ref{tab:parameters}.
Fig.~\ref{LGF-fig} shows the normalized linear growth factor of the density contrast, i.e., $D_+ / a$, as a function of redshift for the 4D EGB and $\Lambda$CDM models. The result indicates that for the low redshifts, specifically for $z<1$, the shape of linear growth factor for the 4D EGB model tends to one obtained from the $\Lambda$CDM model in a way that at the present-time Universe, i.e., $z=0$, the growth factors of both models reach the same value. While the deviation of the 4D EGB growth factor from the corresponding result obtained for the $\Lambda$CDM model increases by increasing the redshift.

In Fig.~\ref{LGFN-fig}, we have indicated the ratio of $D_+/a$ in the 4D EGB model to the corresponding value obtained from the $\Lambda$CDM model as a function of redshift with the amplitudes being normalized to unity at the present-time Universe ($z = 0$) for better comparison between the results obtained from both models. Obviously, the value of this ratio is always lower than one, which informs that the normalized growth factor for the 4D EGB model is always lower than the one obtained for the $\Lambda$CDM model.

It can be inferred from Figs.~(\ref{LGF-fig}) and (\ref{LGFN-fig}) that the amplitudes of the growth factor for the 4D EGB model is smaller than that of the $\Lambda$CDM one, which indicates that the 4D EGB model must have primordial perturbations smaller than the $\Lambda$CDM model to give the same quantity of structures at the present-time Universe.
\begin{table}
\centering
\caption{The best-fit parameters $\gamma$ and $z_{*}$, and the other cosmological parameters.}
\label{tab:parameters}
\begin{tabular}{c|c|c|c|c|c}
\hline
\hline
$\Omega_{c} h^2$ & $\Omega_{b} h^2$ & $H_0$ & $\Omega_{m0}$ & $\gamma$& $z_{*}$ \\
\hline
$0.112$ & $0.0226$ & $70$ & $0.275$ & $3.8652$ & $7.1$ \\
\hline
\hline
\end{tabular}
\end{table}
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-2.pdf}
\caption{The redshift-evolution of the linear growth factor normalized by the scale factor for the 4D EGB and $\Lambda$CDM models. The dot-dashed (black) line indicates this dependency for the 4D EGB model, while the solid (red) line shows the corresponding result for the $\Lambda$CDM model.}
\label{LGF-fig}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-3.pdf}
\caption{Relative comparison between the redshift-evolution of the linear growth factor normalized by the scale factor for the 4D EGB model and the corresponding result obtained from the $\Lambda$CDM model with the amplitudes being normalized to unity at the present-time Universe (i.e., at $z = 0$). The solid (red) line has been depicted for comparison.}
\label{LGFN-fig}
\end{figure}
\section{The Integrated Sachs-Wolfe Effect}\label{sec:iv}
As the CMB photons travel from the last scattering surface to us, they move through gravitational potentials generated by the accumulated matters, e.g., the galaxy clusters. It is well known that the photons become blueshifted when they move into potential wells, while the photons will be redshifted as they move out of the gravitational potential wells. Accordingly, these shifts will be accumulated along the line of sight of the observer. In a Universe with no dark energy or curvature, within the linear regime of fluctuations, the gravitational potentials will not change with time if the expansion of the Universe is dominated by a fluid with a constant EoS. Therefore, during the matter-dominated era, the total shift of photons will cancel out. However, any deviation from the matter domination, i.e., the evolution of the EoS, causes the potentials to vary with time. As a result, an integrated change in the energy of traveling photons occurs as they pass through these evolving potentials. This phenomenon has been observed as anisotropy in the CMB temperature, called the ISW effect.
Moreover, the ISW contribution to the total CMB spectrum is only important on large scales (i.e., low CMB multipoles, $l$), where the possibility of extracting information is limited by cosmic variance. However, there are other ways to detect the ISW effect, which the most practical one is the measurement of the correlation between the CMB and the large scale structure \cite{Crittenden and Turok.}.
The ISW effect is important for two main reasons. First, it is a direct probe of dark energy properties, and second, it may be the source of some of the large-scale anomalies in the CMB temperature data \cite{1407.5623}. In the following sections, we mainly focus on the ISW effect and drive the auto-correlation and cross-correlation angular spectrum for the 4D EGB model, and compare those with the corresponding results obtained from the $\Lambda$CDM model.

\subsection{Gravitational Potential}
While a photon propagates toward a cosmological object from the last scattering surface, the gravitational potential evolve in two periods of time. First, in the early Universe, right after the recombination era due to the presence of the radiation, and second, in the late-time Universe due to the presence of dark energy~\cite{Dodelson2020}. As mentioned before, these evolving potentials lead to CMB temperature anisotropies. Consequently, the evolution of the gravitational potential plays a vital role for studying the ISW effect, namely
\begin{equation}
\bigg( \dfrac{\Delta T}{T_{\rm CMB}} \bigg)_{\rm ISW} = \dfrac{2}{c^2} \int^1_{a_{\rm dec}} \frac{\partial \phi}{\partial a} da,
\label{deltaT/T}
\end{equation}
where $T_{\rm CMB} = 2.725~\rm K$ is the CMB temperature and $a_{\rm dec}$ denotes the time at which the photon decoupling has occurred. It must be emphasized that photon decoupling has occurred at $t_{\rm dec}\simeq 378,000$ years after the big bang, or equivalently at redshift $z_{\rm dec}\simeq 1100$. We have also chosen the comoving coordinate systems, because in this system, cosmological fluids are at rest, and galaxies are spatially fixed \cite{Hobson}. In a flat Universe, the comoving distance between an object and us can be calculated as
\begin{equation}
\chi(a) = c \int_a^1 \frac{da}{a \mathcal{H}(a)},
\end{equation}
where $c$ is the velocity of light in vacuum. With this consideration, Eq.~(\ref{deltaT/T}) can be reformed as
\begin{equation}
\bigg( \dfrac{\Delta T}{T_{\rm CMB}} \bigg)_{\rm ISW} = - \dfrac{2}{c^3} \int_0^{\chi_{\rm H}} d \chi\,a^2 H(a) \frac{\partial \phi}{\partial a}.
\label{deltaT/T-n}
\end{equation}
On the other hand, the gravitational potential is determined by the following perturbed Poisson equation
\begin{equation}
- k^2 \phi = \dfrac{a^4 \Big( 4 \alpha \dot{\mathcal{H}} - \mathcal{Q}(a) \Big)}{4 \mathcal{P}^2(a)} \rho_{\rm m} \delta_{\rm m} (a),
\label{pert-Poisson}
\end{equation}
where $\rho_{\rm m}(a) = \Omega_{\rm m}(a) \rho_{\rm crit}$ is the mean background matter density, $\Omega_{\rm m}(a)$ is the matter density parameter and $\rho_{\rm crit}(a)$ is the critical density and can be specified through Eq.~(\ref{Friedmann-eq}). Furthermore, because the cosmological fluids evolve independent from each other, the matter density follows
\begin{equation}
\Omega_{\rm m} (a) = \dfrac{\Omega_{\rm m0} H_{\rm 0}^2}{a \mathcal{H}^2(a)}.
\end{equation}
Under these assumptions and by inserting Eq.~(\ref{Growth-factor}) into Eq.~(\ref{pert-Poisson}), one can achieve the following relation for the Poisson equation
\begin{equation}
\phi(k,a) = \dfrac{-1}{4 k^2} \mathcal{V}(a) \delta (k,a=1),
\label{pert-Poisson-n}
\end{equation}
where $\mathcal{V}(a)$ is a function of scale factor that is defined as
\begin{equation}
\begin{array}{l}
\mathcal{V}(a) \equiv \dfrac{a^4 (4 \alpha \dot{\mathcal{H}} - \mathcal{Q}(a))}{\mathcal{P}^2(a)} \times \\ \\
\kern 1.5pc \Bigg( \dfrac{2 \kappa^2 (3 \mathcal{H}^2 - a^2 \Lambda)}{a^2} + \dfrac{6 \alpha \mathcal{H}^4}{a^4} \Bigg) \Omega_{\rm m}(a) D_+(a).
\end{array}
\end{equation}
Therefore, substituting Eq.~(\ref{pert-Poisson-n}) into Eq.~(\ref{deltaT/T-n}) leads to the following form of the CMB temperature perturbation caused by the ISW effect
\begin{equation}
\bigg( \dfrac{\Delta T}{T_{\rm CMB}} \bigg)_{\rm ISW} = \dfrac{1}{2 c^3} \int_0^{\chi_H} d \chi ~ a \mathcal{H}(a) \dfrac{d \mathcal{V}}{d a} \dfrac{\delta(k,a=1)}{k^2}.
\end{equation}
In addition to the CMB temperature, we need information from the growth of structures. In order to achieve this purpose, one has to calculate the cross-correlation between the CMB and a foreground density field covering the entire extragalactic sky. On the other hand, the large-scale structure surveys, by counting the abundance of galaxy clusters as a function of redshift, are powerful probes of the growth of structure~\cite{0802.2522}. Moreover, the intrinsic angular galaxy fluctuations can be defined as~\cite{0801.0642}
\begin{equation}
\delta_{\rm g} = \int f(z) \delta_{\rm m}(z) dz ,
\end{equation}
where $f(z)=b(z)\,dN/dz$ is the redshift distribution function of the observed samples, $b(z)$ is the redshift-dependent bias that relates the baryonic matter to dark matter, and $dN/dz$ is the redshift distribution of the survey in the comoving distance. Furthermore, the comoving distance is related to the redshift via $d \chi =-(a c/ \mathcal{H})dz$. Under these assumptions, the intrinsic angular galaxy fluctuations takes the following form
\begin{equation}
\delta_{\rm g} = \int \dfrac{\mathcal{H}(a)}{ac}f(z)D_+(a)\delta_{\rm m}(k,a=1) d\chi .
\end{equation}
By utilizing these functions, we are able to calculate the ISW-auto spectrum~\cite{0803.2239, 1504.02416, 1004.3341}
\begin{equation}\label{ctt-spec}
C^{\rm ISW}_{\rm TT} (l) = \int_0^{\chi_{\rm H}}\dfrac{W^2_{\rm T} (\chi)}{\chi^2} \dfrac{H_0^4 P_{\delta \delta} (k=l/\chi)}{k^4} d \chi ,
\end{equation}
the the ISW-cross spectrum
\begin{equation}\label{ctg-spec}
C_{\rm Tg} (l) = \int_0^{\chi_{\rm H}}\dfrac{W_T (\chi) W_{\rm g} (\chi)}{\chi^2} \dfrac{H_0^2 P_{\rm \delta \delta} (k=l/\chi)}{k^2} d \chi ,
\end{equation}
and the observed galaxies-auto spectrum
\begin{equation}\label{cgg-spec}
C_{\rm gg} (l) = \int_0^{\chi_{\rm H}} \dfrac{W^2_{\rm g} (\chi)}{\chi^2} P_{\rm \delta \delta} (k=l/\chi) d \chi,
\end{equation}
where $P_{\rm \delta \delta}$ is the present matter power spectrum that can be spesified by using the definition of the transfer function~\cite{0803.2239}, $l$ is the multipole order, and $W_{\rm T}$ and $W_{\rm g}$ are the ISW and galaxy window functions, respectively.

In Fig.~\ref{mps-fig}, we have shown the results obtained in Ref.~\cite{2103.12358} for the matter power spectrum as a function of wave number for the 4D EGB model while considering two values of $\beta$ and for the $\Lambda$CDM model. As can be seen from the figure, the matter power spectrum for the 4D EGB model is strengthened with respect to the one extracted from the $\Lambda$CDM model. Also, as expected, the matter power spectrum of the 4D EGB model tends to the corresponding result obtained for the $\Lambda$CDM model as the value of $\beta$ decreases. Nevertheless, for $\beta=10^{-15}$, 4D EGB has a larger power spectrum at all scales.

Other functions that have been introduced in Eqs.~(\ref{ctt-spec})-(\ref{cgg-spec}) are the window functions~\cite{0802.0983, 0803.2239}. The ISW window function in the case of a spatially flat Universe with non-clustering dark energy is
\begin{equation}
W_{\rm T} (\chi) = \dfrac{a}{2 c^3}\mathcal{H}(a) \dfrac{d \mathcal{V}}{d a},
\end{equation}
and the galaxy window function is
\begin{equation}
W_{\rm g} (\chi) = \dfrac{\mathcal{H}(a)}{ac} f(z) D_{+} (a).
\label{Wg-eq}
\end{equation}
The ISW effect depends on the redshift distribution of astronomical surveys that requires a lot of precise works, see, e.g., Ref.~\cite{0801.0642}. In this work, we have employed the latest results of some renowned surveys introduced in Refs.~\cite{1004.3341, 0803.2239, 1502.01595}.
Note that in the above formulas, all scale factors and redshifts must be converted into the comoving distance where the current size of the observable Universe is given approximately by the present-day Hubble distance $d_{H,0} = c H_0^{-1}$.
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-4.pdf}
\caption{The matter power spectrum as a function of wave number for the 4D EGB and $\Lambda$CDM models. The dashed (black) and the dot-dashed (black) lines indicate this dependency for the 4D EGB model while considering $\beta = 10^{-15}$ and $\beta = 10^{-17}$, respectively, whereas the solid (red) line shows the result for the $\Lambda$CDM model.}
\label{mps-fig}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-5.pdf}
\caption{The ISW-auto power spectrum as a function of multipole order $l$ for the 4D EGB and $\Lambda$CDM models. The dashed (black) and the dot-dashed (black) lines indicate this dependency for the 4D EGB model while considering $\beta = 10^{-15}$ and $\beta = 10^{-17}$, respectively, whereas the solid (red) line shows the result for the $\Lambda$CDM model.}
\label{CTT-fig}
\end{figure}
\subsection{Large-Scale Structure Surveys}
In order to measure the ISW effect, one has to find the cross-correlation between CMB maps and survey catalogs. In this section, we briefly review three galaxy catalogs and their analytical redshift distribution functions used in our calculation.

In Fig.~\ref{CTT-fig}, we have shown the ISW-auto power spectrum for the 4D EGB model, and have compared it with the result of the $\Lambda$CDM model. Obviously, the amplitude of the ISW-auto power spectrum for the 4D EGB model is higher than the one obtained from the $\Lambda$CDM model. Also, there is a major difference between the amplitude of the ISW-auto power spectrum of the 4D EGB model and the one extracted from the $\Lambda$CDM model. Since the total CMB temperature-auto power spectrum in the $\Lambda$CDM model is consistent with the observational results from the Planck satellite, a large increase in the ISW-auto power spectrum in the 4D EGB model may lead to the trouble. Although it is difficult to directly detect the amplitude of the ISW-auto power spectrum apart from the total CMB temperature anisotropies, this result may indicate the need for stronger observational constraints on the coupling parameter of the 4D EGB model.

\subsubsection{DUNE Survey}
The Dark Universe Explorer (DUNE) is a wide-field space imager that consists of a 1.2 m telescope and is designed to detect both visible and three near-infrared bands. It is optimized for weak gravitational lensing and the ISW effect as a complementary cosmological probe \cite{0802.0983}. For this survey, the redshift distribution function is defined as
\begin{equation}
f_{\rm DUNE}(z)=b_{\rm eff} \left[ \dfrac{z_*}{\alpha_*} \Gamma(\dfrac{3}{\alpha_*}) \right]^{-1} \bigg( \dfrac{z}{z_*} \bigg)^{2} \exp \left[-\left( \dfrac{z}{z_*} \right)^{\alpha_*}\right].
\end{equation}
where $\Gamma(x)$ is the Gamma funcion, and $b_{\rm eff}$, $z_{*}$ and $\alpha_{*}$ are free parameters to be specified.

In Fig.~\ref{CTg-DUNE-fig}, for the DUNE survey, we have shown the ISW-cross power spectrum as a function of multipole order, $l$, while considering the 4D EGB model, and have compared it with the corresponding result of the $\Lambda$CDM model. As can be seen from the figure, the amplitude of the ISW-cross power spectrum for the 4D EGB model is higher than the one obtained from the $\Lambda$CDM model. The results indicate that the 4D EGB model can amplify the ISW-cross power spectrum, which can be considered as a relative advantage of the 4D EGB model. Moreover, it is clear that the deviation from the $\Lambda$CDM model is directly proportional to the value of $\beta$.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-6.pdf}
\caption{The ISW-cross power spectrum as a function of multipole order $l$ for the 4D EGB and $\Lambda$CDM models. The dashed (black) and the dot-dashed (black) lines indicate this dependency for the 4D EGB model while considering $\beta = 10^{-15}$ and $\beta = 10^{-17}$, respectively, whereas the solid (red) line shows the result for the $\Lambda$CDM model. The results have been depicted for the DUNE sample.}
\label{CTg-DUNE-fig}
\end{figure}
\subsubsection{NVSS Survey}
In addition to visible light, radio astronomy is relatively deep and wide in a way that it provides an excellent opportunity to study large-scale structures and especially the ISW effect. The National Radio Astronomy Observatory (NRAO) Very Large Array (VLA) Sky Survey (NVSS) is a $1.4~\rm GHz$ continuum survey that is covering about $82\%$ of the sky, the entire sky north of $-40~\rm deg$ declination~\cite{NVSS, 1203.3277}. Finding a redshift distribution function for NVSS has difficulties. The common method is by cross-correlation against the other samples whose redshift distributions are known. Albeit this method has some flaws, for instance, there are limited data for $z>2.6$ due to the limitation of the range of the other samples \cite{1407.5623, NVSS}. For this survey, the redshift distribution function is defined as
\begin{equation}
f_{\rm NVSS} (z) = b_{\rm eff}\, \dfrac{\alpha_*^{\alpha_* + 1}}{\Gamma (\alpha_*)} \frac{z^{\alpha_*}}{z_*^{\alpha_* +1}}\exp \left( -\frac{\alpha_* z}{z_*} \right).
\end{equation}

In Fig.~\ref{CTg-NVSS-fig}, for the NVSS survey, we have shown the ISW-cross power spectrum as a function of multipole order, $l$, for the 4D EGB model, and have compared it with the corresponding result of the $\Lambda$CDM model. Clearly, the results exhibit that the 4D EGB model can increase the amplitude of the ISW-cross power spectrum in a way that its deviation from the results of $\Lambda$CDM model changes directly with the value of the dimansionless coupling parameter $\beta$.

\begin{table}
\centering
\caption{The best-fit parameters of the redshift distribution of three different surveies, i.e., DUNE, NVSS and SDSS surveys.}
\label{tab-best-fit}
\begin{tabular}{c|c|c|c|c}
\hline
\hline
Survey~~ & ~~$b_{\rm eff}$~~ & ~~$z_*$ & ~~$\alpha_*$~~ & ~~$m$~~ \\
\hline
DUNE~~ & ~~$1.00$~~ & ~~$0.640$~~ & ~~$1.500$~~ & ~~$ - $~~ \\
\hline
NVSS~~ & ~~$1.98$~~ & ~~$0.790$~~ & ~~$1.180$~~ & ~~$ - $~~ \\
\hline
SDSS~~ & ~~$1.00$~~ & ~~$0.113$~~ & ~~$1.197$~~ & ~~$3.457$~~ \\
\hline
\hline
\end{tabular}
\end{table}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-7.pdf}
\caption{The ISW-cross power spectrum as a function of multipole order $l$ for the 4D EGB and $\Lambda$CDM models. The dashed (black) and the dot-dashed (black) lines indicate this dependency for the 4D EGB model while considering $\beta = 10^{-15}$ and $\beta = 10^{-17}$, respectively, whereas the solid (red) line shows the result for the $\Lambda$CDM model. The results have been depicted for the NVSS sample.}
\label{CTg-NVSS-fig}
\end{figure}
\subsubsection{SDSS Survey}
The Sloan Digital Sky Survey (SDSS) is an imaging and spectroscopic survey that collect galaxies, quasars, and stars. There are different subsamples of SDSS. In this work, we have used the redshift distribution of the main photometric SDSS galaxy sample (SDSS-MphG)
that has been introduced in Ref.~\cite{1004.3341}. For this survey, the redshift distribution function is
\begin{equation}
f_{\rm SDSS} (z) = b_{\rm eff}\, \dfrac{\alpha_*}{\Gamma (\frac{m+1}{\alpha_*})} \frac{z^m}{z_*^{m+1}} \exp \left[ - \left(\dfrac{z}{z_*}\right)^{\alpha_*} \right].
\end{equation}
where $m$ is a free parameter to be determined. We have also presented the best-fit values of the redshift distribution of these three surveys in Table~\ref{tab-best-fit}.

In Fig.~\ref{CTg-SDSS-fig}, for the SDSS survey, we have indicated the ISW-cross power spectrum as a function of multipole order for the 4D EGB model, and have compared it with the corresponding result of the $\Lambda$CDM model. As can be seen, the results are in good agreement with corresponding ones obtained from the DUNE and NVSS surveys.

The obtained findings from these three surveys, confirm that the amplitude of the ISW-cross power spectrum for the 4D EGB model is higher than the one obtained for the $\Lambda$CDM model. Moreover, the results indicate that the ISW-cross power spectrum for the 4D EGB model will tend to the corresponding value extracted from the $\Lambda$CDM model if the dimansionless coupling parameter takes the smaller values.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-8.pdf}
\caption{The ISW-cross power spectrum as a function of multipole order $l$ for the 4D EGB and $\Lambda$CDM models. The dashed (black) and the dot-dashed (black) lines indicate this dependency for the 4D EGB model while considering $\beta = 10^{-15}$ and $\beta = 10^{-17}$, respectively, whereas the solid (red) line shows the result for the $\Lambda$CDM model. The results have been depicted for the SDSS-MphG sample.}
\label{CTg-SDSS-fig}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig-9.pdf}
\caption{The redshift distribution for different surveys. The solid (red) line, the dot-dashed (black) line, and the dashed (black) line indicate the corresponding results for the NVSS, DUNE, and SDSS surveys.}
\label{RD-fig}
\end{figure}

In Fig.~\ref{RD-fig}, we have also depicted the normalized redshift distribution of the mentioned surveys as a function of redshift. As can be seen from the figure, the NVSS survey offers the widest redshift coverage with respect to the other surveys.

\section{Conclusions}\label{sec:v}
The EGB gravity was initially proposed as an alternative for dark energy to explain the expansion of the Universe. But contrary to expectations, the theory does not participate in the gravitational dynamics in 4D, because the GB invariant is a total derivative. In fact, the variation of the GB action with respect to the metric is identically zero in 4D due to the presence of an overall factor $(D-4)$. In this regard, Glavan and Lin proposed the novel 4D EGB gravity by rescaling the coupling constant $\alpha \rightarrow \alpha / (D-4)$ to yield a non-trivial contribution to the gravitational dynamics.

In this work, we have calculated the ISW effect in the 4D EGB gravity. For this purpose, we have initially discussed the theoretical framework of the 4D EGB gravity and different cosmological aspects of this theory. First, we have indicated the behavior of the gravitational potentials in the 4D EGB gravity. The results illustrate that the gravitational potential of the 4D EGB model is lower than the corresponding value obtained from the $\Lambda$CDM model for all values of redshift. It has also been confirmed that the deviation of the gravitational potential of the 4D EGB model from the one obtained from the $\Lambda$CDM model changes directly with the redshift.

Moreover, we have calculated the linear growth factor of the 4D EGB model, and have compared it with the corresponding result obtained for the $\Lambda$CDM model. The results demonstrate that for the low redshifts, specifically for $z<1$, the shape of the linear growth factor for the 4D EGB model tends to the one obtained from the $\Lambda$CDM model in a way that at the present-time Universe, i.e., $z=0$, the growth factors of both models reach the same value. Also, the deviation of the 4D EGB growth factor from the corresponding result obtained for the $\Lambda$CDM model changes directly with the redshift. We have also shown the ratio of $D_+/a$ for the 4D EGB model to the one obtained from the $\Lambda$CDM model as a function of redshift for better comparison. The results represent that the values of the growth factor for the 4D EGB model should be smaller than those of the $\Lambda$CDM model as the growth factor is normalized at the present-time Universe, Hence, it indicates that the 4D EGB model must have primordial perturbations smaller than those in the $\Lambda$CDM model to give the same number of structures at the present-time Universe.

Also, we have shown the matter power spectrum with respect to the wave number for the 4D EGB model while considering two values of $\beta$ and for the $\Lambda$CDM model that have been presented in Ref.~\cite{2103.12358}. The results indicate that the matter power spectrum for the 4D EGB model is strengthened with respect to the one extracted from the $\Lambda$CDM model. Also, it has been confirmed that the matter power spectrum of the 4D EGB model tends to the one obtained for the $\Lambda$CDM model as the value of $\beta$ tends to zero.

Furthermore, we have calculated the ISW-auto power spectrum for the 4D EGB model, and have compared it with the corresponding result of the $\Lambda$CDM model. The results indicate that the amplitude of the ISW-auto power spectrum for the 4D EGB model should be higher than the one obtained from the $\Lambda$CDM model. Moreover, the results display a major difference between the amplitude of the ISW-auto power spectrum for the 4D EGB model and the one extracted from the $\Lambda$CDM model. Since the total CMB temperature-auto power spectrum in the $\Lambda$CDM model is consistent with the observational results from the Planck satellite, a large increase in the ISW-auto power spectrum in the 4D EGB model may not be suitable. This result may mean the need to place stronger observational constraints on the coupling parameter of the 4D EGB model.

Additionally, by employing three different surveys, we have calculated the ISW-cross power spectrum as a function of multipole order while considering the 4D EGB model, and have compared those results with the corresponding findings from the $\Lambda$CDM model. The results exhibit that the amplitude of the ISW-cross power spectrum for the 4D EGB model must be higher than the one obtained from the $\Lambda$CDM model. This means that the 4D EGB model can amplify the ISW-cross power spectrum, which can be considered as a relative advantage of the 4D EGB model. Moreover, it has been confirmed that the deviation of the 4D EGB ISW-cross power spectrum from the one obtained for the $\Lambda$CDM model is directly proportional to the value of the dimensionless coupling parameter $\beta$.

Finally, we have demonstrated the normalized redshift distribution of the DUNE, NVSS, and SDSS surveys as a function of redshift. The results indicate that the NVSS survey yields the widest redshift coverage with respect to the other surveys.

\section*{ACKNOWLEDGEMENT}
M.G. and A.B. warmly thank Abdolali Banihashemi for many useful discussions. M.G. gratefully acknowledges Farbod Hassani for invaluable comments on the numerical calculations. The authors would also like to appreciate Deng Wang who provided us with the 4D EGB matter power spectrum data.