Thermal evolution and stability analysis of phenomenologically emergent dark energy model

The phenomenologically emergent dark energy (PEDE) model is a varying dark energy model with no extra degrees of freedom proposed by Li and Shafieloo (Astrophys J 883(1):L3, 2019) 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 Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda $$\end{document}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.


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, ω (ω = p/ρ, p is the pressure and ρ is the density), that they termed 'dark energy'.The current standard model, the Λ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 Λ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 Λ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 ΛCDM, provided a value of H 0 = 67.4± 0.5 kms −1 Mpc −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± 1.42 kms −1 Mpc −1 [11].The inconsistency between these two measurements is over 4σ.The tension has now increased over time as a result of examining multiple data sets gathered using various techniques [4,12,13].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.
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α data finds H 0 = 71.02+1. 45 −1.37 kms −1 Mpc −1 [1] thereby arXiv:2301.12172v3[gr-qc] 8 Aug 2023 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 Λ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 Λ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σ level in this scenario, and the value of σ 8 is in better agreement with cosmic shear estimates.They have obtained M ν ∼ 0.21 +0.15  −0.14 eV and N ef f = 3.03±0.32,which point to a non-zero neutrino mass with a significance greater than 2σ.However, they obtained the sound horizon as ∼ 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 σ 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 ΛCDM.According to the Bayesian analysis performed on the model by Pan et al. [41], the majority of observational dataset combinations favor Λ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 Λ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 ω < −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 (APEDE) via the Bianchi type-I spacetime metric.When compared to the base ΛCDM, the statistical results for both the PEDE and APEDE 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 Λ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, where H is the Hubble parameter and G is the gravitational constant.ρ γ , ρ m and ρ DE represent the density of radiation, matter and dark energy respectively.The density components can be expressed in terms of their corresponding dimensionless density parameters, Ω i = ρ i /ρ c,0 where ρ c,0 = 3H 2 0 /8πG is the critical density.Eq.1 can then be rewritten as, (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, Ω DE,0 .But, in general the evolution of Ω DE (z) with redshift can be expressed as, where ω DE 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: 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 ρ γ will be negligibly small in late phase.)If matter and dark energy have no interaction, then the conservation equations are as follows: where ω m = 0 and ω DE can be derived as, Substituting the PE dark energy and the equation becomes, Both Ω DE and ω DE follows a symmetrical behavior in logarithmic scale(Fig.1 in [1]).For an early epoch, i.e., when z is large and positive, ω DE converges to a value of − 2 3 ln 10 − 1.As time progresses the equation of state increases, becomes − 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 ω DE < −1 [45] except at the asymptotic limit, z → −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 Λ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 kms −1 Mpc −1 and Ω m0 in the range, 0.1 to 0.5 .The best fit values are then calculated by minimizing the χ 2 function defined as [52,53], where V obs i is the observed value from the datasets, V th i is the corresponding theoretical value obtained from our model and σ i is the error in the measured values.The minimum χ 2 per degrees of freedom or the reduced χ 2 is defined as, χ 2 d.o.f = χ 2 N −Np , 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, Similar to χ 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 ≤ z ≤ 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± 1.74 kms −1 Mpc −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, 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 (α) and the color-luminosity parameter (β) are calibrated to zero for the pantheon sample, hence the observed distance modulus reduces to µ obs i = m−M .For our analysis, we treat M as a nuisance parameter and is allowed to vary in the range, −20 ≤ M ≤ −18.
The luminosity distance of these SN Ia can be calculated from the equation, where z i is the redshift of the SN Ia, c is the speed of light and H(H 0 , Ω m0 , z ′ ) is the Hubble parameter specified for the model.The theoretical distance modulus of SN Ia is then given by, • CMB: The theoretical shift parameter R of the cosmic microwave background is defined as, 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 ± 0.0046 at the redshift z * = 1089.2.By considering this Planck measurement corresponding χ 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], where z A is the redshift corresponding to the acoustic peak.To obtain the χ 2 BAO , we compare the theoretical value with A obs = 0.484 ± 0.016 at the redshift z A = 0.35 obtained from SDSS-BAO distance data [49].
The best estimated parameters for the ΛCDM and PEDE model within 68.3% confidence level are given in Tab.I and Tab.II 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 χ 2 value as well as χ 2 d.o.f are also mentioned in the tables.The χ 2 for the combination of datasets is obtained by adding the individual χ 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 Λ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 ΛCDM except when the Pantheon sample is considered alone.For the Pantheon dataset the χ 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 Ω m0 value in the case of PEDE 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 ΛCDM and PEDE model parameters  −0.618 for PEDE with a χ 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 χ 2 , AIC and BIC values.Similar problem of alleviation of tension by worsening the χ 2 values compared to the ΛCDM model has been pointed out in Ref. [41].The χ 2 value listed in both tables points out the fact that ΛCDM is still favored over PEDE model when the datasets are considered together.

