A thermodynamic point of view on dark energy models

We present a conjugate analysis of two different dark energy models investigating both their agreement with recent data and their thermodynamical properties. The successful match with the data allows to both constrain the model parameters and characterize their kinematical properties. As a novel step, we exploit the strong connection between gravity and thermodynamics to further check models viability by investigating their thermodynamical quantities. In particular, we study whether the cosmological scenario fulfills the generalized second law of thermodynamics and, moreover, we contrast the two model asking whether the evolution of the total entropy is in agreement with the expectation for a closed system. As a general result, we discuss whether thermodynamic constraints can be a valid complementary way to both constrain dark energy models and differentiate among rival scenarios.


Introduction
That the universe is presently undergoing a phase of accelerated expansion is nowadays an observationally unquestionable fact [1][2][3][4][5]. Understanding what is driving this speed up is one of most debated issue of modern cosmology with many contenders on the ground. An apparently simple solution to this fascinating problem may be found in the framework of Einsteinian General Relativity relying on the cosmological constant. Indeed, the ΛCDM model excellently reproduces the most recent observational data [6] assigning to the Λ term the origin of cosmic acceleration, while cold dark matter (CDM) is in charge for assembling the large scale structure. Unfortunately, the remarkable success in fitting the data comes at a non negligible price: the inferred Λ being 120 orders of magnitude larger than the predicted one with no idea of what cancellation mechanism can cure such a spectacular discordance. As a further issue, it is also disturbing that Λ and CDM contribute roughly equally to the today total energy budget which is quite unlikely considering that this only occurs once during the universe evolutionary history. Replacing Λ with an exotic fluid with a varying negative pressure, commonly referred to as dark energy (DE), still does not solve the coincidence problem (see, e.g., [7][8][9] for this and other related issues).
Lacking a solid theoretical basis, in the last decade, it has become popular to rely on phenomenological models [10]. On one hand, DE is modelled as a perfect barotropic fluid with pressure p and energy density ρ related by p = wρ with w the equation of state which becomes the uncertain feature assigning the model. Noting that w = −1 identifies ΛCDM model, the immediate next step is a Taylor expansion in the scale factor a thus leading to w = w 0 + w a (1 − a) which is known as the Chevallier -Polanski -Linder parameterization [26,27]. Such a simple formula is nevertheless able to mimic reasonably well the expansion rate of a large class of theoretically motivated DE models making it the reference choice for dark energy. Indeed, it is now common to quantify the capability of an observable to constrain DE in terms of the expected errors on the (w 0 , w a ) parameters. On the other hand, one can note that, in a Friedman -Robertson -Walker (FRW) framework, the expansion rate is directly determined by the energy density of the fluids composing the universe. As such, one can therefore directly model ρ instead of w with the further attracting possibility to work out models which can mimic both dark energy and dark matter with a single term. This is the case of the Hobbit model proposed by some of us [11] with an energy density scaling with a in such a way that the fluid evolves first as radiation, then as CDM and finally as a cosmological constant.
Despite the dynamical dark energy models capability to interpret cosmic expansion and the wide versatility of such models with respect to other experimental prescriptions, e.g. the perturbations "realm" [12], one should also be aware that the proposed scenarios are not in conflict with the prescriptions descending from fundamental physics. To this end, one can exploit the intriguing connection between gravity and thermodynamics which is well illustrated by the case of black holes (hereafter, BHs). In the semiclassical description of their physics, BHs behave as black bodies emitting thermal radiation and have both a temperature [14] and an entropy [15]. According to the Hawking -Bekenstein formula, the event horizon entropy is proportional to its area through the Planck and Boltzmann constants. Moreover, it has also been shown that BHs in contact with their own radiation satisfy the second law of thermodynamics. These results have been generalized in a cosmological context considering the universe as an isolated thermodynamic system filled with a cosmic fluid and bounded by some horizon [13,16]. Although this connection has been later reinforced [17,18], the deep relation between gravity, thermodynamics and quantum mechanics is far from being fully unveiled. Nevertheless, one can pursue this intriguing link as long as the universe can be considered a closed system with a horizon. In analogy with classical thermodynamics, one should expect that it naturally evolves towards a thermodynamical equilibrium. Therefore, its overall entropy function S(x), being x the relevant variable, should be a positively defined quantity that increases with convex concavity [19]. In practice, any physically viable cosmological model needs to fulfill the two conditions dS(x)/dx ≥ 0 and, in the far future, d 2 S(x)/dx 2 < 0. In particular, one can consider to evaluate the Universe effective entropy amount (including both the horizon contribution and the matter -energy term) and check that its first and second derivative have the correct behaviour. It could then be possible to discriminate among thermodynamically suitable cosmological models by considering the effective contribution of the different components to the entropy evolution [21][22][23][24][25].
In the standard approach to dark energy models, one is typically only concerned with its kinematic and dynamical properties so that the viability of a given parameterization is judged on the basis of how well it reproduces the observed data. We advocate here the need for a next step on. We consider two different dark energy models and constrain their parameters fitting a wide range of up to date observations probing the expansion rate. We then move on asking whether the models also comply with the second law of thermodynamics. Put in other words, we do not halt at the observations level, but push the analysis a step further into the realm of thermodynamics. As we will see, this allows to put stronger constraints on the model parameters and discriminate them at a more fundamental level.
Plan of the paper is as follows. In Section II, we resume the basics of the two cosmological models we consider, while Section III is devoted to matching them with observations. Section IV presents a general discussion on thermodynamical relations of both the horizon and the fluids making the cosmic pie, while Section V demonstrates the potentiality of a thermodynamical analysis applying it to the two models previously matched with the data.
A summary of the results and further prospectives are finally given in Section VI.

