# An 8.0% Determination of the Baryon Fraction in the Intergalactic Medium
from Localized Fast Radio Bursts

Bao Wang Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210023, China  School of Astronomy and Space Sciences, University of Science
and Technology of China, Hefei 230026, China  School of Physics and
Electronics, Hunan Normal University, Changsha 410081, China  Jun-Jie Wei
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023,
China  School of Astronomy and Space Sciences, University of Science and
Technology of China, Hefei 230026, China

###### Abstract

The dispersion measure (DM)–redshift relation of fast radio bursts (FRBs) has
been proposed as a potential new tool for probing intergalactic medium (IGM)
and for studying cosmology. However, the poor knowledge of the baryon fraction
in the IGM ($f_{\mathrm{IGM}}$) and its degeneracy with cosmological
parameters impose restrictions on the cosmological applications of FRBs.
Furthermore, DMs contributed by the IGM ($\mathrm{DM_{IGM}}$) and host galaxy
($\mathrm{DM_{host}}$), important cosmological quantities, cannot be exactly
extracted from observations, which would bring uncontrolled systematic
uncertainties in FRB cosmology. In this work, we use seventeen localized FRBs
to constrain $f_{\mathrm{IGM}}$ and its possible redshift evolution. Other
cosmological probes such as type Ia supernovae, baryon acoustic oscillations,
and cosmic microwave background radiation are combined to break parameter
degeneracy. Taking into account the probability distributions of
$\mathrm{DM_{IGM}}$ and $\mathrm{DM_{host}}$ derived from the the IllustrisTNG
simulation, we obtain a robust measurement of $f_{\mathrm{IGM}}=0.927\pm
0.075$, representing a precision of 8.0%. We find that there is no strong
evidence for the redshift dependence of $f_{\mathrm{IGM}}$ at the current
observational data level. The rapid progress in localizing FRBs will
significantly improve the constraints on $f_{\mathrm{IGM}}$.

Radio transient sources (2008) — Intergalactic medium (813) — Observational
cosmology (1146) — Cosmological parameters (339)

##  1 INTRODUCTION

Fast radio bursts (FRBs) are a class of brief ($\sim$ms) and intense
($\sim$Jy) radio transients with large dispersion measures (DMs), well in
excess of the expected contributions from the Milky Way (Lorimer et al., 2007;
Thornton et al., 2013; Petroff et al., 2016; Platts et al., 2019; Xiao et al.,
2021; Zhang, 2022). Owing to their anomalously high DMs, FRBs are believed to
be of extragalactic or even cosmological origin. To date, more than 600 FRBs
have been detected, and over two dozen of them have been reported to repeat
(CHIME/FRB Collaboration et al., 2021). There are more than 20 FRBs with
definite host galaxies and redshift measurements. These observations suggest
that FRBs are promising tools for studying cosmology. Some proposals include
using localized FRBs to measure the baryon number density of the universe
(Deng & Zhang, 2014; McQuinn, 2014; Macquart et al., 2020; Yang et al., 2022),
constrain the dark energy equation of state (Gao et al., 2014; Zhou et al.,
2014; Walters et al., 2018; Wei et al., 2018; Zhao et al., 2020; Qiu et al.,
2022), constrain cosmic reionization history (Deng & Zhang, 2014; Zheng et
al., 2014; Hashimoto et al., 2021), measure cosmological distance (Yu & Wang,
2017; Kumar & Linder, 2019), or measure the Hubble constant (Hagstotz et al.,
2022; Wu et al., 2022); using strongly lensed FRBs to probe the nature of
compact dark matter (Muñoz et al., 2016; Wang & Wang, 2018), or measure the
Hubble constant and cosmic curvature (Li et al., 2018).

Another cosmological puzzle is the baryon distribution of the universe. While
it is widely believed that more than three-quarters of the baryonic content of
the universe resides in the diffuse intergalactic medium (IGM), with only a
small fraction in galaxies and galaxy clusters (Fukugita et al., 1998; Cen &
Ostriker, 2006), gaining direct observational evidence of the baryon
distribution is challenging. There have been many studies of detecting the
baryon fraction in the IGM, $f_{\mathrm{IGM}}$, through numerical simulations
(Cen & Ostriker, 1999, 2006; Meiksin, 2009) or observations (Fukugita et al.,
1998; Fukugita & Peebles, 2004; Shull et al., 2012; Hill et al., 2016; Muñoz &
Loeb, 2018). For instance, Meiksin (2009) performed numerical simulations and
suggested that $\sim 90$% of the baryons produced by the Big Bang are
contained within the IGM at redshifts of $z\geq 1.5$ (i.e.,
$f_{\mathrm{IGM}}\approx 0.9$). It was observed that $18\pm 4$% of the baryons
exists in galaxies, circumgalactic medium, intercluster medium, and cold
neutral gas at $z\leq 0.4$, or equivalently $f_{\mathrm{IGM}}\approx 0.82$
(Shull et al., 2012). There is an ongoing debate about the value of
$f_{\mathrm{IGM}}$.

The observed DMs of cosmological FRBs are mainly contributed by the IGM. Since
the DM contributed by the IGM ($\mathrm{DM_{IGM}}$) carries important
information on the location of baryons in the late universe, one may combine
the $\mathrm{DM_{IGM}}$ and $z$ information of FRBs to constrain the IGM
baryon fraction $f_{\mathrm{IGM}}$. Indeed, a number of methods have been
proposed to estimate $f_{\mathrm{IGM}}$ by utilizing the DM($z$) data of FRBs
(Li et al., 2019, 2020; Wei et al., 2019; Walters et al., 2019; Qiang & Wei,
2020; Dai & Xia, 2021; Lin et al., 2022; Lemos et al., 2022). However, one
issue that restricts such studies is the strong degeneracy between
cosmological parameters and $f_{\mathrm{IGM}}$. It is hard to determine
$f_{\mathrm{IGM}}$ directly only relying on FRB data. Moreover, there is
another thorny issue that DMs contributed by FRB host galaxies
($\mathrm{DM_{host}}$) and the inhomogeneities in the IGM
($\mathrm{DM_{IGM}}$) cannot be exactly extracted from observations (Macquart
et al., 2020). Previous studies assumed fixed values for them, which would
bring uncontrollable systematic uncertainties in the analysis. A more
plausible approach is to treat them as probability distributions derived from
cosmological simulations (Jaroszynski, 2019; Macquart et al., 2020; Zhang et
al., 2020, 2021; Wu et al., 2022).

In this paper, we present a high-precision measurement of the baryon fraction
in the IGM with seventeen localized FRBs through the $\mathrm{DM_{IGM}}$-$z$
relation. In order to break the degeneracy between cosmological parameters and
$f_{\mathrm{IGM}}$, we combine FRB data with current constraints from type Ia
supernovae (SNe Ia), baryon acoustic oscillations (BAOs), and cosmic microwave
background (CMB) radiation. In our $f_{\mathrm{IGM}}$ estimation, the
reasonable probability distributions of $\mathrm{DM_{host}}$ and
$\mathrm{DM_{IGM}}$ derived from the IllustrisTNG simulation (Zhang et al.,
2020, 2021) are adopted to reduce the systematic errors. Additionally, to
explore the possible evolution of $f_{\mathrm{IGM}}$ with redshift, we also
consider two different parametric models, namely a constant model and a time-
dependent model.

The rest of our paper is organized as follows. In Section 2, we review the FRB
DM measurements and the probability distributions of $\mathrm{DM_{host}}$ and
$\mathrm{DM_{IGM}}$ derived from the IllustrisTNG simulation. In Section 3, we
give an introduction of the compilations of the three other cosmological
probes, including SNe Ia, BAO, and CMB. Monte Carlo Markov Chain (MCMC)
parameter inference results are presented in Section 4. Finally, conclusions
are summarized in Section 5.

##  2 FRB Dispersion Measures

The precise localization of FRBs to their host galaxies provides an ensemble
of DM and $z$ measurements. The DM measurement represents the integrated
column density of free electrons along the line of sight. For an extragalactic
FRB, its observed DM can be separated into the following components:

$\mathrm{DM_{obs}}(z)=\mathrm{DM_{ISM}^{MW}}+\mathrm{DM_{halo}^{MW}}+\mathrm{DM_{IGM}}(z)+\frac{\mathrm{DM_{host}}}{1+z}\;,$
(1)

where $\mathrm{DM_{ISM}^{MW}}$, $\mathrm{DM_{halo}^{MW}}$,
$\mathrm{DM_{IGM}}$, and $\mathrm{DM_{host}}$ represent the DM contributions
from the Milky Way interstellar medium (ISM), the Milky Way halo, the IGM, and
the FRB host galaxy, respectively. The $(1+z)$ factor converts the local
$\mathrm{DM_{host}}$ to the observed one (Deng & Zhang, 2014). Because of the
inhomogeneity of the free electron distribution in the IGM, two sources at the
same redshift but in different sightlines will likely have significant
differences in the measured value of $\mathrm{DM_{IGM}}$. Adopting the flat
$\Lambda$CDM cosmological model, the average value of $\mathrm{DM_{IGM}}$ can
be estimated as (Deng & Zhang, 2014)

$\left\langle\mathrm{DM}_{\mathrm{IGM}}\right\rangle(z)=\frac{3c\Omega_{b}H_{0}}{8\pi
Gm_{p}}\int_{0}^{z}\frac{(1+z^{\prime})f_{\mathrm{IGM}}(z^{\prime})\chi_{e}(z^{\prime})}{\sqrt{\Omega_{m}\left(1+z^{\prime}\right)^{3}+1-\Omega_{m}}}dz^{\prime}\;,$
(2)

where $m_{p}$ is the proton mass, $H_{0}$ is the Hubble constant,
$f_{\mathrm{IGM}}(z)$ is the baryon fraction in the IGM, and $\Omega_{b}$ and
$\Omega_{m}$ are the present-day baryon and matter density parameters. The
free electron number fraction per baryon is
$\chi_{e}(z)=\frac{3}{4}\chi_{e,\mathrm{H}}(z)+\frac{1}{8}\chi_{e,\mathrm{He}}(z)$,
where $\chi_{e,\mathrm{H}}(z)$ and $\chi_{e,\mathrm{He}}(z)$ are the
ionization fractions of hydrogen and helium, respectively. Both hydrogen and
helium are completely ionized at redshifts $z<3$ (Meiksin, 2009; Becker et
al., 2011), allowing one to set $\chi_{e,\mathrm{H}}=\chi_{e,\mathrm{He}}=1$,
which gives $\chi_{e}=7/8$.

The $\mathrm{DM_{IGM}}$ value of a well-localized FRB can be extracted using
$\mathrm{DM_{IGM}}=\mathrm{DM_{obs}}-\mathrm{DM_{ISM}^{MW}}-\mathrm{DM_{halo}^{MW}}-\mathrm{DM_{host}}/(1+z)$.
Here the $\mathrm{DM_{ISM}^{MW}}$ term can be well estimated from the NE2001
model of the ISM free electron distribution (Cordes & Lazio, 2002). The
$\mathrm{DM_{halo}^{MW}}$ term is not well constrained, but is expected to
contribute $50\sim 80$ $\mathrm{pc\;cm^{-3}}$ (Prochaska & Zheng, 2019).
Hereafter we assume that the probability distribution of
$\mathrm{DM_{halo}^{MW}}$ can be described by a Gaussian distribution with
mean $\mu_{\mathrm{halo}}=65$ $\mathrm{pc\;cm^{-3}}$ and standard deviation
$\sigma_{\mathrm{halo}}=15$ $\mathrm{pc\;cm^{-3}}$ (Wu et al., 2022):

$P_{\rm
halo}(\mathrm{DM_{halo}^{MW}})=\frac{1}{\sqrt{2\pi}\sigma_{\mathrm{halo}}}\exp\left[-\frac{\left(\mathrm{DM_{halo}^{MW}}-\mu_{\mathrm{halo}}\right)^{2}}{2\sigma_{\mathrm{halo}}^{2}}\right]\;.$
(3)

Based on the state-of-the-art IllustrisTNG simulation (Springel et al., 2018),
Zhang et al. (2020) selected a large sample of simulated galaxies with similar
properties to observed FRB hosts to derive the distributions of
$\mathrm{DM}_{\mathrm{host}}$ of repeating and non-repeating FRBs. The
distributions of $\mathrm{DM}_{\mathrm{host}}$ can be well described by the
log-normal function (Macquart et al., 2020; Zhang et al., 2020)

$P_{\rm
host}\left(\mathrm{DM}_{\mathrm{host}}\right)=\frac{1}{\sqrt{2\pi}\mathrm{DM}_{\mathrm{host}}\sigma_{\mathrm{host}}}\exp\left[-\frac{\left(\ln\mathrm{DM}_{\mathrm{host}}-\mu_{\mathrm{host}}\right)^{2}}{2\sigma_{\mathrm{host}}^{2}}\right]\;,$
(4)

where $e^{\mu_{\mathrm{host}}}$ and
$e^{2\mu_{\mathrm{host}}+\sigma_{\mathrm{host}}^{2}}(e^{\sigma_{\mathrm{host}}^{2}}-1)$
are the mean and variance of the distribution, respectively. Due to the
diversity of host galaxies, Zhang et al. (2020) computed the
$\mathrm{DM_{host}}$ distributions for repeating FRBs in dwarf galaxies (like
FRB 121102, FRB 180301, FRB 181030, and FRB 190711), repeating FRBs in spiral
galaxies (like FRB 180916 and FRB 201124), and non-repeating FRBs separately.
Here we divide the localized FRBs into these three types according to their
host properties. The evolution of the median of $\mathrm{DM_{host}}$
($\mathrm{\mu_{host}}$) can be fitted by
$e^{\mathrm{\mu_{host}}}(z)=\kappa(1+z)^{\gamma}$, where $\kappa$ and $\gamma$
are given by Zhang et al. (2020). The propagated uncertainty
$\sigma_{\mathrm{host}}$ of $\mathrm{DM_{host}}$ is calculated from the
uncertainties of $\kappa$ and $\gamma$. With this expression of redshift
evolution, we can derive the $\mathrm{DM}_{\mathrm{host}}$ distributions at
any redshift of a localized FRB.

