Constructions of Variable Modified Chaplygin-Jacobi and Chaplygin-Abel Gases from Geometrothermodynamics

Introduction

Observations of Type Ia supernovae (SNe Ia) [1, 2, 3] provided the first compelling evidence that the expansion of the Universe is accelerating. This conclusion has subsequently received support from several independent cosmological probes, including Baryon Acoustic Oscillations (BAO) [4, 5], Cosmic Microwave Background (CMB) measurements [6], and observations from the Wilkinson Microwave Anisotropy Probe (WMAP) [7, 8] and the Planck mission [9]. To account for the observed accelerated expansion of the Universe, cosmological models introduce an exotic component known as dark energy (DE), characterized by a sufficiently negative pressure that produces an effective repulsive gravitational effect. Among the various proposed explanations for DE, the cosmological constant $\Lambda$, commonly interpreted as vacuum energy, remains the simplest framework and forms the basis of the standard $\Lambda$CDM cosmological model. Its contribution is estimated to constitute approximately 70% of the total energy density of the present-day Universe [10, 11, 12]. Although the concordance $\Lambda$CDM model provides an excellent description of a wide range of cosmological observations, the theoretical nature of the cosmological constant remains problematic, particularly due to the fine-tuning and cosmic coincidence issues [13]. Motivated by the theoretical challenges associated with the cosmological constant, a variety of dynamical dark energy scenarios have been proposed in which the energy density evolves with cosmic time [14, 15, 16]. Among these possibilities, quintessence, generally described by a dynamical scalar field, has emerged as one of the widely investigated alternatives to the cosmological constant [17, 18]. In addition to quintessence, several other scalar-field-based frameworks have been explored as possible explanations for the observed late-time acceleration, including different scalar field and tachyonic dark energy models [19, 20, 21]. An alternative approach to explaining the late-time accelerated expansion is to modify the underlying gravitational theory rather than introduce a separate dark energy component. In such frameworks, modifications to the gravitational dynamics can effectively generate the acceleration of the Universe and thereby provide possible alternatives to dark energy [22, 23, 24, 25, 26].

The Chaplygin gas provides another phenomenological framework for describing the accelerated expansion without introducing dark energy as a separate component. This exotic fluid can reproduce the late-time acceleration of the Universe and, in its unified formulation, offers a description in which the dark matter and dark energy sectors emerge from a single cosmic fluid [27, 28, 29]. The Chaplygin gas framework was subsequently extended to the generalized Chaplygin gas (GCG), which retains the unified description of dark matter and dark energy within a single effective fluid [30, 31]. Further generalization led to the modified Chaplygin gas (MCG), whose cosmological evolution can interpolate between a radiation-dominated regime at early times and a $\Lambda$CDM-like behavior at late times [32, 33]. Several further generalizations of the Chaplygin gas framework have been investigated, including the variable Chaplygin gas [34], variable modified Chaplygin gas [35], and viscous Chaplygin gas [36]. These extended formulations provide alternative phenomenological descriptions of dark energy and, in suitable settings, retain the possibility of interpreting dark matter and dark energy as manifestations of a single effective cosmic fluid. Observational data have been extensively employed to place constraints on the parameters governing Chaplygin gas models. In particular, measurements from Type Ia supernovae, the Cosmic Microwave Background (CMB), and WMAP have been used to constrain the parameters of the generalized Chaplygin gas (GCG) [31, 37]. Similar observational analyses based on CMB and WMAP data have also provided constraints on the parameters of the modified Chaplygin gas (MCG) [38]. Furthermore, the behavior and allowed parameter ranges of the MCG have been investigated within various modified theories of gravity [39, 40, 41, 42, 43]. Villanueva [44] formulated a generalized Chaplygin-Jacobi gas by constructing the Hubble parameter from a generalized Chaplygin scalar-field framework and employing Jacobi elliptic functions. Within this formulation, the corresponding inflationary observables and quantities were also derived. Inspired by this approach, we extend the analysis by introducing the equations of state for the modified Chaplygin-Jacobi gas and the modified Chaplygin-Abel gas.

The geometric description of thermodynamic systems through metric structures has been widely studied. Different thermodynamic metrics are used to describe phase transitions through curvature singularities. The Fisher-Rao metric [45], arising from information theory, is expressed as the Hessian of entropy, while Ruppeiner’s metric, based on thermodynamic fluctuation theory, is given by the negative Hessian of entropy [46]. However, these metrics are not invariant under changes of thermodynamic representation or different thermodynamic potentials. For instance, in the energy representation $U=U(E^{a})$, the metric is conformally related to the Weinhold metric, defined as the Hessian of energy [47], with inverse temperature as the conformal factor. The Geometrothermodynamics (GTD) program provides a more general approach based on the invariance of its metrics under Legendre transformations. GTD has been applied to several systems, including ideal gases, van der Waals fluids [48, 49], and black holes. The invariance under changes of representation has also been investigated through the conformal structure of GTD metrics [50]. In this work, we study the geometric properties of thermodynamic systems from a fully invariant perspective, considering invariance under both Legendre transformations and changes of representation. In this framework, dark energy (DE) is modeled as a geometrothermodynamic fluid (GTDF) with a negative equation of state, capable of counteracting gravitational attraction. This DE contribution becomes significant only at late cosmic times, emerging naturally after a transition redshift. The fluid’s pressure is proportional to its volume, so its negative pressure intensifies as the universe expands—an effect rooted in a constant thermodynamic interaction. Within geometrothermodynamics, such interaction is encoded in the constant curvature of the equilibrium manifold. This approach yields an evolving equation of state that asymptotically approaches a cosmological constant, offering a self-consistent late-time dynamics [51]. However, a comprehensive investigation of Chaplygin-type realizations within the GTD framework is still lacking. Furthermore, most existing studies of Chaplygin-gas cosmology consider a limited number of fixed equations of state and are tested against only a restricted set of observational probes. Consequently, it remains important to examine how the inferred cosmological parameters and statistical performance depend on the specific functional form of the model and the choice of observational data, particularly in view of the differences among recent Type Ia supernova compilations such as DES-SN5YR, and Union3.

In this work, we address these issues by investigating two variable Chaplygin-type models within the GTD framework: the Variable Modified Chaplygin-Jacobi Gas (VMCJG) and the Variable Modified Chaplygin-Abel Gas (VMCAG). By deriving their fundamental thermodynamic relations and examining their cosmological evolution, we aim to answer the following key questions:

1. Can the VMCJG and VMCAG models provide a viable description of the late-time expansion history while remaining competitive with $\Lambda$CDM?