Dark energy models
In the framework of General Relativity and a homogenous and isotropic spatially flat universe described by the Robertson -Walker metric, the cosmic background evolution is determined by the Friedman equations, the first one being : where H =ȧ/a is the Hubble expansion rate, a the scale factor (with present day value a 0 = 1), the dot denotes derivative with respect to cosmic time, and the sum is over the fluids contributing to the total energy budget. We will hereafter use the labels (R, M, DE) to refer to the radiation, matter, and dark energy terms. The second Friedmann equation may be replaced by the continuity equation which sets the scale factor dependence of the density. Assuming no interaction between fluids, we can conveniently write where w i (a) is the equation of state (EoS) of the i -th term. Setting w R = 1/3 for radiation and w M = 0 for dust matter, we get with and Ω DE = 1 − Ω R − Ω M because of the flatness condition. The density parameters 1 will then read with f DE (a) = ρ DE (a)/ρ DE (a = 1) as obtained by integrating Eq.(2.2) given an expression for the DE EoS. For later convenience, it is worth noting that the second Friedmann and the continuity equations can be suitably rewritten as follows where prime denotes derivative with respect to a.

The BA parameterization
The ignorance about the nature of dark energy has motivated the search for analytical expressions which can mimic the DE EoS of a vast class of theoretical models. The most popular one is by far the one referred to as the CPL parametrization [26,27]: where w 0 is the present day value and w a sets the first derivative. Note that, while the early universe limit is well defined with w DE (a → 0) = w 0 + w a , the EoS diverges asymptotically in the future (a → ∞). Notwithstanding this problem, the CPL model is routinely used in the literature. However, we will be interested also in the future time evolution of the DE EoS so that we need a model which is as similar as possible to the CPL one, but avoid any singularity both in the early universe and the asymptotic future. This can be accomplished resorting to the parameterization proposed by Barbosa & Alcaniz [28] setting which we will refer to in the following as the BA model. It is immediate to show that the BA EoS asymptotes to w 0 + w a for a → 0 as for the CPL case, but it is not divergent in the future where the EoS goes to w 0 . Integrating Eq.(2.2), we then get : which sets the evolution of the DE for this model.

The Hobbit model
Rather than assuming an expression for the DE EoS, one can also use Eq.(2.2) in the reverse way, i.e. as a straightforward derivation of w DE (a) given an analytical expression for ρ DE (a).
Such an approach has been made popular by the search for unified DE models made out by a single fluid working as both dark matter and DE. As an example, we consider here the Hobbit model [11] with with (α, β) setting the scaling in the late and early universe, respectively, and (s, b) two scaling factors controlling the transition among the different regimes. Inserting (2.10) into (2.2) allows to compute w DE (a) = w Hob (a) which turns out to be It is only a matter of algebra to get the following relevant asymptotic behaviours : where, in the second row, we have used α, β ≥ 0 and s << b < 1. These expressions show that, in the early universe, the Hobbit EoS can mimic both radiation (β = 4) and matter (β = 3) thus eliminating the problem of early DE. At the present time, one can tailor (α, b) so that the EoS is close to the Λ value (w = −1) suggesting that the model is able to give an accelerated expansion. Finally, in the future, the Hobbit fluid reduces to the cosmological constant no matter the values of the model parameters. This nice behaviour motivated some of us to introduce the Hobbit model as a unified DE scenario. However, we will rather consider it as a phenomenological DE model thus also taking a dark matter contribution in Eq.(2.1). We let the data decide whether the matter can be reduced to the baryons contribution only.

