# A viable $f(R)$ gravity model without oscillations in the effective dark
energy

A. Oliveros
Programa de Física, Universidad del Atlántico, Carrera 30 Número 8-49
Puerto Colombia-Atlántico, Colombia<EMAIL_ADDRESS>
###### Abstract

In this study, we propose a reparameterization of a specific viable $f(R)$
gravity model to represent it as a perturbation of the $\Lambda$CDM model. The
$f(R)$ gravity model under consideration includes two parameters, $b$ and $n$,
which control how close the proposed model can be to $\Lambda$CDM, allowing
for arbitrary proximity. Furthermore, it is shown that the Hu-Sawicki (HS)
model is a limiting case of this reparameterized model. Following the existing
literature, we also derive an analytical approximation for the expansion rate
$H(z)$, which shows an excellent agreement between this analytical
approximation and the numerical solution over a wide range of redshifts for
realistic values of the deviation parameter $b$. By appropriately selecting
values for the model parameters, we plot the cosmological parameters
$w_{\rm{DE}}$, $w_{\rm{eff}}$, $\Omega_{\rm{DE}}$, and $H(z)$, as well as the
statefinder quantities $q$, $j$, $s$, and $Om(z)$. We find that their present
values (at $z=0$) are consistent with the observations from Planck 2018 and
the values predicted by the $\Lambda$CDM model. It is important to note that
the examined cosmological and statefinder parameters do not exhibit
significant oscillations of effective dark energy, which could lead to
singular and unphysical solutions at high redshifts. This anomalous behavior
has been avoided here by utilizing the approximate analytical solution for
$H(z)$. Additionally, we conduct a detailed analysis of the evolution of
matter density perturbations within the introduced $f(R)$ gravity model. The
results demonstrate that this viable $f(R)$ gravity model is practically
indistinguishable from the $\Lambda$CDM model at the background level.

Keywords: Modified gravity; Dark energy; $f(R)$ gravity
PACS: 98.80.-k, 04.50. Kd

##  1 Introduction

The discovery of late-time cosmic acceleration in the late 1990s has attracted
a huge attention from the scientific community over the past few decades, and
it has prompted extensive efforts to explain it. The commonly proposed
explanation for this late-time cosmic acceleration is the introduction of a
new energy component in the Universe known as dark energy (DE), characterized
by a negative pressure. However, to date, there is no satisfactory solution to
the DE problem, and its explanation within fundamental theories of physics
remains elusive (for a comprehensive review about this topic see Refs. [1, 2,
3]). Within the wide range of proposals attempting to provide an explanation
for the DE problem, modified gravity theories have positioned themselves as an
excellent alternative. In this context, it is not necessary to consider DE or
new forms of matter to explain the late-time acceleration (see e.g. Refs. [4,
5, 6] for a review). A modified gravity theory that has been successfully
applied in the cosmological and astrophysical context is the $f(R)$ theory of
gravity. In this case, the late-time cosmic acceleration arises as an effect
with geometrical origin [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. However, the selection of a
specific $f(R)$ function is not arbitrary; it must adhere to several
consistency requirements and various constraints that impose conditions for
the cosmological viability of $f(R)$ dark energy models [32]. Despite the
success of $f(R)$ gravity, certain difficulties have arisen during the
analysis of some specific models. For example, some $f(R)$ gravity models
exhibit large frequency oscillations of effective dark energy, which makes
solutions singular and unphysical at a high redshift. Addressing this issue,
Ref. [33] analyzes the oscillatory behavior of the $f(R)$ dark energy around
the phantom divide line, $w_{\rm{DE}}=-1$, both during the matter era and also
in the de Sitter epoch. Similarly, Ref. [34] investigates this characteristic
in viable $f(R)$ gravity models, specifically exploring exponential gravity
and a power form model. In this case, the authors propose improving the models
by introducing a corrective term that effectively stabilizes the oscillations
without compromising the viability properties. In the same way, Ref. [35]
examines a similar behavior within an $f(R)$ Einstein-Gauss-Bonnet gravity
model, while Ref. [36] considers an $f(R)$ Einstein-Gauss-Bonnet model with a
non-minimal coupling between gravity and the kinetic term of a scalar field.
Lastly, Ref. [37] delves into the the late-time dynamics of several $f(R)$
gravity models that exhibit this peculiar behavior.

In the present study, we aim to demonstrate that it is possible to circumvent
this issue within a specific viable $f(R)$ gravity model. Building upon
existing literature, we employ a reparameterization technique to express this
model as a perturbation from the $\Lambda$CDM model. The $f(R)$ gravity model
under consideration includes two parameters, $b$ and $n$, and it is
demonstrated that the closeness of this model to $\Lambda$CDM depends on these
parameters, allowing for arbitrary proximity. Subsequently, we perform an
analytical perturbative expansion in the parameter $b$, yielding an analytic
expression for $H(\Omega_{m0},b;z)$ to all orders in $b$ (for further details,
refer to Ref. [38]). This approach has also been used in other contexts where
an analytical expression for $H(\Omega_{m0},b;z)$ allows for easier
modifications of computational codes, like CLASS. (see Refs. [39] and [40] for
more details). The $f(R)$ gravity model being reparameterized was originally
introduced in Ref. [41], and its late-time cosmological evolution has been
extensively studied in Ref. [42]. Furthermore, by employing the analytic
expression for $H(\Omega_{m0},b;z)$ and appropriately selecting values for the
model parameters, we plot various cosmological parameters, including
$w_{\rm{DE}}$, $w_{\rm{eff}}$, $\Omega_{\rm{DE}}$, and $H(z)$, alongside
statefinder quantities such as $q$, $j$, $s$, and $Om(z)$. Moreover, we
compare their present values (at $z=0$) with observations from Planck 2018 and
the predicted values of the $\Lambda$CDM model. The primary goal is to
demonstrate that the examined cosmological and statefinder parameters do not
exhibit significant oscillations of effective dark energy, which could result
in singular and unphysical solutions at high redshifts. Finally, we conduct a
comprehensive analysis of the evolution of matter density perturbations within
the considered $f(R)$ gravity model.

This paper is organized as follows: In Section 2, we provide a concise
overview of $f(R)$ gravity. Additionally, we carry out the reparameterization
of a particular $f(R)$ gravity model and conduct an analytical perturbative
expansion in the parameter $b$. This expansion leads to an analytic expression
for $H(\Omega_{m0},b;z)$ at all orders in $b$. In Section 3, we plot the
cosmological parameters, namely $w_{\rm{DE}}$, $w_{\rm{eff}}$,
$\Omega_{\rm{DE}}$, and $H(z)$, as well as the statefinder quantities $q$,
$j$, $s$, and $Om(z)$, while considering selected values for the model
parameters. We calculate their present values at $z=0$ and compare them with
observations from Planck 2018, as well as the values predicted by the
$\Lambda$CDM model. Furthermore, we provide a comprehensive analysis of the
evolution of matter density perturbations within the $f(R)$ gravity model
introduced. Our conclusions are presented in Section 4.

##  2 The model

In general, the action for an $f(R)$ gravity model in the presence of matter
components is given by

$S=\int{d^{4}x\sqrt{-g}\left(\frac{f(R)}{2\kappa^{2}}+\mathcal{L}_{M}\right)},$
(2.1)

where $g$ denotes the determinant of the metric tensor $g^{\mu\nu}$,
$\kappa^{2}=8\pi G=1/M_{p}^{2}$, with $G$ being the Newton’s constant and
$M_{p}$ the reduced Planck mass. $\mathcal{L}_{M}$ represents the Lagrangian
density for the matter components (relativistic and non-relativistic perfect
matter fluids). The term $f(R)$ is for now an arbitrary function of the Ricci
scalar $R$. Variation with respect to the metric gives the equation of motion

$f_{R}(R)R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)+(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu})f_{R}(R)=\kappa^{2}T_{\mu\nu}^{(M)},$
(2.2)

where $f_{R}\equiv\frac{df}{dR}$, $\nabla_{\mu}$ is the covariant derivative
associated with the Levi-Civita connection of the metric, and
$\square\equiv\nabla^{\mu}\nabla_{\mu}$. Plus, $T_{\mu\nu}^{(M)}$ is the
matter energy–momentum tensor which is assumed to be a perfect fluid.
Considering the flat Friedman-Robertson-Walker (FRW) metric,

$ds^{2}=-dt^{2}+a(t)^{2}\delta_{ij}dx^{i}dx^{j},$ (2.3)

with $a(t)$ representing the scale factor, the time and spatial components of
Eq. (2.2) are given, respectively, by

$3H^{2}f_{R}=\kappa^{2}(\rho_{\rm{m}}+\rho_{\rm{r}})+\frac{1}{2}(Rf_{R}-f)-3H\dot{f}_{R},$
(2.4)

and

$-2\dot{H}f_{R}=\kappa^{2}\left(\rho_{\rm{m}}+\frac{4}{3}\rho_{\rm{r}}\right)+\ddot{f}_{R}-H\dot{f}_{R},$
(2.5)

where $\rho_{\rm{m}}$ is the matter density and $\rho_{\rm{r}}$ denotes the
density of radiation. The over-dot denotes a derivative with respect to the
cosmic time $t$ and $H\equiv\dot{a}/a$ is the Hubble parameter. If there is no
interaction between non-relativistic matter and radiation, then these
components obey separately the conservation laws

$\dot{\rho}_{\rm{m}}+3H\rho_{\rm{m}}=0,\quad\dot{\rho}_{\rm{r}}+4H\rho_{\rm{r}}=0.$
(2.6)

As usual in the literature, it is possible to rewrite the field equations
(2.4) and (2.5) in the Einstein-Hilbert form:

$3H^{2}=\kappa^{2}\rho,$ (2.7) $-2\dot{H}^{2}=\kappa^{2}(\rho+p),$ (2.8)

where $\rho=\rho_{\rm{m}}+\rho_{\rm{r}}+\rho_{\rm{DE}}$ and
$p=p_{\rm{m}}+p_{\rm{r}}+p_{\rm{DE}}$ correspond to the total effective energy
density and total effective pressure density of the cosmological fluid. In
this case, the dark energy component has a geometric origin, and after a some
manipulation in Eqs. (2.4) and (2.5), we obtain the effective dark energy and
pressure corresponding to $f(R)$-theory given by

$\rho_{\rm{DE}}=\frac{1}{\kappa^{2}}\left[\frac{Rf_{R}-f}{2}+3H^{2}(1-f_{R})-3H\dot{f}_{R}\right],$
(2.9)

and

$p_{\rm{DE}}=\frac{1}{\kappa^{2}}[\ddot{f}_{R}-H\dot{f}_{R}+2\dot{H}(f_{R}-1)-\kappa^{2}\rho_{\rm{DE}}],$
(2.10)

it is easy to show that $\rho_{\rm{DE}}$ and $p_{\rm{DE}}$ defined in this way
satisfy the usual energy conservation equation

$\dot{\rho}_{\rm{DE}}+3H(\rho_{\rm{DE}}+p_{\rm{DE}})=0,$ (2.11)

in this case, we assume that the equation of state parameter for this
effective dark energy satisfies the following relation:

$w_{\rm{DE}}=\frac{p_{\rm{DE}}}{\rho_{\rm{DE}}},$ (2.12)

now, it is well known that the Ricci scalar can be expressed in terms of the
Hubble parameter as

$R=6(2H^{2}+\dot{H}).$ (2.13)

The $f(R)$ gravity model, which plays a central role in this work, is derived
from the following $f(R)$ model:

$f(R)=R-2\,\lambda\,\mu^{2}\,e^{-(\mu^{2}/R)^{n}},$ (2.14)

where $\lambda$ and $n$ are positive real dimensionless parameters, and $\mu$
is a positive real parameter with dimension of $\rm{eV}$. This model was
introduced in Ref. [41], and it behaves very close to $\Lambda$CDM at early
times and satisfy local and cosmological constraints. In Ref. [42] the authors
analyze the late-time evolution of the Universe for this model and perform an
statistical analysis to constrain the free parameters of the model. An
extension of this model (where an $R^{2}$ Starobinsky term is added) able to
explain early time inflation and late time accelerated expansion was studied
in Ref. [43]. Although such an extension, i.e., including an extra $R^{2}$
term in Eq. (2.14), might be of some interest (see [44] for a recent
discussion about this subject), this term does not have relevant contribution
at late-times (dark energy era). Then, in Ref. [45], the author performs a
generalization of this model, introducing a general function of the scalar
curvature in the exponential term. In the literature, other authors have
studied some $f(R)$ gravity models with exponential functions of the scalar
curvature (see for example Refs. [46] and [47]).