To date, more than 20 FRBs have already been localized. Nonetheless, some of
them are not available for our analysis. For example, the DM of FRB 181030 is
only 103.396 $\mathrm{pc\;cm^{-3}}$ (Bhardwaj et al., 2021a), which will be
reduced to a negative value after subtracting DM contributions from the Milky
Way ISM and halo. That is, the integral upper limit
($\mathrm{DM_{E}}\equiv\mathrm{DM_{obs}}-\mathrm{DM_{ISM}^{MW}}-\mathrm{DM_{halo}^{MW}}$)
in the probability of the external DM contribution outside our Galaxy (see
Equation (7)) will become negative. FRB 190520B is co-located with a compact,
persistent radio source and associated with a dwarf host galaxy at a redshift
of 0.241 (Niu et al., 2022). It is a clear outlier from the general trend of
the extragalactic $\mathrm{DM_{E}}$–$z$ relation, with an unprecedented DM
contribution from its host galaxy. Thus its $\mathrm{DM_{host}}$ term can not
be accurately deducted. FRB 200110E is located in a globular cluster in the
direction of the nearby galaxy M81 (Bhardwaj et al., 2021b; Kirsten et al.,
2022). The distance of FRB 200110E is only 3.6 Mpc, and the IGM between the
Milky Way and M81 contributes of the order of $\mathrm{DM_{IGM}}\approx 1$
$\mathrm{pc\;cm^{-3}}$. Thus the cosmological information carried by FRB
200110E is too little. Additionally, the peculiar velocity effect is
significant, which makes it can not be used for cosmological studies. After
excluding these FRBs, we use a sample of 17 FRBs in the redshift range
$0.0337\leq z\leq 0.66$ to constrain $f_{\mathrm{IGM}}(z)$.111After this work
was done, we noticed a new article (Ryder et al., 2022), which has reported
the discovery of a burst, FRB 220610A, in a complex host galaxy system at a
redshift of $z=1.016$. Whereas this burst was not included in our sample.
Table 1 lists the redshifts, $\mathrm{DM_{obs}}$, and $\mathrm{DM_{ISM}^{MW}}$
of our sample. The estimated $\mathrm{DM_{IGM}}$ and measured $z$ values for
the 17 localized FRBs are shown in Figure 1. We have estimated
$\mathrm{DM_{IGM}}$ by subtracting the following from the observed
$\mathrm{DM_{obs}}$ value: $\mathrm{DM_{ISM}^{MW}}$ from the Galactic ISM
model; a median of 65 $\mathrm{pc\;cm^{-3}}$ contributed by
$\mathrm{DM_{halo}^{MW}}$; and a median value of $\mathrm{DM_{host}}$ at
different redshifts estimated from the IllustrisTNG simulation (Zhang et al.,
2020).

Figure 1: The $\mathrm{DM_{IGM}}$-$z$ relation for 17 localized FRBs. The data points are estimations of $\mathrm{DM_{IGM}}$ versus redshift measurements for the 17 localized FRBs. The $\mathrm{DM_{IGM}}$ values are estimated by correcting the observed $\mathrm{DM_{obs}}$ for the contributions from our galaxy and the FRB host galaxy (see the text for details). The solid line shows the model of Equation (2) with the inferred parameters of $\Omega_{m}=0.309\pm 0.006$, $\Omega_{b}h^{2}=0.02245\pm 0.00013$, $H_{0}=67.73\pm 0.44$ $\mathrm{km\;s^{-1}\;Mpc^{-1}}$, and $f_{\mathrm{IGM,0}}=0.927\pm 0.075$. The shaded area represents the uncertainty of the estimated average $\mathrm{DM_{IGM}}$ value at the 95% confidence level due to the uncertainties of the inferred model parameters which we mentioned above. Table 1: Properties of the 17 localized FRBs Name | Redshift | $\mathrm{DM_{obs}}$ | $\mathrm{DM_{ISM}^{MW}}$ | Refs.
---|---|---|---|---
|  | $(\mathrm{pc\;cm^{-3}})$ | $(\mathrm{pc\;cm^{-3}})$ |
FRB 121102 | 0.19273 | $557\pm 2$ | 188.0 | 1
FRB 180301 | 0.3304 | $536\pm 0.2$ | 152.0 | 2
FRB 180916 | 0.0337 | $349.349\pm 0.005$ | 200.0 | 3
FRB 180924 | 0.3214 | $361.42\pm 0.06$ | 40.5 | 4
FRB 181112 | 0.4755 | $589.27\pm 0.03$ | 102.0 | 5
FRB 190102 | 0.291 | $363.6\pm 0.3$ | 57.3 | 6
FRB 190523 | 0.66 | $760.8\pm 0.6$ | 37.0 | 7
FRB 190608 | 0.1178 | $338.7\pm 0.5$ | 37.2 | 8
FRB 190611 | 0.378 | $321.4\pm 0.2$ | 57.83 | 9
FRB 190614 | 0.6 | $959.2\pm 0.5$ | 83.5 | 10
FRB 190711 | 0.522 | $593.1\pm 0.4$ | 56.4 | 9
FRB 190714 | 0.2365 | $504\pm 2$ | 38.0 | 9
FRB 191001 | 0.234 | $506.92\pm 0.04$ | 44.7 | 9
FRB 191228 | 0.2432 | $297.5\pm 0.05$ | 33.0 | 2
FRB 200430 | 0.16 | $380.1\pm 0.4$ | 27.0 | 9
FRB 200906 | 0.3688 | $577.8\pm 0.02$ | 36.0 | 2
FRB 201124 | 0.098 | $413.52\pm 0.05$ | 123.2 | 11

To construct a likelihood function $\mathcal{L}$ from FRB measurements, we
build a model for $\mathrm{DM_{IGM}}$. The model probability distribution for
$\mathrm{DM_{IGM}}$ has been derived from the theoretical treatments of the
IGM and galaxy halos with a standard deviation $\sigma_{\mathrm{DM}}$
dominated by the variance in $\mathrm{DM_{IGM}}$. The $\mathrm{DM_{IGM}}$
distributions derived in both semi-analytic models and cosmological
simulations can be well fitted by a quasi-Gaussian function with a long tail
(McQuinn, 2014; Prochaska & Zheng, 2019; Macquart et al., 2020)

$P_{\mathrm{IGM}}(\Delta)=A\Delta^{-\beta}\exp\left[-\frac{\left(\Delta^{-\epsilon}-C_{0}\right)^{2}}{2\epsilon^{2}\sigma_{\mathrm{DM}}^{2}}\right],\;\;\;\Delta>0\;,$
(5)

where
$\Delta\equiv\mathrm{DM}_{\mathrm{IGM}}/\left\langle\mathrm{DM}_{\mathrm{IGM}}\right\rangle$,
$A$ is a normalization coefficient, and the indices $\epsilon$ and $\beta$ are
related to the inner density profile of gas in halos. Here we take
$\epsilon=3$ and $\beta=3$, as Macquart et al. (2020) did in their treatment.
$C_{0}$ is a free parameter, which can be fitted when the mean
$\langle\Delta\rangle=1$. The motivation for this analytic form (Equation (5))
is that in the limit of small $\sigma_{\mathrm{DM}}$, the $\mathrm{DM_{IGM}}$
distribution should approach a Gaussian owing to the more diffuse halo gas and
the Gaussianity of structure on large scales. Conversely, when the variance is
large, this probability distribution captures the large skew that due to a few
large structures that contribute to the DM of many sightlines. Recently, Zhang
et al. (2021) used the IllustrisTNG simulation to estimate the probability
distributions of $\mathrm{DM_{IGM}}$ at different redshifts realistically.
Following Wu et al. (2022), the best-fit parameters ($A$, $C_{0}$, and
$\sigma_{\mathrm{DM}}$) of the $\mathrm{DM_{IGM}}$ distributions at the
different redshifts that presented by Zhang et al. (2021) are used for our
purpose. The uncertainties of these best-fit parameters may impact our final
$f_{\rm IGM}$ constraint. To investigate whether the uncertainties of these
parameters affect the $\Omega_{b}$ constraint (similar to our $f_{\rm IGM}$
constraint), Yang et al. (2022) derived $\Omega_{b}$ using the best-fit values
of these parameters plus or minus the uncertainties. They found that these
uncertainties have almost no effect on the final $\Omega_{b}$ constraint.
Given the fact that our $f_{\rm IGM}$ constraint is almost the same as that of
Yang et al. (2022), we can come to the same conclusion. Note that since the
$\mathrm{DM_{IGM}}$ distributions are given in discrete redshifts (Zhang et
al., 2021), we extrapolate them to the redshifts of the localized FRBs through
cubic spline interpolation.

Given the model for $\mathrm{DM_{IGM}}$, we estimate the likelihood function
by computing the joint likelihoods of 17 FRBs (Macquart et al., 2020):

$\mathcal{L}_{\rm
FRB}=\prod_{i=1}^{N_{\mathrm{FRB}}}P_{i}\left(\mathrm{DM}_{\mathrm{E},i}\right)\;,$
(6)

where $P_{i}(\mathrm{DM}_{\mathrm{E},i})$ is the probability of the total
observed $\mathrm{DM_{obs}}$ corrected for our galaxy, i.e.,
$\mathrm{DM_{E}}\equiv\mathrm{DM_{obs}}-\mathrm{DM_{ISM}^{MW}}-\mathrm{DM_{halo}^{MW}}=\mathrm{DM_{IGM}}+\mathrm{DM_{host}}/(1+z)$.
For a burst at redshift $z_{i}$, we have

$\displaystyle P_{i}\left(\mathrm{DM}_{\mathrm{E},i}\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\mathrm{DM}_{\mathrm{E},i}}P_{\text{host}}\left(\mathrm{DM}_{\mathrm{host}}\right)$
(7) $\displaystyle\times$ $\displaystyle
P_{\mathrm{IGM}}\left(\mathrm{DM}_{\mathrm{E},i}-\mathrm{DM}_{\mathrm{host},z_{i}}\right)\mathrm{dDM}_{\mathrm{host}}\;,$

where the probability density functions for
$P_{\text{host}}(\mathrm{DM_{host}})$ and $P_{\text{IGM}}(\mathrm{DM_{IGM}})$
are obtained from Equations (4) and (5), respectively. Note that the Milky Way
halo DM distribution ($\mathrm{DM_{halo}^{MW}}$) will be considered as a free
parameter in our analysis. We will marginalize $\mathrm{DM_{halo}^{MW}}$ using
a Gaussian prior of $\mathrm{DM_{halo}^{MW}}=65\pm 15$ $\mathrm{pc\;cm^{-3}}$
over the range of
$[\mu_{\mathrm{halo}}-3\sigma_{\mathrm{halo}},\;\mu_{\mathrm{halo}}+3\sigma_{\mathrm{halo}}]$,
where $\mu_{\mathrm{halo}}=65$ $\mathrm{pc\;cm^{-3}}$ and
$\sigma_{\mathrm{halo}}=15$ $\mathrm{pc\;cm^{-3}}$ (see Equation (3)).

##  3 Other Cosmological Probes

As we discussed in Section 1, in order to break the degeneracy between the IGM
baryon fraction $f_{\mathrm{IGM}}$ and other cosmological parameters
($\Omega_{m}$, $\Omega_{b}$, and $H_{0}$; see Equation (2)), we use up-to-date
cosmological data compilations, including SNe Ia, BAO, and CMB data. The exact
compilations and the details for the likelihoods are described separately in
what follows.

###  3.1 Type Ia Supernovae

The SNe Ia dataset that we use in this work is the Pantheon sample, which
consists of 1048 SNe Ia in the redshift range $0.01<z<2.3$ (Scolnic et al.,
2018). The observed distance modulus of each SN is given as

$\mu_{\rm SN}=m_{\rm corr}-M_{B}\;,$ (8)

where $m_{\rm corr}$ is the corrected apparent magnitude and $M_{B}$ is the
absolute magnitude.

The theoretical distance modulus $\mu_{\rm th}(z)$ is defined as

$\mu_{\rm th}=5\log_{10}\left[\frac{d_{L}(z)}{\mathrm{Mpc}}\right]+25\;,$ (9)

where
$d_{L}(z)=(1+z)\frac{c}{H_{0}}\int_{0}^{z}{\frac{dz^{\prime}}{\sqrt{\Omega_{m}(1+z^{\prime})^{3}+1-\Omega_{m}}}}$
is the luminosity distance in the flat $\Lambda$CDM model. Thus, the
$\chi^{2}_{\rm SN}$ function for the Pantheon data is

$\chi^{2}_{\rm SN}=\Delta\bf{\hat{\mu}}^{\emph{T}}\cdot C_{\rm
SN}^{-1}\cdot\Delta\bf{\hat{\mu}}\;,$ (10)

where $\Delta\hat{\mu}=\hat{\mu}_{\rm SN}-\hat{\mu}_{\rm th}$ is the data
vector, defined by the difference between the SN distance modulus $\mu_{\rm
SN}$ and the theoretical distance modulus $\mu_{\rm th}$, and $\bf{C_{\rm
SN}}$ is the covariance matrix that contains both statistical and systematic
uncertainties of SNe.

###  3.2 Baryon Acoustic Oscillations

Primordial perturbations in the early universe excite acoustic waves in the
plasma, known as BAO. After the recombination period, the propagation of
acoustic waves was frozen. Thus, there is a characteristic scale called the
comoving sound horizon $r_{s}$, which can be approximated as (Aubourg et al.,
2015)