Constraining model parameters
Both the BA and the Hobbit models are characterized by a set of parameters that can be settled in such a way to provide a cosmic scenario as close as possible to the observed one. On the one hand, the BA model may be naively matched to observation forcing (w 0 , w a ) = (−1, 0) so that the concordance ΛCDM scenario is recovered. On the other hand, this is not possible for the Hobbit model which has been contrasted to then available data in [11], but forcing (α, β) = (3,4) to reduce the number of parameters. In order to be more general, we therefore constrain both model parameters by fitting them to a wide dataset. We briefly describe below the data used and the Bayesian fitting procedure adopted and then present the results relevant for the later discussion.

Data and fitting method
As a first observable, we rely on the SNeIa Hubble diagram as traced by the JLA sample [5]. In order to save time without loss of precision, we use the compressed distance modulus sample thus defining the likelihood function as with N µ the number of data,

1)
∆D µ (p) = D obs µ − D th µ (p) the vector with the difference between the observed and predicted distance moduli and C µ the covariance matrix. We remind that the distance modulus is related to the model parameters as the luminosity distance and M SN eIa a nuisance parameter (accounting for difference between the adopted zero-point of the distance scale) which we marginalize over.
The SNeIa Hubble diagram is based on the use of a standard candle. On the contrary, Baryon Acoustic Oscillations (BAO) provide a standard ruler thus offering a further probe of the background expansion. We rely here on different datasets, the lowest redshift one being at z = 0.106 from the 6dFGS [29] team. The likelihood function is simply being z d the drag redshift (approximated as in [30]), Moving to larger redshift, we use the results of the SDSS DR7 team to measure d z (z) at z = (0.20, 0.35) [31]. Since these data are correlated, the likelihood function now is and ∆d z (p) = d obs z − d z (p). The observed values d obs z and the covariance matrix C DR7 are given in [31].
Both the BOSS and WiggleZ collaborations have measured BAO at still larger z observing galaxies over a partially overlapping survey area. In [32], the combined dataset has been used to measureD where r f id s is a fiducial value adopted in the analysis of the single samples. We include these measurements in our analysis defining the likelihood function being ∆D V (p) the usual difference vector. We refer the reader to [32] for the observedD V at z = (0.57, 0.73), the covariance matrix C BW , and the r f id s values. Note that all the BAO likelihood functions also depend on the baryon density parameter Ω b which we marginalize over.
The dataset described above probe the late universe so that it is interesting to add observables related to the early epochs. These are naturally provided by the CMB distance priors which summarize the information content of the Planck13 data [33]. The three quantities of interest are the acoustic scale with z LSS the last scattering surface redshift estimated as in [30], the shift parameter and the physical baryon density ω b = Ω b h 2 . The CMB distance priors likelihood function is then defined as and the covariance matrix C CM B given in [33].
All the data we have considered so far depend on the integrated Hubble expansion rate. It is therefore interesting to complement them with direct measurements of H(z) itself. We use the compilation in [34] which comes from an heterogenous sample of data obtained with different techniques. This ensures us that the covariance matrix is diagonal so that the corresponding likelihood has the typical Gaussian form Finally, to set the zero-point, we also add a Gaussian prior on h (i.e., H 0 in units of 100 km/s/Mpc) setting h obs = 0.706 ± 0.033 as in [35].
In order to perform a combined analysis, we maximize the full likelihood simply defined as the product of the single terms introduced above. To this end, we use a standard Markov Chain Monte Carlo (MCMC) method to sample the parameter space running six chains of equal length and keep adding points until convergence is achieved. The merged chain (after burn in cut and thinning) is then used to infer median and confidence ranges for each parameter.