Following the ideas worked out in Ref. [38], in this work we reparameterize
the model (2.14) in order to express it as a perturbation deviating from the
$\Lambda$CDM Lagrangian. In this sense, choosing $\mu^{2}=b\Lambda$ and
$\lambda=1/b$, the Eq. (2.14) reduce to

$f(R)=R-2\,\Lambda\,e^{-(b\Lambda/R)^{n}},$ (2.15)

and in this form it is evident that this model can be arbitrarily close to
$\Lambda$CDM, depending on the parameters $b$ and $n$. This statement clearly
illustrates that the reason for successfully passing all the observational
tests is primarily due to small perturbations around the $\Lambda$CDM model.
Furthermore

$\lim_{b\rightarrow 0}f(R)=R-2\Lambda,$ (2.16)
$\lim_{b\rightarrow\infty}f(R)=R,$ (2.17)

From this, we can say that this model converges to $\Lambda$CDM for
$b\rightarrow 0$, while for $b\rightarrow\infty$, the model give rise to a
matter dominated universe. Therefore, it is obvious that this model contain
the cosmological constant $\Lambda$ and in this way, the model must satisfy
the solar system tests. In this sense, in Ref. [48] the authors carried out
cosmological constraints on several well-known $f(R)$ models, but also on a
new class of models that are variants of the Hu-Sawicki model (HS) [49] one of
the form

$f(R)=R-\frac{2\Lambda}{1+b\,y(R,\Lambda)},$ (2.18)

which interpolates between the cosmological constant model and a matter
dominated universe for different values of the parameter $b$. It is worth
saying that Eq. (2.18) represents a specific case of the more general form of
the $f(R)$ function [38, 48, 34]:

$f(R)=R-2\Lambda\tilde{y}(R,b);$ (2.19)

is evident that our model given by Eq. (2.15) has the same form as the above.

It is interesting to see that the $f(R)$ model given by Eq. (2.15) contains
the HS model as a limiting case. In Ref. [38], the authors simplify the HS
model through simple algebraic manipulations and they get

$f(R)=R-\frac{2\Lambda}{1+\left(\frac{b\Lambda}{R}\right)^{n}},$ (2.20)

which can be rewritten as

$f(R)=R-2\Lambda\left[1+\left(\frac{b\Lambda}{R}\right)^{n}\right]^{-1},$
(2.21)

and assuming that $(b\Lambda/R)^{n}\ll 1$, then the HS model reduces to

$f(R)=R-2\Lambda\left[1-\left(\frac{b\Lambda}{R}\right)^{n}\right],$ (2.22)

which corresponds to the result obtained by expanding Eq. (2.15) in a Taylor
expansion up to the first-order approximation.

Another significant advantage of the model given by Eq. (2.15) is that we can
obtain an analytic approximation for the expansion rate $H(z)$. Following the
procedure carried out in Ref. [38], we rewrite Eq. (2.4) in terms of
$N=\ln{a}$

$-f_{R}H^{2}(N)+(\Omega_{m0}e^{-3N}+\Omega_{r0}e^{-4N})H_{0}^{2}+\frac{1}{6}(Rf_{R}-f)=f_{RR}H^{2}(N)R^{\prime}(N),$
(2.23)

here prime denotes differentiation with respect to $N$. Additionally, Eq.
(2.13) is given by

$R(N)=6\left[2H^{2}(N)+\frac{1}{2}\frac{dH^{2}(N)}{dN}\right],$ (2.24)

since the model under study here approach to $\Lambda$CDM as $b\rightarrow 0$,
one can express the solution to Eq. (2.23), $H^{2}(N)$, as a Taylor expansion
in the deviation parameter, $b$, as follows

$H^{2}(N)=H^{2}_{\Lambda}(N)+\sum_{i=1}^{M}b^{i}\delta H_{i}^{2}(N),$ (2.25)

where

$\frac{H^{2}_{\Lambda}(N)}{H_{0}^{2}}=\Omega_{m0}\,e^{-3N}+\Omega_{r0}\,e^{-4N}+(1-\Omega_{m0}-\Omega_{r0})=E_{\Lambda}^{2}(z),$
(2.26)

further, as has been demonstrate in Ref.[38], we can consider only two terms
in the above series expansion, in this way Eq. (2.25) reduces to

$H^{2}(N)\approx H^{2}_{\Lambda}(N)+b\,\delta H_{1}^{2}(N)+b^{2}\delta
H_{2}^{2}(N),$ (2.27)

and considering by simplicity $n=1$ in our model, then, $\delta H_{1}^{2}(N)$
and $\delta H_{2}^{2}(N)$, are given by

$\frac{\delta
H_{1}^{2}(N)}{H_{0}^{2}}=\frac{-2H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})^{2}\left[12(H_{\Lambda}^{2}(N))^{2}+\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}+H^{2}_{\Lambda}(N)\left(15\frac{dH_{\Lambda}^{2}(N)}{dN}+2\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\right)\right]}{\left[4H^{2}_{\Lambda}(N)+\frac{dH_{\Lambda}^{2}(N)}{dN}\right]^{3}},$
(2.28) $\displaystyle\frac{\delta H_{2}^{2}(N)}{H_{0}^{2}}=$
$\displaystyle\bigg{[}H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})^{3}\bigg{(}8192(H_{\Lambda}^{2}(N))^{6}+3\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{6}+8H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{4}$
(2.29)
$\displaystyle\times\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}+1024(H_{\Lambda}^{2}(N))^{5}\bigg{(}25\frac{dH_{\Lambda}^{2}(N)}{dN}+3\bigg{(}\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-6(1-\Omega_{m0}-\Omega_{r0})\bigg{)}\bigg{)}$
$\displaystyle+4H_{\Lambda}^{2}(N)\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}\bigg{(}29\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{3}-24(1-\Omega_{m0}-\Omega_{r0})\left(\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\right)^{2}$
$\displaystyle+3\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}\left(\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-198(1-\Omega_{m0}-\Omega_{r0})\right)-8H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\frac{dH_{\Lambda}^{2}(N)}{dN}$
$\displaystyle\times\left(22\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}\right)\bigg{)}+256(H_{\Lambda}^{2}(N))^{4}\bigg{(}83\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}+12\frac{dH_{\Lambda}^{2}(N)}{dN}$
$\displaystyle\times\left(\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-24H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\right)-2H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\bigg{(}15\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-6\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}$
$\displaystyle-\frac{d^{4}H_{\Lambda}^{2}(N)}{dN^{4}}\bigg{)}\bigg{)}+128(H_{\Lambda}^{2}(N))^{3}\bigg{(}61\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{3}+\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}\bigg{(}9\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}$
$\displaystyle-510H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\bigg{)}-3H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\bigg{(}19\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}+4\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}\bigg{)}$
$\displaystyle-2H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\frac{dH_{\Lambda}^{2}(N)}{dN}\bigg{(}147\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}+16\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}-\frac{d^{4}H_{\Lambda}^{2}(N)}{dN^{4}}\bigg{)}\bigg{)}$
$\displaystyle-32(H_{\Lambda}^{2}(N))^{2}\bigg{(}-44\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{4}+\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{3}\bigg{(}-672H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})$
$\displaystyle-6\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\bigg{)}-21H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\left(\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\right)^{3}-3H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\frac{dH_{\Lambda}^{2}(N)}{dN}$
$\displaystyle\times\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}\left(61\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-4\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}\right)-H_{0}^{2}(1-\Omega_{m0}-\Omega_{r0})\left(\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{2}$
$\displaystyle\times\left(637\frac{d^{2}H_{\Lambda}^{2}(N)}{dN^{2}}-34\frac{d^{3}H_{\Lambda}^{2}(N)}{dN^{3}}+\frac{d^{4}H_{\Lambda}^{2}(N)}{dN^{4}}\right)\bigg{)}\bigg{)}\bigg{]}/\left(2\left(4H_{\Lambda}^{2}(N)+\frac{dH_{\Lambda}^{2}(N)}{dN}\right)^{8}\right),$

where, in order to simplify the calculations, we have expressed Eqs. (2.28)
and (2.29) in terms of $H_{\Lambda}^{2}(N)$ and its derivatives. It is
straightforward to show that Eq. (2.28) is equivalent to the one obtained in
Ref. [38] for the HS model. This similarity is not surprising, as demonstrated
earlier, our model closely resembles the HS model at first order in $b$.
Regarding Eq. (2.29), it exhibits the same dependence on $H_{\Lambda}^{2}(N)$
and its derivatives as the equation derived in Ref. [38]. However, in our
case, there are certain numerical coefficients that differ.

Replacing Eqs. (2.28) and (2.29) in Eq. (2.27), and using Eq. (2.26), we
obtain an approximate solution for the Hubble parameter $H(z)$:

$\displaystyle E^{2}(z)=$
$\displaystyle\frac{H^{2}(z)}{H_{0}^{2}}=1-\Omega_{m0}+(1+z)^{3}\Omega_{m0}$
(2.30)
$\displaystyle+\frac{6b(\Omega_{m0}-1)^{2}\left(-4+\Omega_{m0}(9-3\Omega_{m0}+z(3+z(z+3))(1+(3+2z(3+z(z+3)))\Omega_{m0}))\right)}{(4+(-3+z(3+z(z+3)))\Omega_{m0})^{3}}$
$\displaystyle+\frac{b^{2}(\Omega_{m0}-1)^{3}}{(1+z)^{24}\left(\frac{4(1-\Omega_{m0})}{(1+z)^{3}}+\Omega_{m0}\right)^{8}}\bigg{[}5120(\Omega_{m0}-1)^{6}-9216(1+z)^{3}(\Omega_{m0}-1)^{5}\Omega_{m0}$
$\displaystyle-30144(1+z)^{6}(\Omega_{m0}-1)^{4}\Omega_{m0}^{2}+31424(1+z)^{9}(\Omega_{m0}-1)^{3}\Omega_{m0}^{3}-9468(1+z)^{12}$
$\displaystyle\times(\Omega_{m0}-1)^{2}\Omega_{m0}^{4}+4344(1+z)^{15}(\Omega_{m0}-1)\Omega_{m0}^{5}+\frac{37}{2}(1+z)^{18}\Omega_{m0}^{6}\bigg{]},$

for simplicity, we have assumed $\Omega_{r0}=0$ and made the substitution
$N=-\ln{(1+z)}$. In Ref. [50], the authors derived a comparable expression for
the Hubble rate $H(z)$ in the HS model. However, in our study, we have found
that the numerical coefficients differ from those obtained in their work. In
the left panel of Fig. 1 we can see the evolution of $H(z)$ vs. $z$ taking
into account the approximate analytical solution given by Eq. (2.30) and the
numerical solution for $H(z)$ obtained using the formalism described in Ref.
[42]. In the right panel of Fig. 1, the comparison between them is made. In
this case, the comparison is performed using the quantity $\Delta H(z)$, which
is defined by

$\Delta H(z)=\left[\frac{|H_{A}(z)-H_{N}(z)|}{H_{N}(z)}\right]\times 100,$
(2.31)

where $\rm{A}=\rm{Analytical}$ and $\rm{N}=\rm{Numerical}$. From above, we can
say that the difference between these solutions is $\Delta H(z)\sim 0.01\,\%$,
which means that the approximation behaves very well in this scenario.

