# Thermal evolution and stability analysis of phenomenologically emergent dark
energy model

Rosemin John<EMAIL_ADDRESS>Department of Physics, Cochin University
of Science and Technology, Kochi, Kerala 682022, India  Sarath N
<EMAIL_ADDRESS>Indian Institute of Technology, Kanpur 208016, India  Titus
K. Mathew<EMAIL_ADDRESS>Department of Physics, Cochin University of
Science and Technology, Kochi, Kerala 682022, India

# Thermal evolution and stability analysis of phenomenologically emergent dark
energy model

Rosemin John<EMAIL_ADDRESS>Department of Physics, Cochin University
of Science and Technology, Kochi, Kerala 682022, India  Sarath N
<EMAIL_ADDRESS>Indian Institute of Technology, Kanpur 208016, India  Titus
K. Mathew<EMAIL_ADDRESS>Department of Physics, Cochin University of
Science and Technology, Kochi, Kerala 682022, India

###### Abstract

The phenomenologically emergent dark energy (PEDE) model is a varying dark
energy model with no extra degrees of freedom proposed by Li and Shafieloo[1]
to alleviate the Hubble tension. The statistical consistency of the model has
been discussed by many authors. Since the model depicts a phantom dark energy
that increases with redshift, its cosmic evolution, particularly during the
late phase, must be examined. We discover that the model’s Hubble and
deceleration parameters display unusual behaviour in the future, which differs
from $\Lambda$CDM cosmology. We find the model also follows a distinct
evolution in the statefinder plane. The phantom nature of the model leads to
the violation of the null energy condition and a decrease in horizon entropy.
The asymptotic future epoch also seems to be unstable based on our dynamical
system analysis as well as the stability analysis based on dark energy sound
speed.

PEDE model – NEC – GSL – entropy – stability

##  I Introduction

The late acceleration of the universe’s expansion was one of the crucial
cosmological discoveries [2, 3]. This phenomenon couldn’t be explained if the
universe is composed only of matter and radiation, which prompted scientists
to look for a new exotic cosmic component with a negative equation of state,
$\omega$ ($\omega=p/\rho$, $p$ is the pressure and $\rho$ is the density),
that they termed ‘dark energy’. The current standard model, the $\Lambda$CDM
model, incorporates cosmological constant as dark energy, although its precise
nature has not yet been established. Despite being able to explain the late
acceleration in the universe’s expansion, the $\Lambda$CDM model has some
inconsistencies [4, 5]. The cosmological constant problem [6], also known as
the fine tuning problem, stems from the enormous gap between the dark energy
density that was determined from quantum field theory and that obtained from
cosmic observations. The difference is around $10^{121}$ orders of magnitude
[7]. The coincidence problem [8, 9], which concerns the identical densities of
matter and dark energy in the current universe although these components
evolve differently such that the matter density decreases in inverse
proportion to the volume of the universe while the density of the cosmological
constant remains constant throughout evolution, is another significant issue
with $\Lambda$CDM cosmology. There also exists tensions in the parameter
values extracted from cosmic data, among which the recently discovered Hubble
tension has gained the most attention. ‘Hubble tension’ mainly refers to the
difference in the present value of the Hubble parameter extracted from cosmic
microwave background (CMB) data and that from local measurements. The Planck
collaboration [10], from the CMB observations assuming $\Lambda$CDM, provided
a value of $H_{0}=67.4\pm 0.5\
$\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$. The
SH0ES (Supernova H0 for the Equation of State) group of Hubble space
telescope, from the observations of 70 long-period Cepheids in the Large
Magellanic Cloud with other local distance anchors predicted the value of
$H_{0}=74.03\pm 1.42\
$\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$
[11]. The inconsistency between these two measurements is over 4$\sigma$. The
tension has now increased over time as a result of examining multiple data
sets gathered using various techniques [12, 13, 4]. If there are no systematic
errors and since the CMB measurements are model dependent, it is anticipated
that the existing standard cosmological model will be required modifications
to account for this tension.

There have been several attempts to alleviate these problems of the
$\Lambda$CDM model. The easiest and most natural way is to assume a different
form for the dark energy. Such dark energy models include: scalar field
k-essence models [14], quintessence dark energy models with scalar field
having $\omega>-1$ [15, 16], phantom models where dark energy is replaced with
a phantom component which has $\omega<-1$ [17, 18], quintom model of dark
energy with two scalar fields with one being quintessence and the other
phantom [19, 20], interacting dark energy models [21, 22, 23], Chaplygin gas
model [24, 25, 26], holographic dark energy model [27], etc. Modified gravity
theories[28, 29, 30, 31, 32, 33, 34, 35, 36] can be another significant
approach to address the challenges faced by the standard model.

Recently, Li and Shafieloo [1] proposed a new dark energy model called the
Phenomenologically Emergent Dark Energy (PEDE) model, in which dark energy has
a varying density and a phantom nature. It was primarily developed to address
the Hubble tension. The model effectively replaces the cosmological constant
with a hyperbolic tangent function of redshift which causes the dark energy to
emerge as a function of the cosmic time at later times. The absence of
additional parameters makes the model more intriguing. PEDE model using the
Planck 2018 CMB data + Pantheon + BAO + Ly-$\alpha$ data finds
$H_{0}=71.02_{-1.37}^{+1.45}\
$\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$ [1]
thereby resolving the tension on $H_{0}$ within 68% confidence-level. The
higher Hubble constant value is a direct consequence of dominance of the
phantom energy during the late phase. Later, a generalized PEDE model with an
extra degree of freedom was proposed [37], and under the right limiting
conditions, it can be reduced to both the original PEDE model and the
$\Lambda$CDM model. Benaoum et al.[38] have also developed a modified version
of the PEDE model known as Modified Emergent Dark Energy (MEDE), which also
claims to incorporate the $\Lambda$CDM, PEDE model, Chevallier-Polarski-Linder
(CPL) model and other related cosmological models.

The PEDE model was further analyzed statistically in comparison with the
cosmological data by various authors. Yang et al.[39] claim that the model can
still reduce cosmic tensions despite the inclusion of massive neutrino species
and possible additional relativistic degrees of freedom in their revised
analysis of the PEDE scenario. The Hubble constant tension is found to be
below the 2$\sigma$ level in this scenario, and the value of $\sigma_{8}$ is
in better agreement with cosmic shear estimates. They have obtained
$M_{\nu}\sim 0.21_{-0.14}^{+0.15}$ eV and $N_{eff}=3.03\pm 0.32$, which point
to a non-zero neutrino mass with a significance greater than 2$\sigma$.
However, they obtained the sound horizon as $\sim$ 147 Mpc, failing to restore
the lower value that the BAO data prefers. In Ref.[40], the authors have
demonstrated how the PEDE model can solve the coincidence problem along with
Hubble tension. But the $\sigma_{8}$ tension between low-redshift observations
and Planck inferred value, however, cannot be resolved by the PEDE model, and
the analysis with the growth rate dataset leads them to the conclusion that
the PEDE model is less capable of explaining observable evidence in cluster
scales at the perturbation level when compared to the $\Lambda$CDM. According
to the Bayesian analysis performed on the model by Pan et al.[41], the
majority of observational dataset combinations favor $\Lambda$CDM over the
PEDE model. Alestas et al.[42] have also analyzed the model’s viability
through the Akaike information criterion (AIC) and the Bayes factor
considering the full cosmic microwave background temperature anisotropy
spectrum data with other cosmological datasets and they found the model is
strongly disfavored in both criteria. Schöneberg et al.[43] have shown PEDE
model offers a worse fit to the combined datasets (Planck 2018 (including
TTTEEE and lensing) + BAO (including BOSS DR12 + MGS +6dFGS)+Pantheon) than
$\Lambda$CDM. This is due to a significant degradation of the joint fit to BAO
and Pantheon data and is a result of the non-nested nature of the PEDE model.
The model also decreases the angular diameter distance for the CMB, because of
the equation of state $\omega<-1$ that this model exhibits at later times. The
authors also point out that the model has the ability to artificially lower
the tension while severely impairing the fit to the CMB anisotropy data.
Sharma[44] studied the PEDE model and its anisotropic extension
($\mathcal{A}$PEDE) via the Bianchi type-I spacetime metric. When compared to
the base $\Lambda$CDM, the statistical results for both the PEDE and
$\mathcal{A}$PEDE models indicate a good fit with the Hubble+CMB data
combination, however, the addition of BAO or Pantheon with other data
combinations lead to very strong evidence against the PEDE models.

Even while the PEDE model’s cosmological evolution has been extensively
studied, its thermodynamic evolution has not been considered so far.
Therefore, the primary purpose of this research is to investigate the thermal
evolution of the PEDE model. We further perform a stability analysis of this
model intending to compare its feasibility in comparison with the standard
$\Lambda$CDM. In the following part, we briefly present the model before
extracting the parameter values. We begin our analysis by obtaining an
analytical solution for the Hubble and deceleration parameters of the PEDE
model and study their evolution in detail. We then use the statefinder
diagnostic to distinguish the PEDE model from other cosmological models. We
then continue to analyze the status of the laws of thermodynamics. We examine
the Generalized Second Law(GSL) of thermodynamics and the convexity condition
to check whether the model implies maximization of the horizon entropy.
Finally, employing dynamical system analysis and sound speed analysis, we
investigate the stability of the model.

##  II Phenomenologically Emergent Dark Energy