Results
We fit both the BA and Hobbit models to the large dataset described in the previous paragraph. As can be seen from Fig. 1, both models work quite well at reproducing the SNeIa Hubble diagram and the H(z) data. In order to get these plots, we set the model parameters to their best fit values, namely for the BA model, and for the Hobbit case. Note that the radiation density parameter has not been allowed to vary in the fit, but we rather set with ω γ = 2.475 × 10 −5 and N ef f = 3.04 the effective number of massless neutrinos. While the best fit values for the BA model are expected since they are close to the ΛCDM limit (w 0 , w a ) = (−1, 0), the (α, β) values for the Hobbit model are somewhat surprising being quite different from the case (α, β) = (3, 4) which we have identified as particularly interesting. Actually, such a small value for α can be easily explained going back to the approximated formula for w Hob (a = 1). Inserting the best fit (α, b) values quoted above gives w Hob (a = 1) = −0.982 which is very close to the ΛCDM concordance scenario. The best fit β is then set in such a way that w Hob (a) asymptotes to w DE = −0.85 which is still close to the ΛCDM w = −1 value. Put in other words, both the BA and Hobbit best fit models try to mimic as close as possible the concordance ΛCDM scenario which is somewhat expected given the observational success of the cosmological constant.
As a further confirmation of this qualitative discussion, we stress that both models perform equally well to reproduce not only the SNeIa and H(z) data, but also the BAO clustering constraints and the Planck13 priors. As a consequence, the overall likelihood turns out to be almost equal for the two models although this argues in favor of the BA one which has a smaller number of parameters. It is worth noting, however, that the two models have a distinct evolutionary history with w BA (a) remaining almost constant all over the universe expansion due to the very small w a value. On the contrary, w Hob (a) presents a non negligible evolution moving from w Hob = −0.98 today to w Hob = −0.85 in the radiation era.
Although we do not need them in the following, we also report below the constraints (median and 68% confidence values) on the model parameters. For BA we find For the Hobbit model, we have used η b = log (1/b) and η s = log (1/s) as parameters instead of (b, s) since they allow to explore a wider range. We then find We note that both the matter density parameter Ω M and the present day Hubble constant h are well consistent between the two models. While the BA parameters are severely constrained to lie close to their ΛCDM counterparts, some more discussion is needed to interpret the results for the Hobbit model. Indeed, the median (α, β) values are quite different from the best fit ones which are outside the 68% CL (but within the not reported 95% ones). Such an unusual behaviour is a consequence of two effects. First, strong degeneracies among the model parameters can shift the best fit value from the median ones. Indeed, the best fit parameters are chosen to maximize the likelihood, while the median ones characterize the probability distribution function of each parameter after taking into account the fit to the data. The two quantities are different so that it is not unusual to find discrepant values when a large parameter space is constrained by a limited amount of data. Second, most of the data probe the late universe with only Planck priors tracing the early one and no data at all covering the matter dominated region where the Hobbit EoS can change from cosmological constant -like to dust -like. This prevents to put significant constraints on both β and η b thus propagating the uncertainty to other parameters too. Motivated by these considerations, we do not consider troublesome the apparent discrepancy among best fit and median (α, β) values. However, since in the following we will be mainly interested in showing how thermodynamics can discriminate between models fitting equally well the same dataset, we will set both the BA and Hobbit model parameters to their best fit values.