2. How do the two parametrizations differ in their predicted dark-energy dynamics, particularly in terms of the equation-of-state evolution, statefinder diagnostics, and higher-order kinematic parameters?

3. To what extent does the choice of supernova compilation, when combined with BAO data, affect the inferred values of $H_{0}$ and $\Omega_{m0}$, and can the additional flexibility of the Chaplygin-gas equation of state accommodate or reduce such compilation-dependent differences relative to $\Lambda$CDM?

To address these questions, we perform a Markov Chain Monte Carlo likelihood analysis using four complementary data sets: CC+BAO, DESI DR2, DESI DR2+DES-SN5YR, and DESI DR2+Union3. For each combination, we compute the total $\chi^2$ and derive the posterior distributions of the model parameters. To assess the relative viability against $\Lambda$CDM, we evaluate the AIC, BIC, and DIC information criteria, which penalise model complexity while rewarding goodness of fit. Later, using the best-fitting values of the models’ free parameters, we study the evolutions of cosmological and cosmographical parameters and check the classical stability of the models. The structure of this paper is as follows. In § 0.1, we introduce the basic concepts of geometrothermodynamics (GTD) relevant to the present work. § 0.2 focuses on the construction of the entropy for the Modified Chaplygin-Jacobi Gas and Modified Chaplygin-Abel Gas within the GTD formalism, whereas § 0.3 extends the analysis to their variable counterparts, namely the Variable Modified Chaplygin-Jacobi Gas and Variable Modified Chaplygin-Abel Gas. The fundamental field equations describing the cosmological dynamics are presented in § 0.5. In § 1, we describe the methodology adopted for the observational analysis, including the MCMC likelihood approach and the data sets CC+BAO, DESI DR2, DESI DR2+DES-SN5YR, and DESI DR2+Union3. The observational results are presented in § 1.3, which is devoted to the complete cosmological analysis, including the estimation of the model parameters and their comparison using the AIC, BIC, and DIC criteria, while § 1.4 is devoted to the complete cosmological analysis. Finally, our conclusions are given in § 1.5.

Review of Geometrothermodynamics

In 2007, H. Quevedo [52] has developed geometrothermodynamics (GTD) formalism to describe the properties of thermodynamic systems by using differential geometry. In classical equilibrium thermodynamics [53], the thermodynamic system consists of $n$ extensive thermodynamic variables $E^a$ ($a=1,\ldots,n$), their corresponding dual intensive thermodynamic variables $I^a$ and the thermodynamic potential $\Phi$. So define the phase space as the $(2n+1)$-dimensional differential manifold $\mathcal{T}$ with a set of coordinates $\{Z^A\}=\{\Phi,E^a,I^a\}$ where $A=1,2,\ldots,(2n+1)$. This $n$ determines the number of thermodynamic degrees of freedom. The most general metric for GTD is

where the fundamental Gibbs-1 form is

and $\Lambda\in\mathbb{R}$ is constant. The metric $G$ is invariant under the Legendre transformations.

The $n$-dimensional equilibrium submanifold $\mathcal{E}\subset\mathcal{T}$ can be defined by the smooth map $\varphi:\mathcal{E}\to\mathcal{T}$ by $\varphi:\{E^a\}\to\{\Phi(E^a),E^a,I^a(E^a)\}$ so that the condition $\varphi^*(\Theta)=0$ (i.e., $d\Phi=I_a dE^a$ or, $I_a=\frac{\partial\Phi}{\partial E^a}$) can be satisfied where $\varphi^*$ is the pullback of $\varphi$.

The manifold $\mathcal{E}$ is naturally equipped with the induced thermodynamic metric $g$, given by

The equilibrium manifold $\mathcal{E}$ determines an extremal surface in $\mathcal{T}$, so using variational principle, the variation of the volume element is identically zero, i.e., $\delta\int_{\mathcal{E}}\sqrt{\det(g)}\,d^nE=0$. Since the induced metric $g$ depends on $Z^A$, so this variation derives the Nambu-Goto system of differential equations

where $\Box$ denotes the d’Alembert operator and $\Gamma$ denotes the Christoffel symbol for the metric $G_{AB}$. From the above variation, one can obtain a set of differential equations for the thermodynamic potential $\Phi$ whose solutions can be expressed as functions of the extensive variables $\Phi=\Phi(E^a)$ which is known as fundamental equations.

Constructions of Entropy for Modified Chaplygin-Jacobi Gas and Modified Chaplygin-Abel Gas using Geometrothermodynamics

We consider the entropy $S=S(U,V)$ as the thermodynamical potential i.e., $\Phi=S$ and $U,V$ as the extensive variables i.e., $E^a=\{U,V\}$. So $Z^A=\{S,U,V,\frac{1}{T},\frac{P}{T}\}$. The fundamental Gibbs-1 form is

For equilibrium state, the condition $\varphi_S^*(\Theta_S)=0$ reduces to the first law of thermodynamics

such that

In 2015, Villanueva [44] has formulated the generalized Chaplygin-Jacobi gas using Jacobi’s elliptic function. Subsequently, using Jacobi’s elliptic function, we have derived modified Chaplygin-Jacobi (MCJG) gas equation of state

where $A(>0)$, $B(>0)$, $\alpha\in(0,1)$ are constants and $k\in[0,1]$ is the elliptic modulus. Also using Abel’s elliptic functions, we have formulated modified Chaplygin-Abel gas (MCAG) equation of state

where $c,e\in\mathbb{R}$. The above MCJG and MCAG can be expressed in the form

where $A_1=(2k-1)(1+A)-1$, $B_1=kB$, $C_1=\frac{(1-k)(1+A)^2}{B}$ for MCJG and $A_1=(e^2+2c^2)(1+A)-1$, $B_1=c^2B$, $C_1=-\frac{(e^2+c^2)(1+A)^2}{B}$ for MCAG. Since energy density $\rho=U/V$, so the equation (10) reduces to

So equation (6) transforms to

Since left hand side is the exact differential, so right hand side should be the exact differential. Also the coefficients of $dU$ and $dV$ are homogeneous functions of $U$ and $V$. So $1/T$ is the integrating factor and that is equal to $1/(MU+NV)$. So the expression of temperature is

Putting the expression of $T$ in (12) and integrating, we obtain the expression of entropy as in the following form

where $S_0$ is integration constant and $D_1=\frac{\sqrt{(1+A_1)^2+4B_1^2C_1}}{4B_1^2}$.