Figure 1: (_left_) Plot for $H(z)$ vs. $z$ taking into account the numerical
and analytical solutions (in units of $\rm{Km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$).
(_right_) Comparison between the approximate analytical solution given by Eq.
(2.30) and the numerical solution for $H(z)$ obtained using the formalism
described in Ref. [42]. The values for the model parameters used here are:
$\Omega_{m0}=0.3147$, $H_{0}=67.4\,\rm{Km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$ y
$b=0.01$.

##  3 Cosmological analysis in late-time

In this section, we will examine the aforementioned results in order to
conduct a thorough analysis about the feasibility of the model presented in
Eq. (2.15). Additionally, we demonstrate the elimination of oscillations in
the effective dark energy within this particular scenario. In this regard, we
initiate the analysis by examining the behaviour of the cosmological
parameters, namely $w_{\rm{DE}}$, $w_{\rm{eff}}$, and $\Omega_{DE}$, as well
as the statefinder quantities $q$, $j$, $s$, and $Om(z)$ in the late-time
regime. Let’s remember that the cosmological parameters, $w_{\rm{DE}}$,
$w_{\rm{eff}}$, and $\Omega_{\rm{DE}}$ in terms of the redshift $z$, are given
by:

$w_{\rm{DE}}=-1+\frac{1}{3}(1+z)\frac{(\rho_{\rm{DE}}(z))^{\prime}}{\rho_{\rm{DE}}(z)},$
(3.1) $w_{\rm{eff}}=-1+\frac{1}{3}(1+z)\frac{(E^{2}(z))^{\prime}}{E^{2}(z)},$
(3.2) $\Omega_{\rm{DE}}=\frac{1-\Omega_{m0}}{E^{2}(z)};$ (3.3)

similarly, the statefinder quantities $q$, $j$, $s$ and $Om(z)$ in terms of
the redshift $z$, are given by:

$q=-1+\frac{1}{2}(1+z)\frac{(E^{2}(z))^{\prime}}{E^{2}(z)},$ (3.4)
$j=-2-3q+\frac{(1+z)^{2}(E^{2}(z))^{\prime\prime}+(1+z)(E^{2}(z))^{\prime}}{2E^{2}(z)},$
(3.5) $s=\frac{j-1}{3\left(q-\frac{1}{2}\right)},$ (3.6)
$Om(z)=\frac{\frac{E^{2}(z)}{E^{2}(0)}-1}{(1+z)^{3}-1},$ (3.7)

where prime denotes differentiation with respect to $z$. Using Eqs. (2.9) and
(2.30), we can plot the above expressions in terms of the redshift $z$. In the
left panel of Fig. 2 we show the cosmological evolution of $w_{DE}$ as a
function of $z$ and we can see that for $z>8$, $w_{DE}$ is very close to the
$\Lambda$CDM model prediction (i.e. $w_{\Lambda}=-1$). For $0\leq z\leq 8$,
$w_{\rm{DE}}$ present a small variation about -1 and the equation of state
crosses $w_{\rm{DE}}=-1$ at approximately the redshift $z\approx 1.59$. At
early times, we have $1+w_{\rm{DE}}<0$ thus violating the strong energy
conditions (SEC). A similar behavior was found for the HS model in Ref. [39].
In the right panel of Fig. 2, we depict the evolution of $\Omega_{\rm{DE}}$ in
terms of redshift $z$. From this, it is clear that for high redshifts,
$\Omega_{\rm{DE}}$ lays to zero (matter dominated era) and for the distant
future ($z\rightarrow-1$), $\Omega_{\rm{DE}}\rightarrow 1$ (de Sitter phase).

Figure 2: (_left_) Evolution of $w_{\rm{DE}}$ vs. $z$ and (_right_)
$\Omega_{\rm{DE}}$ vs. $z$. In this case, we have used the values for the
model parameters as those used in Fig. 1

Figure 3: (_left_) Evolution of $w_{\rm{eff}}$ as a function of $z$ and $b$.
(_right_) $w_{\rm{eff}}$ vs. $z$.

Now, we depict the cosmological evolution of $w_{\rm{eff}}$ as a function of
the redshift $z$ and the $b$ parameter in Fig. 3. We can see that,
independently on the choice of $b$, $w_{\rm{eff}}$ starts from zero in the
matter dominated era and asymptotically approaches -1 without any appreciable
deviation. In the right panel of Fig. 3, we represent $w_{\rm{eff}}$ in terms
of $z$ considering a fixed value for $b$, in this case $b=0.01$, and is clear
that in this scenario, we do not have the crossing of the phantom divide,
which has been found for other $f(R)$ models due to the oscillatory behaviour
of dark energy [34]. In the Figs. 4 and 5, the statefinder parameters $q$,
$j$, $s$, and $Om(z)$ present a behaviour very close to that obtained using
the $\Lambda$CDM model.

Figure 4: (_left_) Evolution of $q$ vs. $z$ and (_right_) $j$ vs. $z$.

Figure 5: (_left_) Evolution of $s$ vs. $z$ and (_right_) $Om$ vs. $z$.

In order to compare our model against the $\Lambda$CDM model, we use the
dimensionless Hubble parameter, $E(z)=H(z)/H_{0}$, and calculating

$\Delta E(z)=\left[\frac{E(z)}{E(z)_{\Lambda\rm{CDM}}}-1\right]\times 100,$
(3.8)

which, by definition, is zero for the $\Lambda$CDM model (red-dashed line in
the right panel of Fig. 6). Observing Fig. 6 (left panel), it is evident that
certain choices of $b$ result in an evolution of $H(z)$ that closely resembles
the $\Lambda$CDM model. To quantify this similarity, the left panel of Fig. 6
clearly demonstrates that as $b$ decreases, $H(z)$ progressively converges
towards the $\Lambda$CDM model. The largest deviation from $\Lambda$CDM occurs
around $z\simeq 0.90$ and for $b=0.1$, being $\Delta E(z\simeq 0.90)\simeq
0.39\%$. Additionally, regardless of the value of $b$, the approximated
analytical solution for $H(z)$ tends to the $\Lambda$CDM model at higher
redshifts.

Figure 6: (_left_) Evolution of the Hubble parameter, $H(z)$, (in units of
$\rm{km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$) using the approximated analytical
solution given by Eq. (2.30) and the $\Lambda$CDM model as functions of the
red shift. (_right_) Comparison of our model against $\Lambda$CDM using the
dimensionless Hubble parameter through Eq. (3.8).

Table 1 provides a summary of the current ($z=0$) values obtained for the
analyzed cosmological and statefinder parameters. We also include, for
comparison purposes, the predicted values by the $\Lambda$CDM model and the
available measurements for some of the parameters.

Table 1: Current values obtained for the cosmological and statefinder parameters utilizing the approximate analytical solution for $H(z)$ provided by Eq. (2.30), along with fixed values for the free parameters of the $f(R)$ gravity model. The values for the $\Lambda$CDM model are obtained using Eq. (2.26) in Eqs. (3.1)-(3.7). The units for $H_{0}$ are $\rm{km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$.  Parameter | $f(R)$ | Planck 2018 or SNe Ia |  $\Lambda$CDM
---|---|---|---
$q_{0}$ | -0.526 |  $-0.53^{+0.17}_{-0.13}$ (SNe Ia) | -0.528
$j_{0}$ | 0.999 | - | 1.0002
$s_{0}$ | 0.0002 | - | 0.002
$Om_{0}$ | 0.316 | - | 0.315
$w_{\rm{DE}0}$ | -0.998 | $-1.03\pm 0.03$ | -1
$\Omega_{\rm{DE}0}$ | 0.686 | $0.6847\pm 0.0073$ | 0.6852
$w_{\rm{eff}0}$ | -0.684 | - | -0.685
$H_{0}$ | 67.351 | $67.4\pm 0.5$ | 67.4

It is important to note that the cosmological and statefinder parameters
examined earlier do not exhibit significant frequency oscillations of
effective dark energy. Such oscillations would lead to singular and unphysical
solutions at high redshifts. In our study, we have avoided this anomalous
behavior by utilizing the approximate analytical solution for $H(z)$ provided
by Eq. (2.30). It is worth noting that if we were to use the numerical
solution for this model, the oscillations of effective dark energy would be
unavoidable, as depicted in Fig. 7, for instance. In order to compare see Fig.
2 (left panel) and Fig. 4 (right panel). For simplicity, we have omitted the
other plots in this discussion, as they exhibit similar behavior to those
presented in Ref. [42].

Figure 7: (_left_) Evolution of $w_{\rm{DE}}$ vs. $z$ and (_right_) $j$ vs.
$z$. In this case, the numerical solution employed in this study was obtained
by utilizing the same initial conditions and model parameter values as those
utilized in Ref. [42], but taking into account the model given by Eq. (2.15),
along with $b=0.01$ and $n=1$.

###  3.1 Growth rate of the matter density perturbations

To complement the ongoing study, we proceed in this section to analyze the
evolution of matter density perturbations for the $f(R)$ gravity model
introduced in Eq. (2.15). The equation governing the evolution of matter
density perturbations for the $f(R)$ gravity has been derived previously in
the literature, considering the subhorizon approximation ($k^{2}/a^{2}\gg
H^{2}$) [51, 52, 53]. The gauge-invariant matter density perturbation defined
by

$\delta_{\rm{m}}\equiv\frac{\delta\rho_{\rm{m}}}{\rho_{\rm{m}}},$ (3.9)

satisfies the following equation:

$\ddot{\delta}_{\rm{m}}+2H\dot{\delta}_{\rm{m}}-4\pi
G_{\rm{eff}}(a,k)\rho_{\rm{m}}\delta_{\rm{m}}=0,$ (3.10)

with $k$ being the comoving wavenumber and $G_{\rm{eff}}(a,k)$ being the
effective gravitational “constant” given by

$G_{\rm{eff}}(a,k)=\frac{G_{\rm{N}}}{f_{R}}\left[1+\frac{(k^{2}/a^{2})(f_{RR}/f_{R})}{1+3(k^{2}/a^{2})(f_{RR}/f_{R})}\right],$
(3.11)

in the context of the general theory of relativity, it is important to note
that the evolution of matter density perturbations is independent of the
comoving wavenumber $k$. However, in $f(R)$ gravity, this dependence arises in
the effective gravitational constant, as evident in Eq. (3.11). Consequently,
the effective gravitational constant generally exhibits scale dependence. In
Fig. 8 (left panel), we illustrate the evolution of normalized effective
gravitational constant as a function of the redshift $z$ and the scale
dependence on the comoving wavenumber $k$ for our model. In the right panel,
we observe its evolution while considering specific fixed values of $k$. In
both figures, it is evident that for approximately $z>4$, the normalized
effective gravitational constant tends to 1. This means that for high
redshifts, $G_{\rm{eff}}(z,k)$ approaches $G_{\rm{N}}$, which aligns with the
prediction of $\Lambda$CDM. Additionally, for $0\leq z<4$, there is a slight
variation in $G_{\rm{eff}}/G_{\rm{N}}$ compared to the $\Lambda$CDM
prediction.

Figure 8: (_left_) Evolution of $G_{\rm{eff}}/G_{\rm{N}}$ as a function of $z$
and the scale dependence on $k$ and (_right_) $G_{\rm{eff}}/G_{\rm{N}}$ vs $z$
using some fixed values of $k$. Note that in both cases we have replaced $a$
by $1/(1+z)$ in Eq. (3.11)

Furthermore, these results demonstrate a high level of insensitivity to
changes in the comoving wavenumber $k$. However, it should be noted that this
does not imply that $G_{\rm{eff}}(z,k)$ is completely independent of $k$. In
order to illustrate the above, we expand Eq. (3.11) in Taylor series up to
second order in $b$:

$G_{\rm{eff}}(z,k)\approx G_{\rm{N}}+b\,\delta
G_{\rm{eff1}}(z,k)+b^{2}\,\delta G_{\rm{eff2}}(z,k),$ (3.12)

where $\delta G_{\rm{eff1}}(z,k)$ is given by

$\delta
G_{\rm{eff1}}(z,k)=\frac{2G_{\rm{N}}(\Omega_{m0}-1)^{2}[2k^{2}(1+z)^{2}+3H_{0}^{2}(4+(-3+z(3+z(3+z)))\Omega_{m0}]}{3H_{0}^{2}[4+(-3+z(3+z(3+z)))\Omega_{m0}]^{3}},$
(3.13)

and $\delta G_{\rm{eff2}}(z,k)$ is

$\displaystyle\delta G_{\rm{eff2}}(z,k)=$ $\displaystyle
2G_{\rm{N}}(\Omega_{m0}-1)^{3}\bigg{[}-\frac{1}{H_{0}^{4}}8k^{4}(1+z)^{4}(\Omega_{m0}-1)(4+(-3+z(3+z(3+z)))\Omega_{m0})^{2}$
(3.14)
$\displaystyle+3(4+(-3+z(3+z(3+z)))\Omega_{m0})(-2+(3+z(3+z(3+z)))\Omega_{m0})(320$
$\displaystyle+\Omega_{m0}(-592+368z(3+z(3+z))+8(51+23z(3+z(3+z))(-2+z(3$
$\displaystyle+z(3+z))))\Omega_{m0}+(-135+z(3+z(3+z))(3+z(3+z(3+z))(-181+z$
$\displaystyle\times(3+z(3+z)))))\Omega_{m0}^{2}))+\frac{1}{H_{0}^{2}}k^{2}(1+z)^{2}(-2048+\Omega_{m0}(6848-1344z(3+z$
$\displaystyle\times(3+z)-24(345+z(3+z)(3+z))(-166+z(3+z(3+z))))\Omega_{m0}+16(297$
$\displaystyle+z(3+z))(-144++z(3+z(3+z))(105+34z(3+z(3+z)))))\Omega_{m0}^{2}+(-1269$
$\displaystyle+z(3+z(3+z))(-324+z(3+z(3+z))(-1638+z(3+z(3+z))(-532$
$\displaystyle+3z(3+z(3+z))))))\Omega_{m0}^{3}))\bigg{]}/[3(4+(-3+z(3+z)))\Omega_{m0}^{8}].$

In this way, it becomes clear that by employing our model in Eq. (3.11), we
obtain an explicit dependence on the comoving wavenumber, $k$. Figure 9
illustrates the evolution of the second and third terms in Eq. (3.12), which
are given by Eqs. (3.13) and (3.14). From these figures, we can observe that
$b\,\delta G_{\rm{eff1}}/G_{\rm{N}}\ll 1$ and $b^{2}\,\delta
G_{\rm{eff2}}/G_{\rm{N}}\ll 1$. Therefore, the second and third terms in Eq.
(3.12) represent very small corrections to the Newton gravitational constant,
$G_{\rm{N}}$. These results better demonstrate the behavior observed in Fig. 8
(right panel). In principle, the other quantities studied here, such as the
cosmological and statefinder parameters, could be expressed through a series
expansion in the parameter $b$, where the first term would correspond to the
prediction of the $\Lambda$CDM model. However, for the sake of simplicity, we
have omitted the explicit writing of these expressions. In Ref. [39] a series
expansion in the parameter $b$ for $w_{\rm{DE}}$ is present, but for the HS
model.

Figure 9: (_left_) Evolution of $b\,\delta G_{\rm{eff1}}/G_{\rm{N}}$ vs. $z$
and (_right_) $b^{2}\,\delta G_{\rm{eff2}}/G_{\rm{N}}$ vs. $z$, using in both
cases some fixed values of $k$.

In the literature it is a common practice to rewrite Eq. (3.10) in terms of
the growth rate, which is defined by

$f_{\rm{g}}\equiv\frac{d\ln{\delta_{m}}}{d\ln{a}},$ (3.15)

where in this case “g” stands for “growth” and allows us to differentiate of
the $f(R)$ function. Realizing the replacements $dt=-\frac{1}{H(1+z)}dz$ and
$a=\frac{1}{1+z}$ in Eq. (3.10) and after some simple algebraic manipulations,
it take the equivalent form:

$\frac{df_{\rm{g}}(z)}{dz}+\left(\frac{1+z}{2E^{2}(z)}\frac{dE^{2}(z)}{dz}-2-f_{\rm{g}}(z)\right)\frac{f_{\rm{g}}(z)}{1+z}+\frac{3\Omega_{m0}(1+z)^{2}}{2E^{2}(z)}\frac{G_{\rm{eff}}(z,k)}{G_{\rm{N}}}=0,$
(3.16)

which has been expressed in terms of $E^{2}(z)$, since this function is known
in explicit form in our case. To numerically solve Eq. (3.10), we adopt an
initial condition for the growth rate that is consistent with that observed at
very high redshifts, matching that of the $\Lambda$CDM model. In Fig. 10 (left
panel), we illustrate the cosmological evolution of the growth rate
$f_{\rm{g}}$ as a function of redshift $z$ for various fixed values of the
comoving wavenumber $k$. Once again, we observe that its behavior is largely
unaffected by changes in $k$, and it closely resembles the predictions of the
$\Lambda$CDM model. In the right panel of Fig. 10, we compare the growth rate
$f_{\rm{g}}$ obtained using our model with the growth rate associated with the
$\Lambda$CDM model, denoted as $f_{{\rm{g}}\Lambda}$. The growth rate in the
$\Lambda$CDM model, $f_{{\rm{g}}\Lambda}$, is derived from Eq. (3.16) by
substituting $E^{2}(z)$ with $E_{\Lambda}^{2}(z)$, which is given by Eq.
(2.26). Additionally, we consider $G_{\rm{eff}}(z,k)=G_{\rm{N}}$. This
comparison it is carried out by using

$\Delta
f_{\rm{g}}=\left[\frac{|f_{\rm{g}}-f_{{\rm{g}}\Lambda}|}{f_{\rm{g}}}\right]\times
100.$ (3.17)

From Fig. 10 (right panel) is evident that the growth rate $f_{\rm{g}}$
differs slightly from that associated to the $\Lambda$CDM model,
$f_{{\rm{g}}\Lambda}$, only by $0<z<4$. The maximum value reached for $\Delta
f_{\rm{g}}$ is at $z=0$, $\Delta f_{\rm{g}}(z=0)\approx 0.9\%$.

Figure 10: (_left_) Evolution of $f_{\rm{g}}$ vs. $z$ and (_right_) evolution
of the relative difference $\Delta f_{\rm{g}}$ vs. $z$, using in both cases
some fixed values of $k$ and the initial condition $f_{\rm{g}}(z_{i}=50)=1$.

There is another way to represent the growth rate $f_{\rm{g}}$ that has been
widely used in the literature. It takes on the following form:

$f_{\rm{g}}(z)=[\Omega_{m}(z)]^{\gamma(z)},$ (3.18)

where $\gamma$ is known as the growth index. In our case $\Omega_{m}(z)$ is
given by

$\Omega_{m}(z)=\frac{\Omega_{m0}(1+z)^{3}}{E^{2}(z)}.$ (3.19)

Since in this work $E^{2}(z)$ is known, then we can replace Eq. (3.18) along
with Eq. (3.19) in Eq. (3.16), and this way we obtain an ordinary differential
equation (ODE) for the growth index, $\gamma(z)$. To solve numerically this
equation, we consider again that for very high redshifts, the growth index
matches that of the $\Lambda$CDM model (i.e. the growth index for the
$\Lambda$CDM model is $\gamma\approx 6/11$). In Fig. 11 (left panel), we have
plotted the growth index $\gamma$ as a function of redshift $z$. It is evident
that for the model being studied, the growth index is not constant. However,
similar to the other quantities analyzed previously, $\gamma(z)$ is
practically insensitive to changes in $k$. In the right panel of Fig. 11, we
illustrate the cosmological evolution of the growth rate utilizing Eq. (3.18),
while considering the previously obtained growth index, $\gamma(z)$. It is
worth noting that its behavior closely resembles what was observed in the left
panel of Fig. 10.
In Ref. [34], the authors explored various parameterizations for the growth
index $\gamma$ within a cosmological framework where the effective dark energy
originates from two specific $f(R)$ gravity models. They considered specific
ansatz for the growth index $\gamma$, namely $\gamma=\gamma_{0}$,
$\gamma=\gamma_{0}+\gamma_{1}z$, and
$\gamma=\gamma_{0}+\gamma_{1}\frac{z}{1+z}$. They fitted Eq. (3.18) to the
solution of Eq. (3.16) for different values of the comoving wavenumber $k$ in
the two studied $f(R)$ models. However, in our current work, we have not
addressed this aspect, as these parameterizations for the growth index will be
considered in a forthcoming study (currently in progress). In that future
work, we aim to determine the coefficients of the aforementioned ansatz and
the parameters of the model by utilizing a set of observational data to
achieve the best possible fit.

Figure 11: (_left_) Evolution of the growth index $\gamma$ as a function of
redshift $z$, using the initial condition $\gamma(z_{i}=50)=6/11$. (_right_)
$[\Omega_{m}(z)]^{\gamma(z)}$ vs. $z$, using in both cases some fixed values
of $k$.

##  4 Conclusions

In the recent years the $f(R)$ gravity it has positioned as an excellent
alternative to explain diverse phenomena in the cosmological and astrophysical
context. In this sense, in this work we have proposed a reparameterization of
a specific viable $f(R)$ gravity model to represent it as a perturbation from
the $\Lambda$CDM model. The $f(R)$ gravity model under consideration includes
two parameters, $b$ and $n$, and it is demonstrated that the closeness of this
model to $\Lambda$CDM depends on these parameters, allowing for arbitrary
proximity (see Eq. (2.15)). Furthermore, it is shown that the HS model is a
limiting case of this reparameterized model (see Eq. (2.22)). Following the
existing literature, we also derive an analytical approximation for the
expansion rate $H(z)$ (see Eq. (2.30)), which shows an excellent agreement
between this analytical approximation and the numerical solution over a wide
range of redshifts for realistic values of the deviation parameter $b$ (see
Fig. 1). By appropriately selecting values for the model parameters, we have
ploted the cosmological parameters $w_{\rm{DE}}$, $w_{\rm{eff}}$,
$\Omega_{\rm{DE}}$, and $H(z)$, as well as the statefinder quantities $q$,
$j$, $s$, and $Om(z)$ (see Figs. 2, 3, 4, 5, 6). We find that their present
values (at $z=0$) are consistent with the observations from Planck 2018 and
the values predicted by the $\Lambda$CDM model (see Table 1). It is important
to note that the examined cosmological and statefinder parameters do not
exhibit significant oscillations of effective dark energy, which could lead to
singular and unphysical solutions at high redshifts (see Fig. 7). This
anomalous behavior has been avoided by utilizing the approximate analytical
solution for $H(z)$. Additionally, we have performed a detailed analysis of
the evolution of matter density perturbations within the introduced $f(R)$
gravity model, and we have observed that its behavior is largely unaffected by
changes in $k$, and it closely resembles the predictions of the $\Lambda$CDM
model (see Figs. 8, 9, 10, 11). The results demonstrate that this viable
$f(R)$ gravity model is practically indistinguishable from the $\Lambda$CDM
model at the background level.

Acknowledgements
This work was supported by Patrimonio Autónomo-Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MINCIENCIAS-COLOMBIA) Grant No. 110685269447 RC-80740-465-2020,
projects 69723 and 69553.

## References

  * [1] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559-606 (2003) [arXiv:astro-ph/0207347 [astro-ph]].
  * [2] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006) [arXiv:hep-th/0603057 [hep-th]].
  * [3] K. Bamba, S. Capozziello, S. Nojiri and S. D. Odintsov, Astrophys. Space Sci. 342, 155-228 (2012) [arXiv:1205.3421 [gr-qc]].
  * [4] S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59-144 (2011) [arXiv:1011.0544 [gr-qc]].
  * [5] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, Phys. Rept. 513, 1-189 (2012) [arXiv:1106.2476 [astro-ph.CO]].
  * [6] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692, 1-104 (2017) [arXiv:1705.11098 [gr-qc]].
  * [7] J. c. Hwang and H. Noh, Phys. Lett. B 506, 13-19 (2001) [arXiv:astro-ph/0102423 [astro-ph]].
  * [8] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 (2003) [arXiv:hep-th/0307288 [hep-th]].
  * [9] S. Capozziello, V. F. Cardone and A. Troisi, Phys. Rev. D 71, 043503 (2005) [arXiv:astro-ph/0501426 [astro-ph]].
  * [10] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov and S. Zerbini, JCAP 02, 010 (2005) [arXiv:hep-th/0501096 [hep-th]].
  * [11] S. Capozziello, V. F. Cardone and M. Francaviglia, Gen. Rel. Grav. 38, 711-734 (2006) [arXiv:astro-ph/0410135 [astro-ph]].
  * [12] S. Nojiri and S. D. Odintsov, Phys. Rev. D 74, 086005 (2006) [arXiv:hep-th/0608008 [hep-th]].
  * [13] Y. S. Song, W. Hu and I. Sawicki, Phys. Rev. D 75, 044004 (2007) [arXiv:astro-ph/0610532 [astro-ph]].
  * [14] T. Faulkner, M. Tegmark, E. F. Bunn and Y. Mao, Phys. Rev. D 76, 063505 (2007) [arXiv:astro-ph/0612569 [astro-ph]].
  * [15] G. J. Olmo, Phys. Rev. D 75, 023511 (2007) [arXiv:gr-qc/0612047 [gr-qc]].
  * [16] I. Sawicki and W. Hu, Phys. Rev. D 75, 127502 (2007) [arXiv:astro-ph/0702278 [astro-ph]].
  * [17] V. Faraoni, Phys. Rev. D 75, 067302 (2007) [arXiv:gr-qc/0703044 [gr-qc]].
  * [18] R. Bean, D. Bernat, L. Pogosian, A. Silvestri and M. Trodden, Phys. Rev. D 75, 064020 (2007) [arXiv:astro-ph/0611321 [astro-ph]].
  * [19] S. Nojiri and S. D. Odintsov, Phys. Lett. B 657, 238-245 (2007) [arXiv:0707.1941 [hep-th]].
  * [20] S. Capozziello, A. Stabile and A. Troisi, Phys. Rev. D 76, 104019 (2007) [arXiv:0708.0723 [gr-qc]].
  * [21] N. Deruelle, M. Sasaki and Y. Sendouda, Prog. Theor. Phys. 119, 237-251 (2008) [arXiv:0711.1150 [gr-qc]].
  * [22] S. A. Appleby and R. A. Battye, JCAP 05, 019 (2008) [arXiv:0803.1081 [astro-ph]].
  * [23] S. Carloni, P. K. S. Dunsby and A. Troisi, Phys. Rev. D 77, 024024 (2008) [arXiv:0707.0106 [gr-qc]].
  * [24] S. Capozziello, V. F. Cardone and V. Salzano, Phys. Rev. D 78, 063504 (2008) [arXiv:0802.1583 [astro-ph]].
  * [25] P. K. S. Dunsby, E. Elizalde, R. Goswami, S. Odintsov and D. S. Gomez, Phys. Rev. D 82, 023519 (2010) [arXiv:1005.2205 [gr-qc]].
  * [26] S. Capozziello and M. De Laurentis, Annalen Phys. 524, 545-578 (2012)
  * [27] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 99, no.6, 064049 (2019) [arXiv:1901.05363 [gr-qc]].
  * [28] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 99, no.10, 104070 (2019) [arXiv:1905.03496 [gr-qc]].
  * [29] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 101, no.4, 044009 (2020) [arXiv:2001.06830 [gr-qc]].
  * [30] V. K. Oikonomou, Phys. Rev. D 103, no.12, 124028 (2021) [arXiv:2012.01312 [gr-qc]].
  * [31] V. K. Oikonomou, Phys. Rev. D 103, no.4, 044036 (2021) [arXiv:2012.00586 [astro-ph.CO]].
  * [32] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007) [arXiv:gr-qc/0612180 [gr-qc]].
  * [33] E. Elizalde, S. D. Odintsov, L. Sebastiani and S. Zerbini, Eur. Phys. J. C 72, 1843 (2012) [arXiv:1108.6184 [gr-qc]].
  * [34] K. Bamba, A. Lopez-Revelles, R. Myrzakulov, S. D. Odintsov and L. Sebastiani, Class. Quant. Grav. 30, 015008 (2013) [arXiv:1207.1009 [gr-qc]].
  * [35] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, Class. Quant. Grav. 38, no.7, 075009 (2021) [arXiv:2102.02239 [gr-qc]].
  * [36] F. Fronimos, Eur. Phys. J. Plus 136, no.10, 1014 (2021) [arXiv:2110.00353 [gr-qc]].
  * [37] V. K. Oikonomou and I. Giannakoudi, Int. J. Mod. Phys. D 31, no.09, 2250075 (2022) [arXiv:2205.08599 [gr-qc]].
  * [38] S. Basilakos, S. Nesseris and L. Perivolaropoulos, Phys. Rev. D 87, no.12, 123529 (2013) [arXiv:1302.6051 [astro-ph.CO]].
  * [39] R. Arjona, W. Cardona and S. Nesseris, Phys. Rev. D 99, no.4, 043516 (2019) [arXiv:1811.02469 [astro-ph.CO]].
  * [40] W. Cardona, R. Arjona, A. Estrada and S. Nesseris, JCAP 05, 064 (2021) [arXiv:2012.05282 [astro-ph.CO]].
  * [41] L. N. Granda, Eur. Phys. J. C 80, no.6, 539 (2020) [arXiv:2003.09006 [gr-qc]].
  * [42] A. Oliveros and Mario A. Acero, Phys. Dark Univ. 40, 101207 (2023) [arXiv:2302.07022 [gr-qc]].
  * [43] L. Granda, Symmetry 12, no.5, 794 (2020)
  * [44] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 104, no.12, 124065 (2021) [arXiv:2112.06269 [gr-qc]].
  * [45] L. N. Granda, JHEP 12, 205 (2021) [arXiv:2007.13956 [gr-qc]].
  * [46] E. V. Linder, Phys. Rev. D 80, 123528 (2009) [arXiv:0905.2962 [astro-ph.CO]].
  * [47] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, Phys. Rev. D 77, 046009 (2008) [arXiv:0712.4017 [hep-th]].
  * [48] J. Pérez-Romero and S. Nesseris, Phys. Rev. D 97, no.2, 023525 (2018) [arXiv:1710.05634 [astro-ph.CO]].
  * [49] W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007) [arXiv:0705.1158[astro-ph]].
  * [50] J. Sultana, M. K. Yennapureddy, F. Melia and D. Kazanas, Mon. Not. Roy. Astron. Soc. 514, no.4, 5827–5839 (2022) [arXiv:2206.10761[astro-ph]].
  * [51] S. Tsujikawa, Phys. Rev. D 76, 023514 (2007) [arXiv:0705.1032[astro-ph]].
  * [52] S. Tsujikawa, K. Uddin and R. Tavakol, Phys. Rev. D 77, 043007 (2008) [arXiv:0712.0082[astro-ph]].
  * [53] S. Tsujikawa, R. Gannouji, B. Moraes and D. Polarski, Phys. Rev. D 80, 084044 (2009) [arXiv:0908.2669[astro-ph.CO]].