Thermodynamic analysis
It is customary to consider a DE model as well motivated if it passes the matching with observational data with green lights. However, the realm of observationally viable models is always increasing. Lacking an underlying conclusive theory, it is worth addressing the problem from a different point of view. Thermodynamics offers such a distinct look. Indeed, we aim here at investigating whether different models fitting equally well the same data can be distinguished on the basis of their thermodynamical properties. Put in other words, we wonder whether thermodynamics can alleviate the degeneracy surviving the data fitting test.
In a cosmological framework, we actually deal with the Generalized Second Law (GSL) of thermodynamics since, for system with some horizon, the entropy of the horizon itself must be added to those of the fluids within it. In a sense, this arises in a natural way from the consideration that the horizon prevents the observer to see what lies beyond it.
As already seen elsewhere [21], this approach is theoretically well defined in cosmology and translates, within the FRW framework, in the requirements S ′ (a) ≥ 0 and, in the far future, as a consequence of the generalized second law of thermodynamics, S ′′ (a) ≤ 0, where the prime denotes the derivative with respect to a. In such a case, an apparent horizon is naturally defined as well as the relative entropy, which depends on the surface, in analogy with BHs; on the other hand standard thermodynamical definitions of entropy can be applied to the matter -energy components. The appropriate boundary in this context has been shown to be the apparent horizon of FRW universes, which always exists, something that is not generally true for the particle horizon and the future event horizon [13,20].
Let us first compute the different contribution to the total entropy starting from the horizon one. To this end, we take the apparent horizon, i.e., the marginally -trapped surface with vanishing expansion, since -by contrast to other possible choices -the laws of thermodynamics are fulfilled on it [13]. Writing the line element as the apparent horizon is defined by the condition h ab ∂ ar ∂ br = 0, beingr = a(t)r and h ab = diag(−c 2 , a(t) 2 ). Its radius turn out to ber A = c/H(a) so that the area is naively A = 4πr 2 A . The entropy is simply proportional to A so that we get with k B the Boltzmann constant and ℓ P l the Planck length. It is worth noting that S H is defined in terms of the Hubble expansion rate H(a) and hence is determined by the properties (abundance and EoS) of the fluids contributing to the energy budget. Actually, in order to verify thermodynamical prescriptions, rather than the entropy itself, we are interested in its derivatives with respect to the scale factor a.
Restoring physical constants and using the Friedmann and continuity equations, we finally getS where the sum is over the fluids in the universe and hereafter we will use the dimensionless entropies defined as Making use of the Eq.(2.6), we then get the following result for the second derivative of the apparent horizon entropỹ The approach to the thermodynamical equilibrium of the horizon (S ′′ H = 0) is then driven by the EoS of the fluid components, determining their abundance evolution, and its derivative.
To compute the entropy of each fluids within the horizon, we must first determine how their temperature evolves. This readily follows from the Gibbs equation and the condition that dS i is a differential. We then get We can now rely on the Euler equation T s = ρ(1 + w) with s the entropy density to first compute the entropy of the i -th fluid and then its derivative with respect to a. Going to dimensionless quantity, we find A similar derivation can not be used for the matter term since its temperature is identically null all along the cosmic history. We can instead use S M = k B V H n, where V H is the volume of the apparent horizon and n the number density of dust particles. Differentiating with respect to a finally givesS where we have defined with n 0 ∼ 10 −6 m −3 the present day dust number density 2 . We stress that, although with the same dimension, T 0M is not the matter temperature, but only a convenient quantity introduced to giveS ′ M (a) the same form as the radiation and DE terms. The total entropy is the sum of the horizon entropy and those relative to each fluid component. Total first derivative of entropy,S ′ tot , can be conveniently written as . (4.14) having used w R (a) = 1/3. It is instructive to estimate the order of magnitude of the terms entering Eqs.(4.12) -(4.14). To this end, we insert the constants with the correct units in the above definitions to get It is therefore evident that, unless in the very early universe (i.e., a → 0), the matter related terms in Eqs.(4.12) -(4.14) are always negligible thus explaining why we do not care about the exact n 0 value. Finally, it is immediate to get the second derivative of the entropy with respect to the scale factor a. Differentiating Eq.(4.11), we finally get where we do not report the derivatives of the different terms for sake of shortness. Eqs.(4.11) and (4.15) show that the variation of the total entropy depend on the evolution of the density parameter and the EoS of all the fluids contributing to the full energy budget. Actually, the two by far dominant contributions are the horizon and DE ones with the latter playing a crucial role in the fulfillment of the GSL. Indeed, being the only fluid with negative EoS, it is also the only one which can driveS ′ tot (a) towards negative values. It is the balance between the DE and horizon terms which will finally determine whether the entropy decrease or increase during the universe history.

Thermodynamic constraints
As yet said before, given the strong connection between gravity and thermodynamics [15], we suppose that the Universe behaves as an isolated thermodynamic systems so that it approaches to a state of maximum entropy in the long run. As a consequence, in order for a cosmological model to be viable from a thermodynamical point of view, the GSL must be fulfilled so that two following conditions must be passedS Since the entropy can never decrease, Eq.(5.1) must hold over all the universe history from the beginning to its end, i.e. for 0 ≤ a ≤ a max , with a max a very large value. However, both DE models we are considering have been proposed as phenomenological descriptions of this component which we expect to be valid over a limited range in redshift. We will therefore impose Eq.(5.1) over the range (a min , a max ) with a min = (1 + z LSS ) −1 ∼ 0.001 and a max = (1 + z min ) −1 ≃ 1000 where z LSS ∼ 1090 and z min = −0.999 are the redshifts of the last scattering surface redshift and of distant point in the future close enough to the asymptotic limit a → ∞.