$\displaystyle r_{s}(z^{\ast})$ $\displaystyle=$
$\displaystyle\int_{z^{\ast}}^{\infty}\frac{c_{s}(z)}{H(z)}dz$ (11)
$\displaystyle\approx$
$\displaystyle\frac{55.154\exp\left[-72.3\left(\Omega_{v}h^{2}+0.0006\right)^{2}\right]}{\left(\Omega_{b}h^{2}\right)^{0.12807}\left(\Omega_{m}h^{2}\right)^{0.25351}}\;\mathrm{Mpc}\;,$

where $c_{s}(z)$ is the sound speed of the photon-baryon fluid, $z$ is the
redshift at the drag epoch, $\Omega_{v}$ is the present-day neutrino density,
and $h\equiv H_{0}/(100\;\mathrm{km\;s^{-1}\;Mpc^{-1}})$ is the reduced Hubble
constant. Here we use a combination of 11 BAO measurements from Ryan et al.
(2019). Six of these BAO measurements are correlated, in which case $\chi^{2}$
is given by

$\chi^{2}_{\rm
BAO1}=\left(\bf{A_{obs}}-\bf{A_{th}}\right)^{\emph{T}}\cdot\bf{C_{\rm
BAO}^{-1}}\cdot\left(\bf{A_{obs}}-\bf{A_{th}}\right)\;,$ (12)

where $\bf{A_{obs}}$ ($\bf{A_{th}}$) is the vector that contains all of the
six measured (theoretical) values and $\bf{C_{\rm BAO}}$ is the covariance
matrix for the BAO data sets. The other five BAO measurements are
uncorrelated, so

$\chi^{2}_{\rm BAO2}=\sum_{i=1}^{5}\frac{\left[A_{\rm th}(z_{i})-A_{\rm
obs}(z_{i})\right]^{2}}{\sigma_{i}^{2}}\;,$ (13)

where $\sigma_{i}$ is the standard deviation of the $i$-th BAO measurement
$A_{\rm obs}(z_{i})$. The data are combined into a $\chi^{2}$-statistic
$\chi^{2}_{\rm BAO}=\chi^{2}_{\rm BAO1}+\chi^{2}_{\rm BAO2}$.

###  3.3 Cosmic Microwave Background

For the CMB measurements, we use the derived parameters, including the
acoustic scale $l_{A}$, the shift parameter $R$, and $\Omega_{b}h^{2}$ from
the Planck analysis of the CMB (TT, TE, EE $+$ lowE) (Chen et al., 2019;
Planck Collaboration et al., 2020). The acoustic scale is

$l_{A}(z^{\ast})=(1+z^{\ast})\frac{\pi d_{A}(z^{\ast})}{r_{s}(z^{\ast})}\;,$
(14)

where $r_{s}$ is the comoving sound horizon at the recombination and
$d_{A}=d_{L}(1+z)^{-2}$ is the angular diameter distance. The shift parameter
is

$R(z^{\ast})=(1+z^{\ast})\frac{d_{A}(z^{\ast})\sqrt{\Omega_{m}}H_{0}}{c}\;.$
(15)

The redshift at decoupling $z^{\ast}$ is given by

$z^{\ast}=1048\left[1+0.00124\left(\Omega_{b}h^{2}\right)^{-0.738}\right]\left[1+g_{1}\left(\Omega_{m}h^{2}\right)^{g_{2}}\right]\;,$
(16)
$g_{1}=\frac{0.0738\left(\Omega_{b}h^{2}\right)^{-0.238}}{1+39.5\left(\Omega_{b}h^{2}\right)^{0.763}}\;,$
(17) $g_{2}=\frac{0.560}{1+21.1\left(\Omega_{b}h^{2}\right)^{1.81}}\;.$ (18)

By setting $\boldsymbol{x}=(R,\;l_{A},\;\Omega_{b}h^{2})$, the $\chi^{2}_{\rm
CMB}$ value for the CMB data is then

$\chi^{2}_{\rm CMB}=\left(\boldsymbol{x}_{\rm obs}-\boldsymbol{x}_{\rm
th}\right)^{\emph{T}}\cdot\bf{C_{\rm CMB}^{-1}}\cdot\left(\boldsymbol{x}_{\rm
obs}-\boldsymbol{x}_{\rm th}\right),$ (19)

where $\boldsymbol{x}_{\rm obs}=(1.7502,\;301.471,\;0.02236)$ and
$\boldsymbol{x}_{\rm th}$ contain the observed and theoretical values of the
derived parameters, respectively, and $\bf{C_{\rm CMB}}$ is the covariance
matrix for the CMB data. Note that the distance priors derived from the CMB
data are dependent on the specific cosmological model. Here we adopt the
values of $\boldsymbol{x}_{\rm obs}$ and $\bf{C_{\rm CMB}}$ that inferred from
the flat $\Lambda$CDM model.

##  4 Parameter Estimate and Results

To assess how well localized FRBs may help to constrain the evolution of the
IGM baryon fraction $f_{\mathrm{IGM}}(z)$, we consider two different
parametric models. First, a simple constant model,

$f_{\rm{IGM}}=f_{\rm{IGM},0}\;.$ (20)

And second, a time-dependent model given by

$f_{\mathrm{IGM}}(z)=f_{\mathrm{IGM},0}\left(1+\alpha\frac{z}{1+z}\right)\;,$
(21)

where $f_{\rm{IGM},0}$ is the present value of $f_{\rm{IGM}}$ and $\alpha$
quantifies any possible evolution of $f_{\rm{IGM}}$. As massive halos are more
abundant in the late universe, $f_{\rm{IGM}}$ is believed to grow with
redshift (McQuinn, 2014; Prochaska & Zheng, 2019). Therefore, here we require
$\alpha\geq 0$.

The quantities $f_{\mathrm{IGM}}(z)$, $\Omega_{m}$, $\Omega_{b}h^{2}$,
$H_{0}$, $M_{B}$, and $\mathrm{DM_{halo}^{MW}}$ are fitted to the FRB, SN Ia,
BAO, and CMB data simultaneously using the Python MCMC module $emcee$
(Foreman-Mackey et al., 2013). Given the relation
$\mathcal{L}\propto\exp[-\chi^{2}/2]$, the final log-likelihood sampled by
$emcee$ is a sum of the separated likelihoods of FRBs, SNe Ia, BAO, and CMB:

$\ln\mathcal{L}_{\rm tot}=\ln\mathcal{L}_{\rm FRB}-\chi^{2}_{\rm
SN}/2-\chi^{2}_{\rm BAO}/2-\chi^{2}_{\rm CMB}/2\;.$ (22)

In our baseline analysis, we set flat priors on $f_{\mathrm{IGM,0}}\in[0,\;1]$
and $\alpha\in[0,\;2]$, and a Gaussian prior on $\mathrm{DM_{halo}^{MW}}=65\pm
15$ $\mathrm{pc\;cm^{-3}}$ over the 3$\sigma$ range of $[20,\;110]$
$\mathrm{pc\;cm^{-3}}$.

Figure 2: 1D and 2D marginalized posterior distributions with the $1-2\sigma$
contours for the parameters $f_{\mathrm{IGM,0}}$, $\Omega_{m}$,
$\Omega_{b}h^{2}$, $H_{0}$, $M_{B}$, and $\mathrm{DM_{halo}^{MW}}$ in the
constant case of $f_{\rm{IGM}}=f_{\rm{IGM},0}$, constrained by the FRB, SN Ia,
BAO, and CMB data. The 1D marginalized posterior distribution of
$f_{\mathrm{IGM,0}}$ is magnified on the top-right panel. The vertical dashed
line represents the best fit, whereas the shaded areas correspond to the 68%
and 95% confidence regions. Table 2: Constraints on all parameters for two
different parametric models of $f_{\mathrm{IGM}}$.

Model | Constant | Time-dependent
---|---|---
Parameter | Estimation with 68% limits
$f_{\mathrm{IGM,0}}$ | $0.927\pm 0.075$ | $0.837\pm 0.089$
$\alpha$ | – | $<0.882$
$\Omega_{m}$ | $0.309\pm 0.006$ | $0.309\pm 0.006$
$\Omega_{b}h^{2}$ | $0.02245\pm 0.00013$ | $0.02245\pm 0.00013$
$H_{0}/[\mathrm{km\;s^{-1}\;Mpc^{-1}}]$ | $67.78\pm 0.44$ | $67.71\pm 0.42$
$M_{B}$ | $-19.415\pm 0.012$ | $-19.417\pm 0.012$
$\mathrm{DM_{halo}^{MW}}/[\mathrm{pc\;cm^{-3}}]$ | $47\pm 10$ | $49\pm 10$
$-2\ln\mathcal{L}_{\rm max}$ | $1262.855$ | $1263.099$
$\Delta$AIC | – | $2.244$

Note. — References: (1) Chatterjee et al. (2017); (2) Bhandari et al. (2022);
(3) Marcote et al. (2020); (4) Bannister et al. (2019); (5) Prochaska et al.
(2019); (6) Bhandari et al. (2020); (7) Ravi et al. (2019); (8) Chittidi et
al. (2021); (9) Heintz et al. (2020); (10) Law et al. (2020); (11) Ravi et al.
(2022).

For the constant case, there are six free parameters, including the IGM baryon
fraction $f_{\mathrm{IGM,0}}$, the cosmological parameters ($\Omega_{m}$,
$\Omega_{b}h^{2}$, and $H_{0}$), the SN absolute magnitude $M_{B}$, and the
Milky Way halo DM contribution $\mathrm{DM_{halo}^{MW}}$. The 1D marginalized
posterior distributions and 2D plots of the $1-2\sigma$ confidence regions for
these six parameters are displayed in Figure 2. These contours show that, at
the $1\sigma$ confidence level, the inferred parameter values are
$f_{\mathrm{IGM,0}}=0.927\pm 0.075$, $\Omega_{m}=0.309\pm 0.006$,
$\Omega_{b}h^{2}=0.02245\pm 0.00013$, $H_{0}=67.78\pm 0.44$
$\mathrm{km\;s^{-1}\;Mpc^{-1}}$, and $\mathrm{DM_{halo}^{MW}}=47\pm 10$
$\mathrm{pc\;cm^{-3}}$. The corresponding results are summarized in Table 2.
Figure 1 also shows the theoretical curve for
$\langle\mathrm{DM_{IGM}}\rangle$ versus $z$ for the constant case and a model
estimate of the scatter (95% interval) due to the uncertainties of the
inferred parameters. The theoretical curve reflects the trend of the data
well. Moreover, the inferred value of the IGM baryon fraction is compatible
with previous results obtained from observations (e.g., Fukugita et al. 1998;
Fukugita & Peebles 2004; Shull et al. 2012; Hill et al. 2016; Muñoz & Loeb
2018) and simulations (e.g., Cen & Ostriker 1999, 2006). The constraint
accuracy of $f_{\mathrm{IGM,0}}$ is about 8.0%.

Figure 3: Same as Figure 2, but now for the time-dependent case of
$f_{\mathrm{IGM}}(z)=f_{\mathrm{IGM},0}(1+\alpha\frac{z}{1+z})$.

For the time-dependent case, the free parameters are
$\\{f_{\mathrm{IGM,0}},\;\alpha,\;\Omega_{m},\;\Omega_{b}h^{2},\;H_{0},\;M_{B},\;\mathrm{DM_{halo}^{MW}}\\}$.
These seven parameters are constrained to be $f_{\mathrm{IGM,0}}=0.837\pm
0.089$, $\alpha<0.882$, $\Omega_{m}=0.309\pm 0.006$,
$\Omega_{b}h^{2}=0.02245\pm 0.00013$, $H_{0}=67.71\pm 0.42$
$\mathrm{km\;s^{-1}\;Mpc^{-1}}$, and $\mathrm{DM_{halo}^{MW}}=49\pm 10$
$\mathrm{pc\;cm^{-3}}$, which are displayed in Figure 3 and summarized in
Table 2. Note that with the requirement that $f_{\mathrm{IGM}}$ grows with
redshift ($\alpha\geq 0$), only an upper limit on $\alpha$ can be estimated,
which implies that there is no strong evidence for the redshift dependence of
$f_{\mathrm{IGM}}$. This is consistent with the cosmology-insensitive result
obtained from five localized FRBs (Li et al., 2020). The comparison between
columns 2 and 3 of Table 2 suggests that the nuisance parameters
($\Omega_{m}$, $\Omega_{b}h^{2}$, $H_{0}$, $M_{B}$, and
$\mathrm{DM_{halo}^{MW}}$) are almost identical and have little effect on the
adopted parametric model of $f_{\mathrm{IGM}}$.

Because the constant model and the time-dependent model do not have the same
number of free parameters, a comparison of the likelihoods for either being
closer to the correct model must be based on model selection criteria. We use
the Akaike Information Criterion (AIC; Akaike 1974, 1981) to test the
statistical performance of the models, $\mathrm{AIC}=-2\ln\mathcal{L}+2p$,
where $p$ is the number of free parameters. With $\mathrm{AIC}_{1}$ and
$\mathrm{AIC}_{2}$ characterizing models $\mathcal{M}_{1}$ (the constant
model) and $\mathcal{M}_{2}$ (the time-dependent model), respectively, the
difference $\Delta\mathrm{AIC}\equiv\mathrm{AIC}_{2}-\mathrm{AIC}_{1}$
determines the extent to which $\mathcal{M}_{1}$ is favoured over
$\mathcal{M}_{2}$. The evidence of $\mathcal{M}_{1}$ being correct is judged
‘weak’ when the outcome $\Delta\equiv\mathrm{AIC}_{2}-\mathrm{AIC}_{1}$ is in
the range $0<\Delta<2$, ‘positive’ when $2<\Delta<6$, and ‘strong’ when
$\Delta>6$. Therefore, the outcome $\Delta\mathrm{AIC}=2.244$ shows a positive
evidence in favor of the constant model ($\mathcal{M}_{1}$) with respect to
the time-dependent model ($\mathcal{M}_{2}$). Nevertheless, we hold the
opinion that this positive evidence may be due to the relatively low redshifts
in the FRB data. To distinguish between the constant and time-dependent models
better, a larger number of FRBs localized at higher redshifts is required in
the future.