Variable Modified Chaplygin-Jacobi Gas and Variable Modified Chaplygin-Abel Gas from Geometrothermodynamics

Now consider the form of entropy

where $c_i$ ($i=1,2,3,4$) are constants. The induced metric (3) in the space of equilibrium states is given by

where

For this induced metric, we can obtain the non-vanishing scalar curvature. So the internal thermodynamic interaction is possible for the Legendre invariant scalar curvature.

From equation (7), we obtain

and

From equations (20) and (21), we obtain

If $\beta=\alpha$, the above EoS is similar with the MCJG and MCAG equation of states (8) and (9). If we assume $\frac{c_1}{c_2}=(1+\alpha)\sqrt{(2+A)^2-4k}$ and $c_3,c_4=-\frac{1}{2kB}\left[(2k-1)(1+A)-1\pm\sqrt{(2+A)^2-4k}\right]$, then the above expression (22) reduces to

which is termed as “Variable Modified Chaplygin-Jacobi Gas" (VMCJG). Also if we assume $\frac{c_1}{c_2}=(1+\alpha)\sqrt{(e^2(1+A)^2-1)^2-2c_2(1+A)}$ and $c_3,c_4=-\frac{1}{2c^2B}\left[(e^2+2c^2)(1+A)-1\pm\sqrt{(e^2(1+A)^2-1)^2-2c^2(1+A)}\right]$, then the above expression (22) reduces to

which is termed as “Variable Modified Chaplygin-Abel Gas" (VMCAG). In the present geothermodynamic formulation, we assume the thermodynamic volume to be proportional to the physical volume of the FLRW universe, $V\propto a^{-3}$. This choice allows the thermodynamic quantities to be consistently expressed in terms of the cosmological scale factor and its evolution.

More General Equation of State from Geometrothermodynamics

We can extend the form of entropy (15), so the more general form of entropy can be considered as

where $c_i$'s and $d_i$'s are constants. From equation (7), we obtain

and

which turn on the general equation of state (GEoS)

This expression reduces to the VMCJG and VMCAG if $n=2$ and $c_2=-c_1$. For $n=1$, we can get VMCG.

FRW Model

We assume a spatially flat, homogeneous, and isotropic universe described by the Friedmann–Lemaître–Robertson–Walker (FLRW) metric [54],

where $a(t)$ denotes the scale factor. The spatial curvature is set to $k=0$, in agreement with the spatially flat geometry supported by current cosmological observations.

The evolution of the cosmic background is described by the Einstein field equations,

where $R_{\mu\nu}$ and $R$ represent the Ricci tensor and Ricci scalar, respectively, while $T_{\mu\nu}$ denotes the total energy–momentum tensor. We adopt natural units, $8\pi G=c=1$, throughout the analysis.

The cosmic medium is assumed to consist of two non-interacting components: pressureless matter, including dark matter and baryons, with density $\rho_m$ and pressure $p_m=0$, and a dark-energy component characterized by density $\rho_{DE}$ and pressure $p_{DE}$. The total energy–momentum tensor can therefore be written as

where $u^\mu$ is the four-velocity of the fluid satisfying $u^\mu u_\mu=-1$, with $i\in\{m,DE\}$.

For the FLRW background [55], the Einstein equations lead to the Friedmann equation

and the acceleration equation

where $H$ is the Hubble parameter and an overdot denotes differentiation with respect to cosmic time. Equivalently, the evolution of the Hubble parameter may be expressed as

For non-interacting matter and dark-energy components, the conservation equations follow from the Bianchi identity, $\nabla^{\mu}G_{\mu\nu}=0$. The matter component satisfies

where $\rho_{m0}$ is the present-day matter density. The dark-energy component, on the other hand, obeys

The latter equation requires a specific relation between $p_{DE}$ and $\rho_{DE}$ to determine the evolution of the dark-energy density. In the following section, we introduce the modified Chaplygin-type equation of state adopted in this work and derive the corresponding cosmological evolution.

Methodology

Bayesian Framework

The free parameters of the cosmological model are constrained within a Bayesian statistical framework [56, 57]. For a parameter vector $\boldsymbol{\theta}$ and a given dataset $D$, Bayes’ theorem can be expressed as

where $\mathcal{P}(\boldsymbol{\theta}\mid\mathbf{d})$ denotes the posterior probability distribution, $\mathcal{L}(\mathbf{d}\mid\boldsymbol{\theta})$ represents the likelihood of obtaining the observed data for a given set of model parameters, and $\pi(\boldsymbol{\theta})$ is the prior distribution. Assuming Gaussian-distributed observational uncertainties, the likelihood can be written in terms of the corresponding log-likelihood as

where $\Delta\mathbf{X}=\mathbf{X}_{\mathrm{obs}}-\mathbf{X}_{\mathrm{th}}(\boldsymbol{\theta})$ denotes the residual vector between the observed quantities and their corresponding theoretical predictions. Here, $\mathbf{C}$ represents the covariance matrix of the observational data. For datasets for which only independent statistical uncertainties are available, $\mathbf{C}$ reduces to a diagonal matrix with elements $\sigma_i^2$. We impose flat

(uniform) priors on all free parameters, with the allowed ranges chosen according to physically reasonable considerations and summarized in Table I. For statistically independent datasets, the combined likelihood is obtained by adding their individual log-likelihoods:

Model	Parameter	Prior range
VMCJG	$H_0$	$\mathcal{U}(66,73)$
VMCJG	$\Omega_{m0}$	$\mathcal{U}(0.1,0.5)$
VMCJG	$\alpha$	$\mathcal{U}(0,0.6)$
VMCJG	$\beta$	$\mathcal{U}(0.01,1)$
VMCAG	$H_0$	$\mathcal{U}(60,72)$
VMCAG	$\Omega_{m0}$	$\mathcal{U}(0.2,0.35)$
VMCAG	$\alpha$	$\mathcal{U}(0.2,1)$
VMCAG	$\beta$	$\mathcal{U}(0.4,2.6)$

thereby enabling the joint constraints to incorporate complementary information from different cosmological observations.

Markov Chain Monte Carlo Sampling

To sample the posterior distribution described by Eq. (37), we employ the affine-invariant ensemble Markov Chain Monte Carlo (MCMC) algorithm implemented in the emcee package [58]. Unlike conventional single-chain Metropolis–Hastings approaches, the affine-invariant sampler utilizes an ensemble of $N_{\mathrm{walkers}}$ walkers that explore the parameter space simultaneously. The proposal for each walker is constructed using the positions of the remaining walkers, making the method relatively insensitive to the scale and correlations among the model parameters and reducing the need for manually tuned proposal distributions [59, 60].