BA model
Let us first consider the phenomenological BA parameterization for the DE EoS. Plugging in the relevant formulae for the density parameters and temperature scaling, we can investigate whether the thermodynamical constraints are fulfilled or not. To this end, we first set the model parameters (Ω M , w 0 , w a , h) to their best fit values. In such a case we get the results in Fig. 2 where we plot 3S′ tot (a) for different values of τ DE . 3 Since the amplitude ofS ′ tot (a) changes by orders of magnitude with τDE, here and in the following plots, we scale the derivatives with respect to an arbitrary chosen value which we adjust with τDE so that all the curves may be shown in the same plot. We therefore warn the reader to only look at the shape ofS ′ tot (a) and not to directly compare the values. As a remarkable result, we find thatS ′ tot (a) is not positive definite, but may rather become negative for a limited period of time. Put in other words, the functionS ′ tot (a) has at least one zero whose value depends on the present day DE temperature. Such a behaviour can be qualitatively explained as follows. As universe expands, the derivative of the DE entropy decreases until becoming negative, but the positive contribution of the horizon term becomes more and more important thus delaying the violation of the constraint (5.1). Eventually, the horizon term becomes the only one, but now the universe content is dominated by the negative pressure fluid so that the horizon entropy crosses the zero line leading toS ′ tot (a) < 0. The τ DE parameter then controls the transition between the two regimes, while the cosmological parameters (Ω M , w 0 , w a ) determine the sign of the term in square parentheses in Eq.(4.2) which controls the sign ofS ′ H (a). Note that, in the limit (w 0 , w a ) = (−1, 0), the horizon entropy derivative asymptotes to the null value since Ω DE (a >> 1) ≃ 1 and w DE = −1. Such a condition is independent on the value of τ DE so that it is not surprising that all the curves in Fig. 2 asymptotically approach an almost null value given that the fiducial (w 0 , w a ) parameters are indeed close to the critical ones.
One could wonder whether such a result is specific of the fiducial model or can it be cured by changing the model parameters. In Fig. 3, we plotS ′ tot (a) changing one parameter at time and leaving the other fixed to their fiducial values for different τ DE . For all cases we consider, there is always a range whereS ′ tot (a) becomes negative so that the entropy is decreasing in that period which is unphysical. Actually, such a result is consistent with the thermodynamical analysis of DE models presented in e.g., [21]. As we go back in time, the DE becomes subdominant (i.e., there is no early dark energy), but its contribution to the entropy derivative is still non negligible (because of the factor T 0R /T 0M which makes the matter one negligible). Since w BA (a → 0) = w 0 + w a , for fixed w a , the larger is w 0 /w f id 0 , the more the DE fluid behaves as a phantom one driving the dominantS ′ H towards negative values because of the term Ω DE w DE in the square parenthesis. This explains whyS ′ tot (a) has the largest variation with w 0 shown in the central row of Fig. 3, while the similar scaling one would expect for the dependence on w a is suppressed by a fiducial value which is close to zero. Note that, no matter which is the value of τ DE , models with small Ω M have a negative derivative of the entropy over almost the full evolutionary history approaching then null value  in the far future because of the EoS asymptotic to the cosmological constant value. This leads us to argue against them, but they are actually already excluded by the comparison with the data. A somewhat strange trend takes place for w 0 /w 0,f id = 0.5 giving a divergent entropy. It is not easy to trace the origin of this unrealistic behaviour, but we nevertheless note that such values are outside the confidence range allowed by the data. Fig. 4 shows the asymptotic limit of the entropy second derivative as a function of Ω M for given values of the EoS parameters (w 0 , w a ). These results might have been anticipated following the argument in [21] where the authors have shown that models with an eternal acceleration comply the laws of thermodynamics. Indeed, we find thatS ′′ (a max ) is negative for the fiducial BA model since the DE EoS stays approximately constant thanks to the close to zero w a value so that the universe tends to a de Sitter -like state. Moreover, different from the CPL model, the BA EoS has no pathological behaviour in the far future where, on the contrary, the EoS reduces to constant w 0 . Since we have found w 0 ∼ −1, we get a ΛCDM model in the future thus leading to a constant entropy henceS ′′ ∼ 0 which is what Fig. 4 shows considering the values ofS ′′ ref listed in the caption.