Figure 4: Comparison of $f_{\mathrm{IGM,0}}$ constraints with previous results
in the constant (left panel) and time-dependent (right panel) cases. The
shaded areas in both panels represent the 1$\sigma$ and 2$\sigma$ confidence
regions of the corresponding $f_{\mathrm{IGM,0}}$ constraints.

In our analysis, the $\mathrm{DM_{ISM}^{MW}}$ values are estimated from the
Galactic electron density model of NE2001 (Cordes & Lazio, 2002). We also
perform a parallel comparative analysis of the FRB data using the YMW16 model
(Yao et al., 2017). The resulting constraints now turn to be
$f_{\mathrm{IGM,0}}=0.901\pm 0.081$
($f_{\mathrm{IGM,0}}=0.788^{+0.082}_{-0.100}$) for the constant (time-
dependent) case. Comparing these inferred $f_{\mathrm{IGM,0}}$ with those
obtained using the NE2001 model, we see that the adoption of a different
electron distribution model has a minimal influence on the results. Due to the
larger $\mathrm{DM_{ISM}^{MW}}$ contribution at low Galactic latitudes in the
YMW16 Model, the derived $f_{\mathrm{IGM,0}}$ values are slightly smaller than
those in the NE2001 model.

To make a direct comparison with previous works, in Figure 4 we plot some
typical $f_{\mathrm{IGM,0}}$ constraints from different FRB samples, as well
as our constraints from both the NE2001 and YMW16 models. One can see from
Figure 4 that our $f_{\mathrm{IGM,0}}$ constraints are well consistent with
previous results at the $1\sigma$ confidence level.

##  5 Conclusions

The $\mathrm{DM_{IGM}}$-$z$ relation of FRBs has been used for probing the
baryon fraction in the IGM, $f_{\mathrm{IGM}}$. However, such studies have
been restricted by the strong degeneracy between cosmological parameters and
$f_{\mathrm{IGM}}$. Moreover, the DM contribution from the IGM
($\mathrm{DM_{IGM}}$) cannot be effectively distinguished from other DMs
contributed by the Milky Way or host galaxy. In this work, we investigate
precise constraints on $f_{\mathrm{IGM}}$ from the DM measurements of
seventeen localized FRBs. In order to break the parameter degeneracy, we
combine FRB data with three other cosmological probes (including SNe Ia, BAO,
and CMB) to infer cosmological parameters and $f_{\mathrm{IGM}}$
simultaneously. To avoid uncontrollable systematic errors induced by the
$\mathrm{DM_{IGM}}$ deduction, we handle the DM contributions of the host
galaxies and IGM as the probability distributions derived from the the
IllustrisTNG simulation.

Following the analysis method described in Section 2, we explore the possible
redshift dependence of $f_{\mathrm{IGM}}(z)$ considering two different
parametric models, which are expressed as the constant and time-dependent
parameterizations given by Equations (20) and (21). The MCMC analysis is used
to constrain $f_{\mathrm{IGM}}(z)$ and other cosmological parameters. For the
constant model, we infer that $f_{\mathrm{IGM,0}}=0.927\pm 0.075$,
representing a precision of 8.0%. This constraint from FRB observations is
roughly consistent with those obtained from other probes (Fukugita et al.,
1998; Fukugita & Peebles, 2004; Shull et al., 2012; Hill et al., 2016). For
the time-dependent model, whereas only an upper limit on the evolution index
$\alpha$ can be set ($\alpha<0.882$), we can obtain a good limit on the local
$f_{\mathrm{IGM,0}}=0.837\pm 0.089$, which is slightly looser but still
consistent with previous results derived from different methods. According to
the AIC model selection criteria, there is a mild evidence suggesting that the
constant model is preferred over the time-dependent model. However, due to the
fact that the number of current localized FRBs is small and their redshift
measurements are relatively low, we cannot safely exclude the possibility of
an evolving $f_{\mathrm{IGM}}(z)$.

Redshift measurements of a larger sample of FRBs are essential for using this
method presented here to constrain $f_{\mathrm{IGM}}$ and its possible
redshift evolution. Forthcoming radio telescopes such as the Deep Synoptic
Array 2000-dish prototype (Hallinan et al., 2019) and the Square Kilometre
Array (Dewdney et al., 2009), with improved detection sensitivity and
localization capability, will be able to increase the current localized FRB
sample size by orders of magnitude. With the rapid progress in localizing
FRBs, the constraints on $f_{\mathrm{IGM}}$ will be significantly improved,
and the baryon distribution of the universe will be better understood.

We are grateful to the anonymous referee for their helpful comments. This work
is partially supported by the National Key Research and Development Program of
China (2022SKA0130100), the National Natural Science Foundation of China
(grant Nos. 11725314 and 12041306), the Key Research Program of Frontier
Sciences (grant No. ZDBS-LY-7014) of Chinese Academy of Sciences,
International Partnership Program of Chinese Academy of Sciences for Grand
Challenges (114332KYSB20210018), the CAS Project for Young Scientists in Basic
Research (grant No. YSBR-063), the CAS Organizational Scientific Research
Platform for National Major Scientific and Technological Infrastructure:
Cosmic Transients with FAST, the Natural Science Foundation of Jiangsu
Province (grant No. BK20221562), and the Young Elite Scientists Sponsorship
Program of Jiangsu Association for Science and Technology.

## References

  * Akaike (1974) Akaike, H. 1974, IEEE transactions on automatic control, 19, 716
  * Akaike (1981) —. 1981, Journal of econometrics, 16, 3
  * Aubourg et al. (2015) Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516, doi: 10.1103/PhysRevD.92.123516
  * Bannister et al. (2019) Bannister, K. W., Deller, A. T., Phillips, C., et al. 2019, Science, 365, 565, doi: 10.1126/science.aaw5903
  * Becker et al. (2011) Becker, G. D., Bolton, J. S., Haehnelt, M. G., & Sargent, W. L. W. 2011, MNRAS, 410, 1096, doi: 10.1111/j.1365-2966.2010.17507.x
  * Bhandari et al. (2020) Bhandari, S., Sadler, E. M., Prochaska, J. X., et al. 2020, ApJ, 895, L37, doi: 10.3847/2041-8213/ab672e
  * Bhandari et al. (2022) Bhandari, S., Heintz, K. E., Aggarwal, K., et al. 2022, AJ, 163, 69, doi: 10.3847/1538-3881/ac3aec
  * Bhardwaj et al. (2021a) Bhardwaj, M., Kirichenko, A. Y., Michilli, D., et al. 2021a, ApJ, 919, L24, doi: 10.3847/2041-8213/ac223b
  * Bhardwaj et al. (2021b) Bhardwaj, M., Gaensler, B. M., Kaspi, V. M., et al. 2021b, ApJ, 910, L18, doi: 10.3847/2041-8213/abeaa6
  * Cen & Ostriker (1999) Cen, R., & Ostriker, J. P. 1999, ApJ, 514, 1, doi: 10.1086/306949
  * Cen & Ostriker (2006) —. 2006, ApJ, 650, 560, doi: 10.1086/506505
  * Chatterjee et al. (2017) Chatterjee, S., Law, C. J., Wharton, R. S., et al. 2017, Nature, 541, 58, doi: 10.1038/nature20797
  * Chen et al. (2019) Chen, L., Huang, Q.-G., & Wang, K. 2019, J. Cosmology Astropart. Phys, 2019, 028, doi: 10.1088/1475-7516/2019/02/028
  * CHIME/FRB Collaboration et al. (2021) CHIME/FRB Collaboration, Amiri, M., Andersen, B. C., et al. 2021, ApJS, 257, 59, doi: 10.3847/1538-4365/ac33ab
  * Chittidi et al. (2021) Chittidi, J. S., Simha, S., Mannings, A., et al. 2021, ApJ, 922, 173, doi: 10.3847/1538-4357/ac2818
  * Cordes & Lazio (2002) Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0207156
  * Dai & Xia (2021) Dai, J.-P., & Xia, J.-Q. 2021, MNRAS, 503, 4576, doi: 10.1093/mnras/stab785
  * Deng & Zhang (2014) Deng, W., & Zhang, B. 2014, ApJ, 783, L35, doi: 10.1088/2041-8205/783/2/L35
  * Dewdney et al. (2009) Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proceedings, 97, 1482, doi: 10.1109/JPROC.2009.2021005
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fukugita et al. (1998) Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518, doi: 10.1086/306025
  * Fukugita & Peebles (2004) Fukugita, M., & Peebles, P. J. E. 2004, ApJ, 616, 643, doi: 10.1086/425155
  * Gao et al. (2014) Gao, H., Li, Z., & Zhang, B. 2014, ApJ, 788, 189, doi: 10.1088/0004-637X/788/2/189
  * Hagstotz et al. (2022) Hagstotz, S., Reischke, R., & Lilow, R. 2022, MNRAS, 511, 662, doi: 10.1093/mnras/stac077
  * Hallinan et al. (2019) Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 255.  https://arxiv.org/abs/1907.07648
  * Hashimoto et al. (2021) Hashimoto, T., Goto, T., Lu, T.-Y., et al. 2021, MNRAS, 502, 2346, doi: 10.1093/mnras/stab186
  * Heintz et al. (2020) Heintz, K. E., Prochaska, J. X., Simha, S., et al. 2020, ApJ, 903, 152, doi: 10.3847/1538-4357/abb6fb
  * Hill et al. (2016) Hill, J. C., Ferraro, S., Battaglia, N., Liu, J., & Spergel, D. N. 2016, Phys. Rev. Lett., 117, 051301, doi: 10.1103/PhysRevLett.117.051301
  * Jaroszynski (2019) Jaroszynski, M. 2019, MNRAS, 484, 1637, doi: 10.1093/mnras/sty3529
  * Kirsten et al. (2022) Kirsten, F., Marcote, B., Nimmo, K., et al. 2022, Nature, 602, 585, doi: 10.1038/s41586-021-04354-w
  * Kumar & Linder (2019) Kumar, P., & Linder, E. V. 2019, Phys. Rev. D, 100, 083533, doi: 10.1103/PhysRevD.100.083533
  * Law et al. (2020) Law, C. J., Butler, B. J., Prochaska, J. X., et al. 2020, ApJ, 899, 161, doi: 10.3847/1538-4357/aba4ac
  * Lemos et al. (2022) Lemos, T., Gonçalves, R. S., Carvalho, J. C., & Alcaniz, J. S. 2022, arXiv e-prints, arXiv:2205.07926.  https://arxiv.org/abs/2205.07926
  * Li et al. (2019) Li, Z., Gao, H., Wei, J.-J., et al. 2019, ApJ, 876, 146, doi: 10.3847/1538-4357/ab18fe
  * Li et al. (2020) Li, Z., Gao, H., Wei, J. J., et al. 2020, MNRAS, 496, L28, doi: 10.1093/mnrasl/slaa070
  * Li et al. (2018) Li, Z.-X., Gao, H., Ding, X.-H., Wang, G.-J., & Zhang, B. 2018, Nature Communications, 9, 3833, doi: 10.1038/s41467-018-06303-0
  * Lin et al. (2022) Lin, H.-N., Li, X., & Tang, L. 2022, Chinese Physics C, 46, 075102, doi: 10.1088/1674-1137/ac5e92
  * Lorimer et al. (2007) Lorimer, D. R., Bailes, M., McLaughlin, M. A., Narkevic, D. J., & Crawford, F. 2007, Science, 318, 777, doi: 10.1126/science.1147532
  * Macquart et al. (2020) Macquart, J. P., Prochaska, J. X., McQuinn, M., et al. 2020, Nature, 581, 391, doi: 10.1038/s41586-020-2300-2
  * Marcote et al. (2020) Marcote, B., Nimmo, K., Hessels, J. W. T., et al. 2020, Nature, 577, 190, doi: 10.1038/s41586-019-1866-z
  * McQuinn (2014) McQuinn, M. 2014, ApJ, 780, L33, doi: 10.1088/2041-8205/780/2/L33
  * Meiksin (2009) Meiksin, A. A. 2009, Reviews of Modern Physics, 81, 1405, doi: 10.1103/RevModPhys.81.1405
  * Muñoz et al. (2016) Muñoz, J. B., Kovetz, E. D., Dai, L., & Kamionkowski, M. 2016, Phys. Rev. Lett., 117, 091301, doi: 10.1103/PhysRevLett.117.091301
  * Muñoz & Loeb (2018) Muñoz, J. B., & Loeb, A. 2018, Phys. Rev. D, 98, 103518, doi: 10.1103/PhysRevD.98.103518
  * Niu et al. (2022) Niu, C. H., Aggarwal, K., Li, D., et al. 2022, Nature, 606, 873, doi: 10.1038/s41586-022-04755-5
  * Petroff et al. (2016) Petroff, E., Barr, E. D., Jameson, A., et al. 2016, PASA, 33, e045, doi: 10.1017/pasa.2016.35
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Platts et al. (2019) Platts, E., Weltman, A., Walters, A., et al. 2019, Phys. Rep., 821, 1, doi: 10.1016/j.physrep.2019.06.003
  * Prochaska & Zheng (2019) Prochaska, J. X., & Zheng, Y. 2019, MNRAS, 485, 648, doi: 10.1093/mnras/stz261
  * Prochaska et al. (2019) Prochaska, J. X., Macquart, J.-P., McQuinn, M., et al. 2019, Science, 366, 231, doi: 10.1126/science.aay0073
  * Qiang & Wei (2020) Qiang, D.-C., & Wei, H. 2020, J. Cosmology Astropart. Phys, 2020, 023, doi: 10.1088/1475-7516/2020/04/023
  * Qiu et al. (2022) Qiu, X.-W., Zhao, Z.-W., Wang, L.-F., Zhang, J.-F., & Zhang, X. 2022, J. Cosmology Astropart. Phys, 2022, 006, doi: 10.1088/1475-7516/2022/02/006
  * Ravi et al. (2019) Ravi, V., Catha, M., D’Addario, L., et al. 2019, Nature, 572, 352, doi: 10.1038/s41586-019-1389-7
  * Ravi et al. (2022) Ravi, V., Law, C. J., Li, D., et al. 2022, MNRAS, 513, 982, doi: 10.1093/mnras/stac465
  * Ryan et al. (2019) Ryan, J., Chen, Y., & Ratra, B. 2019, MNRAS, 488, 3844, doi: 10.1093/mnras/stz1966
  * Ryder et al. (2022) Ryder, S. D., Bannister, K. W., Bhandari, S., et al. 2022, arXiv e-prints, arXiv:2210.04680.  https://arxiv.org/abs/2210.04680
  * Scolnic et al. (2018) Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101, doi: 10.3847/1538-4357/aab9bb
  * Shull et al. (2012) Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23, doi: 10.1088/0004-637X/759/1/23
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Thornton et al. (2013) Thornton, D., Stappers, B., Bailes, M., et al. 2013, Science, 341, 53, doi: 10.1126/science.1236789
  * Walters et al. (2019) Walters, A., Ma, Y.-Z., Sievers, J., & Weltman, A. 2019, Phys. Rev. D, 100, 103519, doi: 10.1103/PhysRevD.100.103519
  * Walters et al. (2018) Walters, A., Weltman, A., Gaensler, B. M., Ma, Y.-Z., & Witzemann, A. 2018, ApJ, 856, 65, doi: 10.3847/1538-4357/aaaf6b
  * Wang & Wang (2018) Wang, Y. K., & Wang, F. Y. 2018, A&A, 614, A50, doi: 10.1051/0004-6361/201731160
  * Wei et al. (2019) Wei, J.-J., Li, Z., Gao, H., & Wu, X.-F. 2019, J. Cosmology Astropart. Phys, 2019, 039, doi: 10.1088/1475-7516/2019/09/039
  * Wei et al. (2018) Wei, J.-J., Wu, X.-F., & Gao, H. 2018, ApJ, 860, L7, doi: 10.3847/2041-8213/aac8e2
  * Wu et al. (2022) Wu, Q., Zhang, G.-Q., & Wang, F.-Y. 2022, MNRAS, 515, L1, doi: 10.1093/mnrasl/slac022
  * Xiao et al. (2021) Xiao, D., Wang, F., & Dai, Z. 2021, Science China Physics, Mechanics, and Astronomy, 64, 249501, doi: 10.1007/s11433-020-1661-7
  * Yang et al. (2022) Yang, K. B., Wu, Q., & Wang, F. Y. 2022, ApJ, 940, L29, doi: 10.3847/2041-8213/aca145
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Yu & Wang (2017) Yu, H., & Wang, F. Y. 2017, A&A, 606, A3, doi: 10.1051/0004-6361/201731607
  * Zhang (2022) Zhang, B. 2022, arXiv e-prints, arXiv:2212.03972.  https://arxiv.org/abs/2212.03972
  * Zhang et al. (2020) Zhang, G. Q., Yu, H., He, J. H., & Wang, F. Y. 2020, ApJ, 900, 170, doi: 10.3847/1538-4357/abaa4a
  * Zhang et al. (2021) Zhang, Z. J., Yan, K., Li, C. M., Zhang, G. Q., & Wang, F. Y. 2021, ApJ, 906, 49, doi: 10.3847/1538-4357/abceb9
  * Zhao et al. (2020) Zhao, Z.-W., Li, Z.-X., Qi, J.-Z., et al. 2020, ApJ, 903, 83, doi: 10.3847/1538-4357/abb8ce
  * Zheng et al. (2014) Zheng, Z., Ofek, E. O., Kulkarni, S. R., Neill, J. D., & Juric, M. 2014, ApJ, 797, 71, doi: 10.1088/0004-637X/797/1/71
  * Zhou et al. (2014) Zhou, B., Li, X., Wang, T., Fan, Y.-Z., & Wei, D.-M. 2014, Phys. Rev. D, 89, 107303, doi: 10.1103/PhysRevD.89.107303