Within the framework of Einstein’s general relativity, a spatially flat,
homogeneous and isotropic universe following the Friedmann-Lemaître-Robertson-
Walker (FLRW) metric can be explained by the Friedmann equation,

$H^{2}=\frac{8\pi G}{3}(\rho_{\gamma}+\rho_{m}+\rho_{\scaleto{DE}{3.5pt}})$
(1)

where $H$ is the Hubble parameter and $G$ is the gravitational constant.
$\rho_{\gamma},\rho_{m}$ and $\rho_{\scaleto{DE}{3.5pt}}$ represent the
density of radiation, matter and dark energy respectively. The density
components can be expressed in terms of their corresponding dimensionless
density parameters, $\Omega_{i}=\rho_{i}/\rho_{c,0}$ where
$\rho_{c,0}=3H_{0}^{2}/8\pi G$ is the critical density. Eq.1 can then be
rewritten as,

$H^{2}(z)=H_{0}^{2}\left[\Omega_{\gamma
0}(1+z)^{4}+\Omega_{m0}(1+z)^{3}+\Omega_{\scaleto{DE}{3.5pt}}(z)\right].$ (2)

Naught in the subscript denotes the present value of the parameters. In the
standard model, the dark energy density parameter has a constant value,
$\Omega_{\scaleto{DE}{3.5pt},0}$. But, in general the evolution of
$\Omega_{\scaleto{DE}{3.5pt}}(z)$ with redshift can be expressed as,

$\Omega_{\scaleto{DE}{3.5pt}}(z)=\Omega_{\scaleto{DE}{3.5pt},0}\times\exp\Bigg{[}3\int_{0}^{z}\frac{1+\omega_{\scaleto{DE}{3.5pt}}(z^{\prime})}{1+z^{\prime}}dz^{\prime}\Bigg{]}.$
(3)

where $\omega_{\scaleto{DE}{3.5pt}}$ is the equation of state of dark energy.
Therefore, one can construct a dark energy model by assuming a specific form
for either equation of state or density of dark energy. In PEDE model, Li and
Shafieloo[1] have given the following form for dark energy density:

$\Omega_{\scaleto{DE}{3.5pt}}(z)=\Omega_{\scaleto{DE}{3.5pt},0}\times\Big{(}1-\tanh\big{[}\log_{10}(1+z)\big{]}\Big{)}.$
(4)

As the authors claim, in this model, dark energy has no effective presence in
the past as its density will be negligibly small, but emerges at later times
as the universe evolves. Therefore, in this study, we are mainly interested in
the PEDE model’s late phase evolution, during which dark energy becomes the
dominant component. (We disregard the contribution of radiation from here
onward because $\rho_{\gamma}$ will be negligibly small in late phase.)

If matter and dark energy have no interaction, then the conservation equations
are as follows:

$\begin{split}\dot{\rho}_{m}+3H\rho_{m}(1+\omega_{m})=0\\\
\dot{\rho}_{\scaleto{DE}{3.5pt}}+3H\rho_{\scaleto{DE}{3.5pt}}(1+\omega_{\scaleto{DE}{3.5pt}})=0\end{split}$
(5)

where $\omega_{m}=0$ and $\omega_{\scaleto{DE}{3.5pt}}$ can be derived as,

$\omega_{\scaleto{DE}{3.5pt}}(z)=\frac{1}{3}\frac{d\ln\Omega_{\scaleto{DE}{3.5pt}}}{dz}(1+z)-1.$
(6)

Substituting the PE dark energy and the equation becomes,

$\omega_{\scaleto{DE}{3.5pt}}(z)=-\frac{1}{3\ln
10}\times\Big{(}1+\tanh\big{[}\log_{10}(1+z)\big{]}\Big{)}-1$ (7)

Both $\Omega_{\scaleto{DE}{3.5pt}}$ and $\omega_{\scaleto{DE}{3.5pt}}$ follows
a symmetrical behavior in logarithmic scale(Fig.1 in [1]). For an early epoch,
i.e., when $z$ is large and positive, $\omega_{\scaleto{DE}{3pt}}$ converges
to a value of $-\frac{2}{3\ln 10}-1$. As time progresses the equation of state
increases, becomes $-\frac{1}{3\ln 10}-1$ in the present epoch, and
asymptotically reaches $-1$ in the far future. Here the interesting fact to be
noted is that the PEDE is having phantom nature with
$\omega_{\scaleto{DE}{3pt}}<-1$[45] except at the asymptotic limit, $z\to-1.$
This is in clear contrast with the standard model, for which the equation of
state of dark energy stays pegged at -1.

##  III Parameter Estimation

We have estimated and compared the standard $\Lambda$CDM and PEDE model
parameters using the cosmological observation data on : $H(z)$ data [46],
Pantheon sample consisting of type Ia supernovae data [47], cosmic microwave
background (CMB) data [48] and baryon acoustic oscillation (BAO) data [49].
The parameters are obtained by applying Markov Chain Monte Carlo (MCMC) method
using the emcee package as optimizer [50] in the lmfit python library [51]. We
have adopted uniform priors for the parameters as, $H_{0}$ in the range, 65 to
75 $\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$
and $\Omega_{m0}$ in the range, 0.1 to 0.5 . The best fit values are then
calculated by minimizing the $\chi^{2}$ function defined as [52, 53],

$\chi^{2}_{i}=\sum\Big{[}\frac{V_{i}^{obs}-V_{i}^{th}}{\sigma_{i}}\Big{]}^{2}$
(8)

where $V_{i}^{obs}$ is the observed value from the datasets, $V_{i}^{th}$ is
the corresponding theoretical value obtained from our model and $\sigma_{i}$
is the error in the measured values. The minimum $\chi^{2}$ per degrees of
freedom or the reduced $\chi^{2}$ is defined as,
$\chi^{2}_{d.o.f}=\frac{\chi^{2}}{N-N_{p}}$, where $N$ is the number of data
points and $N_{p}$ is the number of variable parameters. We have also included
the Akaike Information Criterion (AIC) and Bayesian Information Criterion
(BIC), which are defined as,

$\begin{split}AIC&=N\ln(\chi^{2}/N)+2N_{p}\\\
BIC&=N\ln(\chi^{2}/N)+\ln(N)N_{p}\end{split}$ (9)

Similar to $\chi^{2}$, models with smaller AIC and BIC values are preferred.
The observations considered for the estimation are detailed below.

  * •

$H(z)$ data: The Hubble data set used in this analysis consists of 52 $H(z)$
measurements in the range $0\leq z\leq 2.36$ out of which 31 data points are
from the differential age (DA) and 20 data points from clustering
measurements[46]. The DA method, introduced by Jimenez and Loeb[54], measures
the Hubble parameter by comparing the ages of passively-evolving galaxies with
similar metallicity and separated by a small redshift interval. The latter is
the clustering of galaxies or quasars which was first proposed in Ref.[55],
using the BAO peak position as a standard ruler in the radial direction. The
remaining data point is the direct measurement result of Hubble constant,
$H_{0}=73.24\pm 1.74$
$\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$[56]
from the HST.

  * •

SN Ia: The type Ia supernovae dataset used for the analysis is the latest
Pantheon sample [47] which comprises the luminosity of 1048 SNe Ia within the
redshift $0.01<z<2.3$. To calculate the observed distance modulus, Scolnic et
al.[47] used the SALT2[57] light curve fitter,

$\mu_{i}^{obs}=m(z)+\alpha X_{1}-\beta C-M$ (10)

where $m$ is the rest frame B-band peak magnitude and $M$ is the absolute
B-band magnitude of a fiducial SN Ia with $X_{1}$ (time stretch of light
curve) and $C$ (supernova color at maximum brightness) set to zero. The
stretch-luminosity parameter ($\alpha$) and the color-luminosity parameter
($\beta$) are calibrated to zero for the pantheon sample, hence the observed
distance modulus reduces to $\mu_{i}^{obs}=m-M$. For our analysis, we treat
$M$ as a nuisance parameter and is allowed to vary in the range, $-20\leq
M\leq-18$.

The luminosity distance of these SN Ia can be calculated from the equation,

$d_{L}(H_{0},\Omega_{m0},z_{i})=c(1+z_{i})\int_{0}^{z_{i}}\frac{dz^{\prime}}{H(H_{0},\Omega_{m0},z^{\prime})}$
(11)

where $z_{i}$ is the redshift of the SN Ia, $c$ is the speed of light and
$H(H_{0},\Omega_{m0},z^{\prime})$ is the Hubble parameter specified for the
model. The theoretical distance modulus of SN Ia is then given by,

$\mu_{i}^{th}(H_{0},\Omega_{m0},z_{i})=5\log_{10}\Big{[}\frac{d_{L}(H_{0},\Omega_{m0},z_{i})}{Mpc}\Big{]}+25.$
(12)

  * •

CMB: The theoretical shift parameter $R$ of the cosmic microwave background is
defined as,

$R(H_{0},\Omega_{m0})=\sqrt{\Omega_{m0}}\int_{0}^{z_{*}}\frac{dz}{h(H_{0},\Omega_{m0},z)}$
(13)

where $z_{*}$ is the redshift at the photon decoupling epoch and $h=H/H_{0}$.
The Planck 2018 [48] observations suggests the distance prior measurement
value $R^{obs}=1.7502\pm 0.0046$ at the redshift $z_{*}=1089.2$. By
considering this Planck measurement corresponding $\chi^{2}$ function can be
defined.

  * •