For each observational dataset, the walkers are initialized in a compact region surrounding a physically plausible starting point, with their initial positions drawn from a restricted portion of the adopted prior volume. Each ensemble is then evolved for $N_{\mathrm{iter}}$ iterations. To remove the influence of the initial conditions, the first $N_{\mathrm{burn}}$ iterations are discarded as the burn-in phase. The remaining samples are subsequently thinned by a factor of $N_{\mathrm{thin}}$, where appropriate, to reduce the effects of autocorrelation. The resulting samples are then flattened to obtain the final posterior distribution used for parameter estimation.

The marginalized constraints on individual parameters and their correlations are represented through one- and two-dimensional posterior distributions. The quoted parameter values correspond to the median of the marginalized posterior distributions, while the associated uncertainties are obtained from the 16th and 84th percentiles, representing the 68% credible intervals. The posterior distributions are visualized using triangle plots generated at the $1\sigma$ and $2\sigma$ confidence levels, corresponding to the 68% and 95% credible regions, respectively. GetDist [61].

For each MCMC analysis, the best-fit parameter vector, $\boldsymbol{\theta}_{\mathrm{best}}$, is identified with the maximum a posteriori (MAP) sample, i.e., the sample having the largest posterior probability within the flattened chain. We adopt the MAP point rather than constructing a parameter vector from the individual marginalized medians because, in the presence of strong parameter degeneracies and correlations, the latter combination need not correspond to an actual point within the high-probability region of the multidimensional posterior distribution.

Data Description

To constrain the model and examine its consistency with observations, we consider several complementary late-time cosmological probes. The analyses are performed both for individual datasets and for suitable combinations of independent datasets.

Cosmic Chronometers and Baryon Acoustic Oscillations (CC+BAO)

We consider a compilation of Hubble parameter measurements, $H(z_i) \pm\sigma_{H,i}$, obtained primarily through the cosmic chronometer (CC) differential-age technique, supplemented by BAO-based measurements [62, 63, 64]. The data cover a broad redshift interval extending approximately over $0 \lesssim z \lesssim2$. For these measurements, the covariance matrix is taken to be diagonal, and the corresponding chi-square statistic is given by $\chi^2_{\mathrm{CC+BAO}} = \sum_i [H_{\mathrm{obs}}(z_i) - H_{\mathrm{th}}(z_i;\theta)]^2/\sigma_{H,i}^2$. The associated likelihood is then obtained from Eq. (38).

DESI DR2 Baryon Acoustic Oscillations

We further employ the Data Release 2 (DR2) BAO measurements from the Dark Energy Spectroscopic Instrument (DESI). The dataset consists of 12 BAO measurements involving the quantities $D_M/r_d$, $D_H/r_d$, and/or $D_V/r_d$, distributed over seven effective redshift bins in the range $0.295 \leq z \leq2.33$. The correlations among the measurements are accounted for through the corresponding Gaussian covariance matrix [65, 66, 67, 68].

The relevant distance measures are defined as

and $D_V = [zD_M^2D_H]^{1/3}$ where $D_M$ denotes the comoving angular-diameter distance. The sound horizon at the drag epoch, $\tau_d$, is fixed to the fiducial value $r_d = 148.3\,\mathrm{Mpc}$.

Type Ia Supernovae

Two independent Type Ia supernova (SN Ia) compilations are considered in our analysis [69, 70, 71]:

• DES-SN5YR: The Dark Energy Survey five-year sample contains 1829 SNe Ia covering a comparable redshift range. Since the complete systematic covariance matrix is not available in our adopted implementation, we employ a diagonal covariance matrix constructed from the reported statistical uncertainties, $\sigma_{\mu,i}$. This treatment is less comprehensive than that adopted for Pantheon+SH0ES and may therefore lead to a modest underestimation of the parameter uncertainties [72, 73, 74].

• Union3: We also consider the Union3 compilation, represented by 22 redshift-binned SN Ia distance-modulus measurements spanning approximately 0.05 $\lesssim z \lesssim$ 2.26. The corresponding 22 $\times$ 22 covariance matrix is incorporated in the analysis, and the absolute-magnitude nuisance parameter is marginalized following the standard procedure [75, 76].

For all three SN Ia datasets, the theoretical distance modulus is calculated according to $\mu(z) = 5\log_{10}[D_L(z)/\mathrm{Mpc}]+25$, with the luminosity distance $D_L(z) = (1+z)D_M(z)$.

Model Comparison Statistics

To quantify the statistical performance of the proposed model and compare it with the standard $\Lambda$CDM scenario for different observational combinations, we employ the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Deviance Information Criterion (DIC) [77, 78, 79]. These quantities are defined as

where $\mathcal{L}_{\max}$ denotes the maximum likelihood, $k$ is the number of free parameters in the model, and $N$ represents the total number of observational data points.

For the DIC, the deviance is defined as $D(\theta) = -2\ln\mathcal{L}(\theta)$, $\overline{D(\theta)}$ is its posterior mean over the MCMC chain, and $p_D = \overline{D(\theta)} - D(\overline{\theta})$ is the effective number of parameters, penalizing model complexity, with $D(\overline{\theta})$ representing the posterior mean of the deviance, which accounts for the complexity of the model.

For all three information criteria, smaller values indicate a statistically preferred model. To facilitate a relative comparison, we define $\Delta\mathrm{AIC}$, $\Delta\mathrm{BIC}$, and $\Delta\mathrm{DIC}$ with respect to the model having the minimum value of the corresponding criterion. In general, differences $\Delta\mathrm{AIC}$ or $\Delta\mathrm{BIC} < 2$ indicate no substantial statistical distinction between the competing models, values in the range 2–6 provide positive evidence against the model with the larger value, while differences exceeding 10 are conventionally interpreted as strong evidence in favor of the model with the lower criterion value.

Observational Analysis