# Linear growth index of matter perturbations in Rastall gravity

Wompherdeiki Khyllep${\ }^{a,b}$, Jibitesh Dutta${\ }^{c,d}$
a Department of Mathematics, North-Eastern Hill University, Shillong,
Meghalaya 793022, India
b Department of Mathematics, St. Anthony’s College, Shillong, Meghalaya
793001, India
c Mathematics Division, Department of Basic Sciences and Social Sciences,
North-Eastern
Hill University, Shillong, Meghalaya 793022, India
d Inter University Centre for Astronomy and Astrophysics, Pune 411 007, India
sjwomkhyllep<EMAIL_ADDRESS>
###### Abstract

Rastall gravity theory shows notable features consistent with physical
observations in comparison to the standard Einstein theory. Recently, there
has been a debate about the equivalence of Rastall gravity and general
relativity. Motivated by this open issue, in the present work, we attempt to
shed some light on this debate by analyzing the evolution of the Rastall based
cosmological model at the background as well as perturbation level. Employing
the dynamical system techniques, we found that at late times, the dynamics of
the model resembles the $\Lambda$CDM model at the background level
irrespective of the choice of Rastall’s parameter. However, at the
perturbation level, we found that the evolution of the growth index heavily
depends on the Rastall’s parameter and displays a significant deviation from
the $\Lambda$CDM model.

##  1 Introduction

The inability of the standard theory of gravity to describe various
observational evidence has led to the need for extending the theory of general
relativity (GR). The structure of extended theories of gravity provides a
useful approach to alleviate the fundamental problems such as dark matter (DM)
and dark energy (DE) associated with the standard model of GR [1, 2, 3].

One of the extended theories of gravity was proposed by Rastall in 1972 [4,
5]. In his theory, Rastall assumed the violation of the conservation of
energy-momentum tensor in curved spacetime without dropping the Bianchi
identities. As the usual particle creation phenomenon violates the energy-
momentum conservation law, therefore, this theory can also be regarded as a
classical formulation of a particle creation process [6, 7, 8, 9]. Rastall
gravity is one of the most promising non-conservative modified theory of
gravity supporting various cosmological and astrophysical observations.
Recently, it got a lot of attention in the physics community [10, 11, 12, 13,
14]. For instance, this theory fits well with observational data related to
the age of the Universe, Hubble parameter and helium nucleosynthesis [15, 16].
Further, this theory can possibly circumvent the entropy and age problems of
standard cosmology [17]. From Mach’s principle perspective, this theory is
more ‘Machian’ than the standard theory of GR [18].

It has been claimed by Visser in [19] that the theory of Rastall gravity is
completely equivalent to GR and further supported by [20] from a
thermodynamical perspective. However, this result has been recently denied by
[21, 22] leaving an open debate that whether this theory is a modified gravity
theory or it is equivalent to GR with an additional modified matter content.
Earlier, the non-equivalence of Rastall gravity from GR has also been pointed
by Smalley in Ref. [23].

In Ref. [11] it has been shown that Rastall cosmology is equivalent to the
$\Lambda$CDM at both background and linear perturbation level except that DE
cluster for the former model. The growth of matter perturbations provides an
efficient approach to predict the matter distribution of the Universe and also
to discriminate various gravitational theories [24, 25]. One of the simple
observational tools used to study the growth history of a model is the so-
called growth index of matter perturbations (denoted by $\gamma$) [26]. The
accurate estimation of the growth index is one of the basic tasks from the
cosmological point of view as it can be used as a tool to test the validity of
GR on extragalactic scales. It is a usual procedure that for each cosmological
model, one requires to analyze its background evolution and the growth index
of matter perturbations. This task helps one to get an overall effect of the
model at the cosmological and astrophysical level. For instance, those DE
models based on GR, the value of the constant growth index can be reduced to
that of the $\Lambda$CDM model (i.e. $\gamma=\frac{6}{11}$) irrespective of
the choice of model parameters. However, in the case of modified gravity based
models, the value shows a significant deviation from $\frac{6}{11}$. Recently,
there are various work in the literature where an analytical form of the
growth index for different models are obtained which include the scalar field
DE [27, 28, 29], DGP[30, 31], $f(R)$ gravity [32, 33, 34], Finsler-Randers
[35], time-varying vacuum models [36], clustered DE [37], holographic DE [38]
and $f(T)$ gravity [39].

At the level of linear perturbation, it has been obtained that the evolution
equation of the growth of matter perturbations for Rastall based cosmological
model coincides with that of the $\Lambda$CDM model. Additionally, in the case
of Rastall model, DE perturbations cluster even if it is in the vacuum energy
form [11]. This might have a strong impact on the value of the growth index,
leading to interesting cosmological signatures of the model which can be
tested with growth rate data. It is therefore interesting to analyze the
growth index of the Rastall model in order to shed some light on the said
debate.

With this motivation the aim of the paper is two-fold. First, we analyze the
dynamics of the Rastall model at the background level using the dynamical
system techniques. The dynamical system analysis is performed in Sec. 3 to
qualitatively extract the evolution described by the basic cosmological
equations presented in Sec. 2. Secondly, we investigate the behavior of this
model at the linear perturbation level in Sec. 4. For this we shall derive the
growth index of the Rastall model and then compare with that of the
$\Lambda$CDM in subsections 4.0.1 and 4.0.2. Finally, the conclusion is given
in Sec. 5.

##  2 Basic cosmological equations of Rastall gravity

In this section, we briefly introduce the basic cosmological equations of
Rastall theory of gravity. In many gravitational theories, the source of
energy-momentum is determined by a vanishing divergence tensor, minimally
coupled to the geometry. However, in this theory, the usual conservation law
of the energy-momentum tensor is not satisfied, instead, it is assumed to
satisfy the following relation [4, 5]

$\displaystyle T^{\mu\nu}_{;\mu}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{8\pi G}R^{;\nu}\,,$ (1)

where semi-colon denotes the covariant derivative, $G$ is the Newton’s
gravitational constant, $T_{\mu\nu}$ is the energy-momentum tensor and $R$ is
the Ricci scalar. Here, $\lambda$ is the coupling constant which measure the
exchange of energy between geometry and matter field. When $\lambda=0$, one
obtain a usual conservation law in GR. The modified Einstein’s field equations
in Rastall gravity framework are given by

$\displaystyle G_{\mu\nu}$ $\displaystyle=$ $\displaystyle 8\pi
G\Big{(}T_{\mu\nu}-\frac{\tilde{\gamma}-1}{2}g_{\mu\nu}T\Big{)}\,,$ (2)
$\displaystyle T^{\mu\nu}_{;\mu}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\gamma}-1}{2}T^{;\nu}\,,$ (3)

where $T$ is the trace of the energy-momentum tensor and $\tilde{\gamma}$ is
the Rastall’s parameter which is related to $\lambda$ as

$\displaystyle\tilde{\gamma}$ $\displaystyle=$
$\displaystyle\frac{1+6\lambda}{1+4\lambda}\,.$ (4)

Clearly, $\tilde{\gamma}=1$ corresponds to the GR case.

We now consider a two perfect fluid system consisting of a pressureless matter
with corresponding energy density $\rho_{m}$ and an exotic component with
energy density $\rho_{x}$ which accounts for acceleration of the universe and
whose equation of state (EoS) is $w_{x}$. It is worth mentioning here that an
extra exotic fluid $\rho_{x}$ behaving as a source of DE is required similar
to the GR case [11]. Further, we assume that the matter component satisfies
the usual conservation law. Therefore, the modified Einstein’s field equations
(2)-(3) become

$\displaystyle G_{\mu\nu}$ $\displaystyle=$ $\displaystyle 8\pi
G\Big{(}T^{m}_{\mu\nu}+T^{x}_{\mu\nu}-\frac{\tilde{\gamma}-1}{2}g_{\mu\nu}(T^{m}+T^{x})\Big{)}\,,$
(5) $\displaystyle T^{\mu\nu}_{;\mu}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\gamma}-1}{2}(T_{x}^{;\nu}+T_{m}^{;\nu})\,,\,\,\,~{}~{}~{}T_{m;\mu}^{\mu\nu}=0\,,$
(6)

where index $m$ denotes the matter component and $x$ denotes the DE component.

At large scale, various observational data favor a homogeneous and isotropic
universe which can be described by a spatially flat Friedmann-Lemaitre-
Robertson-Walker (FLRW) metric given by

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}\,,~{}~{}~{}~{}(\mu,\nu=0,1,2,3)$ (7)

with $g_{\mu\nu}=$diag
$(1,-a^{2}(t),-a^{2}(t)\,r^{2},-a^{2}(t)\,r^{2}\,\sin^{2}\theta)$. Here,
$a(t)$ is a scale factor and $t$ denotes the cosmological time. The above
Einstein field equations (5), (6) with respect to the metric (7) can be
written as

$\displaystyle H^{2}$ $\displaystyle=$ $\displaystyle\frac{4\pi
G}{3}\Big{[}\rho_{x}\,(3-\tilde{\gamma}-3(1-\tilde{\gamma})w_{x})$ (8)
$\displaystyle+(3-\tilde{\gamma})\rho_{m}\Big{]}\,,$
$\displaystyle\dot{H}+H^{2}$ $\displaystyle=$ $\displaystyle\frac{4\pi
G}{3}\Big{[}\rho_{x}\,(3(\tilde{\gamma}-2)w_{x}-\tilde{\gamma})-\tilde{\gamma}\rho_{m}\Big{]},~{}~{}~{}~{}$
(9) $\displaystyle\dot{\rho}_{x}+3H(1+w_{x})\rho_{x}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\gamma}-1}{2}\Big{[}(1-3w_{x})\dot{\rho}_{x}+\dot{\rho}_{m}\Big{]}\,,$
(10) $\displaystyle\dot{\rho}_{m}+3H\rho_{m}$ $\displaystyle=$ $\displaystyle
0\,,$ (11)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter and the overdot denotes
derivative with respect to $t$. On solving Eqs. (10) and (11), one can obtain
the evolution equation of $\rho_{m}$ and $\rho_{x}$ as [13]