Hobbit model
Let us now consider the Hobbit model taking as fiducial model parameters the best fit ones. We obtain theS ′ tot (a) curves in Fig. 5 for different τ DE values. A comparison with the same plot for the fiducial BA model shows a striking difference with the Hobbit proposal, this model being able to fulfill the constrain (5.1) over the full evolutionary history as long as log τ DE > −3.0 eventually approaching a null value in the asymptotic future. Such a nice behaviour can be qualitatively explained looking at the terms entering Eq.(4.11). The only negative is approximately constant and negative as for the BA model, the balance between the negative and positive terms is regulated by the evolution of the energy densities of the different components. On the contrary, for the Hobbit model, both w DE (a) and the coefficient functions W i (a) change in such a way that the net result is a positive derivative. A negative overall value can only be achieved when this balance is not possible anymore which is the case when τ DE is so small to make the function W 1 (a) much larger than the other factors. Moving to the future, the universe becomes dominated by the Hobbit fluid whose EoS is now approaching the asymptotic limit w DE = −1. The universe is therefore driven by a constant Λ term which does not evolve anymore thus giving a constant total entropy so thatS ′ tot (a) vanishes. This qualitative picture is also consistent with what we see in Fig. 6 where we let the model parameter change one at time leaving the others fixed to their fiducial values and exploring three different τ DE values. First, we note that, no matter the model parameter values,S ′ tot (a) stays always positive provided log τ DE > −3.0 so that we can safely conclude that the Hobbit model fulfills the first thermodynamic constraint provided the present day DE temperature is not too smaller than the horizon one T H = H 0 /k B . This is a remarkable result since it suggests that the Hobbit model could be a better physically motivated prescription for the DE evolution with respect to the phenomenological BA model.
Second, we stress that the dependence ofS ′ tot (a) on the model parameters is consistent with the qualitative scenario described above. Indeed, models with smaller Ω M gives larger S ′ tot (a) values since the smaller is Ω M , the larger is the Hobbit fluid contribution both in the past and in the present and near future thus explaining why we find larger (in absolute value) peaks. On the other hand, (α, b) control the width of the transition from positive to negative EoS with smaller values leading to models with a smoother w DE and hence a smallerS ′ tot (a). On the contrary, the (β, s) parameters have only a marginal impact since they determine the shape of the EoS in the very far past (a << 1) so that their impact is hardly appreciable.
It is instructive wondering when the entropy derivative violates the constraint (5.1). Put in other words, we want to determine the zeros (a − , a + ) ofS ′ tot (a) since it isS ′ tot (a) ≤ 0 for a − ≤ a ≤ a + . These are shown in Fig. 7 as function of log τ DE for different model parameters combinations. Note that the curves abruptly stop since for τ DE > τ crit , the functionS ′ tot (a) have no zeros being always non negative and asymptoting to zero only for a → ∞. As can be read from Fig. 7, the critical value depends on the model parameters (Ω M , α, b), while we have checked that (β, s, h) have no impact at all. In particular, we note that the larger is Ω M , the smaller is τ crit in accordance with the qualitative scenario described above. We also note that the larger is Ω M , the narrower is the range (a − , a + ), i.e. the smaller is ∆a ± = a + − a − , so that the better is behaviour of the model from a thermodynamical point of view. Similar considerations argue in favour of models with smaller (α, b) values although they have a minor impact on both τ crit and ∆a ± .
Finally, let us look at the asymptotic value of the second order derivative of the total entropy to see whether the Hobbit model fulfills the constraint (5.2). Fig. 8 showsS ′′ tot (a max ) as a reliable (numerically stable) approximation of this quantity arbitrarily setting τ DE = 1 (but same results are obtained for different values). No matter which are the model parameters values, it turns out thatS ′′ tot (a max ) is always negative and quite close to the null value being of order 10 −7 . This is consistent with the picture of a universe which is dominated in the future by a DE term with w DE ≃ −1, i.e. the Hobbit model reduces to the ΛCDM one in the asymptotic future. In such a regime, the only term contributing to the entropy is the Λ -like one which does not evolve at all. The entropy is therefore constant as a consequence of the universe ending in a thermodynamical equilibrium. Model parameters (Ω M , α, b) control how fast this asymptotic state is achieved.