BAO: The baryon acoustic peak parameter $A$ best describes the distance-
redshift relationship produced by a BAO measurement. The theoretical acoustic
parameter can be defined in terms of the model parameters as [52],

$A(H_{0},\Omega_{m0})=\frac{\sqrt{\Omega_{m0}}}{h(z_{A})^{1/3}}\Big{(}\frac{1}{z_{A}}\int_{0}^{z_{A}}\frac{dz}{h(H_{0},\Omega_{m0},z)}\Big{)}^{2/3}$
(14)

where $z_{A}$ is the redshift corresponding to the acoustic peak. To obtain
the $\chi_{BAO}^{2}$, we compare the theoretical value with $A_{obs}=0.484\pm
0.016$ at the redshift $z_{A}=0.35$ obtained from SDSS-BAO distance data [49].

Table 1: The best estimated values and 68.3% confidence limit for the $\Lambda$CDM model parameters. Datasets | $H_{0}$ | $\Omega_{m0}$ | $M$ | $\chi^{2}$ | $\chi^{2}_{d.o.f}$ | AIC | BIC
---|---|---|---|---|---|---|---
$H(z)$ | $71.592_{-0.991}^{+0.987}$ | $0.248_{-0.014}^{+0.015}$ | … | 29.746 | 0.5949 | -25.045 | -21.142
Pantheon | $69.991_{-3.401}^{+3.416}$ | $0.285_{-0.012}^{+0.013}$ | $-19.357_{-0.108}^{+0.103}$ | 1035.680 | 0.9911 | -6.393 | 8.471
$H(z)$ \+ Pantheon + CMB | $69.276_{-0.572}^{+0.576}$ | $0.289_{-0.006}^{+0.006}$ | $-19.377_{-0.017}^{+0.017}$ | 1074.862 | 0.9789 | -20.453 | -5.441
$H(z)$ \+ Pantheon + CMB + BAO | $69.257_{-0.567}^{+0.570}$ | $0.290_{-0.006}^{+0.006}$ | $-19.378_{-0.017}^{+0.017}$ | 1074.929 | 0.9781 | -21.409 | -6.394
Table 2: The best estimated values and 68.3% confidence limit for the PEDE model parameters. Datasets | $H_{0}$ | $\Omega_{m0}$ | $M$ | $\chi^{2}$ | $\chi^{2}_{d.o.f}$ | AIC | BIC
---|---|---|---|---|---|---|---
$H(z)$ | $73.578_{-1.024}^{+0.855}$ | $0.255_{-0.012}^{+0.014}$ | … | 33.697 | 0.6739 | -18.559 | -14.657
Pantheon | $69.988_{-3.389}^{+3.402}$ | $0.332_{-0.012}^{+0.012}$ | $-19.365_{-0.108}^{+0.103}$ | 1032.984 | 0.9885 | -9.125 | 5.739
$H(z)$ \+ Pantheon + CMB | $71.962_{-0.625}^{+0.625}$ | $0.281_{-0.006}^{+0.006}$ | $-19.330_{-0.018}^{+0.018}$ | 1097.631 | 0.9997 | 2.626 | 17.638
$H(z)$ \+ Pantheon + CMB + BAO | $71.966_{-0.618}^{+0.618}$ | $0.280_{-0.006}^{+0.006}$ | $-19.330_{-0.018}^{+0.017}$ | 1097.642 | 0.9988 | 1.633 | 16.648
Figure 1: The two dimensional posterior contours with one dimensional
posterior distribution (using pygtc open python package [58]) from the
$H(z)$+Pantheon+CMB+BAO datasets for the $\Lambda$CDM and PEDE model
parameters Figure 2: The evolution of Hubble parameter with redshift in
$\Lambda$CDM and PEDE models compared with H(z) data points Figure 3: The
evolution of apparent magnitude of supernovae with redshift in $\Lambda$CDM
and PEDE models compared with Pantheon data points

The best estimated parameters for the $\Lambda$CDM and PEDE model within 68.3%
confidence level are given in Tab.1 and Tab.2 respectively and Fig.1
represents the two dimensional posterior contours and one dimensional
posterior distribution for the best fit parameter values obtained for
$H(z)$+SNIa+CMB+BAO dataset. The minimum $\chi^{2}$ value as well as
$\chi^{2}_{d.o.f}$ are also mentioned in the tables. The $\chi^{2}$ for the
combination of datasets is obtained by adding the individual $\chi^{2}$
minimum values for each dataset. The AIC and BIC values are also displayed in
the table. In Fig.2 the evolution of the Hubble parameter for the PEDE and
$\Lambda$CDM models using the best predicted model parameter values for the
combined dataset are displayed and contrasted with the $H(z)$ data points. The
evolution of apparent magnitude $m$ of supernovae for both the models are also
compared with the Pantheon data points in Fig.3.

It is evident that the PEDE model can give a higher best fit value for the
Hubble constant than $\Lambda$CDM except when the Pantheon sample is
considered alone. For the Pantheon dataset the $\chi^{2}$ values for both
models are similar. The addition of Supernovae data pushes the Hubble constant
in PEDE model to relatively lower values. However, adding the $H(z)$ and CMB
seems to increase the $H_{0}$ value and lower the $\Omega_{m0}$ value in the
case of PEDE model whereas, for $\Lambda$CDM, the reverse is true. For the
combined dataset of $H(z)$\+ Pantheon +CMB, the estimated Hubble constant is
$69.257_{-0.567}^{+0.570}$ for $\Lambda$CDM and $71.966_{-0.618}^{+0.618}$ for
PEDE with a $\chi^{2}$ difference of -22.713. The higher $H_{0}$ value from
the PEDE model can be regarded as an indicator of its ability to resolve the
tension in the Hubble constant. Even though the dataset combination provides a
higher $H_{0}$ value for the PEDE model, it happens at the cost of an increase
in $\chi^{2}$, AIC and BIC values. Similar problem of alleviation of tension
by worsening the $\chi^{2}$ values compared to the $\Lambda$CDM model has been
pointed out in Ref.[41]. The $\chi^{2}$ value listed in both tables points out
the fact that $\Lambda$CDM is still favored over PEDE model when the datasets
are considered together.

##  IV Evolution of Cosmological Parameters

###  IV.1 Hubble parameter

The late phase evolution of the Hubble parameter in PEDE model for a flat
universe can be expressed as,

$H^{2}(z)=H_{0}^{2}\Big{[}\Omega_{m0}(1+z)^{3}+\Omega_{\scaleto{DE}{3.5pt},0}\big{(}1-\tanh\big{[}\log_{10}(1+z)\big{]}\big{)}\Big{]}.$
(15)

where $\Omega_{\scaleto{DE}{3.5pt},0}=1-\Omega_{m0}$. The evolution of the
normalized Hubble parameter, $H(z)/H_{0}$ is plotted in Fig.4, for both the
standard as well as PEDE model with parameter values obtained for the dataset
combination $H(z)$+Pantheon+CMB+BAO in the earlier section.

Figure 4: The evolution of $H/H_{0}$ with $z$

It is evident from the figure that the evolution of the Hubble parameter in
the present model has a substantial difference from the canonical model. The
Hubble parameter in $\Lambda$CDM is a smooth decreasing function of redshift
that reaches a constant value in the future epoch. But in the present PEDE
model, the Hubble parameter decreases until around $z\approx-0.3,$ which
corresponds to a future time and then begins to increase. The dominating
phantom component, which causes a faster acceleration in the universe’s
expansion, is the clear cause of the increase in the Hubble parameter. When
$z\to-1$ the Hubble parameter in PEDE model asymptotically approaches
$1.2H_{0}$. However, this non-monotonic behavior can eventually cause the
violation of the fundamental thermodynamic principle, as we will show in a
later section.

###  IV.2 Deceleration Parameter

As the two quantities, the Hubble parameter and the deceleration parameter,
are the focus of cosmology, we may now turn our attention to the evolution of
deceleration parameter. For the PEDE model, we have observed a peculiar
tendency in the evolution of the Hubble parameter, which also has an impact on
the evolution of the deceleration parameter. The standard definition of the
deceleration parameter is provided by,

$q=-1-\frac{\dot{H}}{H^{2}}.$ (16)

The over dot represents derivative over cosmic time. This equation can
suitably be expressed in terms of a convenient variable, $x=\ln a,$

$q=-1-\frac{1}{2h^{2}}\frac{dh^{2}}{dx}.$ (17)

where $h=H/H_{0}.$ The evolution of deceleration parameter for both the
canonical and PEDE models are plotted in Fig.5.

Figure 5: The evolution of $q$ with $z$

As anticipated, the evolution of the deceleration parameter in PEDE model
shows a marked deviation from the standard behavior. According to the standard
model, the deceleration parameter decreases from its positive values and
reaches -1 asymptotically. For $\Lambda$CDM, the transition redshift and the
value of deceleration parameter for the present epoch (z=0) are found to be
$z_{*}=0.70$ and $q_{0}=-0.56$ respectively. As shown in Fig.5, the
deceleration parameter in PEDE model as well declines from the positive region
and crosses into the negative region corresponding to the accelerating epoch
at a redshift $z_{*}=0.71$. The current value of q is found to be,
$q_{0}=-0.74$. However, below a redshift of $z\approx-0.3$, the $q$ parameter
in PEDE model evolves to values less than -1 and attains a minimum, then
increases and reaches the value -1 asymptotically. This behavior is in line
with that of both the Hubble parameter and dark energy density and is very
unconventional as compared to the standard $\Lambda$CDM cosmology.