IV. EVOLUTION OF COSMOLOGICAL PARAMETERS A. Hubble parameter
The late phase evolution of the Hubble parameter in PEDE model for a flat universe can be expressed as, where Ω DE,0 = 1 − Ω 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.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 Λ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 ≈ −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 → −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.

B. 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, The over dot represents derivative over cosmic time.This equation can suitably be expressed in terms of a convenient variable, x = ln a, where h = H/H 0 .The evolution of deceleration parameter for both the canonical and PEDE models are plotted in Fig. 5.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 Λ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.56respectively.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 ≈ −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 Λ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 µν on physical grounds [59][60][61][62].

• Strong energy condition (SEC)
The SEC states that (T µν − 1 2 T g µν )u µ u ν ≥ 0 for all time like vectors u µ .SEC is satisfied by gravitationally interacting cosmic components.For a perfect fluid, it is equal as, ρ + 3p ≥ 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 µν u µ u ν ≥ 0 for any timelike u µ .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, 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 µ , T µν u µ u ν ≥ 0 and T µν u µ 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, DEC holds true for the standard Λ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 ω DE → −1

• Null energy condition (NEC)
NEC is implied by SEC and WEC, and stated as T µν k µ k ν ≥ 0 for every null vector k µ or, We have analyzed the variation of ρ + p with redshift to check the NEC and the plot is given in Fig. 6.It can be observed that NEC is satisfied throughout the evolution for Λ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 ≲ −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 ρ + p is not greater than zero always, hence the model 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, 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, 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, The evolution of the present model in the statefinder plane is shown in Fig. 7.For reference, note that, for Λ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 Λ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 Λ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 → −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.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], 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, where A H = 4πr 2 H , 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 = c H . Expressed in terms of h, the Eq.27 becomes, 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 ΛCDM model in the same figure.The figure indicates that, in the PEDE model, the entropy increases first, attains a maximum value at z ≈ −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.
This shows that since dh 2 /dx < 0 up to z ≈ −0.3, the rate of entropy, S ′ H > 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 is plotted against redshift (Fig. 9).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.

B. 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], 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 Λ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 < 0, where the derivative can be obtained with respect to the scale factor of the evolution of the universe.Differentiating Eq.29 gives, We found the evolution of S ′′ H which is shown in Fig. 10.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 Λ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 dimen-sionless dynamical variables, 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πG = c = 1, can be expressed as, which makes u + v = 1.From the conservation equations for non-interacting matter and dark energy (Eq.5), The autonomous coupled differential equations can be obtained from the above equations as, 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 → u ′ = u c + δu and v → v ′ = v c + δv, where δu and δ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, The 2 × 2 matrix on the right hand side is called the Jacobian matrix.The Jacobian matrix obtained for the PEDE model is given as, 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.For PEDE model, the eigenvalues corresponding to the critical point (1,0) are obtained as (3 , −3ω DE ).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ω DE , 3(1 + ω DE )).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.
A. Dark energy sound speed The dark energy sound speed is described by the ratio of pressure to density variations as [78], It has been argued that the physical values of v 2 s 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, ω DE ̸ = −1 [78].One can rewrite the equation for v 2 s as, where ṗDE = ω DE ρDE + ωDE ρ DE .Substituting with PE dark energy, the above equation reduces to, When z → −1, v 2 s reaches to a value ≈ −0.7.The Fig. 12 gives the complete evolution of v 2 s .It is evident from the figure that the value of v 2 s 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 2 s < 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), Ω m = H 2 0 Ω m0 /(a 3 H 2 ), Ω γ = H 2 0 Ω γ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 (ω DE < −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±0.618kms −1 Mpc −1 and Ω m0 = 0.280 ± 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 ≈ −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 ≈ 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 ≈ −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 nd law pertaining to the entropy of the universe.After a redshift of z ≈ −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 matterdominated 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 Λ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.
FIG. 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 ΛCDM and PEDE model parameters

FIG. 2 .FIG. 3 .
FIG. 2. The evolution of Hubble parameter with redshift in ΛCDM and PEDE models compared with H(z) data points

FIG. 5 .
FIG. 5.The evolution of q with z

FIG. 9 .
FIG. 9.The evolution of the first derivative of Horizon entropy

FIG. 10 .
FIG. 10.The evolution of the second derivative of Horizon entropy

2 FIG. 12 .
FIG.12.The evolution of squared sound speed with redshift

TABLE I .
The best estimated values and 68.3% confidence limit for the ΛCDM model parameters.

TABLE II .
The best estimated values and 68.3% confidence limit for the PEDE model parameters.