Conclusions
Recent years have witnessed an impressive progress in observational cosmology with the accumulation of wide and high quality datasets probing both the kinematics and the dynamics of the universe. Such a remarkable progress has nevertheless not been rewarded by a similar advance in the understanding of what is driving the cosmic speed up. Given this frustrating situation, it is worth wondering whether the elusive nature of dark energy can be investigated from a different point of view. In an attempt to add further arrows to hit this rewarding target, we have here tackled dark energy adopting a thermodynamical point of view. We have considered the universe as a closed system and, therefore, we have developed a thermodynamical analysis in analogy with black holes thermodynamics. As a matter of facts, each cosmological model must not only match the observed data, but also fulfill the generalized second law (GSL) of thermodynamics and the approach to equilibrium. As a consequence, the overall universe entropy must never decrease and reach a maximum in the long run so that the universe asymptotically attains a state of thermodynamical equilibrium. Mathematically, it is the same as asking that the constraints given by Eqs.(5.1) and (5.2) are fulfilled. We have therefore first evaluated first and second entropy derivatives with respect to the scale factor a for two different DE models, namely the phenomenological BA parameterization (similar to the popular CPL one) and the empirically motivated Hobbit scenario. In practice, after having shown that both are in agreement with present day data (which also allows us to constrain their parameters), we have then checked whether they pass the GSL constraints with green light.
It turns out that this is actually possible for both models, but a critical role is played by the unknown value of the present day DE temperature. As a general result, we find that, for the BA model, the first derivative of the total entropyS ′ tot (a) changes sign at least one time, i.e. the equationS ′ tot (a) = 0 has at least one solution. Depending on τ DE and the model parameters (Ω M , w 0 , w a ) values,S ′ tot (a) can then become positive again or stays negative converging to a close to null limit in the far future. It is worth noting that as far as log τ DE > 0,S ′ tot (a) stays positive over the full evolutionary history of the universe provided model parameters are suitably chosen yet remaining within the 68% CL allowed by the data. On the other hand, the Hobbit model can be reconciled with the GSL quite easily as far as the condition log τ DE > −2.5 is imposed on the present day DE temperature. We indeed find that the functionS ′ tot (a) has no zeros in this case, while the second derivative of the entropy always asymptotes to an almost null value no matter which values are set for τ DE and the model parameters (Ω M , α, β, b, s). This nice thermodynamical behaviour originates from the characteristics of the Hobbit EoS which makes the fluid mimic the cosmological constant Λ in the far future (so that the entropy stays constant and the universe is in a state of thermodynamical equilibrium).
It is somewhat ironic that the thermodynamical analysis points at the Hobbit model as a better physically motivated scenario, but it is the BA EoS which comes out as the preferred one from the fitting procedure because of the lower number of parameters (being the likelihood almost the same for the two models). Such a contradictory result is actually an outcome of how the two models have been conceived. The BA parameterization aims at generalizing the ΛCDM model and it has been designed so that it can efficiently mimic a wide class of DE models including the (unknown) one underlying the data themselves. It is therefore not surprising that it is able to fit the data especially considering that the best fit parameters are indeed quite close to those of the concordance ΛCDM scenario. On the contrary, the Hobbit model has been manufactured by hand so that it mimic standard fluids in the past and the cosmological constant in the recent and future epochs. Since each one of these fluids correctly behaves from a thermodynamical point of view, we expect that the Hobbit one does the same except in the transition regime from matter to the Λ term.
Although the above argument explains why the two independent indicators point in different directions, it is clearly desirable to have both tests select the same model. To this end, one could importance sample the MCMC chains of a given model (not necessarily those we have investigated here) so that the final selected parameters identify configurations that both match the data and fulfill the thermodynamic constraints. Two difficulties are in order here. First, both Eqs.(5.1) and (5.2) are inequalities so that they can not be implemented as a likelihood function to be maximized. On the contrary, for each set of parameters along the chain, one must compute the entropy derivatives and check their sign all over the range (a min , a max ). Only the configurations passing the GSL constraints are then saved thus significantly reducing the parameter space narrowing down the confidence ranges. However, such an approach is straightforward if all the quantities entering Eqs.(4.11) and (4.15) are analytical otherwise one should computeS ′ tot (a) andS ′′ tot (a) numerically which can be quite time consuming (especially if one needs a long chain to achieve MCMC convergence). Actually, this is only a numerical problem which can likely be solved by a suitable algorithm or running the code on fast machines. The most serious hurdle is indeed represented by the ignorance on the DE temperature, i.e., the value of the τ DE parameter. As we have shown for both the BA and Hobbit models, for a given set of model parameters, the number of zeros of the functioñ S ′ tot (a) depend on τ DE . Put in other words, each set along the MCMC chain can be accepted or rejected according to τ DE being smaller or larger than a critical value. As a mandatory step before implementing a joint data -GSL analysis, one must find a way to constrain the present day DE temperature.
Although some significant work is needed for a more efficient implementation, the results in this preliminary analysis show that thermodynamics can offer an additional yet competitive way to compare different DE models. The conjecture that universe evolutions can be depicted from the thermodynamical point of view in analogy with black holes allows to settle two natural constraints that have to be satisfied along the cosmological history. In our opinion, adding more and more data is a rewarding strategy to constrain model parameters but adopting suitable conceptual paradigms can allow to further improve the differentiation process among rival competing approaches.