PASA 2024

# Can an Infinitely Long Object Fit in an Expanding Universe?

Aaron GlanvilleA ∗ Tamara M. DavisA AThe University of Queensland, School of
Mathematics and Physics, QLD 4072, Australia  ∗E-mail<EMAIL_ADDRESS>
###### Abstract

Does space stretch its contents as the universe expands? Usually we say the
answer is no — the stretching of space is not like the stretching of a rubber
sheet that might drag things with it. In this paper we explore a potential
counter example — namely we show that is is impossible to make an arbitrarily
long object in an expanding universe, because it is impossible to hold the
distant end of the object ‘stationary’ with respect to us (as defined in the
Friedmann-Lemaître-Robertson-Walker metric). We show that this does not mean
that expanding space has a force associated with it, rather, some fictitious
forces arise due to our choice of reference frame. By choosing our usual time-
slice (where all comoving observers agree on the age of the universe), we
choose a global frame that does not correspond to the frame of any inertial
observer. As a result, simple relativistic velocity transforms generate an
apparent acceleration, even where no force exists. This effect is similar to
the fictitious forces that arise in describing objects in rotating reference
frames, as in the case of the Coriolis effect.

###### doi:

10.1017/pas.2024.xxx

###### keywords:

Cosmology: observations – Cosmology: theory