##  V Energy conditions

Energy conditions are useful to analyze the physical feasibility as well as
predicting the occurrence of singularities concerning a cosmological model. In
principle, energy conditions are expressing the fact that energy density
measured in a region can’t be negative. The general energy conditions in
cosmology are obtained by limiting the energy momentum tensor, $T_{\mu\nu}$ on
physical grounds [59, 60, 61, 62].

  * •

Strong energy condition (SEC)
The SEC states that $(T_{\mu\nu}-\frac{1}{2}Tg_{\mu\nu})u^{\mu}u^{\nu}\geq 0$
for all time like vectors $u^{\mu}$. SEC is satisfied by gravitationally
interacting cosmic components. For a perfect fluid, it is equal as,

$\rho+3p\geq 0$ (18)

For matter and radiation, the SEC is satisfied in every phase of evolution of
universe. But dark energy, because of having a negative pressure doesn’t hold
this energy condition. The basic reason for the violation of this by dark
energy can be attributed to its anti-gravitating property. This implies that
SEC is satisfied in both radiation and matter dominated eras but violated in
late accelerated epoch. SEC is violated in other scenarios as well,
cosmological inflation is one among them.

  * •

Weak energy condition (WEC)
The energy momentum tensor obeys the inequality $T_{\mu\nu}u^{\mu}u^{\nu}\geq
0$ for any timelike $u^{\mu}$. It is a direct mathematical expression for: the
matter-energy density measured by any observer must not be negative. The
condition can simply be expressed as,

$\rho\geq 0$ (19)

This is the WEC, which is believed to be satisfied by all the cosmic
components throughout the expansion of the universe.

  * •

Dominant energy condition (DEC)
The DEC can be defined as: “for every timelike $u_{\mu}$,
$T^{\mu\nu}u_{\mu}u_{\nu}\geq 0$ and $T^{\mu\nu}u_{\mu}$ is a non-spacelike
vector”[59]. In addition to energy density being non-negative, DEC
accommodates causality as well, i.e, the energy density can’t flow with a
speed faster than that of light and the condition can also be given as,

$\rho\geq|p|$ (20)

DEC holds true for the standard $\Lambda$CDM model. But this condition is
found to be violated in PEDE model during the dark energy dominated late epoch
except at the asymptotic limit where $\omega_{DE}\to-1$

  * •

Null energy condition (NEC)
NEC is implied by SEC and WEC, and stated as $T_{\mu\nu}k^{\mu}k^{\nu}\geq 0$
for every null vector $k^{\mu}$ or,

$\rho+p\geq 0$ (21)

We have analyzed the variation of $\rho+p$ with redshift to check the NEC and
the plot is given in Fig.6.

Figure 6: Violation of NEC shown by PEDE model

It can be observed that NEC is satisfied throughout the evolution for
$\Lambda$CDM model and in the case of PEDE model, it is satisfied during the
initial part of evolution but is violated in the future, when $z\lesssim-0.3$.
The plot indicates that, there occurs a chance for regaining the validity of
NEC in the extreme future. However, we can say that $\rho+p$ is not greater
than zero always, hence the model violates the NEC. The violation of null
energy condition is a general property of phantom cosmology which causes dark
energy density to increase with expansion as suggested by the conservation
equation. The violation of NEC exhibited by phantom cosmologies usually leads
to future singularities. However, in PEDE model, the dark energy density
approaches a constant value rather than increasing boundlessly, there won’t be
any sudden singularities, but still the system may get unstable against
perturbations.

##  VI Statefinder analysis

The presence of a large number of dark energy models necessitates the
comparison of them with the standard model and also among themselves. Such a
comparison is needed to measure the advantage of the model over the standard
model or any other existing models. In 2003, Sahni et al.[63] introduced a
diagnostic tool to facilitate the comparison of cosmological models, using
what is known as statefinder variables, represented as $\\{r,s\\}$. These are
dimensionless geometric variables, constructed from the higher order
derivatives of scale factor with cosmic time [64].

Among statefinder diagnostic pair, $r,$ is called the jerk parameter, which is
defined as,

$r=\frac{\dddot{a}}{aH^{3}}$ (22)

In comparison to the Hubble and deceleration parameters, this parameter
becomes crucial in highlighting the difference between a model and any other
model since it contains the triple derivative of the cosmic scale factor. The
second parameter in the statefinder set $s,$ is constructed by combining the
parameters $q$ and $r$, defined as,

$s=\frac{r-1}{3(q-{1/2})}$ (23)

We can rewrite the equations for $r$ and $s$ in terms of $h.$ As we did
earlier in the case of deceleration parameter, with a change of variable from
$t$ to $x$, the Eqs. 22 and 23 can be rewritten as,

$r=\frac{1}{2h^{2}}\frac{d^{2}h^{2}}{dx^{2}}+\frac{3}{2h^{2}}\frac{dh^{2}}{dx}+1$
(24)
$s=-\left(\dfrac{\dfrac{1}{2h^{2}}\dfrac{d^{2}h^{2}}{dx^{2}}+\dfrac{3}{2h^{2}}\dfrac{dh^{2}}{dx}}{\dfrac{3}{2h^{2}}\dfrac{dh^{2}}{dx}+\dfrac{9}{2}}\right).$
(25)

Figure 7: The evolution of statefinder pair

The evolution of the present model in the statefinder plane is shown in Fig.7.
For reference, note that, for $\Lambda$CDM model, the statefinder pair has a
constant value of $\\{r,s\\}=\\{1,0\\}$ and therefore any departure from this
point signifies the deviation from $\Lambda$CDM model. The $r$ parameter
increases initially in the PEDE model before starting to decline and becoming
closer to 1. On the other hand, the $s$ parameter monotonically rises to 0.
The present values of statefinder parameters for the PEDE model are,
$\\{r_{0},s_{0}\\}=\\{1.469,-0.126\\}$, making the model clearly
distinguishable from the standard model. It can be inferred from Fig.7, the
model approaches the $\Lambda$CDM values of statefinder in the far future,
i.e., when the scale factor goes to infinity. However, during the entire
evolution of the model, until $z\to-1,$ evolution of the parameters satisfies,
$r>1$ and $s<0,$ which implies a Chaplygin gas like behavior with phantom
nature. There are examples of dark energy models in the literature which show
Chaplygin gas behavior with phantom nature [65, 66]. So the predominant
behavior of the model is phantom-like throughout its evolution.

##  VII Thermodynamics of PEDE model

In this section, we will analyze the thermal evolution of the model. We first
consider the status of the generalized second law of thermodynamics and then
the status of the principle of maximization of the entropy.

###  VII.1 Horizon entropy

The second law of thermodynamics constrains the evolution of entropy, as it
states that: “the total entropy of the system should always increase”.The law
also holds true on cosmological scales. The Hubble horizon can be thought of
as the thermodynamic boundary of the universe [67, 68]. By taking account of
horizon entropy, the generalized second law of thermodynamics can be stated as
[69, 70],

$S_{H}^{\prime}+S_{m}^{\prime}\geq 0$ (26)

where $S_{H}$ is the entropy of the Hubble horizon and $S_{m}$ is the entropy
of matter contained within the horizon. The prime denotes derivative with
respect to an appropriate variable like cosmic time or scale factor. Since the
matter entropy is several orders of magnitude less than the horizon
entropy[71], total entropy can be approximated to that of the horizon. The
cosmic horizon entropy has a form similar to black hole event horizon entropy
suggested by Bekenstein [72] and is defined as,

$S_{H}=\frac{A_{H}k_{B}}{4l_{P}^{2}}$ (27)

where $A_{H}=4\pi r_{H}^{2},$ the area of the horizon surface , $k_{B}$ is
Boltzmann constant and $l_{p}$ is the Planck length. $r_{H}$ is the apparent
horizon radius and for a spatially flat universe, $r_{H}=\frac{c}{H}$.
Expressed in terms of $h$, the Eq.27 becomes,

$S_{H}=\frac{\pi c^{5}k_{B}}{\hslash GH_{0}^{2}h^{2}}$ (28)

Substituting the Hubble parameter of the PEDE model in the above equation, the
evolution of the horizon entropy can be studied. Fig.8 presents the plot of
its evolution with redshift. For better comparison, we have shown the entropy
evolution of $\Lambda$CDM model in the same figure.

Figure 8: The evolution of Horizon entropy with redhsift

The figure indicates that, in the PEDE model, the entropy increases first,
attains a maximum value at $z\approx-0.3,$ and then decreases thereafter.
This, violates the generalized $2^{nd}$ law, according to which entropy will
never decrease. To add more clarity, let’s consider the derivative of horizon
entropy with respect to the scale factor. For an increasing entropy evolution,
the derivative must be non-negative.

$\begin{split}S_{H}^{\prime}=&\frac{dS_{H}}{da}\\\ =&-\frac{\pi
c^{5}k_{B}}{\hslash
GH_{0}^{2}}\;\frac{e^{-x}}{h^{4}}\left(\frac{dh^{2}}{dx}\right)\end{split}$
(29)

This shows that since $dh^{2}/dx<0$ up to $z\approx-0.3$, the rate of entropy,
$S_{H}^{\prime}>0,$ which implies an initial increase in entropy. But in
future evolution, where the redshift becomes less than around $-0.3$ the rate
of change of Hubble parameter is such that, $dh^{2}/dx>0,$ corresponding to
which the entropy will decrease. The evolution of $S_{H}^{\prime}$ is plotted
against redshift (Fig.9).