Both reconstructed $H(z)$ curves accurately track the CC+BAO data across the full redshift range ($z < 2.4$), remaining well within the observational uncertainties at all points, as seen from Figure 1. The VMCJG and VMCAG best-fit curves are visually indistinguishable from one another, indicating that the two parametrizations yield nearly identical background expansion histories despite their differing equations of state. This near-degeneracy at the level of $H(z)$ suggests that CC+BAO data alone cannot discriminate between the two models, and that additional probes such as Type Ia supernovae and DESI BAO distance measurements are needed to break this degeneracy. From the getdist plots in Figure 2, both models constrain $H_0$ to $\sim69$–$72\ \mathrm{km\,s^{-1}\,Mpc^{-1}}$ and $\Omega_{m0}$ to $\sim0.24$–$0.30$, consistent with the standard $\Lambda$CDM scenario and lying between the Planck and SH0ES constraints. The posterior distributions also show the expected inverse trend between $H_0$ and $\Omega_{m0}$. For the Chaplygin-Jacobi model, $\alpha$ and $\beta$ are relatively well constrained, with $\alpha\sim0$–$0.45$ and $\beta\sim0.30$–$0.45$. In contrast, the $\alpha$-$\beta$ contours of the Chaplygin-Abel model are strongly elongated along a diagonal direction, indicating a pronounced parameter degeneracy. This feature is also reflected in the bimodal one-dimensional posterior distributions of $\alpha$ and $\beta$. Overall, the Chaplygin-Jacobi parametrization provides comparatively tighter and more well-defined constraints on the model parameters, indicating a better constrained parameter space than the Chaplygin-Abel case.

The best-fit parameters, presented in Table. Table II reveal a clear distinction between the Chaplygin models and $\Lambda$CDM. While the $\Lambda$CDM value of $H_0$ decreases with the inclusion of supernova data, from 69.87 for CC+BAO to 66.94 $\mathrm{km\,s^{-1}\,Mpc^{-1}}$ for DESI DR2+DES-SN5YR, both Chaplygin models consistently favor higher values of $H_0 \sim70$–$72\ \mathrm{km\,s^{-1}\,Mpc^{-1}}$. The difference becomes particularly pronounced for the SN-inclusive datasets, suggesting that the additional freedom in the Chaplygin equations of state allows a higher expansion rate closer to local measurements. A similar trend is found for $\Omega_{m0}$, with both Chaplygin models yielding lower and more stable matter densities than $\Lambda$CDM; the Jacobi model gives the lowest values, $\Omega_{m0} \sim0.237$-0.264. The shape parameters also show a consistent hierarchy, with the Abel model favoring larger $\alpha$ and $\beta$ values ($\alpha\sim0.46$-0.86, $\beta\sim0.87$-1.40) than the Jacobi model ($\alpha\sim0.09$-0.28, $\beta\sim0.35$-0.40). Both parameters increase notably when DES-SN5YR is included, indicating that supernova data play a key role in driving the departure from $\Lambda$CDM, particularly for the Jacobi parametrization, which is also strongly supported by the information criterion analysis.

Model	Dataset	$H_0$ (km s$^{-1}$ Mpc$^{-1}$)	$\Omega_{m0}$	$\alpha$	$\beta$
$\Lambda$CDM	CC+BAO	$69.867^{+0.826}_{-0.903}$	$0.268^{+0.015}_{-0.018}$	–	–
$\Lambda$CDM	DESI DR2	$68.551^{+0.451}_{-0.493}$	$0.296^{+0.008}_{-0.009}$	–	–
$\Lambda$CDM	DESI DR2 + DES-SN5YR	$66.943^{+0.382}_{-0.491}$	$0.329^{+0.010}_{-0.005}$	–	–
$\Lambda$CDM	DESI DR2 + Union3	$68.041^{+0.467}_{-0.555}$	$0.305^{+0.011}_{-0.008}$	–	–
VMCAG	CC+BAO	$70.509^{+0.706}_{-1.253}$	$0.299^{+0.039}_{-0.049}$	$0.487^{+0.313}_{-0.174}$	$0.927^{+0.407}_{-0.215}$
VMCAG	DESI DR2	$71.055^{+0.600}_{-0.480}$	$0.292^{+0.031}_{-0.053}$	$0.457^{+0.294}_{-0.197}$	$0.865^{+0.377}_{-0.259}$
VMCAG	DESI DR2 + DES-SN5YR	$71.790^{+0.154}_{-0.286}$	$0.284^{+0.026}_{-0.028}$	$0.864^{+0.099}_{-0.199}$	$1.377^{+0.136}_{-0.260}$
VMCAG	DESI DR2 + Union3	$71.343^{+0.379}_{-0.633}$	$0.288^{+0.041}_{-0.036}$	$0.553^{+0.272}_{-0.211}$	$0.999^{+0.342}_{-0.294}$
VMCJG	CC+BAO	$69.675^{+1.765}_{-2.176}$	$0.238^{+0.021}_{-0.026}$	$0.090^{+0.130}_{-0.045}$	$0.354^{+0.028}_{-0.027}$
VMCJG	DESI DR2	$69.634^{+1.732}_{-1.132}$	$0.264^{+0.014}_{-0.014}$	$0.128^{+0.105}_{-0.088}$	$0.378^{+0.021}_{-0.025}$
VMCJG	DESI DR2 + DES-SN5YR	$71.396^{+0.636}_{-0.566}$	$0.249^{+0.013}_{-0.012}$	$0.281^{+0.035}_{-0.052}$	$0.399^{+0.012}_{-0.010}$
VMCJG	DESI DR2 + Union3	$70.553^{+0.717}_{-1.121}$	$0.259^{+0.014}_{-0.012}$	$0.200^{+0.065}_{-0.110}$	$0.391^{+0.014}_{-0.023}$

Dataset	Model	$N$	Reduced $\chi^2$	AIC	BIC	DIC
CC+BAO	VMCAG	57	0.573	38.39	46.56	34.91
CC+BAO	VMCJG	57	0.552	37.24	45.42	36.00
CC+BAO	$\Lambda$CDM	57	0.597	36.18	40.27	$35.56^\dagger$
DESI DR2	VMCAG	12	0.800	14.40	16.34	10.81
DESI DR2	VMCJG	12	0.720	13.76	15.70	11.88
DESI DR2	$\Lambda$CDM	12	1.141	14.27	15.24	13.47
DESI DR2 + DES-SN5YR	VMCAG	1841	0.933	1721.87	1743.94	1720.98
DESI DR2 + DES-SN5YR	VMCJG	1841	0.921	1700.03	1722.11	1700.34
DESI DR2 + DES-SN5YR	$\Lambda$CDM	1841	0.933	1718.53	1729.57	1718.29
DESI DR2 + Union3	VMCAG	34	1.050	39.49	45.60	36.92
DESI DR2 + Union3	VMCJG	34	0.946	36.39	42.49	36.25
DESI DR2 + Union3	$\Lambda$CDM	34	1.246	42.63	45.68	42.09

