# Confronting the Chaplygin gas with data: background and perturbed cosmic
dynamics

Renier Hough1, Shambel Sahlu2,3, Heba Sami4, Maye Elmardi1, Anna-Mia Swart4,
and Amare Abebe4,5 1 Centre for Space Research, North-West University,
Potchefstroom, South Africa.  2 Astronomy and Astrophysics Research and
Development Department, Entoto Observatory and Research Centre.  3 Department
of Physics, College of Natural and Computational Science, Wolkite University,
Wolkite, Ethiopia  4 Centre for Space Research, North-West University,
Mahikeng, South Africa.  5 National Institute for Theoretical and
Computational Sciences (NITheCS), South Africa.

###### Abstract

In this paper, we undertake a unified study of background dynamics and
cosmological perturbations in the presence of the Chaplygin gas. This is done
by first constraining the background cosmological parameters of different
Chaplygin gas models with SNIa data, and then feeding these observationally
constrained parameters in the analysis of cosmological perturbations. Based on
the statistical criteria we followed, none of the models has a substantial
observational support but we show that the so-called ‘original’ and
‘generalized’ Chaplygin gas models have some observational support and less
observational support, respectively, whereas the ‘modified’ and ‘modified
generalized’ Chaplygin gas models miss out on the less observational support
category but cannot be ruled out. The so-called ‘generalized cosmic Chaplygin
gas’ model, on the other hand, falls under the no observational support
category of the statistical criterion and can be ruled out. We follow the
$1+3$ covariant formalism of perturbation theory and derive the evolution
equations of the fluctuations in the matter density contrast of the matter-
Chaplygin gas system for the models with some or less statistical support. The
solutions to these coupled systems of equations are then computed in both
short-wavelength and long-wavelength modes.

PACS numbers

04.50.Kd, 98.80.Jk, 98.80.-k, 95.36.+x, 98.80.Cq

Corresponding Author

<EMAIL_ADDRESS>
###### pacs:

04.50.Kd, 98.80.Jk, 98.80.-k, 95.36.+x, 98.80.Cq

††preprint: APS/123-QED

##  I Introduction

The major component of the matter energy density (about 85%) is made from dark
matter, and currently, our Universe is experiencing an accelerated expansion.
The discovery of the accelerating expansion of the Universe and the existence
of dark matter are the primary indicators for the limitation of our knowledge
of physics laws and those two uncharted territories of the Universe have
opened a new era for the development of modern cosmology. In recent times, one
of the most active areas of research in cosmology is trying to understand the
nature of dark matter and dark energy. Different suggestions have been put
forward to understand dark energy. The first suggestion is that the
cosmological constant is the one responsible for the cosmic acceleration
allemandi2005dark ; riess1998observational ; frieman2008dark , and the second
approach is the modification of General Theory of Relativity (GR)
clifton2012modified ; faraoni2005stability . Meanwhile, scholars have been
proposing the Chaplygin Gas (CG) model saadat2013viscous ; dev2003cosmological
to mimic the effects of dark energy at late times and dark matter prior to
that, and can be a possible substitution to our standard model of cosmology.
This model comes in several versions Elmardi:2018ryn  the most common of which
are the original Chaplygin gas (OCG), generalized Chaplygin gas (GCG),
modified generalized Chaplygin gas (MGCG) and generalized cosmic Chaplygin gas
(GCCG) and has mostly been used to study the Friedmann-Lemaître-Robertson-
Walker (FLRW) background expansion history without much emphasis on
observational viability and the implications to large-scale structure
formation. To our knowledge there is no work in the literature that considered
the cosmological perturbations of this fluid model in the $1+3$ covariant
formalism of perturbation theory ellis1989covarianta ; dunsby1991gauge ;
abebe2012covariant . This work attempts to fill this important gap of using
background-constrained CG models in the study of large-scale structure
formation through the covariant perturbations approach.

In the background cosmology, some observational constraints such as
deceleration parameter, the luminosity distance and the observational Hubble
parameter have been studied in the afore-mentioned CG models using two
different datasets to find the best fits for each of the different models, to
be able to make conclusions about the viability of each, as well as getting a
better feel for how these models compare to the concordance model of
cosmology. For our first dataset, we will use Supernovae Type 1A data. This
method is called supernovae cosmology and with the cosmic microwave background
(CMB) radiation cosmology, it is the two most popular ways to determine with
observable data, the values of the cosmological parameters. The class of
supernovae we will be using is when a binary white dwarf (WD) star has
accreted enough hydrogen from the low-mass companion star, to form an outer-
layer, which when compressed enough results in the WD exploding Tayler1994 .
These types of explosions are regarded as standard candles, since the process
is always the same, resulting in their luminosities being relatively similar.
Therefore, when measuring the incoming flux, it can be used to approximate the
distance to the particular supernova, since distance is the only variable when
calculating the flux of a standard candle. We can use the obtained redshift
$(z)$ for each supernova to represent the approximated distance. The reason
for doing this comes down to the fact that we can use the distance modulus
function to measure the expansion of the Universe.

The other aspect of this work focuses on linear cosmological perturbations in
the $1+3$ covariant formalism in multi-fluid systems abebe2012covariant ;
sahlu2020scalar ; ntahompagaze2020multifluid  adapted to include CG fluid
models. After we derive the density perturbation equations and present the
solutions in a step-by step manner, both analytically and numerically, we will
discuss the cosmological implications as far as large-scale structure
formation is concerned on both sub- and super-horizon scales, and see if there
is a breaking of background degeneracy at the level of the perturbations
abebe2012covariant ; abebe2015breaking .

The layout of the manuscript is as follows: In the following section, we
present the background cosmology in CG models. After reviewing the basic
Friedmann equations by unifying the exotic fluid with other components like
radiation and baryonic matter fluids, we present the numerical analysis of
deceleration parameter to examine the current features of the Universe in OCG
and GCG, MGCG and GCCG in Sub-Sec. II.1 \- II.4 respectively. The theoretical
distance modulus are investigated in CG models and computed and compared with
results of existing SNIa data in Sec. III. In Section IV, we also review the
basic spatial gradient variables and derive the linear evolution equations,
applying the scalar and harmonic decomposition techniques and obtain the wave-
number dependent energy density fluctuations for both CG and matter fluids. In
Sub-Sec. IV.3 \- IV.5, we present and discuss the analytical and numerical
solutions of density perturbations by considering baryonic-CG, radiation-CG
dominant Universe for both wave-length ranges and provide some concluding
remarks.

##  II Background cosmology in the CG model

The negative pressure associated with the CG model is related to a positive
energy density by a characteristic equation of state given as
kamenshchik2001alternative

$\displaystyle p=-\frac{A}{\rho^{\alpha}}\;,$ (1)

where $p$ is the pressure and $\rho$ is the energy density, both in a comoving
reference frame with $\rho>0$. $A$ and $\alpha$ are positive constants111The
values for $\alpha$ are given by CG model for $0<\alpha\leq 1$, and the
original OCG at $\alpha=1$ ccgmg .. We assume the energy density of the multi-
fluid universe to consist of CG, baryons and radiation. The conservation
equations for the ith fluid is given as

$\displaystyle\dot{\rho}_{i}+3\frac{\dot{a}}{a}(\rho_{i}+p_{i})=0\;,\quad\mbox{$i=\\{r,b,ch\\}$}.$
(2)

From the Eq. (1) the energy density of the exotic fluid reads

$\displaystyle\rho_{ch}(a)=\left[A+\frac{C_{2}}{a^{3(1+\alpha)}}\right]^{\frac{1}{1+\alpha}}\;,$
(3)

where $C_{2}=e^{C_{1}(1+\alpha)}$ with $C_{1}$ a constant of integration,
which means that $C_{2}$ is a positive constant. We can now look at equation
(3) and determine how the energy density of the CG evolves during different
epochs. In the limiting case where $\frac{C_{2}}{A}\gg a^{3(1+\alpha)}$ for
the early universe, the energy density depends on the scale factor as

$\displaystyle\rho_{ch}(a)=\frac{\sqrt{C_{2}}}{a^{3}}\;.$ (4)

From this, we can conclude that the CG corresponds to baryonic-like matter
which can be interpreted as the dark matter contribution to dust:
$\rho_{ch}\sim a^{-3}$. On the other hand, $\frac{C_{2}}{A}\ll
a^{3(1+\alpha)}$ at late times, the energy density avelino2003lambdacdm
becomes

$\displaystyle\rho_{ch}(a)\approx\pm\sqrt{A}\;,$ (5)

and this equation clearly shows that $\rho_{ch}\sim$ constant, and therefore
the CG corresponds to a cosmological constant quintessence . Based on Eqs. (4)
and (5), the exotic fluid CG act as dark matter and dark energy in the early
and late times 222We assume a flat Universe ($\Omega_{k}=0$, where
$\Omega_{k}$ is the density due to the curvature of space-time). We can thus
write the total energy density as the sum of the different energy density
components:

$\displaystyle\rho_{tot}(a)=\left[A+\frac{C_{2}}{a^{3(1+\alpha)}}\right]^{\frac{1}{1+\alpha}}+\rho_{r0}a^{-4}+\rho_{b0}a^{-3}\;,$
(6)

with $\rho_{b0}$ and $\rho_{r0}$ corresponding to the present-day values of
the baryonic and radiation energy densities, respectively. The total pressure
can also be given in a similar fashion:

$\displaystyle
p_{tot}(a)=-\frac{A}{\left[\rho_{ch}(a)\right]^{\alpha}}+\frac{1}{3}\rho_{r}(a)\;.$
(7)

The fluid equation for the total energy density then becomes

$\displaystyle\dot{\rho}_{tot}+3\frac{\dot{a}}{a}(\rho_{tot}+p_{tot})=0\;.$
(8)

We can also make a change of variables and use redshift, as the former is a
physically measurable quantity. We do this by substituting
$1+z=\frac{a_{0}}{a}$, where $a_{0}$ is the scale factor at the present time
which we can normalize to 1. The Friedmann equation consequently reads as

$\displaystyle 3H^{2}(z)$ $\displaystyle=$
$\displaystyle\left[A+C_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}$
(9) $\displaystyle+\rho_{r0}(1+z)^{4}+\rho_{b0}(1+z)^{3}\;,$

where $H=\dot{a}/a$ is the Hubble parameter. Dividing throughout by
$3H^{2}_{0}$ ($H_{0}$ being the value of the Hubble parameter today) in
equation (9), we obtain

$\displaystyle h^{2}(z)$ $\displaystyle=$
$\displaystyle\left[D_{1}+D_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}$
(10) $\displaystyle+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\;,$

where $h\equiv\frac{H}{H_{0}}$ is the normalized Hubble parameter, and
$D_{1}=\frac{A}{(3H_{0}^{2})^{(1+\alpha)}}$,
$D_{2}=\frac{C_{2}}{(3H_{0}^{2})^{(1+\alpha)}}$,
$\Omega_{r0}\equiv\frac{\rho_{r0}}{3H_{0}^{2}}$ and
$\Omega_{b0}=\frac{\rho_{b0}}{3H_{0}^{2}}$, are fractional (dimensionless)
density parameters. The total fractional density parameter for a flat universe
at any given redshift can be written as

$\displaystyle\Omega_{r}+\Omega_{b}+\Omega_{ch}=1\;,$ (11)

where here:

$\displaystyle\Omega_{r}\equiv\frac{\rho_{r}}{3H^{2}}=\frac{\Omega_{r0}(1+z)^{4}}{\left[D_{1}+D_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}}\;,$
(12)
$\displaystyle\Omega_{b}\equiv\frac{\rho_{b}}{3H^{2}}=\frac{\Omega_{b0}(1+z)^{3}}{\left[D_{1}+D_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}}\;,$
(13)
$\displaystyle\Omega_{ch}\equiv\frac{\rho_{ch}}{3H^{2}}=\frac{\left[D_{1}+D_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}}{\left[D_{1}+D_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}}\;.$
(14)

In the following subsections, we further analyze some cosmological
constraints, namely the deceleration parameter and luminosity distance to
explain cosmic history of the late time due to the presence of the CG fluid.
We then constrain these parameters by comparing the CG models with
$\Lambda$CDM to see which of the CG models best mimics like the $\Lambda$CDM
model.

N.B We will use the values all parameters from our MCMC simulations of the
best fitting parameters values for each tested model including the
$\Lambda$CDM as presented in Table 1.