$\displaystyle\rho_{x}$ $\displaystyle=$
$\displaystyle\rho_{d0}\,a^{-\frac{6(1+w_{x})}{2+(1-\tilde{\gamma})(1-3w_{x})}}+\frac{(1-\tilde{\gamma})\rho_{m0}a^{-3}}{2w_{x}+(\tilde{\gamma}-1)(1-3w_{x})},~{}~{}~{}~{}$
(12) $\displaystyle\rho_{m}$ $\displaystyle=$ $\displaystyle\rho_{m0}a^{-3},$
(13)

where $\rho_{d0}$ is given by

$\displaystyle\rho_{d0}$ $\displaystyle=$
$\displaystyle\rho_{x0}+\frac{\tilde{\gamma}-1}{2w_{x}+(\tilde{\gamma}-1)(1-3w_{x})}\rho_{m0}.$
(14)

Eq. (12) shows that the energy density of an exotic fluid consists of two
components and one of them behaves as DM. In the above equations $\rho_{m0}$,
$\rho_{x0}$ and $\rho_{d0}$ denote the present value of matter density, of an
exotic fluid and of a DE component of an exotic fluid. With the aid of (12),
(13) and (14), the Friedmann equation (8) can be written as

$\displaystyle E^{2}(a)$ $\displaystyle\textstyle=$
$\displaystyle\textstyle\Big{[}\Big{(}\frac{A}{2}\Omega_{m0}+\frac{B}{2}\Omega_{x0}\Big{)}\,\textstyle{a}^{-{\frac{6(1+w_{x})}{B}}}+\frac{A}{2}\Omega_{m0}{a}^{-3}\Big{]}$
(15)
$\displaystyle\textstyle+\,\frac{\left(3-\tilde{\gamma}\right)}{2}\,\Omega_{m0}\,{a}^{-3}\,,$
$\displaystyle\text{where}~{}A$ $\displaystyle=$
$\displaystyle\frac{B(\tilde{\gamma}-1)}{2(w_{x}+1)-B},~{}B=2-(\tilde{\gamma}-1)(1-3w_{x}).$

In the above equation, $E(a)=\frac{H(a)}{H_{0}}$, $\Omega_{m0}=\frac{8\pi
G\rho_{m0}}{3H_{0}^{2}}$, $\Omega_{x0}=\frac{8\pi G\rho_{x0}}{3H_{0}^{2}}$ and
$H_{0}$ is the present value of the Hubble parameter. Note that the expression
within the square bracket of Eq. (15) denotes the energy contribution from an
exotic fluid with one of its components behaving as DM. The term outside the
square bracket denotes the contribution solely from a pressureless matter
fluid. As expected, for $\tilde{\gamma}=1$, the above equation reduces to the
corresponding GR case given by

$\displaystyle E_{GR}^{2}(a)$ $\displaystyle\textstyle=$
$\displaystyle\textstyle\Omega_{x0}\,\textstyle{a}^{-3(1+w_{x})}+\Omega_{m0}\,{a}^{-3}\,,$
(16)

with $x$-component playing the role of DE and $m$-component plays the role of
DM. On differentiating (15) with respect to the scale factor, we obtain

$\frac{d\ln E}{d\ln
a}=-\frac{3\,\Omega_{m}\,w_{x}\,(2\tilde{\gamma}-3)}{3\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3w_{x}+3}\,.$
(17)

Further, from Eqs. (8) and (9), we can define the effective energy density and
pressure as

$\displaystyle p_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{\rho_{x}}{2}(\tilde{\gamma}-1-w_{x}(3\tilde{\gamma}-5))+\frac{\tilde{\gamma}-1}{2}\rho_{m}\,,$
(18) $\displaystyle\rho_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{\rho_{x}}{2}(3-\tilde{\gamma}-3(1-\tilde{\gamma})w_{x})+\frac{3-\tilde{\gamma}}{2}\rho_{m}\,.$
(19)

Then the effective EoS is given by

$\displaystyle w_{\rm eff}=\frac{p_{\rm eff}}{\rho_{\rm
eff}}=-1+\frac{2}{3-\tilde{\gamma}}[1+(3-2\tilde{\gamma})w_{x}\,\Omega_{x}].$
(20)

For accelerated universe, we must have $w_{\rm eff}<-\frac{1}{3}$. In order to
obtain a concrete picture on the background cosmological evolution described
by this model, in the next section, we shall perform a dynamical system
analysis for this model.

##  3 Dynamical system analysis of Rastall model

In this section, we shall perform a complete dynamical system analysis of the
Rastall gravity based model, where we introduce suitable dimensionless
variables to recast the cosmological equations (8)-(11) into an autonomous
dynamical system. In this case, we choose the dimensionless energy density
parameters of matter and of an exotic fluid given by

$\displaystyle\Omega_{m}=\frac{8\pi
G\rho_{m}}{3H^{2}},~{}~{}~{}\Omega_{x}=\frac{8\pi G\rho_{x}}{3H^{2}}\,,$ (21)

respectively, as the dynamical variables of the system which are related by
Eq. (8) as

$\displaystyle\Omega_{x}(3-\tilde{\gamma}-3(1-\tilde{\gamma})w_{x})+(3-\tilde{\gamma})\Omega_{m}=2.$
(22)

It may be noted that depending on the choice of parameters $\tilde{\gamma}$
and $w_{x}$, there is a possibility that either $\Omega_{m}$ or $\Omega_{x}$
or both assume negative values. This feature usually arises in cosmological
models where a phenomenological interaction between different components is
present and it is also related to an uncertainty in defining the relative
energy densities of the model [40]. Since in this Rastall cosmological model
the exotic fluid $\rho_{x}$ and matter fluid $\rho_{m}$ can be interpreted as
effectively interacting quantities (cf. Eq. (10)), it is therefore not
surprising that $\Omega_{m}$ or $\Omega_{x}$ take negative values. Further, it
is important to note here that the condition $\Omega_{x}<0$ is actually
obligated in several $f(R)$ theories in order to alleviate the coincidence
problem and also to match with the current observations [41]. While one could
argue that $\Omega_{m}<0$ is physically unacceptable, the present model
demands this condition (cf. Eq. (22)) and discarding it may lead to incomplete
dynamics. For relevant recent work one can also see [42, 43]. Again, it is
interesting to observe that negative energy density does not always imply
phantom behavior in general (one can see [40], [43]). For instance, in the
present model, even if negative energy density is allowed, it is possible that
both the DE component ($x$-component) and the overall behavior is non-phantom
i.e. $w_{x}>-1$ and $w_{\rm eff}>-1$. This can be seen from Eq. (20) (for
$\tilde{\gamma}=\frac{3}{2}$, $w_{\rm eff}>-1$ even if $\Omega_{x}<0$). As
expected for $\tilde{\gamma}=1$, the relation (22) reduces to that of the GR
case i.e. $\Omega_{m}+\Omega_{x}=1$.

Figure 1: The phase portrait of the system (23) describing the stability of
critical points of the system. Here, $\tilde{\gamma}=0.9$ and $w_{x}=-1$.
Figure 2: The evolution of the relative matter energy density $\Omega_{m}$,
the relative DE density $\Omega_{x}$ and effective EoS $w_{\rm eff}$. Here,
$\tilde{\gamma}=0.9$ and $w_{x}=-1$.

From the above constraint equation (22), it can be seen that one of the
dimensionless density parameters is linearly independent. Then the above
governing cosmological equations (8)-(11) can be recast into a following one-
dimensional system:

$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3}{2(\tilde{\gamma}-3)}\,\Big{(}3w_{x}\Omega_{x}\left(\tilde{\gamma}^{2}-4\Omega_{x}+3\right)+8\tilde{\gamma}w_{x}\,\Omega_{x}$
(23)
$\displaystyle\times\left(\Omega_{x}-1\right)-\Omega_{x}(\tilde{\gamma}^{2}+1)+2\tilde{\gamma}(\Omega_{x}-1)+2\Big{)}\,,$

where prime denotes a derivative with respect to a new time variable $N=\ln
a$. The critical points for this one-dimensional system are
$\Omega_{x}=\frac{1}{8\,{w_{x}}\,\left(2\,\tilde{\gamma}-3\right)}\,\Big{[}\Delta\pm\Big{(}\Delta^{2}-32\,w_{x}(\tilde{\gamma}-1)\,(2\tilde{\gamma}-3)\Big{)}^{\frac{1}{2}}\Big{]}$
where
$\Delta=-{w_{x}}\,\left(2\,\tilde{\gamma}-3\right)^{2}+\left(\tilde{\gamma}-1\right)^{2}+\tilde{\gamma}\,{w_{x}}\,\left(\tilde{\gamma}-4\right)$.
We denote these critical points by $C_{\pm}$. On linearizing Eq. (23) around
critical points $C_{+}$ and $C_{-}$, one obtains small perturbations
$\delta_{+}$, $\delta_{-}$ respectively as

$\displaystyle\delta_{+}^{\prime}$ $\displaystyle=$
$\displaystyle\lambda_{+}\,\delta_{+}\Rightarrow\delta_{+}(N)=\bar{\delta}_{+}\exp^{\lambda_{+}N}$
(24) $\displaystyle\delta_{-}^{\prime}$ $\displaystyle=$
$\displaystyle\lambda_{-}\,\delta_{-}\Rightarrow\delta_{-}(N)=\bar{\delta}_{-}\exp^{\lambda_{-}N}$
(25)

where
$\lambda_{\pm}=\pm\,\frac{3}{2\,(\tilde{\gamma}-3)}\,\Big{(}\Delta^{2}-32\,w_{x}(\tilde{\gamma}-1)\,(2\tilde{\gamma}-3)\Big{)}^{\frac{1}{2}}$.
Here, $\bar{\delta}_{+}$, $\bar{\delta}_{-}$ denote the values of
${\delta}_{+}$, ${\delta}_{-}$ when $N=0$ respectively. Therefore, for
$\tilde{\gamma}>3$ case, the perturbation corresponds to critical point
$C_{+}$ increases exponentially, but for $C_{-}$ the corresponding
perturbation decreases exponentially. The situation is vice versa for
$\tilde{\gamma}<3$ case. Hence, point $C_{+}$ (point $C_{-}$) corresponds to a
stable (unstable) node for $\tilde{\gamma}<3$ and unstable (stable) node for
$\tilde{\gamma}>3$. Therefore, the local stability of critical points of the
system undergoes an abrupt change as we vary $\tilde{\gamma}$ from a value
smaller than 3 to greater than 3. In particular, the two critical points
exchange their stability behavior as $\tilde{\gamma}$ passes through a value
$\tilde{\gamma}=3$. Hence, the system seems to undergo a bifurcation similar
to the so-called transcritical bifurcation except that in the present scenario
the critical points of the system vanish at $\tilde{\gamma}=3$. This occurs
due to the fact that the system (23) blows up at $\tilde{\gamma}=3$ and
therefore, we cannot extract any dynamics at the background level using the
present dynamical variables. Interestingly, for $\tilde{\gamma}=\frac{3}{2}$,
we have checked that the above system contains only one unstable critical
point which corresponds to a radiation dominated universe ($w_{\rm
eff}=\frac{1}{3}$). However, from Eq. (12) the case where $w_{x}=-1$ and
$\tilde{\gamma}=\frac{3}{2}$ is undesirable as it leads to divergence of
$\rho_{x}$. Further, it can be seen that for $\tilde{\gamma}\simeq 1$ and
$w_{x}\simeq-1$, point $C_{+}$ corresponds to $\Omega_{x}\simeq 1$ and point
$C_{-}$ corresponds to $\Omega_{x}\simeq 0$. In general, for
$\tilde{\gamma}<3$ (except $\tilde{\gamma}=\frac{3}{2}$), point $C_{+}$
corresponds to a DE dominated universe and $C_{-}$ corresponds to a matter
dominated universe (as $\Omega_{x}$ for $C_{+}$ dominates that of $C_{-}$ and
$\Omega_{m}$ for $C_{-}$ dominates that of $C_{+}$). Furthermore, for
$\tilde{\gamma}>3$, point $C_{-}$ corresponds to a DE dominated universe and
$C_{+}$ corresponds to a matter dominated universe. Hence, irrespective of
value of $\tilde{\gamma}$ this model describes the late time transition of the
universe from an unaccelerated DM to an accelerated DE domination at the
background level. We have indeed checked that there are no other critical
points that are hidden towards an infinite regime of $\Omega_{x}$. A
particular case of the phase portrait for the system (23) is given in Fig. 1
for $\tilde{\gamma}=0.9$, $w_{x}=-1$. It is also interesting to note that for
$w_{x}=-1$, the late time attractor corresponds to $w_{\rm eff}=-1$ for any
choices of $\tilde{\gamma}$. Thus, the background evolution of this model
coincides exactly with the $\Lambda$CDM model at late times, however the
evolution during the matter domination depends on the value of the parameter
$\tilde{\gamma}$. This behavior is depicted numerically in Fig. 2 by plotting
the evolution of the DE and DM energy density parameters along with an
effective EoS against redshift $z$ (where $1+z=\frac{a_{0}}{a}$ and $a_{0}$ is
the present value of scale factor taken to be 1). In the next section, we
shall analyze the behavior of this model at the linear perturbation level and
compare the evolution with the $\Lambda$CDM model.