Table IV reports $\Delta$AIC, $\Delta$BIC, and $\Delta$DIC, defined as $\mathrm{IC}_{\mathrm{model}}-\mathrm{IC}_{\Lambda\mathrm{CDM}}$; negative values therefore indicate that the given Chaplygin-type model achieves a lower, more favorable information criterion score than $\Lambda$CDM and is statistically preferred by that criterion, while positive values indicate $\Lambda$CDM is preferred. For CC+BAO, $\Lambda$CDM is mildly favored over both Chaplygin models by the AIC and BIC values, aside from a small preference for the Abel model by the DIC ($\Delta\mathrm{DIC}=-0.65$). For DESI DR2 alone, no criterion reaches even positive evidence for $\Lambda$CDM ($+\nu$); the DIC in fact mildly favors both Chaplygin models, though these fall short of strong evidence. The clearest result emerges for DESI DR2 combined with DES-SN5YR: the VMCJG model is very strongly preferred over $\Lambda$CDM, whereas the VMCAG model is disfavored relative to $\Lambda$CDM by all three criteria. For DESI DR2 combined with Union3, both Chaplygin models are preferred over $\Lambda$CDM: the VMCJG model shows strong evidence via the AIC ($\Delta\mathrm{AIC}=-6.24$) and positive evidence via the BIC and DIC, while the VMCAG model shows positive evidence via the AIC and DIC but no significant preference via the BIC ($\Delta\mathrm{BIC}=-0.08$). Taken together, the VMCJG model is consistently preferred over, or statistically indistinguishable from, $\Lambda$CDM across all four dataset combinations, with a decisive preference once supernova data are included, whereas the VMCAG model is comparable to or mildly preferred over $\Lambda$CDM for CC+BAO, DESI DR2, and Union3, but is disfavored specifically for the DES-SN5YR combination.

Dataset	Model	$\Delta$AIC	$\Delta$BIC	$\Delta$DIC
CC+BAO	VMCAG	+2.21	+6.29	−0.65
CC+BAO	VMCJG	+1.06	+5.15	+0.44
CC+BAO	$\Lambda$CDM	0.00	0.00	0.00
DESI DR2	VMCAG	+0.13	+1.10	−2.66
DESI DR2	VMCJG	−0.51	+0.46	−1.59
DESI DR2	$\Lambda$CDM	0.00	0.00	0.00
DESI DR2 + DES-SN5YR	VMCAG	+3.34	+14.37	+2.69
DESI DR2 + DES-SN5YR	VMCJG	−18.50	−7.46	−17.95
DESI DR2 + DES-SN5YR	$\Lambda$CDM	0.00	0.00	0.00
DESI DR2 + Union3	VMCAG	−3.14	−0.08	−5.17
DESI DR2 + Union3	VMCJG	−6.24	−3.19	−5.84
DESI DR2 + Union3	$\Lambda$CDM	0.00	0.00	0.00

Cosmographical Parameters

In this section, we’ll study the nature of the cosmological and cosmographical parameters by taking the best-fit values of the parameters of the VMCJG and VMCAG models [80].

• EoS parameter:

The propagation of VMCJG and VMCAG for the evolution of the universe can be determined by the equation of state parameter $W$ [81].

The two models exhibit a significant difference in the redshift evolution of the dark-energy equation-of-state parameter $w_{DE}(z)$, as plotted in Figure 3. Both models cross the phantom divide, $w=-1$, around $z\sim0.3$-$0.7$, while their present-day values remain near the quintessence regime. Here, $W<-1/3$ shows the dark energy phase of the universe. Finally, $W$ reaches $-1$ as $z\sim-1$, i.e., the universe goes to the $\Lambda$CDM stages due to the contributions of VMCJG and VMCAG. However, their high-redshift behaviors differ markedly: the VMCJG model evolves toward a deep phantom regime, reaching $w_{DE}\sim-2$ to $-2.6$ at $z=2$, whereas the VMCAG model evolves toward $w_{DE}\approx0$, with slightly positive values for some datasets. Thus, the Abel parametrization exhibits a more physically conventional early-time behavior, while the strong high-redshift phantom evolution of the Jacobi model may be difficult to reconcile with the expected matter-dominated early universe. This divergence at $z\gtrsim1$ provides an important discriminating feature that can be tested with future high-redshift observations.

• $W$-$W'$ plane:

The $w$–$w'$ phase-plane trajectories in Fig. 4 reveal distinctly different dynamical behaviors depending on the diagnostic quantity used.