###  II.1 Deceleration parameter in CG models

From the Friedmann equation, we derived the deceleration parameter for OCG and
GCG, and it yields as

$\displaystyle
q(z)=\frac{\left[A+C_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}-3A\left[A+C_{2}(1+z)^{3(1+\alpha)}\right]^{-\frac{\alpha}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+\rho_{b0}(1+z)^{3}}{2\left[A+C_{2}(1+z)^{3(1+\alpha)}\right]^{\frac{1}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+2\rho_{b0}(1+z)^{3}}\;.$
(15)

In the following, we present the numerical results of the above deceleration
parameters for different values of $\alpha$. For the case of $\alpha=1$ and
$\alpha\leq 1$, the results of OCG and GCG model respectively, as presented in
Figs. 1 \- 2.

Figure 1: A plot of the deceleration parameter, $q(z)$, versus the redshift,
$z$ for Eq. (15) in the OCG model, $\alpha=1$. Figure 2: The deceleration
parameter $q(z)$ versus $z$ for Eq. (15) in the GCG model for different values
of $\alpha<1$.

###  II.2 The deceleration parameter in MCG cosmology

Here, we consider the MCG that has an equation of state of the form
Elmardi:2018ryn

$\displaystyle p_{ch}=B\rho_{ch}-\frac{A}{\rho_{ch}^{\alpha}},$ (16)

where $A,~{}B$ and $\alpha$ are constant parameters. For the case of $A=0$,
the well-known equation of state for a perfect fluid is recovered, i.e.
$p=B\rho$. In addition for the case of $B=0$, the equation of state is reduced
to GCG, and then to OCG at $\alpha=1$. The energy density for this model
becomes as

$\displaystyle\rho_{ch}(a)=\left[A\Big{(}\frac{1}{(B+1)}+\mathcal{K}a^{-3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}},$
(17)

where $\mathcal{K}=\frac{C_{3}}{A}$, and $C_{3}$ is a constant of integration.
Therefore, we now have an equation for the energy density in terms of the
scale factor. The total energy density and total pressure are the quantities
that we are concerned with, so we would have to add the baryonic matter and
radiation components by hand (as was done in equations (6) and (7)) to obtain
the normalized Hubble parameter $h(z)$:

$\displaystyle h^{2}(z)$ $\displaystyle=$
$\displaystyle\left[E_{0}\Big{(}\frac{1}{(B+1)}+\mathcal{K}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}$
(18) $\displaystyle+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\;,$

where $E_{0}=\frac{A}{(3H_{0}^{2})^{1+\alpha}}$, and the deceleration
parameter can be written in terms of redshift dependent as

$\displaystyle q(z)$ $\displaystyle=$
$\displaystyle\frac{(3B+1)\left[A\Big{(}\frac{1}{(B+1)}+\mathcal{K}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}-3A\left[A\Big{(}\frac{1}{(B+1)}+\mathcal{K}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{-\frac{\alpha}{1+\alpha}}}{2\left[A\Big{(}\frac{1}{(B+1)}+\mathcal{K}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+2\rho_{b0}(1+z)^{3}}$
(19)
$\displaystyle+\frac{2\rho_{r0}(1+z)^{4}+\rho_{b0}(1+z)^{3}}{2\left[A\Big{(}\frac{1}{(B+1)}+\mathcal{K}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+2\rho_{b0}(1+z)^{3}}$

In the following, we present the numerical results of deceleration parameter
Eq. (19) for different values of of $\alpha$ and $B$ as presented in Figs. 3
\- 4. We consider the max and min values of $B$ from Table 1 for illustrative
purpose.

Figure 3: The deceleration parameter $q(z)$ versus $z$ for MCG mode Eq. (19)
for different values of $\alpha$ and $B=0.0428$. Figure 4: The deceleration
parameter $q(z)$ versus $z$ for MCG mode Eq. (19) for different values of
$\alpha$ and $B=0.2629$.

###  II.3 The deceleration parameter in MGCG cosmology

A more generalized version of the MCG model with some interesting properties
elmardi2016chaplygin ; bouhmadi2015scalar ; morais2015can  is the so-called
modified generalized CG (MGCG) model described by the equation of state

$\displaystyle p_{ch}=B\rho_{ch}-(1+B)\frac{A}{\rho_{ch}^{\alpha}}\;,$ (20)

with A, B, and $\alpha$ as the constant parameters of the model. With the
appropriate choice of these parameters, the model can be shown to have the
equation of state of both the cosmological constant and the MCG as its
limiting cases. The energy density of such a CG gas model is given by

$\displaystyle\rho_{ch}(a)=\left[A+C_{4}a^{-3(\alpha+1)(B+1)}\right]^{\frac{1}{1+\alpha}}\;,$
(21)

where $C_{4}$ is an integration constant. The normalized Hubble and
deceleration parameters of the model are therefore:

$\displaystyle\rho_{ch}(a)=\left[A\Big{(}1+\mathcal{K}_{2}a^{-3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}\;,$
(22)

where $\mathcal{K}_{2}=\frac{C_{4}}{A}\;.$ The normalized Hubble parameter
yield as

$\displaystyle h^{2}(z)$ $\displaystyle=$
$\displaystyle\left[F_{1}\Big{(}1+\mathcal{K}_{2}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}$
(23) $\displaystyle+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\;,$

and the deceleration parameter becomes

$\displaystyle q(z)$ $\displaystyle=$
$\displaystyle\frac{(1+3B)\left[F_{1}\Big{(}1+\mathcal{K}_{2}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}-3(1+B)A\left[F_{3}\Big{(}1+\mathcal{K}_{2}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{-\frac{\alpha}{1+\alpha}}}{2\left[\left(F_{1}\Big{(}1+\mathcal{K}_{2}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right]^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\right)}$
(24)
$\displaystyle+\frac{2\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}}{2\left[\left(F_{1}\Big{(}1+\mathcal{K}_{2}(1+z)^{3(\alpha+1)(B+1)}\Big{)}\right)^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\right]}\;.$

The numerical results of deceleration parameter is presented in Figs. 5 \- 6
with different values of $B$ and $\alpha$.

Figure 5: The decelerating parameter versus cosmological redshift $z$ for MGCG
model for different values of $\alpha$ at B = 0.2683. Figure 6: The
decelerating parameter versus cosmological redshift $z$ for MGCG model for
different values of $\alpha$ at B = 0.0462.

###  II.4 The deceleration parameter in GCCG cosmology

In this section we consider the generalized cosmic Chaplygin gas (GCCG) models
to present the numerical results of the deceleration parameter to explain the
cosmic acceleration. The equation state is given as gonzalez2003you

$\displaystyle
p_{ch}=-\rho^{-\alpha}_{ch}\left[C+\left(\rho_{ch}^{1+\alpha}-C\right)^{-\omega}\right]\;,$
(25)

where $C=\frac{A}{1+\omega}-1$, for $A$ a constant which can take on both
negative and positive numbers, and $0>\omega>-l$, and $l$ here is a positive
definite constant which can take on values larger than unity (Elmardi:2018ryn,
). The continuity equation of the exotic fluid reads as

$\displaystyle\dot{\rho}_{ch}+3\frac{\dot{a}}{a}\left[\rho_{ch}-\rho^{-\alpha}_{ch}\left(C+(\rho_{ch}^{1+\alpha}-C)^{-\omega}\right)\right]=0\;.$
(26)

From this equation, we have a total energy density and total pressure
(including the baryonic matter and radiation components) can be written
respectively as

$\displaystyle\rho_{tot}$ $\displaystyle=$
$\displaystyle\left[C+\left(1+Ba^{-3(\alpha+1)(\omega+1)}\right)^{\frac{1}{1+\omega}}\right]^{\frac{1}{1+\alpha}}$
(27) $\displaystyle+\rho_{r0}a^{-4}+\rho_{b0}a^{-3}$

and the normalized Hubble parameter is given as

$\displaystyle
h^{2}(z)=\frac{1}{3}\left[\frac{C}{H_{0}^{2(1+\alpha)}}+\frac{\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right]^{\frac{1}{1+\omega}}}{3H_{0}^{2(1+\alpha)}}\right)^{\frac{1}{1+\alpha}}+\Omega_{r0}(1+z)^{4}+\Omega_{b0}(1+z)^{3}\;.$
(28)

Then, the deceleration parameter for the GCCG model is

$\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-\frac{3\left[C+\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right)^{\frac{1}{1+\omega}}\right]^{-\frac{\alpha}{1+\alpha}}\left[C+\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right)^{-\frac{\omega}{1+\omega}}\right]}{2\left[C+\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right)^{\frac{1}{1+\omega}}\right]^{\frac{1}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+2\rho_{b0}(1+z)^{3}}$
(29)
$\displaystyle+\frac{\left[C+\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right)^{\frac{1}{1+\omega}}\right]^{\frac{1}{1+\alpha}}+\rho_{b0}(1+z)^{3}+2\rho_{r0}(1+z)^{4}}{2\left[C+\left(1+B(1+z)^{3(\alpha+1)(\omega+1)}\right)^{\frac{1}{1+\omega}}\right]^{\frac{1}{1+\alpha}}+2\rho_{r0}(1+z)^{4}+2\rho_{b0}(1+z)^{3}}\;.$

As we did in the previous sections we present the numerical results of Eq. 29
in the following Figs. 8 for the deceleration parameter in GCCG models. From
the plots, we observe that this model can explain the phantom phase of the
universe as presented in gonzalez2003you  for small values of $\alpha$ and $B$
and which less candidates than the other CG models to explain the accelerating
expansion universe. We use the average values $\omega$ from Table 1 for
numerical purpose.

Figure 7: $q(z)$ versus cosmological redshift $z$ for GCCG model for
$B=0.4170$. Figure 8: $q(z)$ versus cosmological red-shift $z$ for GCCG model
for $B=6209$.

From above all plots Fig. 1 \- 8 , we tried to show the accelerated expansion
within the observationally known range of cosmological redshift $z$. All
results are consistent with what one would expect the expansion history of the
Universe to look like: in the current concordance model, there should be a
decelerating matter-dominated phase followed by a dark-energy-driven late time
acceleration phase except the GCCG model. However the numerical results of
GCCG rule out this model, and its further statistical analysis is investigated
in Sec. III.5.

##  III Supernovae cosmology for CG models

###  III.1 The $\Lambda$CDM model

The distance modulus that can be obtained by combining the different
cosmological distance definitions, found in deza2009encyclopedia , is given by
333This distance modulus is given in terms of Mpc.

$\mu=m-M=25-5\times\log_{10}\left[3000\bar{h}^{-1}(1+z)\int^{z}_{0}\frac{dz^{\prime}}{h(z^{\prime})}\right]\;,$
(30)

where $m$ and $M$ are the apparent and absolute magnitudes for a particular
measured supernova. We also used the definition for the Hubble uncertainty
parameter, which is given as $\bar{h}=\frac{H(z)}{100\frac{km}{s.Mpc}}$, with
$H(z)$ being the Hubble parameter. We obtained the normalized Friedmann
equation $h(z)$, as a result from the definition for the line-of-sight co-
moving distance deza2009encyclopedia .

Now that we have a model to fit against the supernovae data, we will use 359
low- and intermediate-redshift data points obtained from the SDSS-II/SNLS3
Joint Light-curve Analysis (JLA). This will be sufficient to test the
transition phase between the matter dominated epoch (a decelerating expansion)
and the late-time dark energy epoch (an accelerating expansion), which occurs
at $z\approx 0.5$ linder2001understanding ; capozziello2019model . We will be
using the absolute magnitudes for these supernovae in the B-filter that were
calculated in the research papers conley2010supernova ; Neill2009 ; Hicken2009
, and are represented by blue, green and magenta colour dots in the supernovae
best-fit figures444These absolute magnitude values can also be found on NED..

A question that may arise now is: “What about the degeneracy between the
Hubble parameter and the absolute magnitude of the supernovae?”. Strictly
speaking, since there exists a degeneracy between these two parameters, and
usually the absolute magnitude is unknown, we cannot constrain the parameters
using the supernovae data colgain2019hint . In fact, normally you would have
to use either a fixed absolute magnitude value for the supernovae (as done by
colgain2019hint , since it gives a robust enough prediction that the other
free parameters, excluding the Hubble constant, can still be constrained), or
you will need to use Cepheid variable star data, as done by freedman2001final
; riess2009redetermination ; riess20113 , where they determined the error on
the distance modulus created by the uncertainty of the absolute magnitude, and
add this correction to the supernovae’s distance modulus.

However, in our case it is not necessary, due to the fact that we used the
absolute magnitudes as determined by Hicken2009 ; Neill2009 ;
conley2010supernova . In their work, they already incorporated the corrections
to the absolute magnitudes when they calculated it. In fact, the corrections
they made were specifically made to enable them to constrain the Hubble
constant on Supernovae Type 1A data. Therefore, our distance modulus also
incorporates these corrections. Furthermore, we use a Markov Chain Monte Carlo
(MCMC) simulation that takes into account the error on our the distance
modulus, when calculating the best-fitting parameters. These errors include,
the systematics errors that were obtained when acquiring the data (the error
on the apparent magnitude found in the JLA dataset) and the errors on the
absolute magnitudes as given by Hicken2009 ; Neill2009 ; conley2010supernova .
Therefore, we can use our calculated distance modulus as is.

Now that we have a model and the accompanying dataset, we want to compare how
the different CG models fare against the concordance model, namely the
$\Lambda$CDM model. We will start with the best-fitting model for the
$\Lambda$CDM model, to get an idea on what our “true model” will look like. By
assuming $\Omega_{k}=0$ (therefore, a flat universe), as well as making the
substitution $\Omega_{\Lambda}=1-\Omega_{m0}-\Omega_{r0}$, we find the
normalized Friedmann equation for the $\Lambda$CDM model as

$h(z)=\sqrt{\Omega_{r0}\left(1+z\right)^{4}+\Omega_{m0}\left(1+z\right)^{3}+1-\Omega_{m0}-\Omega_{r0}}\;.$
(31)

We can then use our MCMC simulation as mentioned, that were developed in
hough2020constraining  and later used in hough2019constraining ;
hough2020viability , spesifically for the testing of modified gravity models
on supernovae data. This MCMC simulation uses the Metropolis-Hasting algorithm
and for simplicity, uses a Gaussian probability distribution to calculate the
likelihood of each best-fitting parameter value depending on the distance
modulus. This code incorporates the build in Python package EMCEE Hammer
foreman2013emcee  to perform the MCMC simulation. This build-in package uses
different random walkers, each executing the MH algorithm and all starting at
the same initial parameter values, following their own path to determine the
most probable parameter values. This creates a Gaussian distribution based on
each random walker’s ending parameter values. Using the average values for
each probability distribution for each parameter, we will on average obtain
the best-fitting parameter values and its $1\sigma$-deviation
hough2020constraining ; hough2020viability .

Executing this MCMC simulation for the $\Lambda$CDM model, by combining Eqs.
(30) and (31), we find on average the best-fitting parameter value for each
free parameter to be $\bar{h}=0.6965^{+0.0048}_{-0.0048}$ (constrained) for
the Hubble uncertainty parameter, $\Omega_{m0}=0.2668^{+0.0250}_{-0.0241}$
(constrained) for the matter density parameter, and
$\Omega_{r0}=0.00050^{+0.00034}_{-0.00034}$ (unconstrained) for the radiation
density parameter. The MCMC simulation results are shown in Fig. 9.

Figure 9: The MCMC simulation results for the $\Lambda$CDM model’s (Eq. 31)
cosmological free parameters ($\Omega_{m0}$, $\bar{h}$, and $\Omega_{r0}$),
with the “true” values (blue lines: $\Omega_{m0}=0.315$, $\bar{h}=0.674$, and
$\Omega_{r0}=2.47\times 10^{-5}/\bar{h}^{2}$) provided by the Planck2018
collaboration data release aghanim2020planck . We used 100 random walkers and
25 000 iterations.

It must be noted that not all of the parameters are not always fully
constrained. This can be due to the amount of data or the actual simulation
running into the numerical errors that occurs when approximating the values.
There are mainly 3 ways to define how they are constrained. The first is, the
parameter is fully constrained. This occurs when a parameter was able to
create a full Gaussian probability distribution in the available parameter
search space. The second is, the parameter is unreliably constrained. This
occurs when only one side of the Gaussian distribution was found in the search
space. This means that there is a value that is more desirable than the one it
received, since the average is slightly shifted from the peak of the
distribution. Lastly, the parameter is unconstrained, which occurs when a
uniform distribution was obtained for the particular parameter. This means
that according to the MCMC simulation, this parameter does not play a role in
the best-fit (at least on this numerical resolution) for a particular model.
Using this knowledge and the results from the MCMC simulation for the
$\Lambda$CDM model, we can plot the best-fitting model on the data and
obtaining Fig. 10.

Figure 10: The $\Lambda$CDM model’s Eq. (31) best-fitting free parameters for
the Supernovae Type 1A data with cosmological parameter values as
$\bar{h}=0.6965^{+0.0048}_{-0.0048}$ (constrained),
$\Omega_{m0}=0.2668^{+0.0250}_{-0.0241}$ (constrained), and
$\Omega_{r0}=0.00050^{+0.00034}_{-0.00034}$ (unconstrained) as calculated by
the MCMC simulation (top panel). The bottom panel shows the residuals distance
in Mpc between the predicted model values and the data points.

From Fig. 10 on the top panel, as expected, we see that the $\Lambda$CDM model
fit the data extremely well. Even the $1\sigma$-deviation does not really have
an impact on the full range predicted by the $\Lambda$CDM model. Furthermore,
on the bottom panel, where the residuals are shown, we can see that at no
point does the model over-or under-estimate the resulting distance modulus for
each supernovae. We also note that the average off-set the model has, compared
to the data, is $\bar{x}_{res}=-0.0373$ Mpc with a standard deviation of
$\sigma_{res}=0.21477$, showing that this is a thigh relation between the
model and the data points.

As for comparing this results to other results to determine whether or not the
model and simulation does indeed work, we saw that our results are not within
$1\sigma$ from the Hubble constant parameter value that were determined by
Planck2018 collaboration aghanim2020planck . However, this is actually
expected, since it is known that there exists a discrepancy between the
predicted Hubble constant values, due to studies using the CMB data, and the
studies using the supernovae data. It is noted in the research paper
riess20162 , that the supernovae related studies tend to find a higher
predicted value than the CMB related studies. It is still uncertain why this
discrepancy exists. When we compare our results to other supernovae studies,
we do see that we are within $1\sigma$ from their results riess20162 .
Indicating that our MCMC simulation is indeed working as intended.

###  III.2 The various CG model results

Now that we have a way of testing the various different CG models against
observational data, we can go ahead and use each different CG model’s
normalized Friedmann equation in our MCMC simulation. Since all of the MCMC
simulation results relatively look the same and take up a lot of space when
plotting, we will just be providing the necessary results, such as the most
probable best-fitting parameter values, and whether or not the parameter was
constrained. We start with the original CG model. The only difference the
original CG model has compared to the general CG model is $\alpha=1.0$ in Eq.
(10). We can then insert equation (10) into the distance modulus function,
with $\alpha=1.0$, and execute the MCMC simulation. We plot the best-fitting
model on the data in Fig 11.

Figure 11: The original CG model’s Eq. (10), with $\alpha=1.0$, best-fitting
free parameters for the Supernovae Type 1A data with cosmological parameter
values obtained to be $\bar{h}=0.6489^{+0.0712}_{-0.0363}$ (unreliably
constrained), $\Omega_{b0}=0.0447^{+0.0358}_{-0.0303}$ (unreliably
constrained), $\Omega_{r0}=0.00507^{+0.00333}_{-0.00344}$ (unconstrained),
$D_{1}=1.0544^{+0.2861}_{-0.3699}$ (constrained), and
$D_{2}=0.1678^{+0.0656}_{-0.0654}$ (constrained) as calculated by the MCMC
simulation (top panel). The bottom panel shows the residuals distance in Mpc
between the predicted model values and the data points.

From Fig. 11, we see that the original CG model does explain the supernovae
data. This is evident on the bottom panel, since the model does not over-or
under-estimate the supernovae at the low- or intermediate-redshift with the
average off-set being $\bar{x}_{res}=-0.0313$ Mpc, with a standard deviation
on this average of $\sigma_{res}=0.21385$ (which is even smaller than the
$\Lambda$CDM model). However, even though the matter-density and Hubble
uncertainty parameter gave realistic values that are within $1\sigma$ from
their expected values, but they are not fully constrained (leading to the
large error region).

Now that we have seen that the OCG model does indeed show promise in being
able to explain the supernovae data, we can try to fit the general CG model,
where $\alpha$ is no longer fixed, but is a free parameter with boundary
values of $0<\alpha\leq 1.0$. We can then once again insert Eq. (10) into our
distance modulus model, execute the MCMC simulation and plot the results on
the supernovae data in Fig. 12.

Figure 12: The general CG model’s Eq. (10) best-fitting free parameters for
the Supernovae Type 1A data with cosmological parameter values obtained to be
$\bar{h}=0.6575^{+0.0722}_{-0.0424}$ (unreliably constrained),
$\Omega_{b0}=0.0477^{+0.0349}_{-0.0320}$ (unconstrained),
$\Omega_{r0}=0.00050^{+0.00034}_{-0.00034}$ (unconstrained),
$\alpha=0.6318^{+0.2609}_{-0.3633}$ (unreliably constrained),
$D_{1}=0.9160^{+0.2530}_{-0.2567}$ (constrained), and
$D_{2}=0.1952^{+0.0683}_{-0.0661}$ (constrained) as calculated by the MCMC
simulation (top panel). The bottom panel shows the residuals distance in Mpc
between the predicted model values and the data points.

Fig. 12 once again shows that the CG models can explain the supernovae data,
with realistic parameter values for the cosmological parameter, albeit
unreliably constrained, as well as explaining the supernovae data without any
over- or under-estimation on the low- or intermediate-redshift data. Even the
average under-estimation on the model on the data, namely
$\bar{x}_{res}=-0.0282$ Mpc, with a standard deviation on this average of
$\sigma_{res}=0.21455$, which is on average closer to the actual data point
values than even the OCG model (that were closer than the $\Lambda$CDM model),
although the deviation for this model on the data is more in-line with the
$\Lambda$CDM model’s results (with the OCG model faring better in this case).
However, as with the OCG model, we need to be careful when using these
unreliable constraints. At least in the OCG model’s case, both are unreliably
constrained, which shows that it is closer to the actual best-fitting
parameter value, than being unconstrained, such as with the GCG model.

It is though interesting that the $D_{1}$ and $D_{2}$ free parameter values,
are the parameters that the MCMC simulation gave as the dominant parameters
when trying to fit the model. Therefore, it would seem as if $D_{1}$ and
$D_{2}$ can up to a point determine the shape of the function, with the other
parameters only being used to “fine tune” the best-fitting model. This raises
another questions: since $D_{1}$ and $D_{2}$ are related to the free
parameters $A$, $C_{1}$, and $H_{0}$, is it the two free parameters being the
dominant parameters, or is it the $H_{0}$. In the case of the second, would it
then not result in a constrained parameter value for $\bar{h}$, or is this
parameter not needed due to using the Hubble parameter in $D_{1}$ and $D_{2}$,
therefore, with them already constrained, is the Hubble parameter more
reliable than the MCMC simulation predicts?

Moving on to the modified CG model (MCG), we can find the best-fitting
parameters using Eq. 18. We used the boundary conditions $0<B\leq 1/3$ for the
$B$ parameter. The MCMC results for the MCG model is plotted on the supernovae
data and is shown in Fig. 13.

Figure 13: The modified CG model’s Eq .(18) best-fitting free parameters for
the Supernovae Type 1A data with cosmological parameter values obtained to be
$\bar{h}=0.6779^{+0.0772}_{-0.0550}$ (unreliably constrained),
$\Omega_{b0}=0.0484^{+0.0349}_{-0.0333}$ (unconstrained),
$\Omega_{r0}=0.00050^{+0.00033}_{-0.00035}$ (unconstrained),
$\alpha=0.5694^{+0.2939}_{-0.3312}$ (unreliably constrained),
$B=0.1448^{+0.1181}_{-0.1020}$ (unreliably constrained),
$E_{0}=1.0026^{+0.3388}_{-0.2818}$ (constrained), and
$\mathcal{K}=0.1305^{+0.0638}_{-0.0465}$ (constrained) as calculated by the
MCMC simulation (top panel). The bottom panel shows the residuals distance in
Mpc between the predicted model values and the data points.

From Fig. 13, we can see that the modifications to the general CG model holds
up, since it was able to find an average under-estimation on the residual
between the model and the supernovae data points of only
$\bar{x}_{res}=-0.0296$ Mpc, which has a difference of $0.0014$ Mpc further
away relative to the GCG model, but has a tighter relationship with the data,
with a standard deviation of $\sigma_{res}=0.21426$. We also once again note
that this model has no under- or over-estimation on the data at any red-shift
It also worth it to take note that the only parameters that were fully
constrained were the free parameters, $E_{0}$ and $\mathcal{K}$, enhancing our
suspicion that the free parameters determine the general shape of the
function, while the other parameters are only there for the fine adjustments
to make it more accurate.

Moving on to the next CG model, namely the modified generalised CG model
(MGCG) given by eq. 23. The MCMC results for the MGCG model is plotted on the
supernovae data and is shown in Fig. 14.

Figure 14: The modified generalized CG model’s Eq. (23) best-fitting free
parameters for the Supernovae Type 1A data with cosmological parameter values
obtained to be $\bar{h}=0.6722^{+0.0766}_{-0.0514}$ (unreliably constrained),
$\Omega_{b0}=0.0473^{+0.0347}_{-0.0322}$ (unconstrained),
$\Omega_{r0}=0.00049^{+0.00035}_{-0.00033}$ (unconstrained),
$\alpha=0.5992^{+0.2770}_{-0.3571}$ (unreliably constrained),
$B=0.1493^{+0.1190}_{-0.1031}$ (unreliably constrained),
$F_{1}=0.9033^{+0.2787}_{-0.2591}$ (constrained), and
$\mathcal{K}_{2}=0.1460^{+0.0607}_{-0.0456}$ (constrained) as calculated by
the MCMC simulation (top panel). The bottom panel shows the residuals distance
in Mpc between the predicted model values and the data points.

From Fig. 14, it is clear that the MGCG model fits the supernovae data as we
have seen with the other CG models. Furthermore, we once again see that the
cosmological parameters do not play a big role in determining the outcome of
the model, which is not useful for determining physical attributes of the
Universe. We also see that there is no over- or under-estimation of the
observational data at low- or intermediate-redshifts with the average off-set
being $\bar{x}_{res}=-0.0281$ Mpc, with a standard deviation of
$\sigma_{res}=0.21427$ (which is in-line with the results from the MCG model),
therefore confirming that this is the best possible fit with our MCMC
simulation resolution.

We can now go over to our last model’s Hubble parameter equation, namely the
generalized cosmic CG (GCCG) model, given by Eq. 28. At the start of finding
this models best-fitting parameter values, we struggled to find working
boundary conditions. This happened due to the fact that the MCMC simulation
tries different combinations of parameter values to find on average the best-
fitting parameter values. It ended up running into logarithmic errors at the
boundaries for certain combinations. This points to the fact that one small
change in a parameters current value can lead to instabilities. In the end, we
did managed to find plausible parameter values after subjecting the model to
some tight constrains, mostly on the $\omega$ and the $B$-parameters. These
manually inserted constrains are $-0.90<\omega<-0.980$ and $0.01<B<2$. The
MCMC simulation results (with these manually constrained parameter values) is
shown in Fig. 15.

Figure 15: The generalized cosmic CG model’s Eq. (28) best-fitting free
parameters for the Supernovae Type 1A data with cosmological parameter values
obtained to be $\bar{h}=0.7064^{+0.0655}_{-0.0714}$ (unconstrained),
$\Omega_{b0}=0.0559^{+0.0307}_{-0.0359}$ (unconstrained),
$\Omega_{r0}=0.00048^{+0.00036}_{-0.00033}$ (unconstrained),
$\alpha=0.4763^{+0.3042}_{-0.2534}$ (constrained),
$\omega=-0.9710^{+0.0069}_{-0.0052}$ (constrained),
$A=1.1358^{+2.6400}_{-2.6972}$ (unconstrained), and
$B=0.5150^{+0.1059}_{-0.0980}$ (constrained) as calculated by the MCMC
simulation (top panel). The bottom panel shows the residuals distance in Mpc
between the predicted model values and the data points.

From Fig. 15, it is clear that this model is very unstable (since our
1$\sigma$-deviation is very large) as predicted, when we started testing out
various different boundary conditions. Our first instinct, when we saw this
large 1$\sigma$-deviation, was to look at the free parameter A, since we did
find constrained values for $\omega$ and $B$, while A was left unconstrained
with very large errors on it. However, this turned out not to be the problem.
The largest contribution to the large standard deviation was in-fact the
$\omega$ parameter, even though it were constrained with very small errors. It
turns out that this model is so sensitive to the $\omega$ parameter, that any
small changes to this parameter’s value alone can change the entire outcome of
this model. As you can see from Fig. 15, the best-fitting model on average
(Red line), is close to fitting the data with a high accuracy, as seen with
the other Chaplygin gas models. However, if for the time being we accepted
that all of the best-fitting values are correct and fix them to these values,
with the exception of the $\omega$ parameter, we can illustrate how unstable
this model is. Let’s say the $\omega$ parameter has a best-fitting value of
$\omega=-0.9650$. Note that this value is still within the 1$\sigma$-error
from the MCMC results, since we only added $0.006$ to the original best-
fitting value. If we then re-plot this average best-fitting model on the data,
we obtain Fig. 16.

Figure 16: The is an example to show how unstable the GCCG model is. We kept
all of the best-fitting values the same as in Fig. 15, with the exception of
making a minuscule change to the $\omega$ parameter best-fitting value that is
still within in its 1$\sigma$-error received from the MCMC results.

From Fig. 16, it is very clear that even a minuscule change in just one
parameter ended up changing the entire outcome of this model. Therefore,
confirming our initial findings that the parameters in this model can lead to
instabilities forming. This is not a good attribute for a model, since this
model can then be manipulated into any form that is necessary for a given
problem.

However, if we continue with this model’s results from Fig. 15, we do note
that we once again obtained the cosmological parameters to be closely related
to the values found in other supernovae related results riess20162 . Another
value to take note of is the $\omega$ parameter. Even though we had to
manually set the boundary conditions from the $\omega$ parameter to avoid
logarithmic errors. It is interesting that the only range and the eventual
best-fitting value for this parameter does correspond to the accepted dark
energy equation of state of $-1$, as well as being within the
1$\sigma$-deviation of the best-fitting value found by the Planck2018
collaboration on the $\Lambda$CDM model, namely
$\omega_{0}=-1.03^{+0.1}_{-0.1}$ aghanim2020planck . Although they combined
different data sets to obtain a more accurate value. This confirms that the
model was done correctly, since we can reproduce the known parameter values,
even though the model itself still remains unstable.

We can also confirm this from our residuals plot in Fig. 15, since we did not
obtain any over-or under estimation on a specific section of the data,
therefore matching the transition phase going from the matter dominated epoch,
to an accelerating Universe epoch. However, we do have a slight under-
estimation on the distance to the supernovae with an average of
$\bar{x}_{res}=-0.2231$ Mpc. The standard deviation on this average best-fit
over the entire redshift range is $\sigma_{res}=0.21604$, which is comparable
to the previous models. However, this is a larger under-estimation than the
previous tested models, but this is probably due to the instabilities created
by the $\omega$ parameter. We summarise all of the different parameter values
for each model in Tab. 1.

Table 1: The MCMC simulation best fitting parameters values for each tested model, including the $\Lambda$CMD model. The models are listed in the order that they were tested. NB! Due to limitations the following groupings are done: 1) The $\Lambda$CDM’s $\Omega_{m0}$ parameter is grouped in the $\Omega_{b0}$ column. 2) Due to being similar free parameters $D_{1}$, $E_{0}$ and $F_{1}$ are all grouped together, while $D_{2}$, $\mathcal{K}$ and $\mathcal{K}_{2}$ is grouped together. 3) The B-parameter for the MCG and MGCG models are limited to a maximum of $B=\frac{1}{3}$ due to the definitions of these models, while the B-parameter from the GCCG model is only an integration constant, therefore even though their names are the same, they are inherently different. Model | $\bar{h}$ | $\Omega_{b0}/\Omega_{m0}$ | $\Omega_{r0}$ | $D_{1}/E_{0}/F_{1}$ | $D_{2}/\mathcal{K}/\mathcal{K}_{2}$ | $\alpha$ | $B$ | $\omega$
---|---|---|---|---|---|---|---|---
$\Lambda$CDM | $0.6965^{+0.0048}_{-0.0048}$ | $0.2668^{+0.0250}_{-0.0241}$ | $0.00050^{+0.00034}_{-0.00034}$ |  |  |  |  |
Original CG | $0.6489^{+0.0712}_{-0.0363}$ | $0.0447^{+0.0358}_{-0.0303}$ | $0.00507^{+0.00333}_{-0.00344}$ | $1.0544^{+0.2861}_{-0.3699}$ | $0.1678^{+0.0656}_{-0.0654}$ |  |  |
General CG | $0.6575^{+0.0722}_{-0.0424}$ | $0.0477^{+0.0349}_{-0.0320}$ | $0.00050^{+0.00034}_{-0.00034}$ | $0.9160^{+0.2530}_{-0.2567}$ | $0.1952^{+0.0683}_{-0.0661}$ | $0.6318^{+0.2609}_{-0.3633}$ |  |
Modified CG | $0.6779^{+0.0772}_{-0.0550}$ | $0.0484^{+0.0349}_{-0.0333}$ | $0.00050^{+0.00033}_{-0.00034}$ | $1.0026^{+0.3388}_{-0.2818}$ | $0.1305^{+0.0638}_{-0.0465}$ | $0.5694^{+0.2939}_{-0.3312}$ | $0.1448^{+0.1181}_{-0.1020}$ |
Modified GCG | $0.6722^{+0.0766}_{-0.0514}$ | $0.0473^{+0.0347}_{-0.0322}$ | $0.00049^{+0.00035}_{-0.00033}$ | $0.9033^{+0.2787}_{-0.2591}$ | $0.1460^{+0.0607}_{-0.0456}$ | $0.5992^{+0.2770}_{-0.3571}$ | $0.1493^{+0.1190}_{-0.1031}$ |
GCCG | $0.7064^{+0.0655}_{-0.0714}$ | $0.0559^{+0.0307}_{-0.0359}$ | $0.00048^{+0.00036}_{-0.00033}$ | $1.1358^{+2.6400}_{-2.6972}$ |  | $0.4763^{+0.3042}_{-0.2534}$ | $0.5150^{+0.1059}_{-0.0980}$ | $-0.9710^{+0.0069}_{-0.0052}$

###  III.3 Comparing the residuals of the various CG models to the
$\Lambda$CDM’s

Now that we have all of the models that we are going to work with, we can
compare each of their best-fitting residual plots to the $\Lambda$CDM model’s,
to understand how each of them handles the transition phase between the matter
dominate epoch ($z>0.5$) and the accelerating epoch that we are currently
experiencing ($z<0.5$) linder2001understanding ; capozziello2019model . We
visualize this comparisons in Fig. 17

Figure 17: Theoretical residuals comparing the different tested models against
the $\Lambda$CDM model.

From Fig. 17, we learn a few interesting details about the behavior of the
various Chaplygin gas models. Firstly, with the exception of the unstable GCCG
model, all of the Chaplygin gas models have a very similar residuals plot when
compared to the $\Lambda$CDM model, and even then the GCCG model, follows the
same trend just being a mirror image of the other CG models. We also note that
the off-set for the GCCG model is larger compared to the other models, but on
the relative size that these models predict the distance to the various
supernovae, even the GCCG model performs well to a certain extend compared to
the $\Lambda$CDM model’s prediction.

Secondly, we notice that at the transition point between the matter dominated
epoch going over to an accelerating epoch, all five CG models change their
relative behavior, and start matching the trend produced by the $\Lambda$CDM
model. This is quite a distinction for these models, since it means that they
can explain an accelerating Universe with the best-fitting parameter values
that were obtained from observational data. This is an important piece of
information to take note off, since these models do not only have a
theoretical backing to explain an accelerating Universe, but have some support
from observations.

Thirdly, it is interesting to note that as time progress, the four stable CG
models converge onto the $\Lambda$CDM model, with the OCG matching the very
late time distance predictions of the $\Lambda$CDM model almost exactly. All
four of these stable models end up within a distance modulus prediction of
about $\Delta\mu\sim 0.001$ Mpc. Which is relative very small on the scale of
these observations. To put this in perspective, the distance between the Sun
and the center of the Milky Way is estimated at $\sim 0.008$ Mpc
reid1993distance . Therefore, these models can all predict the distance to
supernovae happening in various different galaxies, to a distance less than
the distance to the center of our galaxy relative to the $\Lambda$CDM model’s
prediction. This is quite an admirable result.

Lastly, although some of the models do quite well at late time, none of them
succeeds with following the matter dominated epochs trend set by the
$\Lambda$CDM model. This is worrisome, since a model that can explain both of
the epochs, as well as the transition phase between the two will be more
useful in explaining the various cosmological phenomenons. For instance, if we
were to use higher redshift supernovae data, which is very difficult to obtain
due to the incoming flux decreasing over distance leibundgut2001cosmological ,
these models will struggle more and more relative to the $\Lambda$CDM model.
The reason for this is because they are already diverging away from the
$\Lambda$CDM model at the redshifts that we were able to obtain.

As an interesting extra test that we did, we saw that the all the model’s
relative difference distance modulus begin to plateau at about $z\sim 10$, and
start to match the trend set by the $\Lambda$CDM model once again albeit with
a difference ranging anywhere between $\Delta\mu\sim-0.8$ Mpc for the GCCG
model to $\Delta\mu\sim 0.7$ Mpc for the MCG and MGCG models, with the GCG
model having the smallest relative difference between the $\Lambda$CDM model
and itself with a value of about $\Delta\mu\sim 0.2$ Mpc. Therefore, even
though these models differ in with the values predicted by the $\Lambda$CDM
for earlier times, they do still follow the same trend after $z\sim 10$ and
find predictions with a relative distance modulus difference of less than $\pm
1$ Mpc.

###  III.4 Comparing the deceleration parameters of the various CG models to
the $\Lambda$CDM’s

As done with comparing the of the best-fitting residuals plots, we are able to
compare the behaviors of the deceleration parameter’s from each of the
different tested models against the $\Lambda$CDM plot. As shown in Sec. II,
different combinations of the parameter values can lead to different ways in
which the Universe can be shaped. This is important since the rate at which
the Universe is changing, cannot only tell us about the history of the shape
of the Universe, but it can also enable us to make predictions about the
future fate of the Universe. So using the best-fitting parameter obtained from
the supernovae observations, we will be able to compare how each model
predicts this rate of change compared to the accepted shape, predicted by the
$\Lambda$CDM model also using the same dataset and MCMC simulation. This is
shown in Fig. 18.

Figure 18: The deceleration parameters from each of the different tested
models compared to the $\Lambda$CDM model basd on the best-fitting parameter
values as predicted by the MCMC simulation on the supernovae data.

From Fig. 18, we can once again note a few important points about the
behaviour of these CG models based on their best-fitting parameter values.
Firstly, we can see that the same four Chaplygin gas models as mentioned in
the previous section, obtained deceleration parameter revolutions that are
very similar to one another, in particular the two pairings of OCG/GCG models
and the MCG/MGCG models. The only model that failed in this regard was the
GCCG model. This is probably due to its instabilities that we mentioned
earlier. According to the shape of the GCCG model, the Universe on this time
frame have always been in a state of acceleration and that in the future
($z\sim-1$), the Universe will rapidly accelerate in a very short amount of
time. This might only be an prediction, but based on the fact that this model
is the only model where we obtained an accelerating Universe over the entire
time frame, it might mean that this result is not correct.

Secondly, while focusing on the other more successful models, we see that the
MGCG model, matches on average the $\Lambda$CDM model’s trend the best,
however, it must be noted that at the present day ($z=0$), its the OCG and GCG
models that came the closest to the rate at which the size of the Universe
expands as predicted by the $\Lambda$CDM model, while the MCG and MGCG models
predicted that the Universe is expanding at a rate of about $\sim 25\%$ faster
than the $\Lambda$CDM model.

Thirdly, we note that these four models and the $\Lambda$CDM model predicts
the transitioning period between the matter dominated epoch and the
accelerating epoch to be between $0.4<z<0.8$, with the MGC model being the
closest to the $z\sim 0.5$ value found in the literature
linder2001understanding ; capozziello2019model . Remember that this is based
on our results from the MCMC simulations just to test the viability of these
models, which can differ from the major experimental papers that constrains
these parameters by combining various different dataset. This is just to
observe whether or not these models do follow the accepted trends when tested
with observational data.

###  III.5 Statistical analysis

From a statistical analysis stand point, we can determine the small details
that might have been missed in the visual comparisons to the $\Lambda$CDM
model. The statistical analysis test that we will be using is the Akaike
information criterion (AIC) and Bayesian/Schwarz information criterion (BIC)
selections. These information criteria evaluate the plausibility of an
alternative model explaining the data compared to an “accepted/true model”. In
our case the $\Lambda$CDM model will be considered as the “true model”. These
selection take into account the best-fitting likelihood function values of the
various models, but then also punishes the model based on the amount of free
parameters it uses. The reason for this is due to the fact that with more free
parameters, it is easier to manipulate a model into any shape you need,
therefore making it easier to fit the model to any dataset. This means that a
model that can explain the data to the same accuracy with less free parameters
that can be manipulated, would be more likely to be statistically accepted
than the models with lots of free parameters. The values for each model’s AIC
and BIC selections are calculated using the following equations:

$\begin{split}\bullet\quad AIC&=\chi^{2}+2K,\\\ \bullet\quad
BIC&=\chi^{2}+K\log(n),\end{split}$ (32)

where $\chi^{2}$ is calculated using the model’s Gaussian likelihood function
$\mathcal{L}(\hat{\theta}|data)$ value, $K$ is the amount of free parameters
for the particular model, while $n$ is the amount of data points in our
dataset. Since the value of $\chi^{2}$ can be any positive values based on the
number of data points, the calculated values for the AIC and BIC can by very
random. Therefore, we will have to use the difference in AIC and BIC values of
each model compared to the “true model’s AIC and BIC values. This statistical
test is also used in other similar research papers nunes2017new . To make
sense of the information criteria difference values between the CG models (aka
the alternative models) and the $\Lambda$CDM model, we will use the Jeffrey’s
scale in order to make conclusions about the viability of the various CG
models. It should be noted that this is not an exclusive scale and should be
handled with care nesseris2013jeffreys . The Jeffrey’s scale ranges are:

$\begin{split}\bullet\quad&\Delta IC\leq 2-\textrm{substantial observational
support},\\\ \bullet\quad&4\leq\Delta IC\leq 7-\textrm{less observational
support},\\\ \bullet\quad&\Delta IC\geq 10-\textrm{no observational
support}.\\\ \end{split}$ (33)

Lastly, since we will already have the calculated $\chi^{2}$-values, we will
also be able to calculate a reduced $\chi^{2}_{red.}$-value for each tested
model. This will allow us to mathematically compare the over-all best-fits,
namely the goodness of fit, for each model on the supernovae data. All of the
calculated statistical analysis values for each model, including the
$\Lambda$CDM model is shown Tab. 2.

Table 2: The best-fit for each tested model, including the $\Lambda$CMD model. The models are listed in the order from the largest likelihood function value $\mathcal{L}(\hat{\theta}|data)$ to the smallest likelihood of being viable. The reduced $\chi^{2}$-values are given as an indication of the goodness of fit for a particular model. The AIC and BIC values are shown, as well as the $\Delta IC$ for each information criterion. The $\Lambda$CDM model is chosen as the “true model”. Model | $\mathcal{L}(\hat{\theta}|data)$ | $\chi^{2}$ | Red. $\chi^{2}$ | $AIC$ | $|\Delta AIC|$ | $BIC$ | $|\Delta BIC|$
---|---|---|---|---|---|---|---
Original CG | -120.1716 | 240.3432 | 0.6789 | 250.3432 | 2.9485 | 269.7598 | 10.7152
Modified CG | -120.3719 | 240.7438 | 0.6839 | 254.7438 | 7.3491 | 281.9271 | 22.8825
Modified Generalised CG | -120.5274 | 241.0548 | 0.6848 | 255.0548 | 7.6601 | 282.2381 | 23.1935
General CG | -120.6764 | 241.3527 | 0.6837 | 253.3527 | 5.9580 | 276.6527 | 17.6081
$\Lambda$CDM | -120.6973 | 241.3947 | 0.6781 | 247.3947 | 0 | 259.0446 | 0
Generalized Cosmic CG | -281.7253 | 563.4506 | 1.6007 | 577.4506 | 330.0559 | 604.6338 | 345.5892

From Tab. 2, we can take note of a few interesting pieces of information.
Firstly, we found four different CG models, that obtained a better likelihood
function value than the $\Lambda$CDM model based on a Gaussian probability
distribution, with the OCG model obtaining the largest likelihood function
value. However, in the reduced $\chi^{2}$-values, where the number of
parameters are taken into account when determining the goodness of fit, the
$\Lambda$CDM model has the best value with only the OCG model managing to
match this accuracy. The MCG, GCG and MGCG models has better fits to the data
than the $\Lambda$CDM model, but due to their increasing amount of free
parameters, they did not obtain a goodness of fit, that matches the
$\Lambda$CDM and the OCG model’s results. Now that we have determined that
there are four possible alternative models, with the Generalized Cosmic CG
model being statistically the only model to obtain a likelihood value that is
less than the $\Lambda$CDM model, we can move towards the information criteria
to determine their viability of being an alternative model to explain the
accelerated expansion of the Universe. In terms of the AIC test, we see that
the OCG and the GCG models actually obtained some observational support, with
the OCG model just missing out on the substantial observational support
category, but is still higher the boundary for less observational support.
Therefore, it can be concluded that the OCG model has some observational
support according to the AIC criterion, while the GCG has less observational
support. The MCG and the MGCG models misses out on the less observational
support, but is not yet in the no observational support, so we cannot rule it
out.

In terms of the BIC criteria, we did not obtain one model to have some
observational support category, but the OCG model were the closest of being in
one of the categories. Therefore, statistically, based on the likelihood, the
goodness of fit, the AIC and BIC criteria, the OCG model is the most likely to
be an alternative model to the $\Lambda$CDM model, with the GCG, MGC and the
MGCG models not being ruled out, but will have to be tested on other datasets
before being accepted or rejected. Of these three boundary case models, the
MGCG model has the best chance of being a viable option due to it following
the trends produced by the $\Lambda$CDM model on average the closest, as well
as obtaining the most likelihood- and $\chi^{2}$-values between all the
different models.

Lastly, the GCCG model were statistically not viable to be an alternative
model, and when taking into account Figs. 17 and 18, it would have come at no
supersize, since it was visibily the weakest performing model, based on the
best-fitting MCMC results. Furthermore, due to the MCMC model obtaining a
best-fitting OCG model that can be considered to be an alternative model,
pending further testing, we know that the MCMC simulation were indeed
calibrated correctly, and that the error in the GCCG model is dependent on the
$\omega$ parameter for which the model is very sensitive to changes. It would
probably be wise to fix this parameter in the future based on other
experiments and then determine a best-fitting function for the GCCG model on
the remaining parameters.

##  IV Linear perturbations in Chaplygin gas cosmology

After we explore the background cosmology of CG models, here we also extended
our knowledge at a linear cosmological perturbation level to clearly see the
contributions of CG for the growth matter density fluctuations. Based on the
Jeffrey’s scale ranges of our statistical results 33, we consider that OCG and
GCG which are substantial observational support and less observational support
of CG models. The density contrast is studied in GR, $\Lambda$CDM and
baryoinc-Chaplygin gas approaches. Based on the statistical analysis as
presented in III.5, GMCG and GCCG model are not favored with observation and
not considered in this section to study the density contrasts. However MCG is
considered to study at least this models is failed at the background and see
how it behave at the perturbative level. To do this, we consider the non-
interacting non-relativistic fluid $\rho_{b}$ and relativistic fluid
$\rho_{\rm{r}}$, $\rho_{m}\equiv\rho_{r}+\rho_{b}$ with exotic fluid
$\rho_{ch}$ in our Universe. In the background FLRW cosmological, we define
spatial gradients of gauge-invariant variables such as

$\displaystyle
D^{m}_{a}\equiv\frac{a}{\rho_{m}}\tilde{\nabla}_{a}\rho_{m}\;,\quad
Z_{a}\equiv a\tilde{\nabla}_{a}\theta\;,$ $\displaystyle
D^{ch}_{a}\equiv\frac{a}{\rho_{ch}}\tilde{\nabla}_{a}\rho_{ch}\;,\quad
D^{t}_{a}\equiv\frac{a}{\rho_{t}}\tilde{\nabla}_{a}\rho_{t}\;,$

where the energy density $D^{m}_{a}$, volume expansion of the fluid $Z_{a}$ as
follows dunsby1992covariant ; abebe2012covariant ; carloni2006gauge ;
ananda2008detailed , $D^{ch}_{a}$ is the energy density for CG and $D^{t}_{a}$
is the spatial energy density for total fluids [Matter fluid +CG]. Those four
gradient variables are key to examine the evolution equation for matter
density fluctuations. For a perfect multi-fluid system, the following
conservation equations, considered in dunsby1991gauge ; dunsby1992covariant
hold:

$\displaystyle\dot{\rho}_{i}=-\theta(\rho_{i}+p_{i})=-\theta(1+w_{i})\rho_{i}\;,$
(35)
$\displaystyle\dot{\rho}_{t}=-\theta(\rho_{t}+p_{t})=-\theta(1+w_{t})\rho_{t}\;,$
(36) $\displaystyle(\rho_{t}+p_{t})\dot{u}_{a}+\tilde{\nabla}_{a}p_{t}=0\;.$
(37)

From these conservation equation, the 4-acceleration in the energy frame of
the total fluid can be given as

$\displaystyle\dot{u}_{a}$
$\displaystyle=-\frac{\tilde{\nabla}_{a}p_{t}}{\rho_{t}+p_{t}}$
$\displaystyle=-\frac{1}{a\left(1+w_{t}\right)\rho_{t}}\left[w_{m}\rho_{m}D^{m}_{a}+w_{ch}\rho_{ch}D^{ch}+a\rho_{ch}\tilde{\nabla}_{a}w_{ch}\right]$

with $w_{m}$ and $w_{ch}$ the equation of state parameters for matter fluid
and Chaplygin gas, and $w_{t}$ is the effective equation of state parameter
for the total fluid. Another key equation for a general fluid is the so-called
Raychaudhuri equation and it can be expressed as

$\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{1}{2}(\rho_{t}+3p_{t})+\tilde{\nabla}^{a}\dot{u}_{a}\;.$
(39)

The above-defined spatial gradient variables evolve according to:

$\displaystyle\dot{D}^{m}_{a}-\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{w_{m}\rho_{m}}{\rho_{t}}\theta
D^{m}_{a}-\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{w_{ch}\rho_{ch}}{\rho_{t}}\theta
D^{ch}_{a}-\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{\rho_{ch}}{\rho_{t}}a\theta\tilde{\nabla}_{a}w_{ch}+(1+w_{m})Z_{a}=0\;,$
(40)
$\displaystyle\dot{D}^{ch}_{a}-\left(\frac{1+w_{ch}}{1+w_{t}}\right)\frac{w_{ch}\rho_{ch}}{\rho_{t}}\theta
D^{ch}_{a}-\left(\frac{1+w_{ch}}{1+w_{t}}\right)\frac{w_{m}\rho_{m}}{\rho_{t}}\theta
D^{m}_{a}+\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{\rho_{m}}{\rho_{t}}a\theta\tilde{\nabla}_{a}w_{ch}+(1+w_{ch})Z_{a}=0\;,$
(41) $\displaystyle\dot{D}^{t}_{a}-w_{t}\theta D^{t}_{a}+(1+w_{t})Z_{a}=0\;,$
(42) $\displaystyle\dot{Z}_{a}+\frac{2}{3}\theta
Z_{a}+\left[\frac{1}{2}\left(1+3w_{m}\right)\rho_{m}-\frac{w_{m}\rho_{m}}{(1+w_{t})\rho_{t}}\left(\frac{1}{3}\theta^{2}+\frac{1}{2}(1+3w_{t})\rho_{t}\right)+\frac{w_{m}\rho_{m}}{(1+w_{t})\rho_{t}}\tilde{\nabla}^{2}\right]D^{m}_{a}$
$\displaystyle\quad~{}+\left[\frac{1}{2}\left(1+3w_{ch}\right)\rho_{ch}-\frac{w_{ch}\rho_{ch}}{(1+w_{t})\rho_{t}}\left(\frac{1}{3}\theta^{2}+\frac{1}{2}(1+3w_{t})\rho_{t}\right)+\frac{w_{ch}\rho_{ch}}{(1+w_{t})\rho_{t}}\tilde{\nabla}^{2}\right]D^{ch}_{a}$
$\displaystyle\quad~{}+\frac{\rho_{ch}}{(1+w_{t})\rho_{t}}\left[\rho_{t}-\frac{1}{3}\theta^{2}+\tilde{\nabla}^{2}\right]a\tilde{\nabla}_{a}w_{ch}=0\;.$
(43)

Note that the above equations (40) - (43) are the generalized evolution
equations of the perturbations for matter-Chaplygin gas models. To simplify
these evolution equations, we will consider OCG, GCG and MCG models by
specifying the gradient of their equation of state parameters as
$a\tilde{\nabla}_{a}w_{ch}$, and present evolution equation for each models
accordingly to analyze the growth of density contrasts with cosmic-time.

###  IV.1 The evolutions for GCG model

From the GCG equation of state form (1), we have

$a\tilde{\nabla}_{a}w_{ch}=-w_{ch}(\alpha+1)D^{ch}_{a}\;.$ (44)

Then, the evolution equations for GCG models are yield as:

$\displaystyle\dot{D}^{m}_{a}+(1+w_{m})Z_{a}-\frac{\theta
w_{m}(1+w_{m})\rho_{m}}{\rho_{t}(1+w_{t})}D^{m}_{a}+\frac{\theta\alpha
w_{ch}(1+w_{m})\rho_{ch}}{\rho_{t}(1+w_{t})}D^{ch}_{a}=0\;,$ (45)
$\displaystyle\dot{D}^{ch}_{a}+(1+w_{ch})Z_{a}-\frac{\theta
w_{m}(1+w_{ch})\rho_{m}}{\rho_{t}(1+w_{t})}D^{m}_{a}-\frac{\theta
w_{ch}}{(1+w_{t})\rho_{t}}\Big{(}(1+w_{m})(1+\alpha)\rho_{m}+(1+w_{ch})\rho_{ch}\Big{)}D^{ch}_{a}=0\;,$
(46) $\displaystyle\dot{D}^{t}_{a}+(1+w_{t})Z_{a}-\theta w_{t}D^{t}_{a}=0\;,$
(47) $\displaystyle\dot{Z}_{a}+\frac{2}{3}\theta
Z_{a}+\Big{\\{}\frac{1}{2}(1+3w_{m})\rho_{m}+\frac{w_{m}\rho_{m}}{\rho_{t}(1+w_{t})}\Big{(}-\frac{1}{3}\theta^{2}-\frac{1}{2}\rho_{t}(1+3w_{t})+\tilde{\nabla}^{2}\Big{)}\Big{\\}}D^{m}_{a}$
$\displaystyle~{}~{}~{}+\Big{\\{}\frac{1}{2}(1-3w_{ch}\alpha)\rho_{ch}-\frac{w_{ch}\rho_{ch}\alpha}{\rho_{t}(1+w_{t})}\Big{(}-\frac{1}{3}\theta^{2}-\frac{1}{2}\rho_{t}(1+3w_{t})+\tilde{\nabla}^{2}\Big{)}\Big{\\}}D^{ch}_{a}=0\;.$
(48)

By applying the scalar and harmonic decomposition of the vector gradient
variables defined in Eq. (IV) as presented in abebe2012covariant ;
sahlu2020scalar , the evolution of the $kth$ harmonic of the scalar
perturbations:

$\displaystyle\dot{\Delta}^{k}_{m}+(1+w_{m})Z_{k}-\frac{\theta
w_{m}(1+w_{m})\rho_{m}}{\rho_{t}(1+w_{t})}\Delta^{k}_{m}+\frac{\theta\alpha
w_{ch}(1+w_{m})\rho_{ch}}{\rho_{t}(1+w_{t})}\Delta^{k}_{ch}=0\;,$ (49)
$\displaystyle\dot{\Delta}^{k}_{ch}+(1+w_{ch})Z_{k}-\frac{\theta
w_{m}(1+w_{ch})\rho_{m}}{\rho_{t}(1+w_{t})}\Delta^{k}_{m}-\frac{\theta
w_{ch}}{(1+w_{t})\rho_{t}}\Big{[}(1+w_{m})(1+\alpha)\rho_{m}+(1+w_{ch})\rho_{ch}\Big{]}\Delta^{k}_{ch}=0\;,$
(50) $\displaystyle\dot{\Delta}^{k}_{t}+(1+w_{t})Z_{k}-\theta
w_{t}\Delta^{k}_{t}=0\;,$ (51) $\displaystyle\dot{Z}_{k}+\frac{2}{3}\theta
Z_{k}+\Big{\\{}\frac{1}{2}(1+3w_{m})\rho_{m}+\frac{w_{m}\rho_{m}}{\rho_{t}(1+w_{t})}\Big{[}-\frac{1}{3}\theta^{2}-\frac{1}{2}\rho_{t}(1+3w_{t})-\frac{k^{2}}{a^{2}}\Big{]}\Big{\\}}\Delta^{k}_{m}$
$\displaystyle~{}~{}~{}~{}+\Big{\\{}\frac{1}{2}(1-3w_{ch}\alpha)\rho_{ch}-\frac{w_{ch}\rho_{ch}\alpha}{\rho_{t}(1+w_{t})}\Big{[}-\frac{1}{3}\theta^{2}-\frac{1}{2}\rho_{t}(1+3w_{t})-\frac{k^{2}}{a^{2}}\Big{]}\Big{\\}}\Delta^{k}_{ch}=0\;.$
(52)

In redshift space as presented in sahlu2020scalar ; ntahompagaze2020multifluid
; sami2021covariant , the above equations can be recast as555From here on-
wards, we drop the index $k$ indicating the harmonic to avoid overcrowding of
notations.:

$\displaystyle\Delta^{\prime}_{m}+\frac{3w_{m}(1+w_{m})\Omega_{m}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{m}-\frac{3\alpha
w_{ch}(1+w_{m})\Omega_{ch}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{ch}-\frac{(1+w_{m})}{H(1+z)}Z=0\;,$
(53)
$\displaystyle\Delta^{\prime}_{ch}+\frac{3w_{m}(1+w_{ch})\Omega_{m}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{m}+\frac{3w_{ch}}{(1+w_{t})\Omega_{t}(1+z)}\Big{[}(1+w_{m})(1+\alpha)\Omega_{m}+(1+w_{ch})\Omega_{ch}\Big{]}\Delta_{ch}$
(54) $\displaystyle\quad\quad-\frac{(1+w_{ch})}{H(1+z)}Z=0\;,$
$\displaystyle\Delta^{\prime}_{t}+\frac{3w_{t}}{(1+z)}\Delta_{t}-\frac{(1+w_{t})}{H(1+z)}Z=0\;,$
(55) $\displaystyle
Z^{\prime}-\frac{2}{(1+z)}Z-\frac{3H}{(1+z)}\Big{\\{}\frac{1}{2}(1+3w_{m})\Omega_{m}+\frac{w_{m}\Omega_{m}}{\Omega_{t}(1+w_{t})}\Big{[}-1-\frac{1}{2}\Omega_{t}(1+3w_{t})-\frac{k^{2}}{3H^{2}a^{2}}\Big{]}\Big{\\}}\Delta_{m}$
$\displaystyle~{}~{}~{}-\frac{3H}{(1+z)}\Big{\\{}\frac{1}{2}(1-3w_{ch}\alpha)\Omega_{ch}-\frac{w_{ch}\Omega_{ch}\alpha}{\Omega_{t}(1+w_{t})}\Big{[}-1-\frac{1}{2}\Omega_{t}(1+3w_{t})-\frac{k^{2}}{3H^{2}a^{2}}\Big{]}\Big{\\}}\Delta_{ch}=0\;.$
(56)

To rewrite the equations in a dimensionally uniform way, we define

$h\equiv H/H_{0}\;,{\mathcal{Z}}\equiv Z/H_{0}\;,\kappa\equiv k/H_{0},$ (57)

so that our new evolution equations for GCG model are given as:

$\displaystyle\Delta^{\prime}_{m}+\frac{3w_{m}(1+w_{m})\Omega_{m}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{m}-\frac{3\alpha
w_{ch}(1+w_{m})\Omega_{ch}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{ch}-\frac{(1+w_{m})}{(1+z)h}\mathcal{Z}=0\;,$
(58)
$\displaystyle\Delta^{\prime}_{ch}+\frac{3w_{ch}}{(1+z)(1+w_{t})\Omega_{t}}\Big{(}(1+w_{m})(1+\alpha)\Omega_{m}+(1+w_{ch})\Omega_{ch}\Big{)}\Delta_{ch}+\frac{3w_{m}(1+w_{ch})\Omega_{m}}{\Omega_{t}(1+w_{t})(1+z)}\Delta_{m}$
$\displaystyle\quad\quad-\frac{(1+w_{ch})}{(1+z)h}\mathcal{Z}=0\;,$ (59)
$\displaystyle(1+z)h\Delta^{\prime}_{t}+3w_{t}h\Delta_{t}-(1+w_{t})\mathcal{Z}=0\;,$
(60)
$\displaystyle(1+z)\mathcal{Z}^{\prime}-2\mathcal{Z}-3h\Big{\\{}\frac{1}{2}(1+3w_{m})\Omega_{m}-\frac{w_{m}\Omega_{m}}{\Omega_{t}(1+w_{t})}\Big{[}1+\frac{1}{2}\Omega_{t}(1+3w_{t})+(1+z)^{2}\frac{\kappa^{2}}{h^{2}}\Big{]}\Big{\\}}\Delta_{m}$
$\displaystyle\quad\quad\quad\quad-3h\Big{\\{}\frac{1}{2}(1-3w_{ch}\alpha)\Omega_{ch}+\frac{w_{ch}\Omega_{ch}\alpha}{\Omega_{t}(1+w_{t})}\Big{[}1+\frac{1}{2}\Omega_{t}(1+3w_{t})+(1+z)^{2}\frac{\kappa^{2}}{h^{2}}\Big{]}\Big{\\}}\Delta_{ch}=0\;.$
(61)

For the case of $\alpha=1$, the evolution equations of OCG can be found from
these equations (58) - (61). The numerical results of density contrast in OCG
and GCG models will present in Sec. IV.5, to investigate the contribution of
Chaplygin gas like fluid for the formation of large-scale.

###  IV.2 The evolutions for MCG model

With the equation of state (16) for the MCG model, $a\tilde{\nabla}_{a}w_{ch}$
is yeild as

$\displaystyle a\tilde{\nabla}_{a}w_{ch}=(\alpha+1)(B-w_{ch})D^{ch}_{a}\;.$
(62)

By applying the same manners as presented in GCG model the evolution equation
for MCG read as

$\displaystyle{\Delta}^{\prime}_{m}+\left(\frac{3(1+w_{m})}{(1+w_{t})}\right)\frac{w_{m}\Omega_{m}}{(1+z)\Omega_{t}}\Delta_{m}+\left(\frac{3(1+w_{m})}{(1+w_{t})}\right)\frac{\Omega_{ch}}{(1+z)\Omega_{t}}\bigg{(}B+\alpha(B-w_{ch})\bigg{)}\Delta^{ch}-\frac{(1+w_{m})}{h(1+z)}\mathcal{Z}=0\;,$
(63)
$\displaystyle{\Delta}^{\prime}_{ch}+\frac{3}{(1+w_{t})(1+z)\Omega_{t}}\Big{(}(1+w_{m})(1+\alpha)(B-w_{ch})\Omega_{m}+w_{ch}(1+w_{ch})\Omega_{ch}\Big{)}\Delta^{ch}+\frac{3(1+w_{ch})w_{m}\Omega_{m}}{(1+z)(1+w_{t})\Omega_{t}}\Delta^{m}$
$\displaystyle\quad\quad-\frac{(1+w_{ch})}{h(1+z)}\mathcal{Z}=0\;,$ (64)
$\displaystyle{\Delta}^{\prime}_{t}+3\left(\frac{1}{1+z}\right)w_{t}\Delta_{t}-\left(\frac{1}{1+z}\right)\frac{(1+w_{t})}{h}\mathcal{Z}=0\;,$
(65)
$\displaystyle\mathcal{Z}^{\prime}-\frac{2}{(1+z)}\mathcal{Z}-\frac{3h}{(1+z)}\left[\frac{1}{2}\left(1+3w_{m}\right)\Omega_{m}+\frac{w_{m}\Omega_{m}}{(1+w_{t})\Omega_{t}}\left(-1-\frac{1}{2}(1+3w_{t})\Omega_{t}-\frac{k^{2}}{3a^{2}h^{2}}\right)\right]\Delta^{m}$
(66)
$\displaystyle~{}\quad-\frac{3h}{(1+z)}\Bigg{[}\frac{1}{2}\Omega_{ch}\Big{(}1+3B+3\alpha(B-w_{ch})\Big{)}-\frac{\Omega_{ch}\Big{(}B+\alpha(B-w_{ch})\Big{)}}{(1+w_{t})\Omega_{t}}\left(1+\frac{1}{2}(1+3w_{t})\Omega_{t}+\frac{k^{2}}{3a^{2}h^{2}}\right)\Bigg{]}\Delta^{ch}=0\;.$

N.B For the case of $B=0$ above equations 63 \- IV.2 reduced to Eqs 58 \- 61.
By admitting the definition of Eq. (57), the evolution equations for MGC model
yield as

$\displaystyle\Delta^{\prime}_{m}+\frac{3}{1+z}\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{w_{m}\Omega_{m}}{\Omega_{t}}\Delta_{m}+\frac{3}{1+z}\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{\Omega_{ch}}{\Omega_{t}}\left[B+(B-w_{ch})\alpha\right]\Delta_{ch}-\frac{1}{h(1+z)}(1+w_{m}){\mathcal{Z}}=0\;,$
(67)
$\displaystyle\Delta^{\prime}_{ch}-\frac{3}{1+z}\Bigg{\\{}\left[(1+\alpha)(B-w_{ch})\right]\left(\frac{1+w_{m}}{1+w_{t}}\right)\frac{\Omega_{m}}{\Omega_{t}}-w_{ch}\left(\frac{1+w_{ch}}{1+w_{t}}\right)\frac{\Omega_{ch}}{\Omega_{t}}\Bigg{\\}}\Delta_{ch}+\frac{3}{1+z}\left(\frac{1+w_{ch}}{1+w_{t}}\right)\frac{w_{m}\Omega_{m}}{\Omega_{t}}\Delta_{m}$
$\displaystyle\quad\quad-\frac{1}{h(1+z)}(1+w_{ch}){\mathcal{Z}}=0\;,$ (68)
$\displaystyle\Delta^{\prime}_{t}+\frac{3}{1+z}w_{t}\Delta_{t}-\frac{1}{h(1+z)}(1+w_{t}){\mathcal{Z}}=0\;,$
(69)
$\displaystyle{\mathcal{Z}}^{\prime}-\frac{2}{1+z}{\mathcal{Z}}-\frac{1}{h(1+z)}\left[\frac{3}{2}\left(1+3w_{m}\right)\Omega_{m}-\frac{3w_{m}\Omega_{m}}{(1+w_{t})\Omega_{t}}\left(h^{3}+\frac{1}{2}(1+3w_{t})\Omega_{t}\right)-\frac{w_{m}\Omega_{m}}{(1+w_{t})\Omega_{t}}(1+z)^{2}\kappa^{2}\right]\Delta_{m}$
$\displaystyle\quad+\frac{1}{h(1+z)}\left\\{\frac{3\Omega_{ch}}{2(1+w_{t})\Omega_{t}}\left[3\left((B-w_{ch})(w_{t}+\frac{1}{3})\alpha+(w_{t}+\frac{1}{3})B-(1+w_{t})((w_{ch}+\frac{1}{3})\right)\Omega_{t}+2\left[B+(B-w_{ch})\alpha\right]h^{3}\right]\right.$
$\displaystyle\quad\left.+\frac{\left[B+(B-w_{ch})\alpha\right]\Omega_{ch}}{(1+w_{t})\Omega_{t}}(1+z)^{2}\kappa^{2}\right\\}\Delta_{ch}=0\;.$
(70)

The density contrast in MCG models are presented in Sec. IV.5, to investigate
the contribution of Chaplygin gas like fluid for the formation of large-scale.

###  IV.3 The growth of density contrast in GR limits

From Eqs. (59), in the absence of Chaplygin gas like fluids the linear
evolution equation is reduced to the well known results in GR as

$\displaystyle\Delta^{\prime}_{m}+\frac{3w_{m}\Omega_{m}}{1+z}\Delta_{m}-\frac{(1+w_{m})}{(1+z)h}\mathcal{Z}=0\;,$
(71)
$\displaystyle\mathcal{Z}^{\prime}-\frac{2}{1+z}\mathcal{Z}-\frac{3h}{1+z}\Big{\\{}\frac{1}{2}(1+3w_{m})\Omega_{m}$
(72)
$\displaystyle-\frac{w_{m}}{(1+w_{m})}\Big{[}1+\frac{1}{2}\Omega_{m}(1+3w_{m})+(1+z)^{2}\frac{\kappa^{2}}{h^{2}}\Big{]}\Big{\\}}\Delta_{m}=0\;.$

The normalized energy density can be defined as presented in sahlu2020scalar

$\delta_{[m,ch]}(z)=\frac{\Delta_{[m,ch]}(z)}{\Delta(z_{in})}\;,$ (73)

where $\Delta_{in}$ is the initial value of $\Delta_{[m,ch]}(z)$ at
$z_{in}=1100$, since the variation of CMB temperature detected observational
in the order of $10^{-3}$ smoot1992structure  at $z\approx 1100$. The
numerical results of equation (IV.3) is presented in Fig. 19. From the plot we
clearly see the contribution of baryonic matter for the formation large-scale
structure which is expected theoretical and observation ally aspects.

Figure 19: The density contrasts of baryonic matter with red-shift $z$ in
standard GR.

###  IV.4 The growth of density contrast in $\Lambda$CDM limits

For $\Lambda$CDM limits, the linear evolution equations are given as

$\displaystyle\Delta^{\prime}_{m}+\frac{3w_{m}\Omega_{m}}{1+z}\Delta_{m}-\frac{(1+w_{m})}{(1+z)h}\mathcal{Z}=0\;,$
(74)
$\displaystyle\mathcal{Z}^{\prime}-\frac{2}{(1+z)}\mathcal{Z}-\frac{3h}{(1+z)}\Big{(}\frac{1}{2}\Omega_{m}(1+3\omega_{m})$
$\displaystyle-\frac{w_{m}\Omega_{m}}{\left(1+w_{t}\right)\Omega_{t}}\Big{(}1+\frac{1}{2}\Omega_{m}(1+3w_{m})-$
$\displaystyle\frac{1}{2}\Omega_{\Lambda}-\frac{(1+z)^{2}k^{2}}{3h^{2}}\Big{)}\Big{)}\Delta^{m}=0\;.$
(75)

The numerical results of equation (20) is presented in Fig. 20

Figure 20: The density contrasts of baryonic matter with red-shift $z$ in
$\Lambda$CDM.

###  IV.5 The growth of density contrast in baryonic matter-Chaplygin gas

In this subsection, we explain the contribution of these coupling system
(i.e., baryonic and Chaplygin gas) for the formation of large scale structure.
The long wavelength mode where $\kappa\ll 1$, and short-wavelength mode within
the horizon, where $\kappa\gg 1$ are considered in the whole system to present
the numerical results of density contrasts in OCG, GCG, MCG. From Table. 1 ,
we use the best fitting parameters values for each models. In the following
Figures. 21 \- 23 we present the behavior of the density contrast in short-
and long-wavelength modes for baryonic fluid, total and Chaplygin gas or
exotic fluids with cosmological red-shift.

Figure 21: Left: The density contrasts of baryonic-Chaplygin gas with $z$ for
short-wavelength mode at $\kappa=10$ for OCG model. Right: The density
contrasts of baryonic matter-Chaplygin gas with $z$ for long-wavelength mode
at $\kappa=0$ for OCG model.

Figure 22: Left: The density contrasts of baryonic-Chaplygin gas with $z$ for
short-wavelength mode at $\kappa=10$ for GCG model. Right: The density
contrasts of baryonic matter - Chaplygin gas with $z$ for long-wavelength mode
at $\kappa=0$ for GCG model. We use $\alpha=0.6318$ from Table 1 for
illustrative purpose.

Figure 23: The growth of density contrasts in the baryonic matter-Chaplygin
gas mixture for the short-wavelength mode at $\kappa=10$ for MCG model (left)
and long-wavelength mode at $\kappa=0$ (right). We use $B=0.1448$ and
$\alpha=0.5694$ from Table 1 for illustrative purpose.

##  V Conclusions

In this manuscript, the background cosmological parameters have been studied
and compared against supernova cosmological data for different Chaplygin gas
models, namely: OCG, GCG, MCG, GMCG and GCCG. The comparison of the
deceleration parameters of the various CG models with $\Lambda$CDM has been
done to identify which model can best explain the late-time accelerating
universe. From the results in Fig. 18 the GCCG model can be ruled out, and our
AIC and BIC criteria analyses also supported this as presented in Table 2. We
also considered the comparison of the residuals of the various CG models and
the $\Lambda$CDM as well, and from the results we observe that the GCCG model
is unrealistic compared to the other considered CG models in the paper. Our
statistical analysis test indicated that OCG, GCG and MCG are alternative
models are relatively “accepted/true model” by using the AIC and BIC selection
criteria.
The linear cosmological perturbations have also been investigated in a matter-
Chaplygin gas multi-fluid setting around the FLRW background using $1+3$
covariant formalism. Based on our statistical analysis for the background
cosmology, the three Chaplygin gas models that performed the best, namely OCG,
GCG and MCG were considered for the study of the formation of large-scale
structure. From the plots 21 and 22, we observe the contribution of baryonic,
Chaplygin gas and total fluid for the formation of large-scale. The results
show that whereas the density contrast of the baryonic matter is growing with
cosmic time (decreases with redshift), the density contrast of the Chaplygin
gas decays with cosmic time. It can be clearly seen that the Chaplygin gas
supports the formation of large-scale structures (as it behalves like dark
matter) in the early universe, but discourages it at late times as it behaves
like dark energy.

In general, both the OCG and GCG provide cosmological scenarios that are more
or less consistent with the $\Lambda$CDM model, both at background and
perturbation levels, but with the former having a better observational support
based on the statistical criteria described herein. From our statistical
analysis the MCG model has $\Delta IC>7$ which means, although not ruled out,
it does not have any significant observational support. We went a step further
and investigated the perturbations of this model, and found that even at this
level the growth rates are not in line with expected $\Lambda$CDM results. The
MGCG and GCCG models that were rejected by the background statistical criteria
were not included in our analysis of the perturbations.

##  VI Acknowledgments

SS acknowledges that this work is based on the research supported by Entoto
Observatory and Research Center, with grant number ESSTI-AA-IP-004-1021. SS is
grateful for the North-West University, Department of Physics, for hosting me
during the initial preparation of the work. HS gratefully acknowledges the
financial support from the Mwalimu Nyerere African Union scholarship and the
National Research Foundation (NRF) free-standing scholarship with grant number
112544. RH also acknowledges the NRF for the financial support for this work
with a grant number 120850. AA acknowledges that this work is based on the
research supported in part by the NRF of South Africa with grant number
112131.

## References

  * (1) Allemandi, G., Borowiec, A., Francaviglia, M. & Odintsov, S. D.  Dark energy dominance and cosmic acceleration in first-order formalism.  _Physical Review D_ 72, 063505 (2005).
  * (2) Riess, A. G. _et al._ Observational evidence from supernovae for an accelerating universe and a cosmological constant.  _The Astronomical Journal_ 116, 1009 (1998).
  * (3) Frieman, J. A., Turner, M. S. & Huterer, D.  Dark energy and the accelerating universe.  _Annu. Rev. Astron. Astrophys._ 46, 385–432 (2008).
  * (4) Clifton, T., Ferreira, P. G., Padilla, A. & Skordis, C.  Modified gravity and cosmology.  _Physics reports_ 513, 1–189 (2012).
  * (5) Faraoni, V. & Nadeau, S.  Stability of modified gravity models.  _Physical Review D_ 72, 124005 (2005).
  * (6) Saadat, H. & Pourhassan, B.  Viscous varying generalized chaplygin gas with cosmological constant and space curvature.  _International Journal of Theoretical Physics_ 52, 3712–3720 (2013).
  * (7) Dev, A., Alcaniz, J. & Jain, D.  Cosmological consequences of a chaplygin gas dark energy.  _Physical Review D_ 67, 023515 (2003).
  * (8) Elmardi, M. Y.  _Topics in Relativistic Cosmology_.  Ph.D. thesis (2018).
  * (9) Ellis, G. F., Hwang, J. & Bruni, M.  Covariant and gauge-independent perfect-fluid robertson-walker perturbations.  _Physical Review D_ 40, 1819 (1989).
  * (10) Dunsby, P. K.  Gauge invariant perturbations in multi-component fluid cosmologies.  _Class. Quant. Grav._ 8, 1785 (1991).
  * (11) Abebe, A., Abdelwahab, M., De la Cruz-Dombriz, Á. & Dunsby, P. K.  Covariant gauge-invariant perturbations in multifluid ${f(R)}$ gravity.  _Class. Quant. Grav._ 29, 135011 (2012).
  * (12) Tayler, R. J.  _The stars: their structure and evolution_ (Cambridge University Press, 1994).
  * (13) Sahlu, S., Ntahompagaze, J., Abebe, A., Cruz-Dombriz, Á. & Mota, D.  Scalar perturbations in ${f(T)}$ gravity using the $1+3$ covariant approach.  _European Physical Journal C_ 80, 1–19 (2020).
  * (14) Ntahompagaze, J., Sahlu, S., Abebe, A. & Mbonye, M. R.  On multifluid perturbations in scalar–tensor cosmology.  _International Journal of Modern Physics D_ 29, 2050120 (2020).
  * (15) Abebe, A.  Breaking the cosmological background degeneracy by two-fluid perturbations in $f(R)$ gravity.  _International Journal of Modern Physics D_ 24, 1550053 (2015).
  * (16) Kamenshchik, A., Moschella, U. & Pasquier, V.  An alternative to quintessence.  _Physics Letters B_ 511, 265–268 (2001).
  * (17) Elmardi, M. & Abebe, A.  Cosmological Chaplygin gas as modified gravity.  _Journal of Physics: Conference Series_ 905, 012015 (2017).
  * (18) Avelino, P., Beca, L., De Carvalho, J. & Martins, C.  The $\lambda$cdm limit of the generalized Chaplygin gas scenario.  _Journal of Cosmology and Astroparticle Physics_ 2003, 002 (2003).
  * (19) Kamenshcik, A., Moschella, U. & Pasquier, V.  An alternative to quintessence (2008).  ArXiv:gr-qc/0103004v2.
  * (20) Elmardi, M., Abebe, A. & Tekola, A.  Chaplygin-gas solutions of ${f(R)}$ gravity.  _International Journal of Geometric Methods in Modern Physics_ 13, 1650120 (2016).
  * (21) Bouhmadi-López, M., Brilenkov, M., Brilenkov, R., Morais, J. & Zhuk, A.  Scalar perturbations in the late universe: viability of the chaplygin gas models.  _Journal of Cosmology and Astroparticle Physics_ 2015, 037 (2015).
  * (22) Morais, J., Bouhmadi-López, M. & Capozziello, S.  Can ${f(R)}$ gravity contribute to (dark) radiation?  _Journal of Cosmology and Astroparticle Physics_ 2015, 041 (2015).
  * (23) González-Díaz, P. F.  You need not be afraid of phantom energy.  _Physical Review D_ 68, 021303 (2003).
  * (24) Deza, M. M. & Deza, E.  Encyclopedia of distances.  In _Encyclopedia of distances_ , 1–583 (Springer, 2009).
  * (25) Linder, E. V.  Understanding the optimal redshift range for the supernovae hubble diagram.  _arXiv preprint astro-ph/0108280_ (2001).
  * (26) Capozziello, S. & Sen, A. A.  Model-independent constraints on dark energy evolution from low-redshift observations.  _Monthly Notices of the Royal Astronomical Society_ 484, 4484–4494 (2019).
  * (27) Conley, A. _et al._ Supernova constraints and systematic uncertainties from the first three years of the supernova legacy survey.  _The Astrophysical Journal Supplement Series_ 192, 1 (2010).
  * (28) Neill, J. D. _et al._ The Local Hosts of Type Ia Supernovae.  _Astrophys. J._ 707, 1449–1465 (2009).  eprint 0911.0690.
  * (29) Hicken, M. _et al._ CfA3: 185 Type Ia Supernova Light Curves from the CfA.  _Astrophys. J._ 700, 331–357 (2009).  eprint 0901.4787.
  * (30) Colgáin, E. Ó.  A hint of matter underdensity at low z?  _Journal of Cosmology and Astroparticle Physics_ 2019, 006 (2019).
  * (31) Freedman, W. L. _et al._ Final results from the hubble space telescope key project to measure the hubble constant.  _The Astrophysical Journal_ 553, 47 (2001).
  * (32) Riess, A. G. _et al._ A redetermination of the hubble constant with the hubble space telescope from a differential distance ladder.  _The Astrophysical Journal_ 699, 539 (2009).
  * (33) Riess, A. G. _et al._ A 3% solution: determination of the hubble constant with the hubble space telescope and wide field camera 3.  _The Astrophysical Journal_ 730, 119 (2011).
  * (34) Hough, R.  _Constraining modified gravity models with cosmological data_.  Ph.D. thesis, North-West University (South Africa) (2020).
  * (35) Hough, R., Abebe, A. & Ferreira, S.  Constraining $f(R)$-gravity models with recent cosmological data.  In _SAIP Conference Proceedings_ (2019).
  * (36) Hough, R., Abebe, A. & Ferreira, S.  Viability tests of ${f(R)}$-gravity models with supernovae type 1a data.  _The European Physical Journal C_ 80, 1–15 (2020).
  * (37) Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J.  emcee: the mcmc hammer.  _Publications of the Astronomical Society of the Pacific_ 125, 306 (2013).
  * (38) Aghanim, N. _et al._ Planck 2018 results-vi. cosmological parameters.  _Astronomy & Astrophysics_ 641, A6 (2020).
  * (39) Riess, A. G. _et al._ A 2.4% determination of the local value of the hubble constant.  _The Astrophysical Journal_ 826, 56 (2016).
  * (40) Reid, M. J.  The distance to the center of the galaxy.  _Annual review of astronomy and astrophysics_ 31, 345–372 (1993).
  * (41) Leibundgut, B. & Sollerman, J.  A cosmological surprise: the universe accelerates.  _Europhysics news_ 32, 121–125 (2001).
  * (42) Nunes, R. C., Pan, S., Saridakis, E. N. & Abreu, E. M.  New observational constraints on ${f(R)}$ gravity from cosmic chronometers.  _Journal of Cosmology and Astroparticle Physics_ 2017, 005 (2017).
  * (43) Nesseris, S. & Garcia-Bellido, J.  Is the jeffreys’ scale a reliable tool for bayesian model comparison in cosmology?  _Journal of Cosmology and Astroparticle Physics_ 2013, 036 (2013).
  * (44) Dunsby, P. K., Bruni, M. & Ellis, G. F.  Covariant perturbations in a multifluid cosmological medium.  _The Astrophysical Journal_ 395, 54–74 (1992).
  * (45) Carloni, S., Dunsby, P. K. & Rubano, C.  Gauge invariant perturbations of scalar-tensor cosmologies: The vacuum case.  _Physical Review D_ 74, 123513 (2006).
  * (46) Ananda, K. N., Carloni, S. & Dunsby, P. K.  A detailed analysis of structure growth in ${f(R)}$ theories of gravity.  _Class. Quant. Grav._ 26, 235018 (2008).
  * (47) Sami, H., Sahlu, S., Abebe, A. & Dunsby, P. K.  Covariant density and velocity perturbations of the quasi-newtonian cosmological model in ${f(T)}$ gravity.  _arXiv preprint arXiv:2105.00646_ (2021).
  * (48) Smoot, G. F. _et al._ Structure in the cobe differential microwave radiometer first-year maps.  _The Astrophysical Journal_ 396, L1–L5 (1992).