Figure 9: The evolution of the first derivative of Horizon entropy

The generalized second law of thermodynamics is violated as indicated by the
drop in horizon entropy and its negative derivative in the later stages of
evolution. It is to be noted that, the redshift at which the entropy starts
decreasing is the same as that of NEC violation.

###  VII.2 Maximization of entropy

We plan to examine the entropy maximization in order to provide further
insight into the thermal stability of the PEDE model. As already known any
isolated system evolves to a state of equilibrium, at which the entropy is
maximum. The condition, the convexity condition, to be satisfied by a system,
to achieve a final maximum entropy state is that [73],

$S^{\prime\prime}<0\hskip 14.22636pt;\hskip 14.22636pt\text{at least in the
long run}.$ (30)

That is the second derivative of the total entropy with respect to any
convenient cosmological variable, must be less than zero, at least in the end
stage of the evolution. This otherwise suggests that the entropy must be a
convex function of the said variable. The standard $\Lambda$CDM model is found
to satisfy this condition at its end de Sitter epoch [74]. To check the
condition for the present model, second derivative of total entropy is
required. With the earlier assumption of total entropy can be approximated to
horizon entropy, the condition becomes $S_{H}^{\prime\prime}<0,$ where the
derivative can be obtained with respect to the scale factor of the evolution
of the universe. Differentiating Eq.29 gives,

$\begin{split}S_{H}^{\prime\prime}&=\frac{dS_{H}^{\prime}}{da}\\\ &=-\frac{\pi
c^{5}k_{B}}{\hslash
GH_{0}^{2}}\;\frac{e^{-2x}}{h^{4}}\left[\frac{d^{2}h^{2}}{dx^{2}}-\frac{dh^{2}}{dx}-\frac{2}{h^{2}}\left(\frac{dh^{2}}{dx}\right)^{2}\right]\end{split}$
(31)

We found the evolution of $S_{H}^{\prime\prime}$ which is shown in Fig.10.

Figure 10: The evolution of the second derivative of Horizon entropy

It is evident from the figure that the second derivative of the entropy in
PEDE scenario increases first, attains a maximum, and then decreases
approaching zero from above. This is a clear violation of the convexity
condition. This implies that the entropy evolution of PEDE model doesn’t
behave like that of an ordinary thermodynamic system[73]. This boundlessness
in entropy can cause dynamical instabilities in the future evolution of the
universe.

To sum up, the thermal evolution of the PEDE model shows two major issues.
First, it violates the generalized second law of thermodynamics and second,
the entropy is not bounded. These may doom the plausibility of this dark
energy model in comparison with the standard $\Lambda$CDM model.

##  VIII Dynamical system analysis

In the previous section, we analyzed the thermal evolution of the PEDE model
of dark energy and it reveals that the model is breaking the conventional
thermodynamics. One can analyze the dynamical system behavior of the model to
get what it actually implies at the end stage of the universe. Even though we
may get a glimpse of this information from the Hubble parameter evolution, we
are performing here the dynamical system analysis to analyze exactly the
dynamical behavior of the model in the asymptotic limits [75, 76].

To facilitate the analysis, we define the simple dimensionless dynamical
variables,

$u=\frac{\rho_{m}}{3H^{2}}\hskip 5.69046pt,\hskip
5.69046ptv=\frac{\rho_{\scaleto{DE}{3.5pt}}}{3H^{2}}$ (32)

We require a few basic equations to analyze the dynamical system’s behavior
using the differential equations on $u$ and $v$. The first Friedmann equation
(Eq.1),taking $8\pi G=c=1$, can be expressed as,

$3H^{2}=\rho_{m}+\rho_{\scaleto{DE}{3.5pt}}$ (33)

which makes $u+v=1.$ From the conservation equations for non-interacting
matter and dark energy (Eq.5),

$\begin{split}\frac{d\rho_{m}}{dx}&=-3\rho_{m}\\\
\frac{d\rho_{\scaleto{DE}{3.5pt}}}{dx}&=-3\rho_{\scaleto{DE}{3.5pt}}(1+\omega_{\scaleto{DE}{3.5pt}})\end{split}$
(34)

The autonomous coupled differential equations can be obtained from the above
equations as,

$\begin{split}\frac{du}{dx}&=-3u+3u\Big{[}u+v(1+\omega_{\scaleto{DE}{3.5pt}})\Big{]}=f(u,v)\\\
\frac{dv}{dx}&=-3v(1+\omega_{\scaleto{DE}{3.5pt}})+3v\Big{[}u+v(1+\omega_{\scaleto{DE}{3.5pt}})\Big{]}=g(u,v)\end{split}$
(35)

The properties of the critical points, which indicate the equilibrium points,
are of relevance here. To find the critical points, we set $f(u,v)=0$ and
$g(u,v)=0$. The resulting points are $(u_{c},v_{c})=(1,0),(0,0),(0,1)$. It can
be shown that the point (1,0) corresponds to matter dominated epoch while
(0,1) corresponds to the dark energy dominated asymptotic future epoch. The
(0,0) point is a trivial solution and that implies a universe with no matter
and dark energy which is called an empty universe or Milne universe. The point
(0,0) is not taken into more consideration as it yields a solution that does
not have significance in the present context.

To study the stability of the other critical points, the method of linear
perturbation is applied i.e., $u\to u^{\prime}=u_{c}+\delta u$ and $v\to
v^{\prime}=v_{c}+\delta v,$ where $\delta u$ and $\delta v$ are
infinitesimally small deviations from the equilibrium values. Linearizing the
system of differential Eq.35, with respect to the perturbation, we can obtain
a matrix equation of the form,

$\begin{bmatrix}\delta u^{\prime}\\\ \delta
v^{\prime}\end{bmatrix}=\begin{bmatrix}\left(\dfrac{\partial f}{\partial
u}\right)_{c}&\left(\dfrac{\partial f}{\partial v}\right)_{c}\\\
\left(\dfrac{\partial g}{\partial u}\right)_{c}&\left(\dfrac{\partial
g}{\partial v}\right)_{c}\end{bmatrix}\begin{bmatrix}\delta u\\\ \delta
v\end{bmatrix}$ (36)

The $2\times 2$ matrix on the right hand side is called the Jacobian matrix.
The Jacobian matrix obtained for the PEDE model is given as,

$\begin{bmatrix}-3+6u+3v(1+\omega_{\scaleto{DE}{3.5pt}})&3u(1+\omega_{\scaleto{DE}{3.5pt}})\\\
3v&3u+(1+\omega_{\scaleto{DE}{3.5pt}})(6v-3)\end{bmatrix}$ (37)

The eigenvalues can be found by diagonalizing the matrix and then the
stability of the system depends upon the sign of eigenvalues. If both the
eigenvalues evaluated for the critical point are negative, the system is
stable. The phase space trajectories, irrespective of their origin will
converge to that critical point. The system is said to be unstable if both the
eigenvalues are positive, and hence the trajectories will diverge away from
the critical point. If one of the eigenvalues is positive and the other is
negative, the system is said to be in saddle at that critical point.

Figure 11: The phase space diagram corresponding to
$\omega_{\scaleto{DE}{3.5pt}}(z=-1)=-1$

For PEDE model, the eigenvalues corresponding to the critical point (1,0) are
obtained as ($3$ , $-3\omega_{\scaleto{DE}{3.5pt}}$). Since the equation of
state of dark energy is found to be less than -1, it turns out that, both the
eigenvalues will become positive, making the system unstable in the matter
dominated epoch. The instability drives to further expansion of the universe.
The second critical point (0,1), refers to the future state and has
eigenvalues,
$(3\omega_{\scaleto{DE}{3.5pt}},3(1+\omega_{\scaleto{DE}{3.5pt}}))$. In the
asymptotic state, which corresponds to the de Sitter epoch, the equation of
state of dark energy becomes $-1.$ It will lead to a set of eigenvalues
$(-3,0$), with one being negative while the other is zero. The nature of this
critical point is not trivial. Let us now look at the phase space trajectory,
given in Fig.11 for these critical points. It suggests that the critical point
(1,0) indeed is a source point, as the trajectories are diverging from it.
While regarding the second equilibrium point, (0,1), it is seen that the
trajectories appear to converge to a line, rather than to a point, parallel to
the $v$ axis corresponding to the dark energy in the dynamical system plane.
To be specific, the primary feature of the ‘stable point’ is that, the
surrounding trajectories must converge to that point. Since that feature is
not strictly obeyed here, there is no guarantee that the critical point is a
stable one. At the same time, in an infinitesimal region around the point
(0,1) the trajectories can be considered converging. Let us analyze a bit
further regarding the nature of the critical point (0,1). There are papers in
the literature that argues, if one of the eigenvalues of a critical point is
zero, the nature of that can be judged from the sign of the other eigenvalue
[77]. However, it is to be noted from the thermodynamic analysis, that the
entropy is not maximizing in the asymptotic future corresponding to the
critical point (0,1). The decreasing entropy as observed in Fig.8 makes this
state of the system thermodynamically unstable in the future epoch and
correspondingly the critical point (0,1) can said to be thermally unstable.
But the dynamical system analysis based on the phase space diagram is
suggesting that the system is at least locally stable. That is we come across
a situation where the critical point (0,1) is thermodynamically unstable but
dynamically it seems to be apparently stable at least in the local sense. To
clarify this apparent paradoxical situation further we will analyze the
behavior of sound speed in the following subsection.