##  4 Analysis of growth index of matter perturbations

Here, we shall study the growth of linear matter fluctuations of the present
model in the matter dominated epoch. We first present the evolution equation
of the matter perturbations of first order within the sub-horizon scale, then
we study the evolution of growth index of matter perturbations for this model.
Following [11, 12], we consider the perturbed metric with synchronous gauge
given by

$\displaystyle\bar{g}_{\mu\nu}$ $\displaystyle=$ $\displaystyle
g_{\mu\nu}+h_{\mu\nu},~{}~{}~{}h_{\mu 0}=0\,,$ (26)

where $\bar{g}_{\mu\nu}$ is the total metric and $h_{\mu\nu}$ is a small
fluctuation of a background metric $g_{\mu\nu}$ given by Eq. (7). For analysis
of growth index, we consider the following physical quantities to be perturbed
viz. the energy density and pressure of an $x$-fluid ($\rho_{x}$, $p_{x}$),
the energy density of a matter ($\rho_{m}$), the four-velocity of an $x$-fluid
and matter ($u_{x}$, $u_{m}$) as

$\displaystyle\bar{\rho}_{x}$ $\displaystyle=$
$\displaystyle\rho_{x}+\delta\rho_{x},~{}~{}~{}\bar{\rho}_{m}=\rho_{m}+\delta\rho_{m},~{}~{}~{}\bar{p}_{x}=p_{x}+\delta
p_{x},~{}~{}~{}$ $\displaystyle\bar{u}_{x}$ $\displaystyle=$ $\displaystyle
u_{x}+\delta u_{x},~{}~{}~{}\bar{u}_{m}=u_{m}+\delta u_{m}.$ (27)

The basic differential equation describing the evolution of the matter
overdensities $\delta_{m}$ defined as $\frac{\delta\rho_{m}}{\rho_{m}}$ is
given by

$\displaystyle\ddot{\delta}_{m}+2\nu H\dot{\delta}_{m}-4\pi
G\mu\rho_{m}\delta_{m}=0\,.$ (28)

In the above equation, the quantities $\nu$ and $\mu$ parametrizes the
deviation of any gravitational theory from GR. In the case of DE models within
the GR framework, we have $\mu=\nu=1$, however, in the case of modified
gravity models, we have $\nu=1$ and $\mu\neq 1$. On using the perturbed metric
(26) in this model, we get

$\displaystyle\nu$ $\displaystyle=$ $\displaystyle
1~{}~{}\text{and}~{}~{}\mu=\tilde{\gamma}+\left(\tilde{\gamma}+3(2-\tilde{\gamma})\,w_{x}\right)\frac{\delta_{x}}{\delta_{m}}\,\frac{\rho_{x}}{\rho_{m}},$
(29)

where $\delta_{x}=\frac{\delta\rho_{x}}{\rho_{x}}$. For a detailed and
complete perturbed equations of this model, one can see [11, 12]. In order to
obtain the quantity $\mu$, we first need to determine the functional form of
$\frac{\delta_{x}}{\delta_{m}}$. For this we shall assume an adiabatic initial
condition of density perturbation given by [44, 45, 46, 47]

$\frac{\delta\rho_{x}}{\dot{\rho}_{x}}=\frac{\delta\rho_{m}}{\dot{\rho}_{m}}.$
(30)

The above adiabatic condition is defined without a prior choice of the
underlying gravitational theory [45], [46]. Further, we note that adiabatic
condition is in agreement with observations and simple inflationary models
usually predict adiabatic perturbations [48]. Using Eqs. (10), (11), the
condition (30) can be written as

$\frac{\delta_{x}}{\delta_{m}}=\frac{\frac{\tilde{\gamma}-1}{2}\frac{\Omega_{m}}{\Omega_{x}}+(1+w_{x})}{1+\frac{1-\tilde{\gamma}}{2}(1-3w_{x})}.$
(31)

Interestingly, for $\tilde{\gamma}=1$, we get the corresponding result of GR
[44] i.e.

$\displaystyle\frac{\delta_{x}}{\delta_{m}}=1+w_{x}\,.$ (32)

On employing (31), we see that the evolution equation of matter perturbations
(28) for $w_{x}=-1$, coincides with that of the $\Lambda$CDM model
irrespective of the choice of $\tilde{\gamma}$. This is in agreement with the
result obtained in [11]. The only difference in this model in comparison to
the $\Lambda$CDM is that an extra relation which relates DE and DM
perturbations arises which is given by [11]

$\delta\rho_{x}=\frac{\tilde{\gamma}-1}{2(3-2\tilde{\gamma})}\delta\rho_{m}.$
(33)

Therefore, DE clusters even though it is in the vacuum energy form. This
suggest that the parameter $\tilde{\gamma}$ might somehow affect the behavior
of the growth index. It worth mentioning here that the condition (31) reduces
to Eq. (33) for $w_{x}=-1$ and which is indeed consistent with a full set of
linear perturbation equations of this model presented in [11]. Note that in
Eq. (33), the case $\tilde{\gamma}=\frac{3}{2}$ is not feasible as discussed
earlier. On employing the favorable choice of parameters values
$\tilde{\gamma}\simeq 1.41$, $w_{x}\simeq-1$ and $\Omega_{m}\simeq 0.296$ [12]
on Eq. (31), we obtain $\frac{\delta_{x}}{\delta_{m}}\simeq 0.08$ which is in
agreement with those predicted by previous studies for clustered DE models
[49]. Now, we shall focus on the investigation of the growth index of matter
perturbations denoted by $\gamma$. In this regard, we consider a growth rate
of clustering $f$ given by [26]

$f(a)=\frac{d\ln\delta_{m}}{d\ln a}\simeq\Omega^{\gamma}_{m}(a)\,,$ (34)

where

$\Omega_{m}(a)=\frac{\Omega_{m0}\,a^{-3}}{E^{2}(a)}.$ (35)

The growth index $\gamma$ is used to distinguish a modified gravity theory
from the standard GR theory on cosmological scales. On differentiating Eq.
(35), one can obtain the equation

$\frac{d\Omega_{m}}{da}=-3\frac{\Omega_{m}(a)}{a}\Big{(}1+\frac{2}{3}\,\frac{d\ln
E}{d\ln a}\Big{)}.$ (36)

Using the first equality of (34), the Eq. (28) can be written as

$a\frac{df}{da}+\left(2\nu+\frac{d\ln E}{d\ln
a}\right)\,f+f^{2}=\frac{3\mu\,\Omega_{m}}{2}.$ (37)

In what follows, we analyze the behavior of growth index $\gamma$ for two
cases separately.

Figure 3: (a) The region in $(w_{x},\tilde{\gamma})$ parameter space where the
asymptotic value of growth index $\gamma_{\infty}$ is greater than that of the
$\Lambda$CDM model. The solid circle corresponds to the $\Lambda$CDM model.
(b) The evolution of $\gamma_{\infty}$ for Rastall’s model as a function of
$\tilde{\gamma}$ (in dashed curve). The solid line corresponds to the
$\Lambda$CDM model. (c) The relative difference
$\Delta(\gamma_{\infty})=[1-\gamma_{\infty}/\gamma_{\infty}^{(\Lambda)}]\%$ of
the asymptotic value of growth index for the Rastall’s model with respect to
the $\Lambda$CDM model. In (b) and (c), we have taken $w_{x}=-1$.

###  4.0.1 Constant growth index

In this subsection, we consider the simplest choice of the growth index i.e.
the asymptotic value of growth index denoted by $\gamma_{\infty}$. An
analytical approach for finding the asymptotic growth index was developed by
Steigerwald et al. [50]. An application of an asymptotic value of growth index
to various class of DE models was discussed in [36]. Within an analytical
method developed by Steigerwald et al., the asymptotic value of growth index
$\gamma_{\infty}$ is given by

$\gamma_{\infty}=\frac{3(M_{0}+M_{1})-2(H_{1}+N_{1})}{2+2X_{1}+3M_{0}}\,,$
(38)

where the above quantities are defined as

$\displaystyle M_{0}$ $\displaystyle=$
$\displaystyle\mu|_{\omega=\Theta},~{}~{}~{}M_{1}=\frac{d\mu}{d\omega}\Big{|}_{\omega=\Theta},~{}~{}~{}N_{1}=\frac{d\nu}{d\omega}\Big{|}_{\omega=\Theta},$
$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle-\frac{X_{1}}{2}=\frac{d(\ln E/d\ln
a)}{d\omega}\Big{|}_{\omega=\Theta}\,,$ (39)

with $\omega=\ln\Omega_{m}(a)$ and
$\Theta=\ln\left(\frac{2}{3-\tilde{\gamma}}\right)$. It is worth mentioning
that the growth index is determined within the matter dominated epoch which
corresponds to $\Omega_{x}=0$. According to relation (22), this corresponds to
$\Omega_{m}=\frac{2}{3-\tilde{\gamma}}$. After some algebraic calculations,
based, on Eqs. (17), (29), (31) we get

$\displaystyle M_{0}$ $\displaystyle=$
$\displaystyle\frac{2\,(3\,\tilde{\gamma}\,w_{x}+\tilde{\gamma}-3\,w_{x})}{3\,\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3\,w_{x}+3}\,,~{}~{}~{}N_{1}=0\,,$
$\displaystyle M_{1}$ $\displaystyle=$
$\displaystyle\frac{-2(1+w_{x})(\tilde{\gamma}-3)(3\tilde{\gamma}w_{x}-\tilde{\gamma}-6w_{x})}{\left(3\,\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3\,w_{x}+3\right)^{2}}\,,$
$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle\frac{-6\,w_{x}\left(2\,\tilde{\gamma}-3\right)}{\left(3-\tilde{\gamma}\right)\left(3\,\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3\,w_{x}+3\right)}\,.$
(40)

Thus, the value of the asymptotic growth index $\gamma_{\infty}$ is given by

$\textstyle\gamma_{\infty}=\frac{3\,w_{x}\left(2\,\tilde{\gamma}-3\right)\left(3\,{\tilde{\gamma}}^{2}w_{x}-{\tilde{\gamma}}^{2}-12\,\tilde{\gamma}\,w_{x}+14\,\tilde{\gamma}-3\,w_{x}-33\right)}{\left(12\,{\tilde{\gamma}}^{2}w_{x}+2\,{\tilde{\gamma}}^{2}-72\,\tilde{\gamma}\,w_{x}-3\,\tilde{\gamma}+72\,w_{x}-9\right)\left(3\,\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3\,w_{x}+3\right)}\,.$
(41)

For $w_{x}=-1$, we get
$\gamma_{\infty}=\frac{3(\tilde{\gamma}-5)}{5\tilde{\gamma}-27}$, which
interestingly depends on $\tilde{\gamma}$ even though the evolution equation
for matter perturbations (i.e. Eq. (28)) coincides with the corresponding
equation for the standard $\Lambda$CDM model. As expected, the above formula
reduces to its standard value for $\Lambda$CDM model i.e
$\gamma^{(\Lambda)}_{\infty}=\frac{6}{11}$ for $\tilde{\gamma}=1$ and
$w_{x}=-1$. In Fig. 3, we plot the region in the $(w_{x},\tilde{\gamma})$
parameter space where the asymptotic value of growth index $\gamma_{\infty}$
is greater than that of the $\Lambda$CDM model [i.e.
$\gamma_{\infty}^{(\Lambda)}=\frac{6}{11}$]. In Fig. 3, we plot the asymptotic
value of growth index as a function of a parameter $\tilde{\gamma}$ for
$w_{x}=-1$. It can be concluded that when $w_{x}=-1$, the value of
$\gamma_{\infty}$ for this model is greater than that of $\Lambda$CDM for
$\tilde{\gamma}<1$ or $\tilde{\gamma}>\frac{27}{5}$. The relative deviation of
$\gamma_{\infty}$ for this model with respect to the $\Lambda$CDM model as a
function of $\tilde{\gamma}$ is given in Fig. 3.

Further, by using the favorable choice of parameters values
$\tilde{\gamma}=1.41$, $w_{x}=-1$ and $\Omega_{m0}=0.296$ [12], we obtain
$\gamma_{\infty}=0.5304$. Therefore, we see that the asymptotic value of
growth index $\gamma_{\infty}$ of this model is somewhat smaller than that of
the $\Lambda$CDM model by approximately $3\%$. This is expected as the
agglomerations of DE perturbations usually lowers the value of the growth
index [51]. It has been also observed that a redshift dependent
parametrization of $\gamma$ can provide a more accurate fit for the growth
rate of matter clustering in comparison to a constant parametrization.
Further, the variation of the growth index may contain crucial information
about the underlying gravitational theory. In what follows, we, therefore,
analyze the evolution of a time-dependent growth index.

Figure 4: The evolution of $\gamma(z)$ as a function of redshift $z$ for
Rastall’s model (red and blue) and the $\Lambda$CDM model (black) with
$\tilde{\gamma}=0.89$ in panel (a) and $\tilde{\gamma}=1.41$ in panel (b). In
panel (c) we plot the relative difference
$\Delta(\gamma_{\infty})=[1-\gamma(z)/\gamma^{(\Lambda)}(z)]\%$ of the varying
growth index for the Rastall’s model with respect to the $\Lambda$CDM model
for $\tilde{\gamma}=1.41$ and $0.89$.

###  4.0.2 Varying growth index

Another possibility is the choice of the growth index $\gamma$ which changes
with redshift. In this case, we utilize the method of Polarski and Gannouji
[52]. On substituting Eq. (34) in Eq. (37) and employing Eq. (36) we get

$a\ln(\Omega_{m})\frac{d\gamma}{da}+\Omega^{\gamma}_{m}-3\Big{(}\gamma-\frac{1}{2}\Big{)}\Big{(}1+\frac{2}{3}\frac{d\ln
E}{d\ln a}\Big{)}+\frac{1}{2}-\frac{3}{2}\,\mu\,\Omega^{1-\gamma}_{m}=0\,.$
(42)