In the first case, all four dataset combinations approach the $\Lambda$CDM fixed point $(-1, 0)$ and then move into the phantom region ($w_{\mathrm{DE}} < -1$), following a pronounced freezing-type path ($w'_{\mathrm{DE}} < 0$) down to a minimum of about $-1.5$ to $-1.65$ near $w_{\mathrm{DE}} \sim-1.7$ to $-2.0$, before turning back toward less negative $w'$ as $w$ becomes strongly phantom ($\lesssim-2.5$) [82, 83]. Among these, the DESI DR2 trajectory shows the most extreme excursion, while the CC+BAO curve stays closest to $\Lambda$CDM, and the DESI DR2 + DES-SN5YR and DESI DR2 + Union3 paths nearly overlap throughout.

In contrast, the second case shows all trajectories confined to $w' > 0$, indicating pure thawing behavior with no phantom crossing. Here $w_{\mathrm{DE}}$ ranges only from about $-0.9$ to $+0.1$, staying well away from the $\Lambda$CDM point. The curves converge to a common minimum of $w' \approx0.14$ near $w \approx-0.6$, then diverge again, with the DESI DR2 + DES-SN5YR combination reaching the highest $w'$ values at both ends.

Overall, the two panels demonstrate that whether the reconstructed dark-energy dynamics allows phantom crossing and freezing evolution or remains confined to standard thawing quintessence depends critically on the reconstruction quantity used. The choice of dataset primarily influences the amplitude of the features, not their qualitative nature.

fig:4: Evolution of $\omega- \omega'$ analysis for the (a) Variable Modified Chaplygin-Jacobi Gas (b) Variable Modified Chaplygin-Abel Gas models for the considered observational datasets.

• Deceleration parameter:

fig:5: Evolution of the deceleration parameter $q(z)$ for the (a) Variable Modified Chaplygin-Jacobi Gas (b) Variable Modified Chaplygin-Abel Gas models for the considered observational datasets.

In order to discuss the deceleration or acceleration phase of the universe, we need to study the deceleration parameter, which is given as $q = -\frac{\ddot{a}}{aH^2}$. For the deceleration phase of the universe, $q > 0$, and for the acceleration phase of the universe, $q < 0$ [84, 85, 86].

From Fig. 5, we can see that both models show the expected transition from a decelerating matter-dominated phase at high redshift ($q > 0$) to late-time acceleration ($q < 0$), with a broadly consistent transition redshift of $z_t \sim0.6$-1.0. The main difference lies in the strength of acceleration: the VMCAG model reaches more negative values, with $q \approx-1.1$ for the SN-inclusive datasets, briefly entering the phantom-like regime ($q < -1$), while the VMCJG model remains moderate, with $q \gtrsim-0.85$. At high redshift, both models approach $q \sim0.3$-$0.6$, close to the matter-dominated value $q = 0.5$. The DESI DR2-only Jacobi fit shows a non-monotonic feature with a double zero-crossing around $z \sim0.3$-$0.9$, likely reflecting the weaker constraints on $(\alpha,\beta)$ from BAO data alone, consistent with the broader posterior contours discussed earlier.

• Cosmographic parameters:

There are another parameters named as cosmographical parameters for understanding the past and future evolution of the universe [87, 88, 89]. The cosmography is the study of the scale factor by expanding it by Taylor series with respect to the cosmic time. The cosmographical parameter like jerk ($J$), snap ($S$) parameters are [90]

So the cosmographical parameters are higher-order derivatives of the deceleration parameter $q$.

Figure 6 demonstrates that the jerk parameter $j(z)$ provides an additional distinction between the two parametrizations. The VMCJG model exhibits an oscillatory evolution around the $\Lambda$CDM value $j = 1$, remaining mostly below unity at low redshift, reaching $j \sim1.3$-$1.7$ near $z \sim1.2$-$1.5$, and approaching $j \sim1.25$-$1.3$ by $z = 2$. The DESI DR2-only case shows a pronounced dip near $j \approx0.5$, reflecting the weaker constraints from BAO data alone. In contrast, the VMCJG model follows a smoother evolution: after a minimum of $j \sim0.2$-$0.3$ near $z \sim0$-$0.2$, it rises above $j = 1$ around $z \sim1$-$1.3$ and reaches $j \sim3$-$5.4$ at $z = 2$. Hence, although the VMCAG model exhibits a more conventional $j(z)$ evolution, it shows a substantially stronger high-redshift deviation from $\Lambda$CDM in higher-order kinematics, while the VMCJG model displays a more oscillatory and dataset-dependent behavior.

The snap parameter $s(z)$ further differentiates the two parametrizations at fourth order, as evaluated from Figure 7. The VMCJG model shows a strongly non-monotonic evolution, with a pronounced minimum of $s \sim-1.8$ to $-8.5$ around $z \sim0.7$-$1.0$, followed by partial recovery toward $z = 2$. Additional extrema in the DESI DR2-only reflect the dataset-dependent features seen in the lower-order diagnostics. In contrast, the VMCAG model remains relatively smooth, staying near zero at $z \lesssim1$, but then declines sharply to very large negative values, reaching $s \sim-9$ to $-83$ by $z = 2$. Thus, the VMCJG model exhibits intermediate-redshift oscillatory behavior, whereas the Abel model develops a much stronger high-redshift divergence, indicating a substantially different departure from $\Lambda$CDM despite its more conventional $w_{DE}(z)$ evolution.

• Statefinder Diagnostic: The statefinder diagnostic provides a useful geometrical tool for distinguishing different dark-energy models through higher-order kinematic quantities [91, 92, 93]. It is defined by the dimensionless parameters $(r,S)$, where the statefinder pair is given by

where $q$ is the deceleration parameter and an overdot denotes differentiation w.r.t. cosmic time. The $\Lambda$CDM model corresponds to the fixed point $(r,s)=(1,0)$, making the statefinder plane particularly useful for quantifying deviations from the concordance model.

From Figure 8, in both $r$-$S$ plane diagrams, all four dataset combinations pass through or near the $\Lambda$CDM fixed point $(S,r)=(0,1)$, confirming consistency with the standard model at some epoch. From there, the curves bend into the upper-left region $(S<0, r>1)$, associated with Chaplygin gas-like behavior, then loop back through the fixed point and continue into the lower-right quadrant $(S>0, r<1)$, characteristic of quintessence. This indicates that each trajectory samples both regimes over its evolution—early phases resemble Chaplygin gas, while later stages shift toward quintessence.

The DESI DR2 and DESI DR2+Union3 curves show the most pronounced curvature, peaking at $r\approx1.55$-1.6, whereas the DESI DR2+DES-SN5YR combination remains more damped, with a peak around $r\approx1.35$-1.47. All trajectories also develop a small loop or kink well into the quintessence region $(S\approx0.1$-0.4, $r\approx0.3$-0.8), indicating rapid, non-monotonic changes in the deceleration parameter’s derivatives at late times; the CC+BAO curve extends furthest in $S$ $(\approx0.4)$ with the least pronounced loop.

Despite these differences in amplitude and loop extent, the qualitative agreement across all datasets supports a dynamical dark energy history that transitions between Chaplygin-gas like and quintessence-like phases near the present epoch, rather than a static $\Lambda$CDM cosmology.

• Stability:

To study the stability against small perturbations, we can write the squared speed of sound in the form $c_s^2=\frac{dp}{d\rho}$. If $c_s^2>0$, we can say that the model is classically stable while $c_s^2<0$ represents a classically unstable.

The squared sound speed $c_s^2(z)$ provides a crucial test of the perturbative stability of the two models, as seen from Figure 9. The VMCJG model is generally stable at low redshift, with $c_s^2 > 0$, but develops a negative sound speed at higher redshift, reaching $c_s^2 \sim-2$ to $-2.5$ by $z = 2$ for most datasets. In contrast, the VMCAG model becomes unstable in the low-redshift regime, with $c_s^2 \sim-0.3$ to $-0.4$ over approximately $z \sim-0.3$ to $0.6$–$1.0$, before recovering at higher redshift. Both models also exhibit superluminal values ($c_s^2 > 1$) at the extremes of the redshift range. Thus, the Jacobi parametrization appears more viable over the observationally constrained low-redshift regime, whereas the Abel model faces a stability issue at late times. The superluminal behavior in both cases remains an important theoretical caveat.

Discussions and Conclusions

In this work, we explore the potential applications of thermodynamic geometry to obtain new physically relevant insights into cosmology. In particular, the central idea is that fluids possessing distinctive geometric properties may also exhibit unique cosmological behavior.

We constrained the free parameters of both models using four complementary observational data sets: CC+BAO, DESI DR2, DESI DR2+DES5Y, and DESI DR2+Union3. Through a Markov Chain Monte Carlo likelihood analysis, we obtained posterior distributions for the model parameters. The reconstructed parameter values reveal a clear physical distinction between the two Chaplygin-gas realizations and their relation to the standard $\Lambda$CDM baseline. Across all four dataset combinations, both the VMCAG and VMCJG models recover a higher Hubble constant $H_0$ than $\Lambda$CDM, with the discrepancy being most pronounced for the DESI DR2 + DES-SN5YR combination. In this case, $\Lambda$CDM attains its lowest inferred value, $H_0 = 66.943^{+0.382}_{-0.491}$, whereas both Chaplygin models remain at $H_0 \gtrsim71$. This behavior indicates that the additional degrees of freedom in the Chaplygin-gas equation of state absorb part of the low-$H_0$/high-$\Omega_{\mathrm{m}0}$ pull that the DES-SN5YR compilation exerts on $\Lambda$CDM, thereby partially decoupling the inferred Hubble constant from the well-known $\Omega_{\mathrm{m}0}$-$H_0$ anti-correlation driven by supernova calibration systematics. Consistent with this interpretation, the matter density parameter $\Omega_{\mathrm{m}0}$ in the VMCJG scenario is systematically lower than the corresponding $\Lambda$CDM value for every dataset, with differences of approximately 0.03-0.08. This suggests that a fraction of the total energy budget attributed to pressureless matter in $\Lambda$CDM is instead being reinterpreted as an emergent unified dark-fluid contribution within this model.

The information criteria analysis reveals a data-dependent preference for the Chaplygin models. For DESI DR2 alone, all three models are nearly degenerate. However, for DESI DR2 + DES-SN5YR, the VMCJG model is strongly favoured over $\Lambda$CDM ($\Delta\mathrm{AIC} = -18.50$, $\Delta\mathrm{DIC} = -17.95$), while the VMCAG model is disfavoured ($\Delta\mathrm{AIC} = +3.34$). This indicates that the Jacobi parametrization effectively absorbs the low-$H_0$/high-$\Omega_{\mathrm{m}0}$ tension introduced by the supernova data. For DESI DR2 + Union3, both Chaplygin models are favored, with the Jacobi model again performing best ($\Delta\mathrm{AIC} = -6.24$). Overall, the VMCJG model emerges as the most competitive alternative to $\Lambda$CDM, though the preference is dataset-dependent and warrants confirmation with future surveys.

The cosmological analysis further distinguishes the two parametrizations. The equation-of-state parameter $w_{\mathrm{DE}}(z)$ shows that both models cross the phantom divide $w = -1$ around $z \sim0.3$–$0.7$, with present-day values in the quintessence regime. However, their high-redshift behaviors differ markedly: the VMCJG model evolves towards a deep phantom regime reaching $w_{\mathrm{DE}} \sim-2$ to $-2.6$ at $z = 2$, while the VMCJG behavior model evolves towards $w_{\mathrm{DE}} \approx0$. The $w$-$w^{\prime}$ phase-plane analysis confirms this distinction, with the VMCJG model exhibiting freezing-type behaviour ($w^{\prime} < 0$) with phantom crossing, whereas the VMCAG model remains confined to pure thawing quintessence ($w^{\prime} > 0$) with no phantom crossing. The deceleration parameter $q(z)$ shows that both models exhibit the expected transition from deceleration to acceleration at $z_t \sim0.6$-1.0. The VMCAG model reaches more negative values ($q \approx-1.1$), briefly entering the phantom-like regime, while the Jacobi model remains moderate ($q \gtrsim-0.85$). The jerk $j(z)$ and snap $s(z)$ parameters further differentiate the two models: the VMCJG model exhibits oscillatory behaviour around $\Lambda$CDM with moderate high-redshift deviations, while the VMCAG model shows smoother evolution at low redshift but develops substantially stronger high-redshift divergences, reaching $j \sim3$-5.4 and $s \sim-9$ to $-83$ at $z = 2$. The statefinder $(r,s)$ diagnostic confirms these trends, with the VMCJG model remaining closer to the $\Lambda$CDM fixed point while the VMCAG model exhibits larger departures. The statefinder diagnostic reveals that both models trace out trajectories that pass through the $\Lambda$CDM fixed point $(s,r) = (0,1)$ and sample both Chaplygin gas-like ($s < 0$, $r > 1$) and quintessence-like ($s > 0$, $r < 1$) regimes over their evolution, indicating a dynamical transition rather than a static cosmological constant. The DESI DR2 and DESI DR2+Union3 combinations exhibit the most pronounced curvature, peaking at $r \sim1.55$-1.6, while the DESI DR2+DES-SN5YR combination remains more damped. All trajectories develop small loops or kinks in the quintessence region, signaling rapid non-monotonic changes in the deceleration parameter’s derivatives at late times.

The perturbative stability analysis, however, reveals important caveats: the Jacobi model is stable at low redshift but becomes unstable at high redshift ($c_s^2 < 0$ for $z \gtrsim2$), while the Abel model exhibits low-redshift instabilities ($c_s^2 \sim-0.3$ to $-0.4$ for $z \sim-0.3$ to 0.6-1.0) before recovering at higher redshift. Both models also show superluminal behavior ($c_s^2 > 1$) at extreme redshifts. These findings suggest that while the background dynamics of both Chaplygin parametrizations remain competitive with $\Lambda$CDM, their perturbative stability requires further investigation before they can be considered fully viable alternatives.

In conclusion, we have demonstrated that both Chaplygin-gas parametrizations derived from geometrothermodynamics offer viable dynamical alternatives to $\Lambda$CDM, with the VMCJG model emerging as the most competitive candidate, particularly for the DESI DR2 + DES-SN5YR dataset, and that the extra equation-of-state freedom in both unified dark-fluid descriptions absorbs much of the compilation dependent shift in $H_0$ and $\Omega_{m0}$ seen in the standard model — pointing to unified dark-fluid geometrothermodynamics as a physically motivated and observationally robust alternative to $\Lambda$CDM. However, the perturbative stability analysis reveals important theoretical considerations that point to future work. While the background dynamics are viable, the perturbative sector requires further investigation. Future work will focus on extending this analysis to include CMB data, constructing a full effective field theory description to ensure perturbative stability, and applying the GTD framework to other generalized dark-energy scenarios.

Acknowledgment

AK is thankful to IIEST, Shibpur, India, for providing Institute Fellowship (SRF).

Data Availability Statement

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

References