##  1 Introduction

The expansion of the universe is the cornerstone of modern cosmology
(Lemaître, 1927; Hubble & Humason, 1931). However, the description of
cosmological recession in the framework of expanding space, and indeed the
interpretation of expanding space itself, has been subject to significant
debate in the literature (Davis et al., 2003; Whiting, 2004; Barnes et al.,
2006; Grøn & Elgarøy, 2007; Chodorowski, 2007; Francis et al., 2007; Peacock,
2008; Cook & Burns, 2009; MacLaurin, 2015). Such discussion has centred around
two fundamental questions: “What does it physically mean to say space is
expanding?” and “Does this interpretation of expanding space provide a unique
description of the observed cosmological recession?”

Through this paper, we focus on illuminating some of the unique features of
expanding space by asking: “Is there a limit to the size of an object that can
be constructed in an expanding universe?” This question is interesting, as in
non-expanding space there is nothing that prevents us from making an
arbitrarily large object. Even in the case of expanding space, we usually
emphasise the “stretching” of space does not stretch or rip its contents.
Indeed, even separating space from its contents goes against the general
relativistic perspective that “Space tells matter how to move; matter tells
space how to curve.” (Misner et al., 1973). However, in the standard
Friedmann-Lemaître-Robertson-Walker (FLRW) metric, space can expand faster
than the speed of light (Davis & Lineweaver, 2004). The Hubble-Lemaître law
holds, allowing us to define a distance known as the Hubble radius
($D_{H}=c/H$) beyond which cosmological recession is superluminal. Objects
beyond this radius would require inward peculiar velocities ($v_{\rm p}$)
faster than the speed of light to fix their position with respect to us. Since
superluminal peculiar velocities are not possible, it would seem that
something is indeed stopping us from making a non-expanding object larger than
the Hubble sphere. In the absence of any force or tension, what feature of
expanding space prevents us from constructing an arbitrarily large object?

This paper continues the tradition of Davis et al. (2003), who explored the
physical consequences of expanding space using the radial dynamics of test
particles released in expanding universes. Their paper motivated the use of
the tethered galaxy system as a thought experiment to explore what it means to
say “space is expanding”. Through this thought experiment, a massless test
galaxy is defined with a fixed proper distance from some observer and zero
proper velocity, analogous to a galaxy that is physically tethered to some
central observer. This galaxy is then “released” and allowed to move freely.
Using this system Davis et al. (2003) showed that in the limit of non-
relativistic velocities, any induced radial motion of untethered galaxies is a
result of the acceleration of the expansion of the universe, but not the
expansion itself. Galaxies released from rest in accelerating cosmologies will
recede outward, and galaxies released in decelerating (but expanding)
cosmologies will move inward, often shooting past the observer and “rejoining
the Hubble flow”111For a detailed discussion of what it means to “rejoin the
Hubble flow” see Barnes et al. (2006). Here we are referring in particular to
their definitions 1 and 3. in the opposite direction. Davis et al. (2003)
argued such dynamics are consistent with the general relativistic description
of expanding space; this expansion alone is not associated with any force on
free particles, but accelerating expansion is.

In the case of a relativistic velocity of recession, velocity corrections
result in subtle but important modifications to the radial dynamics of these
tethered galaxies. With the inclusion of these relativistic effects, a
tethered galaxy released from rest in an expanding, non-accelerating cosmology
will still experience an apparent acceleration. It will accelerate inward
against the direction of the Hubble flow (Whiting, 2004; Grøn & Elgarøy,
2007). Although these relativistic dynamics are well accepted in the
literature, the physical interpretation of these results has undergone
significant debate. Papers such as Whiting (2004) and Peacock (2008) argue
that the inward motion of tethered galaxies released from rest in expanding
universes is qualitatively inconsistent with the picture provided through
expanding space. However, Barnes et al. (2006), Grøn & Elgarøy (2007), and
Francis et al. (2007) use the radial motion of free particles from the
relativistic geodesic equation to argue that expanding space provides a
natural description of the motion of a free particle, as long as relativistic
effects are carefully considered.

To shed light on the physical consequences of expanding space, we build upon
the body of work provided through the tethered galaxy experiment to explore
whether a physical tether could be used to hold a galaxy stationary at the
Hubble radius (the distance at which $v_{\text{rec}}=c$). We use the simplest
case of a universe that neither accelerates, nor decelerates – the empty
universe. Through this paper, we interpret this thought experiment through the
lens of several different reference frames, in order to separate events and
effects from the coordinates that describe them. As we shall see, the
definition of “stationary” itself is ambiguous, and what has zero velocity in
one coordinate system is not necessarily stationary in another.

