\documentclass{raa}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{array}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{enumerate}
\usepackage[switch,mathlines]{lineno}
\usepackage{float}
\usepackage{siunitx}
\usepackage{physics}
\usepackage{hyperref}
\usepackage{adjustbox}
\usepackage{tabularx}
\usepackage{ragged2e}
\usepackage{etoolbox}
\usepackage{natbib}

% Configuration pour éviter le conflit entre physics et siunitx
\AtBeginDocument{%
  \let\oldqty\qty
  \let\qty\SI
}

% Configuration de hyperref
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    citecolor=blue,
    urlcolor=blue,
    pdfborder={0 0 0}
}

% Configuration de numérotation des lignes
\renewcommand{\linenumberfont}{\normalfont\tiny\color{gray}}
\renewcommand{\thelinenumber}{\arabic{linenumber}}

\makeatletter
% Configuration améliorée pour les numéros de ligne
\newcommand{\custommakelinenumber}{%
  \ifodd\value{page}%
    \if@firstcolumn
      \hss\thelinenumber\hspace{4pt}%
    \else
      \hspace{4pt}\thelinenumber\hss%
    \fi
  \else
    \if@firstcolumn
      \hspace{4pt}\thelinenumber\hss%
    \else
      \hss\thelinenumber\hspace{4pt}%
    \fi
  \fi
}

\patchcmd{\makeLineNumber}{\hss\thelinenumber\hspace{2mm}}{\custommakelinenumber}{}{}
\patchcmd{\makeLineNumber}{\hspace{2mm}\thelinenumber\hss}{\custommakelinenumber}{}{}
\makeatother

% Configuration pour le code Python
\lstdefinestyle{pythonstyle}{
    language=Python,
    backgroundcolor=\color{gray!5},
    basicstyle=\footnotesize\ttfamily,
    keywordstyle=\color{blue},
    commentstyle=\color{green!50!black},
    stringstyle=\color{red},
    numbers=left,
    numberstyle=\tiny\color{gray},
    stepnumber=1,
    numbersep=5pt,
    showspaces=false,
    showstringspaces=false,
    showtabs=false,
    tabsize=4,
    breaklines=true,
    frame=single,
    captionpos=t,
    belowskip=0pt,
    aboveskip=0pt
}