###  VIII.1 Dark energy sound speed

The dark energy sound speed is described by the ratio of pressure to density
variations as[78],

$v_{s}^{2}=\dfrac{dp_{\scaleto{DE}{3.5pt}}}{d\rho_{\scaleto{DE}{3.5pt}}}$ (38)

It has been argued that the physical values of $v_{s}^{2}$ ought to be in the
range from 0 to 1. Outside of this region, one encounters gradient/tachyonic
instabilities and/or instabilities linked to superluminal propagation[79]. A
real value of sound speed shows a regular propagating mode for density
perturbation while an imaginary value corresponds to an irregular wave
equation implying a classical instability of a given perturbation. Since a
cosmological constant suffers no spatial perturbations and cosmology with dark
energy in the form of a cosmological constant thus has no sensitivity to the
dark energy sound speed, it is vital to keep in mind that the study is only
applicable when the equation of state of dark energy,
$\omega_{\scaleto{DE}{4pt}}\neq-1$[78]. One can rewrite the equation for
$v_{s}^{2}$ as,

$v_{s}^{2}=\dfrac{\dot{p}_{\scaleto{DE}{3.5pt}}}{\dot{\rho}_{\scaleto{DE}{3.5pt}}}$
(39)

where
$\dot{p}_{\scaleto{DE}{3.5pt}}=\omega_{\scaleto{DE}{3.5pt}}\dot{\rho}_{\scaleto{DE}{3.5pt}}+\dot{\omega}_{\scaleto{DE}{3.5pt}}\rho_{\scaleto{DE}{3.5pt}}$.
Substituting with PE dark energy, the above equation reduces to,

$v_{s}^{2}=-\Big{(}\dfrac{2\tanh\big{[}\log_{10}(1+z)\big{]}}{3\ln
10}+1\Big{)}$ (40)

When $z\rightarrow-1$, $v_{s}^{2}$ reaches to a value $\approx-0.7$. The
Fig.12 gives the complete evolution of $v_{s}^{2}$.

Figure 12: The evolution of squared sound speed with redshift

It is evident from the figure that the value of $v_{s}^{2}$ is negative
including the far future epoch marking instabilities at the perturbation
level. The perturbation has become an irregular wave equation. This negative
squared speed shows an exponentially growing mode for a density perturbation.
In other words, an increasing density perturbation induces a lowering
pressure, supporting the emergence of instability[78, 80, 81]. Moreover, the
increment of instability is inversely proportional to the wavelength of the
perturbations[82], and therefore the PEDE model for which $v_{s}^{2}<0$ is
asymptotically unstable. Previously when we analyzed the stability conditions
corresponding to the critical point (0,1), we found that the corresponding
state seems to be stable only locally. The thermal instability of the same
state compelled us to investigate further. In this regard, the acoustic speed
evolution now shows that this asymptotic state is not stable. Hence taking
account of the thermal evolution and evolution of the sound speed we have to
conclude that the critical point (0,1) may be unstable.

Before concluding the section we would like to mention a recent work in which
the authors have performed the dynamical stability analysis of the very same
model. Hernández-Almada et al.[83] have studied the dynamical stability of the
PEDE model of dark energy, in a three-dimensional plane with variables,
$T=H_{0}/(H_{0}+H),\,\Omega_{m}=H_{0}^{2}\Omega_{m0}/(a^{3}H^{2}),\,\Omega_{\gamma}=H_{0}^{2}\Omega_{\gamma
0}/(a^{4}H^{2})$. Interestingly they obtained the eigenvalues corresponding to
the future epoch as all negative, implying a stable nature and from the 3D
phase diagram (Fig.7 in [83]), the trajectories appear to be converging to a
point in T axis. It is also to be noted that, the authors[83], predict an
effect of the PE dark energy in the very early stage of the evolution of the
universe, such that it can even drive the early acceleration equivalent to
inflation. This argument is in fact contrary to the claim of originators of
this very model, that the PEDE dark energy will emerge as the universe expand
such that it can’t have any commendable effect on the early stage of
evolution. Moreover, they haven’t studied the thermal evolution of the model.
However, it is essentially difficult to comprehend the dynamical stability of
a system that is thermally unstable.

##  IX Conclusions

In this work, we explore the background evolution of the cosmological model,
with the phenomenologically emergent dark energy. We analyzed the evolution of
basic parameters and the status of the model using statefinder diagnostic. We
also examined the thermodynamic evolution and dynamical system behavior of the
model.

Major observations derived from our study are given below.

  * •

Unlike the standard model, the PEDE model offers variable dark energy that
grows as the universe expands. The evolution of the equation of state shows
that the dark energy exhibits phantom behavior
($\omega_{\scaleto{DE}{4pt}}<-1$) till the asymptotic future epoch. Due to the
phantom energy, the Hubble parameter gets shifted to higher values for the
present epoch thereby alleviating the Hubble tension. The parameter values
obtained for the PEDE model from our analysis are: $H_{0}=71.966\pm 0.618\
$\mathrm{k}\mathrm{m}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$ and
$\Omega_{m0}=0.280\pm 0.006$ for the data combination $H(z)$+Pantheon+CMB+BAO.

  * •

The evolution of Hubble and deceleration parameters are found to be much
different from the standard model. Hubble parameter decreases initially, as
seen in Fig.4, but after $z\approx-0.3$, it appears to be increasing. The
increase in $H$ during the late phase is due to the dominating phantom fluid
which is not expected in the standard cosmological scenario.

The evolution of the deceleration parameter also reflects the divergent
behavior of the Hubble parameter and equation of state (see Fig.5). Like the
standard model, the PEDE model predicts an earlier decelerated epoch dominated
by matter and a later accelerated period dominated by dark energy. The
redshift of the transition from decelerated to accelerated epoch is determined
as $z\approx 0.7$. However, owing to the growth of phantom dark energy density
in PEDE model, the value of $q$ in the future decreases from -1 at
$z\approx-0.3$ and eventually returns to -1 at the very end.

  * •

When we investigated the model’s adherence to the null energy condition in
Fig.6, we observed that it is violated at the same redshift that the
deceleration parameter initially passes the -1 value and the Hubble parameter
begins to rise. The NEC’s violation suggests that dark energy density is
increasing in scale factor powers, which could make the universe unstable. The
statefinder analysis for the model is then performed. The anomaly is also
evident here. In contrast to the standard model, which has values of
$\\{r,s\\}=\\{1,0\\}$ throughout the evolution, the PEDE model has $r$ values
that are more than 1 and $s$ values that are less than zero during the late
phase except the very end. Statefinder analysis reassures the presence of
phantom fluid.

  * •

We analyzed the thermodynamic evolution of the model. The model doesn’t meet
the entropy criteria according to the generalized $2^{\text{nd}}$ law
pertaining to the entropy of the universe. After a redshift of $z\approx-0.3$,
the entropy is observed to be declining, and its first order derivative turns
negative. The model also fails to meet the convexity requirement since the
second derivative approaches zero from above, indicating that the entropy is
not maximizing revealing an unstable asymptotic future epoch.

  * •

We performed the phase space analysis to extract the stability conditions of
the model. The matter-dominated epoch appeared to be unstable, as both
eigenvalues were positive, implying that the universe would expand further.
However, it is found that one of the eigenvalues for the asymptotic de Sitter
epoch is zero and the other is negative. The paths appeared to be converging
on a line as shown in Fig.11. We draw the conclusion that the endpoint may be
unstable based on our findings of thermodynamic analysis. The stability
analysis using dark energy sound speed supports the prediction of an unstable
future epoch.

At the outset, our general conclusion is that, even though the PEDE model may
alleviate the Hubble tension, its general behavior is not all promising. The
phantom behavior of the model leads to NEC violation and thermodynamically
unstable future epoch.

At this juncture, we would like to mention a different phenomenological model,
which is different from the present case. This model consists of two dark
energy components, one of which is Einstein’s cosmological constant and the
other a variable dark energy component that would have emerged from the
remnants of the scalar field responsible for inflation after the majority of
the scalar field had transformed into matter[84, 85, 86]. The authors claim
that the model can alleviate the Hubble tension by the action of the
additional dark energy component at the stages after recombination. In the
presence of the variable dark energy component, the Hubble constant decreases
with time more slowly than in the $\Lambda$CDM model. The model doesn’t leads
to any phantom behaviour and doesn’t seem to violate the entropy maximization
in the asymptotic future epoch. However, the effect of this model on CMB power
spectrum and on other datasets need to be studied in detail. Within the realm
of dynamical dark energy proposals, the model seems to provide a plausible
explanation for the phenomenological evolution of dark energy content of the
universe.

## Acknowledgments

One of the authors Sarath N is thankful to UGC, Govt. of India, for providing
financial support through Junior Research fellowship. The authors thank
Dheepika M. and Manosh T. M. for discussions.

## Data Availability

The pantheon data underlying this article are available at
https://github.com/dscolnic/Pantheon. The H(z) data is from Ref.[46]. The CMB
data is available in Ref.[48] and BAO data is from Ref.[49]