Through Section II, we provide an overview of the global properties of the
empty universe. The empty universe does not experience any acceleration or
deceleration, and comoving particles recede with a constant velocity. We use
this feature to reaffirm that expanding space, when disentangled from
acceleration, is not associated with any force. In Section III we demonstrate
the relativistically corrected tethered galaxy system of Davis et al. (2003)
predicts an inward acceleration in the empty universe, and that this inward
motion, as a relativistic velocity effect, needs no force, or spacetime
curvature, to initiate it. In Section IV we extend this tethered galaxy system
to connect a central observer to a galaxy which lies on the Hubble radius. In
order to fix the position of a galaxy at the far end of our tether, we must
supply an inward peculiar velocity equal to $c$. Our inability to hold a
galaxy beyond the Hubble radius stationary is seemingly in contrast with the
conclusion that expanding space is not associated with any force, developed in
Section II.

We resolve this apparent contradiction by expressing this “Hubble tether”
system in the inertial frame of the central observer. By doing this, we show a
tether that is stationary at a constant time in FLRW coordinates must be
defined using points of varying inertial velocity. A tether along a constant
FLRW time-slice that is longer than the Hubble radius is impossible to
construct, as it requires $v_{\rm p}>c$ even though there is no Newtonian-like
force from expanding space. We express the inertial velocity profile of this
tether, and show it restricts the length of objects to within the Hubble
sphere, even in the limit of an empty universe.

##  2 The Canvas of Empty Space

###  2.1 General FLRW metric and expansion

The expansion of the universe is most commonly described within the framework
of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, in which galaxies
tracing the flow of expansion maintain a constant comoving coordinate
($\chi$). For a constant comoving coordinate, the amount these galaxies drift
apart is proportional to scalefactor ($R$), such that their distance from the
origin is $D=R\chi$ (Weinberg, 2008).

More fully, the standard FLRW metric is defined,