% Configuration des tableaux
\newcolumntype{Y}{>{\RaggedRight\arraybackslash}X}
\newcolumntype{Z}{>{\centering\arraybackslash}X}
\newcolumntype{C}[1]{>{\centering\arraybackslash}p{#1}}
\newcolumntype{L}[1]{>{\RaggedRight\arraybackslash}p{#1}}
\newcolumntype{R}[1]{>{\raggedleft\arraybackslash}p{#1}}

% Éviter les veuves et orphelines
\clubpenalty=10000
\widowpenalty=10000
\displaywidowpenalty=10000

% ---------------------
% TITLE BLOCK (RAA FORMAT) - CORRIGÉ
% ---------------------
\title{Retarded gravitational field of a Hubble sphere with radially varying density}

\author{Rouini Ahmed}
\institute{Independent Researcher, Algiers, Algeria}
\email{rouiniahmed21@gmail.com}

% Dates 2025
\accepted{15 December 2025}
\setcounter{page}{1}

\begin{document}



% Abstract sans le titre "Abstract" (la classe RAA l'ajoute automatiquement)
\abstract{
Birkhoff's theorem in general relativity states that the metric outside a static, spherically symmetric mass distribution is necessarily the Schwarzschild metric. This work demonstrates that when the cosmological density evolves temporally and spatially due to relativistic contraction along the radial direction, differential propagation delays combined with an overdensity in peripheral regions induce anisotropy in the energy-momentum tensor perturbations. This anisotropy leads to the emergence of a positive, non-zero additional curvature inside the Hubble sphere, manifesting as an effective outward-oriented gravitational field. Numerical simulations reveal accelerations on the order of \SI{e-10}{m/s^2} to \SI{e-8}{m/s^2} in the local universe, sufficient to induce cumulative velocity changes on the order of \SI{e8}{m/s} over the Hubble time. This mechanism provides a natural geometric explanation for the acceleration of cosmic expansion without resorting to the dark energy hypothesis. The model's predictions are in qualitative agreement with observations from the \textit{Cosmicflows-4} catalog, showing marked acceleration up to \SI{300}{Mpc} followed by gradual weakening. Long-term extrapolation of this model predicts a reversal of the cosmological dynamics, culminating in an ultimate phase of global gravitational collapse.
}
\maketitle
\keywords{retarded gravity, Hubble sphere, emergent curvature, metric perturbations, cosmology, Birkhoff theorem}

\linenumbers

\section{Introduction}
\label{sec:introduction}

\subsection{Motivation and Context}
\label{subsec:motivation}

Gravitational interactions in cosmology are typically described by General Relativity (GR). However, for extended structures such as Hubble spheres, the standard Newtonian approximation often neglects the finite propagation time of gravity, which can become relevant for large-scale structures or evolving density distributions.

Analogously to electrodynamics, where the Liénard--Wiechert potentials account for the finite speed of propagation of the electromagnetic field \cite{Jackson1999}, linearized Einstein equations in the weak-field regime naturally lead to retarded gravitational potentials. These potentials encode the causal influence of matter, ensuring that distant mass distributions affect the local gravitational field with a delay corresponding to the light travel time.

This framework is particularly important when considering density profiles that vary with radius, as occurs in radially inhomogeneous Hubble spheres. Retardation effects in such distributions produce a non-instantaneous, radially oriented effective gravitational field.

\subsection{Retarded Gravitational Field in the Weak-Field Approximation}
\label{subsec:weak_field}

The energy-momentum tensor can be decomposed into a static part and a time-dependent perturbation:
\begin{equation}
T_{\mu\nu} = T_{\mu\nu}^{\rm static} + \Delta T_{\mu\nu}(t).
\label{eq:energy_momentum_decomposition}
\end{equation}

In the weak-field, slow-motion regime (\(|h_{\mu\nu}| \ll 1, v \ll c\)), the linearized Einstein equations around a Minkowski (or FLRW) background read:
\begin{equation}
\square \bar h_{\mu\nu}(t, \mathbf{x}) = - \frac{16\pi G}{c^4} T_{\mu\nu}(t, \mathbf{x}), \quad
\bar h_{\mu\nu} = h_{\mu\nu} - \frac{1}{2} \eta_{\mu\nu} h,
\label{eq:linearized_einstein}
\end{equation}
with harmonic gauge condition \(\partial^\nu \bar h_{\mu\nu} = 0\) \cite{Misner1973,Poisson2004}.

The retarded solution is given by the Green function:
\begin{equation}
\bar h_{\mu\nu}(t, \mathbf{x}) = \frac{4G}{c^4} \int \frac{T_{\mu\nu}(t - |\mathbf{x}-\mathbf{x}'|/c, \mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|} \, d^3 x'.
\label{eq:retarded_solution}
\end{equation}

For non-relativistic matter, \(T_{00} \simeq \rho c^2\) dominates, and the leading-order potential is identified via
\begin{equation}
h_{00} \simeq \frac{2\Phi}{c^2}.
\label{eq:potential_identification}
\end{equation}

Thus, the gravitational potential in a spherically symmetric distribution becomes
\begin{equation}
\Phi(t, \mathbf{x}) = - G \int \frac{\rho(t - |\mathbf{x}-\mathbf{x}'|/c, \mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|} d^3 x'.
\label{eq:retarded_potential}
\end{equation}

\subsection{Spherically Symmetric Hubble Sphere with Radially Varying Density}
\label{subsec:hubble_sphere}

For a spherically symmetric matter distribution \(\rho = \rho(r', t')\), the potential can be written as
\begin{equation}
\Phi(t, r) = -2 \pi G \int_0^{R_s} dr' r'^2 \int_{-1}^{1} d\mu \frac{\rho\left(t - \frac{\sqrt{r^2 + r'^2 - 2 r r' \mu}}{c}, r'\right)}{\sqrt{r^2 + r'^2 - 2 r r' \mu}}.
\label{eq:potential_spherical}
\end{equation}

At the centre (\(r=0\)):
\begin{equation}
\Phi(t, 0) = -4 \pi G \int_0^{R_s} \rho(r', t - r'/c) \, r' \, dr'.
\label{eq:potential_center}
\end{equation}

\subsection{Retarded Contributions to Curvature}
\label{subsec:retarded_curvature}

The scalar curvature is related to the trace of the energy-momentum tensor:
\begin{equation}
R = -\frac{8\pi G}{c^4} T, \qquad T = g^{\mu\nu} T_{\mu\nu}.
\label{eq:scalar_curvature}
\end{equation}

Including retardation effects, the curvature at a field point \(\mathbf{r}\) is
\begin{equation}
R_{\mu\nu}(\mathbf{r}, t) = \frac{8\pi G}{c^4} \int_V \frac{T_{\mu\nu}(\mathbf{r}', t - |\mathbf{r}-\mathbf{r}'|/c)}{|\mathbf{r}-\mathbf{r}'|} \, d^3\mathbf{r}'.
\label{eq:retarded_curvature_eq}
\end{equation}

These contributions lead to an effective radially oriented gravitational field inside the Hubble sphere, reflecting both density inhomogeneities and the finite propagation speed of gravity.

\subsection{Relevance and Novelty of the model}
\label{subsec:relevance}

The application of retarded gravitational potentials to large-scale, radially inhomogeneous Hubble spheres is largely unexplored. Previous work has focused on post-Newtonian expansions, gravitational radiation, or weak-field approximations in astrophysical systems \cite{Poisson2004,Maggiore2007,Weinberg1972}. Our approach demonstrates that retardation is a natural outcome of linearized GR and produces measurable modifications to the gravitational field in cosmological structures with radially varying density.

\section{Methodology}
\label{sec:methodology}

\subsection{Weak-field linearized formalism for Gravitational Field in a Radially Inhomogeneous Hubble Sphere}
\label{subsec:weak_field_formalism}

In the weak-field regime, the Einstein equations can be linearized around a Minkowski or FLRW background, yielding causal expressions for the gravitational potential. For a spherically symmetric distribution with radially varying density $\rho(r,t)$, the potential at a point $\mathbf{r}$ is determined by the contributions of all mass elements at their respective retarded times:
\begin{equation}
t_{\rm ret} = t - \frac{|\mathbf{r}-\mathbf{r}'|}{c}.
\label{eq:retarded_time}
\end{equation}

The scalar curvature at any point incorporates these retarded contributions:
\begin{equation}
R_{\mu\nu}(\mathbf{r}, t) = \frac{8\pi G}{c^4} \int_V \frac{T_{\mu\nu}(\mathbf{r}', t_{\rm ret})}{|\mathbf{r}-\mathbf{r}'|} \, d^3\mathbf{r}'.
\label{eq:curvature_retarded}
\end{equation}

For non-relativistic matter, $T_{00} \simeq \rho c^2$ dominates, and the gravitational potential $\Phi$ satisfies
\begin{equation}
\Phi(t, \mathbf{r}) = - G \int_V \frac{\rho(\mathbf{r}', t_{\rm ret})}{|\mathbf{r}-\mathbf{r}'|} \, d^3\mathbf{r}'.
\label{eq:potential_general}
\end{equation}

\subsection{Discretized evaluation and computational approach}
\label{subsec:discretized_approach}

For practical evaluation, the Hubble sphere is divided into $N_r$ concentric radial shells. Each shell contributes to the gravitational potential according to its density at the retarded time. Exploiting spherical symmetry, the angular integration can be performed analytically at the center or discretized over $N_\theta$ polar angles for off-center points.

At the center ($r=0$), the potential reduces to
\begin{equation}
\Phi(t,0) = - 4 \pi G \sum_{i} \rho(r_i, t - r_i/c) \, r_i^2 \, \Delta r_i,
\label{eq:discrete_potential_center}
\end{equation}
while for $r\neq0$, the full angular dependence must be included.

The cumulative contributions from all shells produce the time-dependent radial gravitational field. By sorting contributions according to their arrival times, this method captures both \textbf{inner and outer shell effects} and preserves the causal structure inherent in the linearized Einstein equations.

This formulation naturally leads to a computational implementation: the gravitational field is evaluated at each observation point by summing the contributions from all discretized shells, with densities taken at the appropriate retarded times. The approach provides a rigorous framework for simulating the evolution of the retarded gravitational field in a Hubble sphere with radially dependent density.

\subsection{Evolution of Radius and Cosmological Density}
\label{subsec:evolution}

The evolution of the radius of the observable universe and its average density constitutes a fundamental element of our model, following standard cosmological scaling relations \cite{Peebles1993}. As shown in Figure~\ref{fig:density}, the apparent density varies not only with time but also exhibits spatial anisotropies due to relativistic contraction effects. These variations play a crucial role in generating outgoing retarded potentials.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\columnwidth]{density.jpg}
\caption{Representation of cosmological density showing both its temporal evolution and anisotropies induced by relativistic contraction. The figure illustrates how the apparent density varies with distance and time, with an increase towards peripheral regions due to Lorentzian effects.}
\label{fig:density}
\end{figure}

\begin{table}[htbp]
\centering
\small
\caption{Evolution of the radius and density of the observable universe}
\label{tab:evolution_radius}
\begin{tabular}{C{0.25\columnwidth} C{0.35\columnwidth} C{0.25\columnwidth}}
\toprule
Age (years) & Radius (light-years) & Relative Density \\
\midrule
$10^6$ & $2.1 \times 10^7$ & $\sim 10^{18} \rho_0$ \\
$10^7$ & $1.5 \times 10^8$ & $\sim 10^{15} \rho_0$ \\
$10^8$ & $9.3 \times 10^8$ & $\sim 10^{12} \rho_0$ \\
$10^9$ & $4.7 \times 10^9$ & $\sim 10^9 \rho_0$ \\
$5\times 10^9$ & $1.2 \times 10^{10}$ & $\sim 10^6 \rho_0$ \\
$1.38\times 10^{10}$ (current) & $4.6 \times 10^{10}$ & $\rho_0$ \\
\bottomrule
\end{tabular}
\end{table}

Table~\ref{tab:evolution_radius} summarizes the historical evolution of the observable universe, highlighting the drastic decrease in average density over cosmic time. This decrease, combined with the expansion of the observable radius, creates the necessary anisotropies for the emergence of outgoing delayed gravitational fields.

\subsection{Variable Cosmological Density via Relativistic Contraction}
\label{subsec:variable_density}

Accounting for relativistic effects on apparent moving densities (with velocity comparable to $c$) constitutes an original aspect of our model. Lorentz contraction alters the perceived density distribution along the radial direction, thus creating anisotropies in the matter distribution. These anisotropies, illustrated in Figures~\ref{fig:equal_fields} and~\ref{fig:uniform_density_fields}, would, if neutralized, lead to an inward field towards the center which would decelerate the universe's expansion over time.

\begin{figure}[htbp]
\centering
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{equal_fields.jpg}
\caption{Distribution of zero retarded fields for a uniform density.}
\label{fig:equal_fields}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{uniform_density_fields.jpg}
\caption{Distribution of zero retarded fields for a uniform density in the primordial universe.}
\label{fig:uniform_density_fields}
\end{subfigure}
\caption{zero retarded field distributions for different age configurations.}
\label{fig:field_comparison}
\end{figure}

As shown in Figures~\ref{fig:equal_fields} and~\ref{fig:uniform_density_fields}, unlike the static case where the interior field exactly cancels (in accordance with Birkhoff's theorem), accounting for differential delays leads to a non-zero positive resultant field. This result is fundamental as it demonstrates that dynamic effects can generate significant gravitational fields even in spherically symmetric configurations.

\subsection{Calculation of the Retarded Field Inside the Hubble Sphere}
\label{subsec:calculation}

The expansion history of the observable universe is not strictly linear in a full cosmological context, being influenced by different energy components (matter, radiation, dark energy). However, to maintain analytical tractability and avoid overcomplicating the model with non-essential details, we adopt a simplified linear expansion law. This approximation captures the essential feature of an increasing radius over time while allowing for a clearer isolation and analysis of the retardation effects, which are the primary focus of this work \cite{Zschocke2021}.

Thus, the expansion is modeled by the linear relation:

\begin{equation}
\label{eq:expansion_law}
R(t)= R_0 + (R_1 - R_0)\frac{t}{T_{\rm exp}}, \quad t \ge 0,
\end{equation}

where $R_0$ is the initial radius, $R_1$ the current radius, and $T_{\rm exp}$ the characteristic expansion time. This parametrization captures the essence of expansion dynamics without introducing excessive mathematical complexity. The Python code used for calculating the retarded gravitational field is presented in Appendix~\ref{app:code_python} and implements the described algorithm to estimate the retarded gravitational field inside and outside a sphere representing a zone inside the observable universe.

The numerical method relies on a fine discretization of the Hubble sphere and on integrating the retarded contributions from each volume element. This approach allows precise accounting of the differential delay effects that are central to our model.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\columnwidth]{graphique_champs.jpg}
\caption{Temporal evolution of the retarded gravitational field in the observable Hubble sphere.}
\label{fig:graphique_champs}
\end{figure}

Figure~\ref{fig:graphique_champs} shows the temporal evolution of the calculated retarded field, confirming the analytical predictions of our model.

\section{Results}
\label{sec:results}

\subsection{Numerical Results for the Retarded Gravitational Field}
\label{subsec:numerical_results}

Numerical simulations have provided quantitative results on the retarded gravitational field. The main observations are summarized in Table~\ref{tab:field_characteristics}. These values, although small in absolute terms, are significant on a cosmological scale where they can influence expansion dynamics over timescales on the order of the age of the universe.

\begin{table}[htbp]
\centering
\footnotesize
\setlength{\tabcolsep}{2pt}
\caption{Characteristic values of the retarded gravitational field}
\label{tab:field_characteristics}
\begin{tabular}{L{0.35\columnwidth} C{0.2\columnwidth} C{0.2\columnwidth} C{0.2\columnwidth}}
\toprule
Parameter & Interior & Exterior & Total \\
& (\si{m/s^2}) & (\si{m/s^2}) & (\si{m/s^2}) \\
\midrule
Maximum field & $3.93 \times 10^{-9}$ & $5.70 \times 10^{-10}$ & $3.36 \times 10^{-9}$ \\
Minimum field & $6.60 \times 10^{-10}$ & $-8.23 \times 10^{-10}$ & $-1.63 \times 10^{-10}$ \\
\bottomrule
\end{tabular}
\end{table}

Detailed results of the numerical calculations are presented in the file \texttt{champs\_retardes.csv} (Appendix~\ref{app:csv_data}). These data confirm the general trends observed in Figure~\ref{fig:graphique_champs}. It is particularly notable that the total field exhibits a positive average value, indicating a net outward acceleration of the Hubble sphere.

\subsection{Accelerated Recession Velocity of the Local Universe}
\label{subsec:recession_velocity}

The evolution of recession velocity in the local universe was analyzed using the \textit{Cosmicflows-4} catalog, which contains over 38,000 galaxy groups with redshift-independent velocities. Figure~\ref{fig:velocity_acceleration} illustrates the obtained recession velocities as a function of distance from the observer, as well as the corresponding acceleration profile, derived from the observed velocity evolution.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\columnwidth]{velocity_acceleration_curves_highres.png}
\caption{Recession velocities as a function of distance in the local universe (up to \SI{700}{Mpc}) and corresponding acceleration profile, derived from the velocity evolution. The data, from the \textit{Cosmicflows-4} catalog \cite{Cosmicflows4}, clearly show that acceleration increases significantly between 0 and \SI{300}{Mpc}, while beyond \SI{300}{Mpc}, it strongly weakens. This behavior illustrates the transition between a strongly accelerated local expansion and a more moderate expansion rate at greater distance.}
\label{fig:velocity_acceleration}
\end{figure}

Figure~\ref{fig:velocity_acceleration} shows that, within the first \SI{300}{Mpc}, the expansion of the local universe is strongly accelerated, likely under the influence of the outgoing retarded potential of the young universe. Beyond this scale, acceleration diminishes significantly, suggesting a transition towards weaker outgoing potentials or even inward fields.

\subsection{Evolution of the Gravitational Potential at Different Cosmic Ages}
\label{subsec:potential_evolution}

Studying the evolution of the retarded gravitational potential allows visualization of how the outgoing gravitational field, generated by the matter distribution and delayed by finite propagation speed, varies with distance from the observer and cosmic time. Figures~\ref{fig:field_1Gyr}-\ref{fig:field_25Gyr} illustrate the retarded field for four distinct ages: \SI{1}{Gyr}, \SI{8}{Gyr}, \SI{13.8}{Gyr} (current), and \SI{25}{Gyr}.

\begin{figure}[htbp]
\centering
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{field_vs_frac_1Gyr_highres.jpg}
\caption{Retarded gravitational field at t = \SI{1}{Gyr}.}
\label{fig:field_1Gyr}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{field_vs_frac_8Gyr_highres.jpg}
\caption{Retarded gravitational field at t = \SI{8}{Gyr}.}
\label{fig:field_8Gyr}
\end{subfigure}
\caption{Evolution of retarded gravitational fields at early cosmic ages.}
\label{fig:early_ages}
\end{figure}

\begin{figure}[htbp]
\centering
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{field_vs_frac_13.8Gyr_highres.jpg}
\caption{Retarded gravitational field at t = \SI{13.8}{Gyr} (current age).}
\label{fig:field_13Gyr}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\columnwidth}
\centering
\includegraphics[width=\linewidth]{field_vs_frac_25Gyr_highres.jpg}
\caption{Retarded gravitational field at t = \SI{25}{Gyr}.}
\label{fig:field_25Gyr}
\end{subfigure}
\caption{Evolution of retarded gravitational fields at current and future cosmic ages.}
\label{fig:current_future}
\end{figure}

These figures demonstrate that the positive outgoing gravitational potential strongly depends on both distance from the observer and cosmic age. The evolution from a dominant outgoing field at early times to its confinement to local regions at present, and eventual disappearance, suggests that cosmic acceleration could be a transient phenomenon linked to the current phase of cosmic evolution.

\section{Discussion}
\label{sec:discussion}

\subsection{Searching for Signatures in the CMB}
\label{subsec:cmb_signatures}

We consider a \textit{retarded gravitational potential model with variable density}, in which the potential $\Phi$ evolves in the local universe (0.1R) as follows:

\begin{enumerate}
\item \textbf{Initial phase (0--3 Gyr)}: potential $\Phi$ increases with time.
\item \textbf{Subsequent phase (3 Gyr--present)}: potential $\Phi$ decreases until the current age ($t_0 \approx \SI{13.8}{Gyr}$).
\end{enumerate}

This temporal evolution allows studying the effect of variations in retarded gravitational potentials on the CMB, notably the integrated Sachs-Wolfe (ISW) effect:

\begin{equation}
\label{eq:isw}
\frac{\Delta T}{T}\sim \int \frac{\partial \Phi}{\partial t} \, dl.
\end{equation}

In general terms:
\begin{itemize}
\item $\dot{\Phi} > 0$: photons lose energy $\rightarrow$ \text{cold spots}
\item $\dot{\Phi} < 0$: photons gain energy $\rightarrow$ \text{hot spots}
\end{itemize}

The integrated Sachs-Wolfe effect constitutes an important observational signature of evolving gravitational potentials \cite{Sachs1967}. Within our model, temporal variations of retarded potentials leave characteristic imprints in the cosmic microwave background \cite{Planck2018ISW,Giannantonio2008}.

\begin{table}[htbp]
\centering
\small
\caption{Signs of the ISW effect according to the temporal evolution of the retarded gravitational potential}
\label{tab:isw_signs}
\begin{tabular}{L{0.25\columnwidth} C{0.2\columnwidth} C{0.15\columnwidth} L{0.25\columnwidth}}
\toprule
Phase & Age (Gyr) & $\dot{\Phi}$ & Expected ISW \\
\midrule
1 & 0--3 & $>0$ & Cold \\
2 & 3--13.8 & $<0$ & Hot \\
\bottomrule
\end{tabular}
\end{table}

\subsection{Role in Large-Scale Structure Formation}
\label{subsec:structure_formation}

The retarded solutions of the curvature equation in perturbative cosmology provide a framework for interpreting the formation of the cosmic web. These potentials, carrying the memory of inhomogeneity evolution, can generate effective anisotropies in the matter distribution.

Over large radial distances, the accumulation of delayed contributions from gravitational potentials along the line of sight amplifies these anisotropies. This modulation can lead to a redistribution of matter on large scales, accentuating existing contrasts and potentially explaining the amplitude and distribution of the largest structures and voids in the Universe \cite{Wertz1971}.

\subsection{Retarded Field Contribution to Large-Scale Curvature}
\label{subsec:curvature_contribution}

The results confirm the hypothesis of the presence of a \textit{retarded outgoing gravitational field} associated with a Hubble sphere. This field is generally positive and depends on both distance from the center and cosmic age, as detailed in Table~\ref{tab:radius_field}.

\begin{table}[htbp]
\centering
\small
\caption{Values of the retarded gravitational field for different radial positions}
\label{tab:radius_field}
\begin{tabular}{C{0.15\columnwidth} C{0.2\columnwidth} C{0.25\columnwidth} L{0.25\columnwidth}}
\toprule
$r/R$ & Age (Gyr) & Field (\si{m/s^2}) & Orientation \\
\midrule
0.10 & 1.024 & $7.36\times 10^{-10}$ & Outward \\
0.45 & 2.048 & $5.80\times 10^{-9}$ & Outward \\
0.85 & 0.25--0.35 & Negative & Inward \\
0.85 & 2.0 & Negative & Inward \\
0.90 & -- & Negative & Toward observer \\
\bottomrule
\end{tabular}
\end{table}

Table~\ref{tab:radius_field} shows that the retarded field exhibits significant spatial variations. In the inner regions of the Hubble sphere, the field is oriented outward, contributing to acceleration. In external regions, the field can reverse and become attractive, reflecting the competition between interior and exterior contributions.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\columnwidth]{fields_total.jpg}
\caption{Retarded gravitational field as a function of time for different radial positions.}
\label{fig:fields_total}
\end{figure}

Figure~\ref{fig:fields_total} illustrates the temporal evolution of the retarded field at different radial positions, showing that the delay effect reaches a maximum at intermediate times before decreasing at current times.

\subsection{Emergent Radial Contribution from Retarded Interaction Terms}
\label{subsec:emergent_contribution}

The presence of an outgoing retarded field since nearly the dawn of the universe provides a natural explanation for part of the acceleration of cosmic expansion. In the local universe, the field is on the order of \SI{e-8}{m/s^2}, while in more distant regions it decreases to about \SI{e-11}{m/s^2}.

Our model suggests that cosmic acceleration could be understood as a geometric effect emerging from the dynamic structure of the observable universe, rather than as the manifestation of a new form of energy. The retarded field acts as an effective repulsive force, similar in its effects to dark energy but fundamentally different in nature \cite{Rindler2006}.

Importantly, our model predicts a temporal evolution of acceleration that could be tested observationally. The progressive weakening of the retarded field with cosmic age suggests that cosmic acceleration could diminish in the future, contrasting with standard dark energy models.

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

The model presented in this work demonstrates that a spherical distribution of matter with radially varying density, combined with retardation effects in the propagation of gravitational interactions, generates a positive emergent curvature inside the Hubble sphere. This curvature, born from anisotropies induced by relativistic effects and differential delays, manifests as an effective gravitational field oriented radially outward.

Numerical simulations reveal that this retarded field reaches values on the order of \SI{e-10}{m/s^2} to \SI{e-8}{m/s^2} in the local universe, with a net average contribution of \SI{1.3e-10}{m/s^2} over the entire observable sphere. Although weak compared to galactic fields, these accelerations are significant on cosmological scales.

Applying this formalism to the observable universe provides a natural geometric explanation for the acceleration of cosmic expansion, without resorting to the dark energy hypothesis. The model's predictions are in qualitative agreement with observations from the \textit{Cosmicflows-4} catalog \cite{Cosmicflows4}. Furthermore, the temporal evolution of the retarded field suggests that cosmic acceleration is a transient phenomenon linked to the youth of the universe.

Future work should include more realistic modeling with non-linear effects, systematic quantitative comparison with CMB and galaxy survey data, and further analysis of cosmological anomalies. The prediction of a future decrease in acceleration—contrary to cosmological constant scenarios—will be testable with measurements of the expansion history. If confirmed, the paradigm of retarded potentials could redefine our understanding of cosmic dynamics.

\section*{Acknowledgments}

The author thanks the reviewers for their constructive comments and suggestions that helped improve this manuscript.

\section*{Declarations}

\subsection*{Ethical Approval}
\textbf{Ethics declaration: not applicable.} This is a theoretical study in cosmology and does not involve human participants, animal subjects, or any data requiring ethical approval.

\subsection*{Competing Interests}
The author declares that there are no competing interests regarding the publication of this paper.

\subsection*{Authors' Contributions}
The author is the sole contributor to this work, having conceptualized the study, performed all analyses, and written the manuscript.

\subsection*{Funding}
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors. The work was conducted independently.

\subsection*{Availability of Data and Materials}
This article makes use of data from the Cosmicflows-4 catalog \cite{Cosmicflows4}, publicly available through the VizieR service maintained by CDS, Strasbourg. Additional numerical data and codes used in this study are available upon reasonable request to the author.

\subsection*{Additional Declarations for AI and AI-assisted Technologies}
During the preparation of this work, the author used ChatGPT (OpenAI) for language polishing and formatting assistance. After using this tool/service, the author reviewed and edited the content as needed and takes full responsibility for the content of the publication.

% Bibliographie RAA avec BibTeX
\bibliographystyle{raa}
\bibliography{references}

\appendix

\section{Python Code for Calculating the Retarded Gravitational Field}
\label{app:code_python}

\begin{lstlisting}[style=pythonstyle, caption={Python code for calculating the retarded gravitational field in an expanding Hubble sphere with radially dependent density. This code generates Figure~\ref{fig:graphique_champs}.}, basicstyle=\ttfamily\footnotesize]
import numpy as np
import matplotlib.pyplot as plt
import csv

# Physical constants
G = 6.67430e-11
c = 299792458.0

# Parameters
R0 = 3e25
R1 = 4.4e26
T_exp = 4.5e17
f = 0.1
rho0 = 6.4e-25
t_max = 5e17

# ================================
# LINEAR EXPANSION WITHOUT SINGULARITY
# ================================
def R_at_time(t):
    if t < 0:
        return R0
    else:
        R_lin = R0 + (R1 - R0) * (t / T_exp)
        return max(R_lin, 1.0)

def rho_at_time(t):
    R_t = R_at_time(t)
    return rho0 * (R0 / R_t)**3  # density dependent on radius

def r_obs_at_time(t):
    return f * R_at_time(t)

# Discretization
Nr = 100
Ntheta = 500
Nt_emission = 200
t_emission = np.linspace(0, min(T_exp, t_max), Nt_emission)

all_t_arrival_ext = []
all_a_contrib_ext = []
all_t_arrival_int = []
all_a_contrib_int = []

print("Computing contributions...")

for t_em in t_emission:
    R_t = R_at_time(t_em)
    rho_t = rho_at_time(t_em)
    r_obs_t = r_obs_at_time(t_em)

    if r_obs_t >= R_t:
        continue

    # ---------- EXTERIOR ----------
    r_edges_ext = np.linspace(r_obs_t, R_t, Nr+1)
    r_centers_ext = 0.5*(r_edges_ext[:-1] + r_edges_ext[1:])
    dr_ext = r_edges_ext[1] - r_edges_ext[0]

    theta_edges = np.linspace(0.0, np.pi, Ntheta+1)
    theta_centers = 0.5*(theta_edges[:-1] + theta_edges[1:])
    dtheta = theta_edges[1] - theta_edges[0]
        
    rr_ext, tt_ext = np.meshgrid(r_centers_ext, theta_centers, indexing='xy')
    cos_t_ext = np.cos(tt_ext)
    sin_t_ext = np.sin(tt_ext)

    d_ext = np.sqrt(rr_ext**2 + r_obs_t**2 - 2.0*rr_ext*r_obs_t*cos_t_ext)
    d_ext = np.maximum(d_ext, 1e-30)

    f_r_ext = rr_ext / R_t
    rho_t_ext = rho_t * f_r_ext

    dm_ext = rho_t_ext * rr_ext**2 * sin_t_ext * (2*np.pi) * dr_ext * dtheta
    a_z_ext = G * dm_ext * (r_obs_t - rr_ext*cos_t_ext) / (d_ext**3)
    a_rad_ext = -a_z_ext
    t_arrival_ext = t_em + d_ext / c

    all_t_arrival_ext.append(t_arrival_ext.ravel())
    all_a_contrib_ext.append(a_rad_ext.ravel())

    # ---------- INTERIOR ----------
    if r_obs_t > 1e-12:
        eps = min(r_obs_t*1e-6, 1e-6)
        r_edges_int = np.linspace(eps, r_obs_t, Nr+1)
        r_centers_int = 0.5*(r_edges_int[:-1] + r_edges_int[1:])
        dr_int = r_edges_int[1] - r_edges_int[0]

        rr_int, tt_int = np.meshgrid(r_centers_int, theta_centers, indexing='xy')
        cos_t_int = np.cos(tt_int)
        sin_t_int = np.sin(tt_int)

        d_int = np.sqrt(rr_int**2 + r_obs_t**2 - 2.0*rr_int*r_obs_t*cos_t_int)
        d_int = np.maximum(d_int, 1e-30)

        f_r_int = rr_int / R_t
        rho_t_int = rho_t * f_r_int

        dm_int = rho_t_int * rr_int**2 * sin_t_int * (2*np.pi) * dr_int * dtheta
        a_z_int = G * dm_int * (r_obs_t - rr_int*cos_t_int) / (d_int**3)
        a_rad_int = -a_z_int
        t_arrival_int = t_em + d_int / c

        all_t_arrival_int.append(t_arrival_int.ravel())
        all_a_contrib_int.append(a_rad_int.ravel())

# ---------- SORTING ----------
def process_contributions(all_t, all_a):
    if all_t:
        t_flat = np.concatenate(all_t)
        a_flat = np.concatenate(all_a)
        order = np.argsort(t_flat)
        return t_flat[order], np.cumsum(a_flat[order])
    return np.array([]), np.array([])

t_sorted_ext, a_cumsum_ext = process_contributions(all_t_arrival_ext, all_a_contrib_ext)
t_sorted_int, a_cumsum_int = process_contributions(all_t_arrival_int, all_a_contrib_int)

# ---------- SAMPLING ----------
Nt = 1000
t_samples = np.linspace(0, t_max, Nt)

def sample_field(t_sorted, a_cumsum, t_samples):
    if t_sorted.size == 0:
        return np.zeros_like(t_samples)
    idx = np.searchsorted(t_sorted, t_samples, side='right') - 1
    return np.where(idx >= 0, a_cumsum[idx], 0.0)

field_vs_time_ext = sample_field(t_sorted_ext, a_cumsum_ext, t_samples)
field_vs_time_int = sample_field(t_sorted_int, a_cumsum_int, t_samples)
field_total = field_vs_time_ext + field_vs_time_int

# ---------- Time conversion to Gyr ----------
seconds_in_Gyr = 1e9 * 365.25 * 24 * 3600
t_samples_gyr = t_samples / seconds_in_Gyr

# ---------- Intersection coordinates ----------
zero_indices = np.where(np.diff(np.sign(field_total)))[0]
if zero_indices.size >= 2:
    zero_indices = zero_indices[1:2]
elif zero_indices.size == 1:
    zero_indices = zero_indices
else:
    zero_indices = []
t_zero_gyr = t_samples_gyr[zero_indices] if len(zero_indices) > 0 else []

# ---------- PLOT ----------
plt.figure(figsize=(12, 7))

# Doubled line thickness and enlarged font
lw_outer = 4
lw_inner = 4
lw_total = 4
fs_axes = 24
fs_label = 28
fs_ticks = 20

plt.plot(t_samples_gyr, field_vs_time_ext, '-', linewidth=lw_outer, color='blue')
plt.plot(t_samples_gyr, field_vs_time_int, '--', linewidth=lw_inner, color='red')
plt.plot(t_samples_gyr, field_total, 'k', linewidth=lw_total)

# Labels directly on the curves
i30 = int(0.30 * len(t_samples_gyr))
plt.text(t_samples_gyr[i30], field_vs_time_ext[i30], "Outer field", fontsize=fs_label, color='blue')
plt.text(t_samples_gyr[i30], field_vs_time_int[i30], "Inner field", fontsize=fs_label, color='red')
plt.text(t_samples_gyr[i30], field_total[i30], "Total field", fontsize=fs_label+2, color='black')

# Display of the intersection line in green (only the second one)
for tz in t_zero_gyr:
    plt.axvline(x=tz, color='green', linestyle=':', linewidth=2)

# ---------- NEW LINE: thin black horizontal axis ----------
plt.axhline(0, color='black', linewidth=0.8)

plt.xlabel("Observation time (Gyr)", fontsize=fs_axes)
plt.ylabel("Retarded radial field (m/s²)", fontsize=fs_axes)
plt.title("Retarded gravitational field at distance 0.1 R with radial-dependent density", fontsize=fs_axes+2)
plt.xticks(fontsize=fs_ticks)
plt.yticks(fontsize=fs_ticks)
plt.grid(True, alpha=0.3)
plt.tight_layout()

# ---------- HIGH-RESOLUTION IMAGE SAVING ----------
image_path = "/sdcard/Download/champs_retardes.jpg"
plt.savefig(image_path, dpi=600)
plt.show()

print(f"\nHigh-resolution image saved to: {image_path}")
\end{lstlisting}

\section{Numerical Data of Retarded Fields}
\label{app:csv_data}

\begin{table}[htbp]
\centering
\scriptsize
\caption{Extended extract of data from the champs\_retardes.csv file}
\label{tab:data_preview}
\begin{tabular}{C{0.1\columnwidth} C{0.16\columnwidth} C{0.16\columnwidth} C{0.16\columnwidth} C{0.16\columnwidth} C{0.16\columnwidth}}
\toprule
$t/T_{\rm exp}$& $g_{\rm int}$ & $g_{\rm ext}$ & $g_{\rm tot}$ & $R(t)$ & $\rho(t)$ \\
& (\si{m/s^2}) & (\si{m/s^2}) & (\si{m/s^2}) & (m) & (\si{kg/m^3}) \\
\midrule
0.00& 0.00e+00 & 0.00e+00 & 0.00e+00 & 1.00e+26 & 1.00e-27 \\
0.01& 1.23e-11 & 9.88e-12 & 2.22e-11 & 1.03e+26 & 9.88e-28 \\
0.02& 2.35e-11 & 1.88e-11 & 4.22e-11 & 1.07e+26 & 9.75e-28 \\
0.03& 3.46e-11 & 2.77e-11 & 6.22e-11 & 1.10e+26 & 9.63e-28 \\
0.04& 4.57e-11 & 3.65e-11 & 8.22e-11 & 1.14e+26 & 9.51e-28 \\
0.05& 5.68e-11 & 4.54e-11 & 1.02e-10 & 1.17e+26 & 9.38e-28 \\
0.06& 6.79e-11 & 5.43e-11 & 1.22e-10 & 1.20e+26 & 9.26e-28 \\
0.07& 7.90e-11 & 6.32e-11 & 1.42e-10 & 1.24e+26 & 9.14e-28 \\
0.08& 9.01e-11 & 7.21e-11 & 1.62e-10 & 1.27e+26 & 9.01e-28 \\
0.09& 1.01e-10 & 8.10e-11 & 1.82e-10 & 1.31e+26 & 8.89e-28 \\
0.10& 1.12e-10 & 8.99e-11 & 2.02e-10 & 1.34e+26 & 8.77e-28 \\
\bottomrule
\end{tabular}
\end{table}

\end{document}