## References

  * [1] X. Li, A. Shafieloo, The Astrophysical Journal 883(1), L3 (2019).  DOI 10.3847/2041-8213/ab3e09.  URL https://doi.org/10.3847%2F2041-8213%2Fab3e09
  * [2] A.G. Riess, A.V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P.M. Garnavich, R.L. Gilliland, C.J. Hogan, S. Jha, R.P. Kirshner, et al., The Astronomical Journal 116(3), 1009 (1998)
  * [3] S. Perlmutter, G. Aldering, G. Goldhaber, R. Knop, P. Nugent, P.G. Castro, S. Deustua, S. Fabbro, A. Goobar, D.E. Groom, et al., The Astrophysical Journal 517(2), 565 (1999)
  * [4] L. Perivolaropoulos, F. Skara, New Astronomy Reviews 95, 101659 (2022).  DOI 10.1016/j.newar.2022.101659.  URL https://doi.org/10.1016%2Fj.newar.2022.101659
  * [5] P. Bull, Y. Akrami, J. Adamek, T. Baker, E. Bellini, J.B. Jimenez, E. Bentivegna, S. Camera, S. Clesse, J.H. Davis, et al., Physics of the Dark Universe 12, 56 (2016)
  * [6] S. Weinberg, Reviews of modern physics 61(1), 1 (1989)
  * [7] E.J. Copeland, M. Sami, S. Tsujikawa, International Journal of Modern Physics D 15(11), 1753 (2006)
  * [8] N. Sivanandam, Physical Review D 87(8), 083514 (2013)
  * [9] H.E. Velten, R. Vom Marttens, W. Zimdahl, The European Physical Journal C 74, 1 (2014)
  * [10] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, et al., Astronomy & Astrophysics 641, A6 (2020)
  * [11] A.G. Riess, S. Casertano, W. Yuan, L.M. Macri, D. Scolnic, The Astrophysical Journal 876(1), 85 (2019).  DOI 10.3847/1538-4357/ab1422.  URL https://doi.org/10.3847%2F1538-4357%2Fab1422
  * [12] E.D. Valentino, O. Mena, S. Pan, L. Visinelli, Weiqiang, A. Melchiorri, D.F. Mota, A.G. Riess, J. Silk, Classical and Quantum Gravity 38(15), 153001 (2021).  DOI 10.1088/1361-6382/ac086d.  URL https://doi.org/10.1088/1361-6382/ac086d
  * [13] S.D. Odintsov, D.S.C. Gómez, G.S. Sharov, Nuclear Physics B 966, 115377 (2021)
  * [14] C. Armendariz-Picon, V. Mukhanov, P.J. Steinhardt, Phys. Rev. Lett. 85, 4438 (2000).  DOI 10.1103/PhysRevLett.85.4438.  URL https://link.aps.org/doi/10.1103/PhysRevLett.85.4438
  * [15] S. Tsujikawa, Classical and Quantum Gravity 30(21), 214003 (2013)
  * [16] S. DeDeo, R. Caldwell, P.J. Steinhardt, Physical Review D 67(10), 103509 (2003)
  * [17] R. Caldwell, Physics Letters B 545(1-2), 23 (2002).  DOI 10.1016/s0370-2693(02)02589-3.  URL https://doi.org/10.1016%2Fs0370-2693%2802%2902589-3
  * [18] R.R. Caldwell, M. Kamionkowski, N.N. Weinberg, Physical review letters 91(7), 071301 (2003)
  * [19] Z.K. Guo, Y.S. Piao, X. Zhang, Y.Z. Zhang, Physics Letters B 608(3-4), 177 (2005).  DOI 10.1016/j.physletb.2005.01.017.  URL https://doi.org/10.1016%2Fj.physletb.2005.01.017
  * [20] Y.F. Cai, E.N. Saridakis, M.R. Setare, J.Q. Xia, Physics Reports 493(1), 1 (2010)
  * [21] S. Kumar, R.C. Nunes, Phys. Rev. D 94, 123511 (2016).  DOI 10.1103/PhysRevD.94.123511.  URL https://link.aps.org/doi/10.1103/PhysRevD.94.123511
  * [22] S. Pan, G.S. Sharov, W. Yang, Phys. Rev. D 101, 103533 (2020).  DOI 10.1103/PhysRevD.101.103533.  URL https://link.aps.org/doi/10.1103/PhysRevD.101.103533
  * [23] E.D. Valentino, A. Melchiorri, O. Mena, S. Vagnozzi, Physics of the Dark Universe 30, 100666 (2020).  DOI 10.1016/j.dark.2020.100666.  URL https://doi.org/10.1016%2Fj.dark.2020.100666
  * [24] P. Avelino, L. Beca, J. De Carvalho, C. Martins, Journal of Cosmology and Astroparticle Physics 2003, 002 (2003).  DOI 10.1088/1475-7516/2003/09/002
  * [25] L. Amendola, F. Finelli, C. Burigana, D. Carturan, Journal of Cosmology and Astroparticle Physics 2003(07), 005 (2003).  DOI 10.1088/1475-7516/2003/07/005.  URL https://doi.org/10.1088%2F1475-7516%2F2003%2F07%2F005
  * [26] H.B. Sandvik, M. Tegmark, M. Zaldarriaga, I. Waga, Phys. Rev. D 69, 123524 (2004).  DOI 10.1103/PhysRevD.69.123524.  URL https://link.aps.org/doi/10.1103/PhysRevD.69.123524
  * [27] M. Li, Physics Letters B 603(1-2), 1 (2004)
  * [28] T.P. Sotiriou, V. Faraoni, Rev. Mod. Phys. 82, 451 (2010).  DOI 10.1103/RevModPhys.82.451.  URL https://link.aps.org/doi/10.1103/RevModPhys.82.451
  * [29] A.D. Felice, S. Tsujikawa, Living Reviews in Relativity 13(1) (2010).  DOI 10.12942/lrr-2010-3.  URL https://doi.org/10.12942%2Flrr-2010-3
  * [30] S. Nojiri, S. Odintsov, V. Oikonomou, Physics Reports 692, 1 (2017).  DOI 10.1016/j.physrep.2017.06.001.  URL https://doi.org/10.1016%2Fj.physrep.2017.06.001
  * [31] R.C. Nunes, Journal of Cosmology and Astroparticle Physics 2018(05), 052 (2018)
  * [32] M. Rossi, M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A.A. Starobinsky, C. Umilta, Physical Review D 100(10), 103524 (2019)
  * [33] M. Braglia, M. Ballardini, W.T. Emond, F. Finelli, A.E. Gümrükçüoğlu, K. Koyama, D. Paoletti, Physical Review D 102(2), 023529 (2020)
  * [34] M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A.A. Starobinsky, C. Umiltà, Journal of Cosmology and Astroparticle Physics 2020(10), 044 (2020).  DOI 10.1088/1475-7516/2020/10/044.  URL https://doi.org/10.1088%2F1475-7516%2F2020%2F10%2F044
  * [35] T. Adi, E.D. Kovetz, Physical Review D 103(2), 023530 (2021)
  * [36] S.F. Yan, P. Zhang, J.W. Chen, X.Z. Zhang, Y.F. Cai, E.N. Saridakis, Physical Review D 101(12) (2020).  DOI 10.1103/physrevd.101.121301.  URL https://doi.org/10.1103%2Fphysrevd.101.121301
  * [37] X. Li, A. Shafieloo, The Astrophysical Journal 902(1), 58 (2020).  DOI 10.3847/1538-4357/abb3d0.  URL https://doi.org/10.3847/1538-4357/abb3d0
  * [38] H. Benaoum, W. Yang, S. Pan, E. Di Valentino, International Journal of Modern Physics D 31(03), 2250015 (2022)
  * [39] W. Yang, E. Di Valentino, S. Pan, O. Mena, Physics of the Dark Universe 31, 100762 (2021)
  * [40] M. Rezaei, T. Naderi, M. Malekjani, A. Mehrabi, The European Physical Journal C 80(5) (2020).  DOI 10.1140/epjc/s10052-020-7942-6.  URL https://doi.org/10.1140%2Fepjc%2Fs10052-020-7942-6
  * [41] S. Pan, W. Yang, E. Di Valentino, A. Shafieloo, S. Chakraborty, Journal of Cosmology and Astroparticle Physics 2020(06), 062 (2020)
  * [42] G. Alestas, D. Camarena, E. Di Valentino, L. Kazantzidis, V. Marra, S. Nesseris, L. Perivolaropoulos, Phys. Rev. D 105, 063538 (2022).  DOI 10.1103/PhysRevD.105.063538.  URL https://link.aps.org/doi/10.1103/PhysRevD.105.063538
  * [43] N. Schöneberg, G.F. Abellán, A.P. Sánchez, S.J. Witte, V. Poulin, J. Lesgourgues.  The $h_{0}$ olympics: A fair ranking of proposed models (2021).  DOI 10.48550/ARXIV.2107.10291.  URL https://arxiv.org/abs/2107.10291
  * [44] S. Sharma, Physics of the Dark Universe 30, 100717 (2020).  DOI https://doi.org/10.1016/j.dark.2020.100717.  URL https://www.sciencedirect.com/science/article/pii/S2212686420304301
  * [45] R. Caldwell, Physics Letters B 545(1), 23 (2002).  DOI https://doi.org/10.1016/S0370-2693(02)02589-3.  URL https://www.sciencedirect.com/science/article/pii/S0370269302025893
  * [46] J.J. Geng, R.Y. Guo, A.Z. Wang, J.F. Zhang, X. Zhang, Communications in Theoretical Physics 70(4), 445 (2018).  DOI 10.1088/0253-6102/70/4/445.  URL https://doi.org/10.1088%2F0253-6102%2F70%2F4%2F445
  * [47] D.M. Scolnic, D. Jones, A. Rest, Y. Pan, R. Chornock, R. Foley, M. Huber, R. Kessler, G. Narayan, A. Riess, et al., The Astrophysical Journal 859(2), 101 (2018)
  * [48] L. Chen, Q.G. Huang, K. Wang, Journal of Cosmology and Astroparticle Physics 2019(02), 028 (2019)
  * [49] C. Blake, E.A. Kazin, F. Beutler, T.M. Davis, D. Parkinson, S. Brough, M. Colless, C. Contreras, W. Couch, S. Croom, D. Croton, M.J. Drinkwater, K. Forster, D. Gilbank, M. Gladders, K. Glazebrook, B. Jelliffe, R.J. Jurek, I.h. Li, B. Madore, D.C. Martin, K. Pimbblet, G.B. Poole, M. Pracy, R. Sharp, E. Wisnioski, D. Woods, T.K. Wyder, H.K.C. Yee, Monthly Notices of the Royal Astronomical Society 418(3), 1707 (2011).  DOI 10.1111/j.1365-2966.2011.19592.x.  URL https://doi.org/10.1111/j.1365-2966.2011.19592.x
  * [50] D. Foreman-Mackey, D.W. Hogg, D. Lang, J. Goodman, Publications of the Astronomical Society of the Pacific 125(925), 306 (2013).  DOI 10.1086/670067.  URL https://dx.doi.org/10.1086/670067
  * [51] M. Newville, R. Otten, A. Nelson, A. Ingargiola, T. Stensitzki, D. Allan, A. Fox, F. Carter, Michał, D. Pustakhod, lneuhaus, S. Weigand, R. Osborn, Glenn, C. Deil, Mark, A.L.R. Hansen, G. Pasquevich, L. Foks, N. Zobrist, O. Frost, A. Beelen, Stuermer, kwertyops, A. Polloreno, S. Caldwell, A. Almarza, A. Persaud, B. Gamari, B.F. Maier.  lmfit/lmfit-py: (2021).  DOI 10.5281/zenodo.4516644.  URL https://doi.org/10.5281/zenodo.4516644
  * [52] S. N, T.K. Mathew, Monthly Notices of the Royal Astronomical Society 510(4), 5553 (2021).  DOI 10.1093/mnras/stab3773.  URL https://doi.org/10.1093%2Fmnras%2Fstab3773
  * [53] M. Dheepika, T.K. Mathew, The European Physical Journal C 82(5) (2022).  DOI 10.1140/epjc/s10052-022-10365-2.  URL https://doi.org/10.1140%2Fepjc%2Fs10052-022-10365-2
  * [54] R. Jimenez, A. Loeb, The Astrophysical Journal 573(1), 37 (2002).  DOI 10.1086/340549.  URL https://doi.org/10.1086%2F340549
  * [55] E. Gaztañaga, A. Cabrè, L. Hui, Monthly Notices of the Royal Astronomical Society 399(3), 1663 (2009).  DOI 10.1111/j.1365-2966.2009.15405.x.  URL https://doi.org/10.1111%2Fj.1365-2966.2009.15405.x
  * [56] A.G. Riess, L.M. Macri, S.L. Hoffmann, D. Scolnic, S. Casertano, A.V. Filippenko, B.E. Tucker, M.J. Reid, D.O. Jones, J.M. Silverman, R. Chornock, P. Challis, W. Yuan, P.J. Brown, R.J. Foley, The Astrophysical Journal 826(1), 56 (2016).  DOI 10.3847/0004-637x/826/1/56.  URL https://doi.org/10.3847%2F0004-637x%2F826%2F1%2F56
  * [57] J. Guy, M. Sullivan, A. Conley, N. Regnault, P. Astier, C. Balland, S. Basa, R.G. Carlberg, D. Fouchez, D. Hardin, I.M. Hook, D.A. Howell, R. Pain, N. Palanque-Delabrouille, K.M. Perrett, C.J. Pritchet, J. Rich, V. Ruhlmann-Kleider, D. Balam, S. Baumont, R.S. Ellis, S. Fabbro, H.K. Fakhouri, N. Fourmanoit, S. Gonzá lez-Gaitán, M.L. Graham, E. Hsiao, T. Kronborg, C. Lidman, A.M. Mourao, S. Perlmutter, P. Ripoche, N. Suzuki, E.S. Walker, Astronomy &Astrophysics 523, A7 (2010).  DOI 10.1051/0004-6361/201014468.  URL https://doi.org/10.1051%2F0004-6361%2F201014468
  * [58] S. Bocquet, F.W. Carter, Journal of Open Source Software 1(6), 46 (2016).  DOI 10.21105/joss.00046.  URL https://doi.org/10.21105/joss.00046
  * [59] S.W. Hawking, G.F.R. Ellis, _The large scale structure of space-time_ , vol. 1 (Cambridge university press, 1973)
  * [60] S.M. Carroll, M. Hoffman, M. Trodden, Physical Review D 68(2), 023509 (2003)
  * [61] J. Santos, J. Alcaniz, N. Pires, M.J. Reboucas, Physical Review D 75(8), 083523 (2007)
  * [62] M. Bouhmadi-Ló pez, A. Errahmani, P. Martín-Moruno, T. Ouali, Y. Tavakoli, International Journal of Modern Physics D 24(10), 1550078 (2015).  DOI 10.1142/s0218271815500789.  URL https://doi.org/10.1142%2Fs0218271815500789
  * [63] V. Sahni, T.D. Saini, A.A. Starobinsky, U. Alam, Journal of Experimental and Theoretical Physics Letters 77(5), 201 (2003)
  * [64] U. Alam, V. Sahni, T. Deep Saini, A. Starobinsky, Monthly Notices of the Royal Astronomical Society 344(4), 1057 (2003)
  * [65] Y.B. Wu, S. Li, M.H. Fu, J. He, General Relativity and Gravitation 39, 653 (2007)
  * [66] H.M. Sadjadi, Physics Letters B 687(2-3), 114 (2010)
  * [67] P.C.W. Davies, Classical and Quantum Gravity 4(6), L225 (1987).  DOI 10.1088/0264-9381/4/6/006.  URL https://doi.org/10.1088/0264-9381/4/6/006
  * [68] T.M. Davis, P.C.W. Davies, C.H. Lineweaver, Classical and Quantum Gravity 20(13), 2753 (2003).  DOI 10.1088/0264-9381/20/13/322.  URL https://doi.org/10.1088%2F0264-9381%2F20%2F13%2F322
  * [69] M. Setare, Physics Letters B 641(2), 130 (2006)
  * [70] N.J. Mohan, P. Krishna, A. Sasidharan, T.K. Mathew, Classical and Quantum Gravity 37(7), 075007 (2020)
  * [71] C.A. Egan, C.H. Lineweaver, The Astrophysical Journal 710(2), 1825 (2010).  DOI 10.1088/0004-637x/710/2/1825.  URL https://doi.org/10.1088/0004-637x/710/2/1825
  * [72] J.D. Bekenstein, Phys. Rev. D 7, 2333 (1973).  DOI 10.1103/PhysRevD.7.2333.  URL https://link.aps.org/doi/10.1103/PhysRevD.7.2333
  * [73] D. Pavon, N. Radicella, General Relativity and Gravitation 45 (2012).  DOI 10.1007/s10714-012-1457-x
  * [74] P. Krishna, T.K. Mathew, Physical Review D 96(6), 063513 (2017)
  * [75] J. Wainwright, Classical and Quantum Gravity 6(10), 1409 (1989)
  * [76] T.K. Mathew, et al., arXiv preprint arXiv:2204.12097 (2022)
  * [77] N. Banerjee, N. Roy, General Relativity and Gravitation 47(8), 1 (2015)
  * [78] U.K. Sharma, V.C. Dubey.  Interacting rényi holographic dark energy with parametrization on the interaction term (2020).  DOI 10.48550/ARXIV.2001.02368.  URL https://arxiv.org/abs/2001.02368
  * [79] S. Vagnozzi, L. Visinelli, O. Mena, D.F. Mota, Monthly Notices of the Royal Astronomical Society 493(1), 1139 (2020).  DOI 10.1093/mnras/staa311.  URL https://doi.org/10.1093%2Fmnras%2Fstaa311
  * [80] K.Y. Kim, H.W. Lee, Y.S. Myung, Physics Letters B 660(3), 118 (2008).  DOI https://doi.org/10.1016/j.physletb.2007.12.045.  URL https://www.sciencedirect.com/science/article/pii/S0370269308000038
  * [81] Y.S. Myung, Physics Letters B 652(5-6), 223 (2007).  DOI 10.1016/j.physletb.2007.07.033.  URL https://doi.org/10.1016%2Fj.physletb.2007.07.033
  * [82] C. Armendariz-Picon, T. Damour, V. Mukhanov, Physics Letters B 458(2-3), 209 (1999)
  * [83] A. Hernández-Almada, G. Leon, J. Magaña, M.A. García-Aspeitia, V. Motta, Monthly Notices of the Royal Astronomical Society 497(2), 1590 (2020)
  * [84] G. Bisnovatyi-Kogan, arXiv preprint arXiv:2002.05602 (2020)
  * [85] G.S. Bisnovatyi-Kogan, Universe 7(11), 412 (2021)
  * [86] G. Bisnovatyi-Kogan, A. Nikishin, Astronomy Reports 67(2), 115 (2023)