$ds^{2}=-c^{2}dt^{2}+R(t)^{2}\left[d\chi^{2}+S^{2}_{k}(\chi)(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\right],$ (1)

for interval $ds$, speed of light $c$, cosmological time $t$, comoving
separation $d\chi$, and scalefactor $R$ (with dimensions of distance),
expressed in spherical coordinates (Peacock, 1998). We are only interested in
radial distances, so $d\theta=d\phi=0$, leaving only $d\chi$. Thus we can see
that the distance defined by $D=R\chi$ represents the distance between the
origin and a comoving coordinate $\chi$ at a constant ‘cosmological time’
($t={\rm constant}$ or $dt=0$; where $t$ in FLRW also corresponds to
wristwatch time for comoving observers).

Differentiating this distance with respect to cosmological time gives two
components to the total velocity, $v_{\rm t}$,

$v_{\rm t}=\frac{dR}{dt}\chi+R\frac{d\chi}{dt}=v_{\rm r}+v_{\rm p}.$ (2)

Here the recession velocity, $v_{\rm r}$, comes entirely from the change in
scalefactor $R$, and the peculiar velocity, $v_{\rm p}$, comes entirely from a
change in comoving coordinate. Note that when you define the Hubble parameter
to be $H\equiv\dot{R}/R$ (where the overdot represents differentiation with
respect to FLRW time $t$), you recover the Hubble-Lemaître law, $v_{\rm
r}=HD$. To achieve a total velocity of zero, we require the peculiar velocity
to balance the recession velocity, $v_{\rm p}=-v_{\rm r}$.

None of the above requires general relativity, it is all determined by the
FLRW metric, which is the general metric for a homogeneous, isotropic
universe. However, in order to derive how scalefactor changes with time, and
how objects move through comoving coordinates, we do need a theory of gravity.
The Friedmann equation is the solution to general relativity with which we can
calculate the evolution of the universe,

$\left(\dfrac{\dot{a}}{a}\right)^{2}=H_{0}^{2}\left[\Omega_{\rm
K}a^{-2}+\Omega_{\rm M}a^{-3}+\Omega_{\rm R}a^{-4}+\Omega_{\Lambda}\right],$
(3)

where for convenience we have normalised the scalefactor to its value at the
present day (subscript zero), $a\equiv R/R_{0}$, so $D=R_{0}a\chi$. Each
$\Omega$ represents a different component of the universe, with $\Omega_{\rm
M},\;\Omega_{\rm R},\;{\rm and}\;\Omega_{\Lambda}$ representing respectively
matter, radiation, and cosmological constant; and $\Omega_{\rm
K}=1-\Omega_{\rm M}-\Omega_{\rm R}-\Omega_{\Lambda}$ relating to the curvature
of the universe.

###  2.2 Metrics of the Empty universe

####  2.2.1 Empty FLRW Metric

In exploring what it means to say “space is expanding,” it is important to
isolate the physical consequences of expanding space from those of
accelerating (or decelerating) expansion. In order to disentangle these
effects we make use of the cosmology where no such accelerations are
experienced: the empty universe. If apparent accelerations arise in the empty
universe, we know we must attribute them to whatever it means to say space is
expanding, as there are no other accelerating features available. We therefore
use this section to provide a brief overview of the global properties of the
empty universe.

Using the Friedmann equation with $\Omega_{\rm M}=\Omega_{\rm
R}=\Omega_{\Lambda}=0$, we find $\dot{a}=H_{0}$; as expected, the expansion
rate of the empty universe does not change with time. So the scalefactor of
this empty universe follows,

$a(t)=H_{0}t.$ (4)

For a galaxy at a comoving coordinate $\chi$, we may use Eq. 4 to calculate
its proper distance as a function of time, yielding,

$D(t)=H_{0}tR_{0}\chi.$ (5)

Comoving test particles in the expanding, empty universe recede with a
constant velocity $H_{0}R_{0}\chi$ throughout expansion. Such comoving
particles are not accelerated by any force or tension from uniformly expanding
space. This provides the backdrop of the empty universe in which we define the
Hubble tether system.

####  2.2.2 Milne Universe

Figure 1: Spacetime diagrams of the empty universe expressed in normalised
FLRW coordinates (with $c=1$) on the left, and using Minkowski coordinates as
defined by a Milne observer on the right. The same lines appear in both plots,
the only difference being the coordinate system expressing them. A constant
FLRW time slice is shown in blue, while the inertial observer’s constant time
slice is in red. Each point along the blue line corresponds to an observer
that has seen the same time pass since the beginning of the universe as the
observer at the origin. Each point along the red line corresponds to a
comoving particle with the same inertial time as our central observer at
$D=0,\ t=1$. We populate this cosmology with equally spaced comoving particles
in FLRW coordinates (black), which recede with a constant velocity. These
particles occupy the entire area of this spacetime diagram in FLRW
coordinates, or the future light cone in Minkowski space. Note that the
density of particles is homogeneous along a constant time slice in our FLRW
representation, but increases with radius along a constant time slice in our
Milne representation. Test particles that lie on the Hubble radius recede with
a constant FLRW velocity of $c$, denoted by the dotted line; in the Milne
metric these same test particles have sub-luminal velocities of $v_{\rm
M}=\tanh(1)\approx 0.762$. In both spacetime diagrams, we also include the
past light cone of our observer at $D=0,\ t=1$ (denoted by the yellow line).

An interesting alternative, but completely equivalent, description of the
empty universe is provided through the Milne model. The Milne universe recasts
the empty FLRW universe into the coordinates of a single, inertial observer,
defined to be the origin of a Minkowski spacetime. By expressing the dynamics
of particles from the inertial reference frame of our Milne observer (a
globally inertial reference frame), we aim to resolve confusion resulting from
how distances are defined in the spatially curved frame of the empty FLRW
universe.222Note that even though the Milne universe can not be directly
compared to the general FLRW cosmologies (with non-zero mass and energy), the
empty universe it corresponds to is a convenient cosmology to test our
assumptions of expanding space when disentangled from acceleration. The
physical interpretation we assign to the features of expanding space in FLRW
co-ordinates must be able to explain the observations of a Milne observer. If
our assumptions do not hold in this simple case (with no accelerating
features), then we can not expect them to hold in the more general case of
cosmologies with accelerating/decelerating expansion. The standard Minkowski
metric is defined,

$ds^{2}=-c^{2}dT^{2}+dl^{2}+l^{2}d\Phi^{2},$ (6)

for radial distance $l$, inertial time $T$, and 3D angular coordinates $\Phi$.
In the Milne universe, massless test particles spread from the big bang
($l=0,\ T=0$) with a constant velocity ($v_{\rm M}=l/T$) through this pre-
existing Minkowski space. Any test particle can be chosen as our “privileged”
inertial observer, and any choice will yield the same distribution of test
particles, meaning the Milne universe has no unique centre (Rindler, 1977).
While the Milne and empty FLRW universe describe the same cosmology, the key
difference between these models lies in their choice of spacetime coordinates.
The empty FLRW universe is explicitly defined such that the density of
particles along a constant time-slice is homogeneous, which is satisfied when
every particle agrees on their local, wristwatch time. Since the Milne and
empty FLRW universe describe the same cosmology, the Milne universe shares
this pattern of test particles. As seen in Figure 1 however, this distribution
of test particles is no longer homogeneous when measured along a constant
Milne time-slice. Time in the Milne universe is defined using the inertial
time of the privileged (but arbitrarily chosen) observer at the origin.
Objects that are further away (and as such, have been receding with a greater
velocity) undergo special relativistic length contraction from the reference
frame of this inertial observer. This length contraction results in a density
profile which increases with distance along our inertial Milne time-slice. The
pattern of test particles in both models is exactly the same, we have only
changed the coordinates used to represent them. We may directly transform the
coordinates of the empty FLRW universe into the coordinates of our inertial
Milne observer through the relation,

$l=ct\sinh\left(\chi\right);\quad T=t\cosh\left(\chi\right).$ (7)

The empty FLRW and Milne universe are interchangeable models that describe the
same pattern of test particles within the same, empty cosmology. However, as
we shall develop in the next section, their differing choices of time-slice
yield different interpretations of recession and expansion that we wish to
explore.

###  2.3 The distribution of test particles in empty-FLRW vs Milne universes

Physical processes are coordinate independent, but the language we use to
describe these processes is often rooted in the coordinate system we use. In
order to highlight how using an inertial reference frame affects the
description of the empty universe, consider a series of comoving test
particles defined with an equal comoving separation in the FLRW metric. Along
a constant FLRW time-slice (the slice where all test particles agree on the
age of the universe), one would record a homogeneous test particle density
extending over all space (Figure 1, left).

By contrast, when this is converted to the Milne universe, we find a density
that increases with radius. As developed earlier, when we cast our empty
universe in the reference frame of a privileged, inertial observer, special
relativistic length contraction results in an increasing test particle density
with distance (Figure 1, right). This effect causes all test particles (and by
extension, the infinite FLRW universe) to be contained within the finite
distance of the future light cone of the big bang in the Milne universe.

Our choice of time-slice has particularly significant consequences when
considering the recession of a galaxy on the Hubble radius. In the FLRW metric
objects can recede with velocities exceeding $c$ without any relativistic
violation. This superluminal recession occurs outside the inertial reference
frame of any observer, resulting in no special relativistic contradiction
(Murdoch, 1977; Davis & Lineweaver, 2004). As such, the more distant the
galaxy, the faster it can recede without limit. However, the Milne universe
(cast on a flat, Minkowski spacetime) is subject to special relativistic
constraints, meaning recession greater than $c$ is not possible in this
coordinate system.

The fact these two coordinate systems demonstrate such differing
interpretations of recession for the same universe highlights the fact that
velocity is a coordinate dependent property. A coordinate-independent feature
would be the four-velocity, but that is not the velocity that appears in
Hubble’s law, and more broadly it is not the standard velocity invoked to
interpret cosmological recession. Importantly, dimensionless observables have
to remain true regardless of coordinates, such as redshift.

By contrasting the description of cosmological recession in the empty universe
expressed in the framework of the Milne and FLRW coordinate systems, we aim to
isolate the properties of expansion that are fundamental physical
characteristics, from those that are a consequence of our coordinate
representation.

##  3 The Tethered Galaxy System

With the backdrop of the empty universe now firmly established, we use this
section to explore the tethered galaxy system as a framework for understanding
the features of expanding space. In this section, we closely follow the
tethered galaxy system outlined in Davis et al. (2003), replicating their
results in the limit of non-relativistic velocities. We then apply a
relativistic correction provided (but not directly used) in their paper.

In generating these results, we make use of the fact that (in the limit of
non-relativistic velocities) momentum with respect to the local comoving frame
decays as $1/a$ (Misner et al., 1973). We demonstrate that such scalefactor
arguments (when relativistic effects are included) display the same behaviour
as the full general relativistic derivation, which you can find in Grøn &
Elgarøy (2007).

###  3.1 Non-Relativistic Limit

In an expanding universe, distant galaxies undergo cosmological recession,
$v_{\rm r}=HD$. In the tethered galaxy system, we imagine that, as an initial
condition, an inward peculiar velocity is supplied, which is of an equal
magnitude to this recessional velocity $v_{{\rm p},0}=-v_{{\rm r},0}$. As
such, the proper distance $D$ of the galaxy from the observer is initially
fixed and the initial total velocity $v_{{\rm t},0}=0$ (as if it had been
physically tethered out of the Hubble flow). Equivalently,

$R_{0}\dot{\chi}_{0}=-\dot{R}_{0}\chi_{0}.$ (8)

The peculiar momentum of this tethered galaxy ($p=mv_{\rm p}$) will decay with
respect to the comoving background frame as $1/a$ (Misner et al., 1973). In
the non-relativistic limit, the peculiar velocity tethering this galaxy decays
by this same factor,

$v_{\rm p}=\frac{v_{{\rm p},0}}{a}.$ (9)

Combining this result with Equation 8 gives,

$R\dot{\chi}=\frac{-\dot{R}_{0}\chi_{0}}{a}.$ (10)

Recalling $R=R_{0}a$ and integrating from $\chi_{0}$ to $\chi$ and $t_{0}$ to
$t$ gives,

$\displaystyle D$ $\displaystyle=$ $\displaystyle
R\chi_{0}\left(1-H_{0}\int_{t_{0}}^{t}\frac{dt}{a^{2}}\right),$ (11)
$\displaystyle=$ $\displaystyle
R\chi_{0}\left(1-H_{0}\int_{1}^{a}\frac{da}{\dot{a}a^{2}}\right).$ (12)

One can compute this integral numerically for any arbitrary cosmology by using
Friedmann’s equation to supply $\dot{a}$. We calculate this for the same set
of cosmologies described in Davis et al. (2003), directly replicating their
results as the solid lines in Figure. 2. The empty universe is shown by the
blue line, and you can see that in the non-relativistic limit a galaxy that
starts at rest ($v_{\rm t}=0$) in the empty universe remains at rest.

Figure 2:  The dynamics of a tethered galaxy in four different cosmologies as
specified in the legend. The non-relativistic solution is given as solid lines
(replicating Fig. 2 of Davis et al., 2003), while the relativistic solutions
are given as dashed lines. In the non-relativistic case the tethered galaxy
moves inward for cosmologies with decelerating expansion, outward for
cosmologies with accelerating expansion, and stays at a constant distance for
cosmologies with constant expansion. In contrast, the relativistic solution
shows a small inward co-ordinate acceleration in the universe with constant
expansion. We clarify that the trajectories of these untethered galaxies
correspond to geodesic motion. This is a spurious acceleration arising because
of our coordinate choice, and is not supplied by any force. This acceleration
is a co-ordinate acceleration (i.e, $\frac{d^{2}D}{dt^{2}}\neq 0$), and does
not correspond to any relativistic 4-acceleration. This release distance
($100$ Mpc) is similar in magnitude to the BAO scale ($\sim 150$ Mpc), however
since the BAO feature corresponds to a comoving standard ruler, rather than a
fixed physical size, one would not expect peculiar velocity decays to have an
effect.

###  3.2 Relativistic Limit

As the peculiar velocity of tethering enters relativistic regimes, velocity
corrections result in subtle but important modifications to the dynamics of
untethered galaxies (Whiting, 2004; Grøn & Elgarøy, 2007). For a relativistic
peculiar velocity, we must use the relativistically adjusted momentum formula,

$p=\gamma mv_{\rm p},$ (13)

where $\gamma=(1-v_{\rm p}^{2}/c^{2})^{-1/2}$. Since $p\propto 1/a$, we
obtain,

$v_{\rm p}=\frac{\gamma_{0}v_{{\rm p},0}}{\sqrt{a^{2}+\gamma_{0}^{2}v_{{\rm
p},0}^{2}/c^{2}}}.$ (14)

Figure 3 shows this results in a continuum of behaviours between photons
(whose peculiar velocities do not decay) to non-relativistic massive particles
(whose peculiar velocities decay as $1/a$). In between these extremes,
relativistic massive particles have velocities that decay more slowly than
$1/a$. Interestingly, while the wavelengths ($\lambda$) of photons are
redshifted by the expanding universe, so are the de Broglie wavelengths of
massive particles. In both cases $\lambda\propto a$, however for massive
particles the redshifting of their de Broglie wavelength manifests itself as a
decrease in peculiar velocity, while for photons it just results in a longer
wavelength.

Using this relativistically adjusted peculiar velocity function to replace
Equation 9 in the derivation of Equation 12 yields,

$D=R\chi_{0}\left(1-\gamma_{0}H_{0}\int_{1}^{a}\frac{da}{\dot{a}a\sqrt{a^{2}+\gamma_{0}^{2}v_{{\rm
p},0}^{2}/c^{2}}}\right).$ (15)

When applied to the same cosmologies as in the previous section this generates
the dashed lines showing the relativistic trajectories of untethered galaxies
in Fig. 2.

It is clear that in the full relativistic solution, we find that tethered
galaxies released from rest in the expanding, empty universe experience an
apparent inward acceleration. Interestingly, the tether does not get stretched
apart by expanding space as we might initially assume, but rather collapses.
This inward acceleration has been described by others in the literature,
notably Whiting (2004) and Grøn & Elgarøy (2007).

###  3.3 Why peculiar velocities decay

There is a simple physical explanation of why peculiar velocities decay.
Imagine a vertical worldline in the FLRW spacetime diagram of Figure 1. (Such
a worldline is only physically possible within the Hubble sphere, although
that detail is not important for this argument.) A vertical worldline means
the object is stationary in these coordinates ($v_{\rm t}=0$), and therefore
it must have an inward peculiar velocity equal and opposite to the outward
recession velocity ($v_{\rm p}=-v_{\rm r}$). As time goes on, the object gets
overtaken by comoving galaxies that are moving more and more slowly (the solid
black lines become more vertical). So to stay at a constant distance, the
object’s peculiar velocity has to decrease. This example shows that peculiar
velocity decay is just a bookkeeping exercise, arising because we have split
our velocity into two components. The total velocity $v_{\rm t}$ is zero and
stays zero — no accelerations have been applied.

Figure 3: How peculiar velocities decay as the universe expands, showing the
continuum between relativistic and non-relativistic behaviour. The non-
relativistic solution has $v_{\rm p}=v_{\rm p,0}/a$ and is shown as a dotted
line. The relativistic solution is given by Eq. 14 and is shown as solid lines
for 7 different initial velocities as labelled. The initial scalefactor is
$a=0.02$, and the $v_{\rm init}=0.9c$ and $0.6c$ cases begin from the same
scalefactor as the dotted line for easy comparison.

The surprising thing is that in the relativistic case, the peculiar velocity
does not decay fast enough to keep the object at $v_{\rm t}=0$. When we apply
the relativistic correction, the inward peculiar velocity decays at a reduced
rate, generating an inferred inward acceleration as in Figure 2. This inward
acceleration of tethered galaxies released in the empty universe is just a
result of relativistic addition of velocities, and does not require any
Newtonian-like force.

Figure 4: Evolution of Hubble tether system whose reference frame is defined
to coincide with a constant FLRW timeslice from $t=1$ to $t=1.2$. As a
massless test galaxy is tethered out of the Hubble radius, the FLRW distance
of the galaxy from the central observer is fixed. This causes the tethered
galaxy to follow the green, vertical worldline centred around $D=1$.
Similarly, all of the elements of our tether have fixed FLRW positions, and
follow parallel vertical worldlines. After following these vertical
worldlines, these points define a tether of the same length at the later time
$t=1.2$ (given by the later blue line). In the Minkowski coordinates of the
Milne universe the Hubble tether follows a hyperbola. In order to stay the
same length in FLRW space, it must shrink in Minkowski space. The green lines
are sloped inwards, indicating the velocities through this inertial frame for
the different parts of the tether. When the end is on the Hubble sphere, the
inward velocity needed for that end to stay at a constant FLRW distance is
$c$. This is why you cannot tether a particle beyond the Hubble sphere.

##  4 The “Hubble Tether” System

###  4.1 Constructing a “Hubble Tether”

Now we take the relativistic limit to the extreme. In order to evaluate the
physical consequences of expanding space in maximally relativistic regimes, we
consider a massless test galaxy that lies on the Hubble radius of some central
observer, which is held in place by a tether spanning this distance. We
imagine the job of designing this “Hubble tether” system has been delegated to
a student of cosmology, who defines the tether to lie along a constant FLRW
time-slice (given by the horizontal blue line of the left panel of Figure 1).
To achieve this, they require every point along this tether be defined with
massless test particles, that are at rest in FLRW coordinates (i.e.
$v_{\text{t}}=0$), and stay at rest.333The term “at rest” relative to the
Hubble expansion is often used to refer to comoving galaxies. Here we use it
to mean at rest with respect to the FLRW coordinates. This means each point
has a peculiar velocity equal but opposite to the recessional velocity at that
point, and they are adjusted as necessary to stay that way, yielding the
spacetime diagram of Figure 4 (left).444For a more observer-centred approach
to constructing this system, one could imagine each element along our tether
uses radar ranging to fix their local positions (following the treatment of
Lewis et al., 2008). In the limit of small distances the empty FLRW universe
is spatially flat, meaning each element could accelerate/decelerate as
necessary to remain stationary with respect to their adjacent elements,
defining our stationary tether. Note that this is different to defining a
“rigid” rod. We choose to analyse this system without referring to rigidity,
due to the technical challenges associated with defining such a property in a
relativistic context (Ch. 16 and Ch. 2 of Panofsky & Phillips, 1964; Landau et
al., 1971, respectively). In particular, rigidity may seem to imply there are
forces internal to the rod holding it in place, raising the issue of faster-
than-light communication to transmit such forces. Our “rod” or “tether” is
simply defined to be a series of adjacent particles, of fixed separation in
the frame defined by FLRW co-ordinates. Whether or not these particles move
away from or toward each other (in a Minkowski frame) when allowed to move
freely will determine whether a rod or tether would need internal forces in
order to hold itself together or stop itself collapsing.

What does such a system look like from the reference frame of an inertial
observer? When we cast this system in the inertial co-ordinates of our Milne
observer (an observer at rest along $D=l=0$), the Hubble tether system
(defined over a constant FLRW time-slice) now stretches across a hyperbolic
slice. Using the Milne-FLRW conversion of Equation 7, we can express the time
evolved Hubble tether system in the inertial coordinates of the central
observer, yielding Figure 4 (right).

By casting our time evolved Hubble tether system in the coordinates of an
inertial observer, we find the fixed elements that define a constant tether in
FLRW coordinates must be constructed using a range of inertial velocities.
While the observer end of our tether follows a stationary vertical worldline
in the Milne universe (about $l=0$), the points extending along this tether
require accelerating kinematic profiles to maintain their FLRW position. It is
this requirement that the tether be at rest against the frame defined by the
standard FLRW coordinate system, which does not correspond to the inertial
reference frame of any observer, that causes this tether to collapse by
construction.

In contrast, if we were to define the tether at a constant length in the
arguably more natural reference frame of an inertial observer, we arrive at
the spacetime diagram of Figure 5. Note that when the elements of our tether
are defined with constant separations in this inertial reference frame, our
object does fulfil the condition of Born rigidity. When we cast our inertial
tether back into FLRW coordinates, we find the far end is embedded in the past
(which has a higher test particle density). Additionally, in order to remain
at a constant position with respect to our inertial observer, the elements
defining this tether must move outward in FLRW coordinates. This shows that a
tether which is at rest in the inertial Minkowski frame, would appear to
lengthen in the FLRW frame, even though no stretching is occuring! Conversely,
this shows that a tether which is “at rest” in FLRW coordinates in not truly
at rest in the inertial frame of any observer.

Figure 5: Alternative view of a “stationary” tether, this time defined to be
stationary in the coordinates of the Milne model. The inertial tether, when
expressed in FLRW coordinates, actually has its distant end embedded in the
past— in a universe with a higher test particle density. Additionally, in
order to remain at a fixed distance as measured by our inertial observer, the
elements defining this tether move along curved worldlines in FLRW co-
ordinates. As a result, our student assessing the tether design would infer
the particles defining this tether are moving apart.

###  4.2 Inertial Velocities along tether

While the previous section allows us to describe the behaviour of our constant
tether in the empty universe, a key part of our question still remains— What,
if anything, is significant about a Hubble tether? Is there anything
significant about the Hubble sphere in the description of our tethered galaxy
system? To answer this we return to our constant FLRW tether, and explore the
inertial velocities of the elements along this length (with an FLRW velocity
of 0) from the central observer to the Hubble sphere. To develop this profile,
we constructed two lines of constant length in FLRW coordinates, separated by
some infinitesimal time $dt$. This corresponds to the time evolved Hubble
tether system of Figure 4, in the limit of an infinitesimal time separation.
Using the FLRW-Milne relation of Equation 7, we solved for the initial
($l_{\textrm{init}},T_{\textrm{init}}$) and time separated
($l_{\textrm{final}},T_{\textrm{final}}$) Milne coordinates of the elements
along these horizontal FLRW lines. This allows for a numerical evaluation of
the Milne velocities required to yield $v=0$ in FLRW coordinates with the
simple relation:

$\frac{dl}{dT}=\frac{(l_{\textrm{final}}-l_{\textrm{init}})}{(T_{\textrm{final}}-T_{\textrm{init}})}.$
(16)

The inertial velocity profile of particles along the length of our Hubble
tether is given as in Figure 6.

Figure 6: Profile of inertial velocities required to construct Hubble tether
system along $t=1$, as a function of both Milne and FLRW distance. As one
follows a hyperbolic time slice of constant FLRW $t$, the peculiar velocities
required to tether these points with a fixed FLRW distance increases
dramatically. At the end of our Hubble tether (along the dotted line of
$l\approx 1.175$), an inertial velocity of precisely $c$ is required to tether
this system in FLRW coordinates.

We can see the Hubble sphere represents a fundamental limit to how far we can
define such a chain of particles, regardless of any rigidity constraints. In
order for a tether that extends to the Hubble radius to maintain a constant
length in FLRW coordinates, we must supply an inward inertial velocity equal
to $c$, beyond which is relativistically forbidden. Imagine we define a series
of observers, each of which have rockets to change their inertial velocity. As
we attempt to extend this chain of observers along a constant FLRW time-slice,
the inertial velocities required to fix their FLRW position increases
dramatically. As we extend this chain to a point on the Hubble radius of the
central observer, we must supply an inertial velocity of $c$ to fix our FLRW
position. The worldline of a galaxy on the Hubble radius is subluminal in
Milne coordinates (as seen in Figure 1), but fixing the position of this
particle in FLRW coordinates beyond the Hubble sphere requires superluminal
motion through this inertial frame, which is why such a tether cannot be
constructed.

##  5 Discussion and Conclusion

With the discovery of the accelerating expansion of the universe, the need for
clarity in how we interpret this expansion has become increasingly important.
The description of the expanding universe, and the interpretation provided
through expanding space, fundamentally frame our approach to open problems in
cosmology. By exploring the dynamics of a particle tethered out of the Hubble
flow in the empty universe, we show such a particle experiences an inward
acceleration. While such a result has been presented in the literature
previously, our treatment highlights how this is well described as a
relativistic velocity decay. Such dynamics do not require any force or
acceleration from expanding space.

The properties of expanding space appear to come into question again when
considering the tethering of a galaxy on the Hubble radius (with a recessional
velocity of $c$) in the empty universe. We find such a tether cannot be
constructed, seemingly in contrast to the conclusion that expanding space does
not stretch the matter it contains. This apparent contradiction is resolved
when considering the effect of our choice of time-slice. By choosing the
standard FLRW time-slice, we choose a global frame that does not correspond to
the frame of any inertial observer. Simple relativistic transforms introduce
accelerations where no such force exists, similar to the fictitious forces
which arise in rotating reference frames. By enforcing the requirement that
each element of our tether is at rest in the frame defined by the standard
FLRW coordinate system (which does not correspond to any inertial reference
frame), we require they accelerate in the reference frame of an inertial
observer. It is this requirement that causes our tether to collapse, and not
any force from the expansion of space.

It is interesting that we cannot construct a Hubble tether which is stationary
with respect to the global frame defined by the FLRW coordinate system because
objects at $D>D_{H}$ cannot be held at a fixed FLRW distance – they must move
away, seemingly ‘tearing’ any attempt to make a tether in this regime. Whereas
within the Hubble sphere once you set up the initial condition of $v_{t}=0$
the tether does not rip apart, it collapses. That behaviour is true for an
empty universe, but in general it differs depending on what type of universe
you live in. Assuming a perfect fluid with equation of state $w=P/\rho$, Grøn
& Elgarøy (2007) determined the initial displacement necessary to observe
inward acceleration, as a function of $w$,

$\chi_{0}^{2}>-\dfrac{1}{2}(1+3\omega)\left(\dfrac{c}{H_{0}}\right)^{2}.$ (17)

If $w=-1$ this corresponds to exactly the Hubble sphere, so you do not get
inward acceleration. However, in the empty universe we have been considering,
for which $w=-1/3$ you get inward acceleration from any $\chi_{0}>0$.
Decelerating universes ($w>-1/3$) also obviously have inward acceleration
regardless of the distance.

To return to the question we posed in the abstract, does this mean we cannot
construct an infinitely long tether? If we define this tether along a constant
inertial time-slice, there is no reason it could not extend to the edge of the
Milne universe (the future light cone of the big bang). In such a case, the
tether would follow the red inertial time-slice of Figure 1, right. As this
constant tether extended further and further, the far end would extend further
and further back in cosmological time ($t$). In order for this ruler to extend
to the edge of the Milne universe (where particles have an inertial
recessional velocity of $c$), the far end would be embedded at $t=0$ in the
FLRW universe.

While questions regarding the fundamental nature of expanding space may seem
esoteric, it is interesting to consider the potential physical implications of
relativistic peculiar velocity decays within our own universe. Super-clusters
form filaments that can span up to 10 Gly, occupying up to 5% of the
observable universe (Horvath et al., 2013). Given the existence of structure
at such large scales, it is interesting to consider whether analogues of
proper tethers exist on cosmological scales. If systems with fixed proper
distances do indeed exist on such scales, one would expect the relativistic
decay of peculiar velocities would provide an increased inward “kick” to
support a slightly higher growth rate for such structures than calculated in a
purely Newtonian simulation.

###### Acknowledgements.

We thank Colin MacLaurin, Geraint Lewis, and Angela Ng for useful discussions.
TMD is the recipient of an Australian Research Council Australian Laureate
Fellowship (project number FL180100168) funded by the Australian Government.
AG is the recipient of an Australian Government Research Training Program
(RTP) Scholarship.

##  Appendix A Conversion from Milne to empty FLRW universe

Here we demonstrate the Milne and FLRW universe are identical. Given a
constant velocity of recession ($v_{\rm M}=\frac{l}{T}$) in conjunction with
the conversion of Equation 7, we may express our special relativistic $\gamma$
in terms of our comoving space and time ($\chi,t$) coordinates,

$\gamma=\dfrac{1}{\sqrt{1-v_{m}^{2}/c^{2}}},$ (18)

$\gamma=\dfrac{1}{\sqrt{1-\left(\frac{ct\sinh(\chi)}{t\cosh(\chi)}\right)^{2}/c^{2}}},$
(19)

$\gamma=\cosh{\chi}.$ (20)

Substituting the time relation of Equation 7 and Equation 20 into the Lorentz
boosted Minkowski time coordinate ($T=\gamma T^{\prime}$) yields,

$T^{\prime}=\dfrac{t\cosh(\chi)}{\cosh{(\chi)}},$ (21)

$T^{\prime}=t.$ (22)

Next, we take the derivative of $l^{\prime}$,

$dl^{\prime}=ct\cosh\left(\chi\right)d\chi+c\sinh(\chi)dt.$ (23)

Along a constant time-slice ($dt=0$);

$dl^{\prime}=ct\cosh\left(\chi\right)d\chi.$ (24)

Substitution the space relation of Equation 7, along with Equtions 22 and 24
into our standard Lorentz boosted Minkowski metric yields

$ds^{2}=-c^{2}dt^{2}+c^{2}t^{2}\left(d\chi^{2}+\sinh^{2}(\chi)\right).$ (25)

This provides the form of the standard FLRW metric (Equation 1) with a
curvature parameter $k=-1$, $R(t)=ct$, along a radial line of sight
($d\theta=d\phi=0$).

## References

  * darkbluergb0,0,0.597656
  * Barnes et al. (2006) Barnes L. A., Francis M. J., James J. B., Lewis G. F., 2006, darkblueMNRAS, 373, 382
  * Chodorowski (2007) Chodorowski M. J., 2007, darkblueMNRAS, 378, 239
  * Cook & Burns (2009) Cook R. J., Burns M. S., 2009, darkblueAmerican Journal of Physics, 77, 59
  * Davis & Lineweaver (2004) Davis T. M., Lineweaver C. H., 2004, darkbluePASA, 21, 97
  * Davis et al. (2003) Davis T. M., Lineweaver C. H., Webb J. K., 2003, darkblueAmerican Journal of Physics, 71, 358
  * Francis et al. (2007) Francis M. J., Barnes L. A., James J. B., Lewis G. F., 2007, darkbluePASA, 24, 95
  * Grøn & Elgarøy (2007) Grøn Ø., Elgarøy Ø., 2007, darkblueAmerican Journal of Physics, 75, 151
  * Horvath et al. (2013) Horvath I., Hakkila J., Bagoly Z., 2013, The largest structure of the Universe, defined by Gamma-Ray Bursts
  * Hubble & Humason (1931) Hubble E., Humason M. L., 1931, darkblueApJ, 74, 43
  * Landau et al. (1971) Landau L. D., Lifshitz E. M., Hamermesh M., 1971, The classical theory of fields. Pergamon
  * Lemaître (1927) Lemaître G., 1927, Annales de la Société Scientifique de Bruxelles, 47, 49
  * Lewis et al. (2008) Lewis G. F., Francis M. J., Barnes L. A., Kwan J., James J. B., 2008, darkblueMNRAS, 388, 960
  * MacLaurin (2015) MacLaurin C., 2015, Expanding Space, Redshifts, and Rigidity- Conceptual Issues in Cosmology
  * Misner et al. (1973) Misner C. W., Thorne K. S., Wheeler J. A., 1973, Gravitation. W. H. Freeman and Company
  * Murdoch (1977) Murdoch H. S., 1977, QJRAS, 18, 242
  * Panofsky & Phillips (1964) Panofsky W. K. H., Phillips M., 1964, Classical electricity and magnetism. Addison-Wesley
  * Peacock (1998) Peacock J., 1998, Cosmological Physics. Cambridge University Press
  * Peacock (2008) Peacock J. A., 2008, A diatribe on expanding space (arXiv:0809.4573)
  * Rindler (1977) Rindler W., 1977, Essential Relativity: Special, General, and Cosmological. Springer-Verland
  * Weinberg (2008) Weinberg S., 2008, Cosmology. Oxford University
  * Whiting (2004) Whiting A. B., 2004, The Observatory, 124, 174