Evaluating the above equation at the present time which corresponds to a
redshift $z=0$ (i.e. $a=1$) we get

$\displaystyle-\gamma^{\prime}(1)\ln\Omega_{m0}+\Omega^{\gamma(1)}_{0}-3\Big{(}\gamma(1)-\frac{1}{2}\Big{)}\,\Big{(}1+\frac{2}{3}\,\frac{d\ln
E}{d\ln a}\Big{)}_{a=1}$
$\displaystyle+\frac{1}{2}-\frac{3}{2}\,\mu_{0}\,\Omega^{1-\gamma(1)}_{m0}=0\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(43)

where

$\displaystyle\mu_{0}$ $\displaystyle=$
$\displaystyle\mu(1)\simeq\tilde{\gamma}+\left(\tilde{\gamma}+3(2-\tilde{\gamma})\,w_{x}\right)$
(44)
$\displaystyle\Bigg{(}\frac{\frac{\tilde{\gamma}-1}{2}\frac{\Omega_{m0}}{\Omega_{x0}}+(1+w_{x})}{1+\frac{1-\tilde{\gamma}}{2}(1-3w_{x})}\Bigg{)}\Bigg{(}\frac{\Omega_{x0}}{\Omega_{m0}}\Bigg{)}\,,$
$\displaystyle\frac{d\ln E}{d\ln a}\Bigg{|}_{a=1}$ $\displaystyle=$
$\displaystyle-\frac{3\,\Omega_{m0}\,w_{x}\,(2\tilde{\gamma}-3)}{3\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3w_{x}+3}\,.$
(45)

In this work, we consider an approximated solution of Eq. (42) expressed as a
first order Taylor’s expansion about the present epoch i.e. $a(z)=1$ given by

$\gamma(a)=\gamma_{0}+\gamma_{1}(1-a).$ (46)

For large redshift ($z\gg 1$) i.e. $a(z)\simeq 0$, we denote the asymptotic
growth index by $\gamma_{\infty}$, so that
$\gamma_{\infty}\simeq\gamma_{0}+\gamma_{1}$ where $\gamma(1)=\gamma_{0}$.
Using Eqs. (4.0.2), (45) and (46), one can expressed $\gamma_{1}$ in terms of
$\gamma_{0}$ as

$\displaystyle\gamma_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\ln\Omega_{m0}}\Bigg{[}\Omega_{m0}^{\gamma_{0}}-3\Big{(}\gamma_{0}-\frac{1}{2}\Big{)}\Big{(}1-\frac{2\,\Omega_{m0}\,w_{x}\,(2\tilde{\gamma}-3)}{3\tilde{\gamma}\,w_{x}-\tilde{\gamma}-3w_{x}+3}\Big{)}$
(47)
$\displaystyle-\frac{3}{2}\mu_{0}\Omega^{1-\gamma_{0}}_{m0}+\frac{1}{2}\Bigg{]}.$

Finally, substituting $\gamma_{0}=\gamma_{\infty}-\gamma_{1}$, in Eq. (47) and
employing the expression of $\gamma_{\infty}$ from Eq. (41), we can express
$\gamma_{0}$ and $\gamma_{1}$ in terms of $\Omega_{m0}$, $\tilde{\gamma}$ and
$w_{x}$. Note that for various class of models within GR, we have
$|\gamma_{1}(z=0)|\lesssim 0.02$. However, this value can be greater in models
outside the GR framework [53]. Here, if we use parameters values
$\tilde{\gamma}=0.89$, $w_{x}=-1$ and $\Omega_{m0}=0.296$, then we obtain
$(\gamma_{0},\gamma_{1})=(0.571,-0.024)$. However, if we use parameters values
$\tilde{\gamma}=1.41$, $w_{x}=-1$ and $\Omega_{m0}=0.296$, then we obtain
$(\gamma_{0},\gamma_{1})=(0.458,-0.054)$. Also, for $\Lambda$CDM model where
$\tilde{\gamma}=1$, $w_{x}=-1$ and $\Omega_{m0}=0.3$, we have
$(\gamma_{0},\gamma_{1})=(0.556,-0.011)$. In Figs. 4 and 4, we present a
comparison of a varying growth index of this model and the $\Lambda$CDM model
for $(\tilde{\gamma},w_{x})=(0.89,-1)$ and $(1.41,-1)$ respectively. It can be
seen that for $\tilde{\gamma}=0.89$, the evolution of growth index is greater
than that of $\Lambda$CDM model for smaller redshift and approaches the
$\Lambda$CDM for large redshift. However, for a favourable choice of model
parameters i.e. $\tilde{\gamma}=1.41$, $w_{x}\simeq-1$, we see that the growth
index is smaller than that of the $\Lambda$CDM model throughout the
evolutionary history of the universe and ultimately approach to that of the
$\Lambda$CDM (see Fig. 4). From Fig. 4 it can be also seen that at the present
epoch, the growth index of this model is smaller than the $\Lambda$CDM by
about $18\%$. Hence, it is interesting to note here that for $w_{x}=-1$, even
though evolution of matter perturbations coincides with that of the
$\Lambda$CDM, however the evolution of the growth index varies from the
$\Lambda$CDM for different values of $\tilde{\gamma}$. Therefore, one need to
impose tight constraints on the parameter $\tilde{\gamma}$, in order to test
possible deviations of this theory from GR.

##  5 Conclusion

In the present work, motivated by an open debate on the comparison between
Rastall gravity (with the vacuum form of DE) and GR, we analyze the behaviour
of the Rastall’s cosmological model at the background and linear perturbation
level. At the background level, this model shows no deviation from the
$\Lambda$CDM model at late time irrespective of the value of Rastall’s
parameter $\tilde{\gamma}$. However, at the level of linear perturbation, we
observed that even though the evolution equation of the matter overdensities
for this model coincides with that of the $\Lambda$CDM, the constant value of
growth index of matter perturbations shows significant deviation from the
$\Lambda$CDM. Furthermore, we obtained that the varying growth index is
smaller than the corresponding case of the $\Lambda$CDM for
$\tilde{\gamma}>1$. In particular, for best-fit parameters values, we obtained
that the asymptotic value of the growth index is about $3\%$ smaller than the
$\Lambda$CDM, but the varying growth index is about $18\%$ smaller than the
$\Lambda$CDM at the present epoch. Therefore, the growth index value is
outside the range of GR based DE models. Although we do not perform an
exhaustive study of the Rastall model, the results of the present work, indeed
shed some light on the non-equivalence of Rastall gravity and GR. Further, the
smaller value of the growth index for this model is also expected due to the
presence of DE perturbations. The influence of DE perturbations on the
behaviour of the growth index is crucial, particularly in the light of future
generation surveys to test the existence of DE perturbations in the observed
Universe [54]. Moreover, the agglomeration of DE perturbations may influence
the early structure formation at the level of galaxies and cluster of galaxies
[11]. The results obtained here also signify the possible deviation of the
present model from the $\Lambda$CDM at the nonlinear regime of cosmic
perturbations as claimed in Ref. [11]. A deeper analysis along this line is a
subject of future work to further test the viability of Rastall theory in
comparison to GR with upcoming precise observational data.

## Acknowledgments

JD was supported by the Core research grant of SERB, Department of Science and
Technology India (File No.CRG/2018/001035) and the Associate program of IUCAA.
The authors thank the referee for constructive suggestions which lead to the
improvement of the work.

## References

  * [1] T. Clifton, P. G Ferreira, A. Padilla, C. Skordis, Phys. Rept. 513 (2012) 1-189.
  * [2] C. M. Will, Living. Rev. Relativ. 17 (2014) 4.
  * [3] S. Nojiri, S. D. Odintsov, V. K. Oikonomou, Phys. Rept. 692 (2017) 1.
  * [4] P. Rastall, Phys. Rev. D 6 (1972) 3357.
  * [5] P. Rastall, Can. J. Phys. 54 (1976) 66.
  * [6] L. Parker,  Phys. Rev. Lett. 21 (1968) 562.
  * [7] L. Parker, Phys. Rev. D 183 (1969) 1057.
  * [8] L. Parker, Phys. Rev. D 3 (1971) 346.
  * [9] G. W. Gibbons, S. W. Hawking, Phys. Rev. D 15 (1977) 2738.
  * [10] C. E. M. Batista, J. C. Fabris, M. H. Daouda, Nuovo Cim. B 125 (2010) 957.
  * [11] C. E. M. Batista, M. H. Daouda, J. C. Fabris, O. F. Piattella, D. C. Rodrigues, Phys. Rev. D 85 (2012) 084008.
  * [12] C. E. M. Batista, J. C. Fabris, O. F. Piattella, A. M. Velasquez-Toribio, Eur. Phys. J. C 73 (2013) 2425.
  * [13] K. A. Bronnikov, J. C. Fabris, O. F. Piattella, D. C. Rodrigues, Eur. Phys. J. C 77 (2017) 409.
  * [14] Y. Heydarzade, F. Darabi, Phys. Lett. B 771 (2017) 365.
  * [15] A. S. Al-Rawaf, O. M. Taha, Phys. Lett. B 366 (1996) 69.
  * [16] A. S. Al-Rawaf, Int. J. Mod. Phys. D 14 (2005) 1941.
  * [17] J.C. Fabris, R. Kerner, J. Tossa, Int. J. Mod. Phys. D 9 (2000) 111.
  * [18] V. Majernik, L. Richterek, [gr-qc/0610070] (2006).
  * [19] M. Visser, Phys. Lett. B 782 (2018) 83.
  * [20] D. Das, S. Dutta, S. Chakraborty, Eur. Phys. J. C 78 (2018) 810.
  * [21] F. Darabi, H. Moradpour, I. Licata, Y. Heydarzade, C. Corda, Eur. Phys. J. C 78 (2018) 25.
  * [22] S. Hansraj, A. Banerjee, P. Channuie, Ann. of Phys. 400 (2019) 320-345.
  * [23] L. L. Smalley, J. Phys. A Math. Gen. 16 (1983) 2179.
  * [24] T. Okumura, et al., Publ. Astron. Soc. Jpn 68 (2016) 38.
  * [25] M. Ishak, Living Rev. Rel. 22 (2019) 1.
  * [26] P. J. E. Peebles, Principles of Physical Cosmology, Princeton University, Princeton, NJ (1993).
  * [27] V. Silveira, I. Waga, Phys. Rev. D 50 (1994) 4890.
  * [28] L. M. Wang, P. J. Steinhardt, Astrophys. J. 508 (1998) 483.
  * [29] S. Nesseris, L. Perivolaropoulos, Phys. Rev. D 77 (2008) 023504.
  * [30] A. Lue, R. Scoccimarro, G. D. Starkman, Phys. Rev. D 69 (2004) 124015.
  * [31] Y. Gong, Phys. Rev. D 78 (2008) 123010.
  * [32] R. Gannouji, B. Moraes, D. Polarski, J. Cosmol. Astropart. Phys. 02 (2009) 034.
  * [33] S. Tsujikawa, R. Gannouji, B. Moraes, D. Polarski, Phys. Rev. D 80 (2009) 084044.
  * [34] L. Boubekeur, E. Giusarma, O. Mena, H. Ramrez, Phys. Rev. D 90 (2014) 103512.
  * [35] S. Basilakos, P. Stavrinos, Phys. Rev. D 87 (2013) 043506.
  * [36] S. Basilakos, J. Sola, Phys. Rev. D 92 (2015) 123501.
  * [37] A. Mehrabi, S. Basilakos, F. Pace, Mon. Not. R. Astron. Soc. 452 (2015) 2930.
  * [38] A. Mehrabi, S. Basilakos, M. Malekjani, Z. Davari, Phys. Rev. D 92 (2015) 123513.
  * [39] S. Basilakos, Phys. Rev. D 93 (2016) 083007.
  * [40] M. Quartin, M. O. Calvao, S. E. Joras, R. R. Reis, I. Waga, J. Cosmol. Astropart. Phys. 0805 (2008) 007.
  * [41] J. D. Barrow, T. Clifton, Phys. Rev. D 73 (2006) 103520.
  * [42] H. Zonunmawia, W. Khyllep, N. Roy, J. Dutta, N. Tamanini, Phys. Rev. D 96 (2017) 083527.
  * [43] J. Dutta, W. Khyllep, E. N. Saridakis, N. Tamanini, S. Vagnozzi, J. Cosmol. Astropart. Phys. 1802 (2018) 041.
  * [44] G. Ballesteros, A. Riotto, Phys. Lett. B 668 (2008) 171.
  * [45] A. R. Liddle, D. H. Lyth, Cosmological inflation and large-scale structure, Univ. Pr., 2002. Cambridge, UK; 400 p.
  * [46] D. Baumann, Cosmology (2018),
http://cosmology.amsterdam/education/cosmology/

  * [47] L. Amendola and S. Tsujikawa, Dark Energy Theory and Observations, Cambridge University Press, Cambridge, England (2010).
  * [48] C. B. Netterfield et al. [Boomerang Collaboration], Astrophys. J. 571 (2002) 604.
  * [49] S. Basilakos, Mon. Not. Roy. Astron. Soc. 449 (2015) 2151.
  * [50] H. Steigerwald, J. Bel, C. Marinoni, J. Cosmol. Astropart. Phys. 05 (2014) 042.
  * [51] R. C. Batista, Phys. Rev. D 89 (2014) 123508.
  * [52] D. Polarski, R. Gannouji, Phys. Lett. B 660 (2008) 439.
  * [53] R. Gannouji, D. Polarski, J. Cosmol. Astropart. Phys. 0805 (2008) 018.
  * [54] D. Sapone, E. Majerotto, M. Kunz, B. Garilli, Phys. Rev. D 88 (2013) 043503.