© The Authors 2026

1. Introduction

Cosmology is experiencing a phase of heated debate, as due to the impressive amount of new data points, certifying unexpected dark-energy evolution within the cosmic background (Abdul Karim et al. 2025)¹ and seeking alternative dark-matter particle candidates following the failure of the weakly interactive massive particle (WIMP) hypothesis². In this mare magnum of challenges for modern cosmology, the perplexing Hubble tension refers to the discrepancy between local measurements and values inferred from early Universe observations of the present-day expansion rate of the Universe, i.e., the Hubble constant, H₀. More precisely, this tension arises when comparing local measurements by the SH0ES collaboration that, utilizing the Cepheid-calibrated cosmic distance ladder, reports H₀^R = (73.04 ± 1.04) km/s/Mpc (Riess et al. 2022) or cosmographic methods (Aviles et al. 2012; Dunsby & Luongo 2016; Luongo & Muccino 2020; Luongo & Quevedo 2014; Bamba et al. 2012), and values obtained from the cosmic microwave background (CMB) within the ΛCDM model, leading to H₀^P = (67.36 ± 0.54) km/s/Mpc (Planck Collaboration VI 2020). This tension at a 4.1σ confidence level³ appears quite unexpected and may hint at new physics (Cortês & Liddle 2024; Hu et al. 2024; Hu & Wang 2023; Gomez-Valent & Solà Peracaula 2024; Yang et al. 2019; Pan et al. 2020; Dainotti et al. 2021, 2022; De Simone et al. 2025; Dainotti et al. 2025), i.e., new fields (Niedermann & Sloth 2020; Chatrchyan et al. 2025), new dark interactions (Blinov et al. 2020; Nygaard et al. 2023; Mirpoorian et al. 2025; Odintsov et al. 2021), or issues related to experimental procedures (Pedrotti et al. 2025; Lee et al. 2025; Banik & Kalaitzidis 2025). For a review, see, for example, Vagnozzi (2023).

To address the persistent H₀ tension, we focused on two main avenues. The first is relaxing the assumption that Milky Way dust properties are well-defined at low redshifts, thereby allowing for spatial or temporal variations in dust composition and extinction laws, at z ∼ 0. Alternatively, assuming that type Ia supernovae (SNe Ia) remain reliable, high-precision, and standardizable candles, a second possibility to resolve the Hubble tension may lie in early Universe physics, potentially requiring modifications to the pre-recombination dynamics or extensions to the standard ΛCDM framework, at redshifts before the CMB.

Approaches to alleviating the H₀ tension are broadly classified into late- and early-time scenarios. In practical terms, provided that the sound horizon angular scale, θ_★ = r_★/[(1 + z_★)D_A(z_★)],  at the recombination⁴ redshift, z_★, is measured with a precision of 0.03% (Planck Collaboration VI 2020), late-time solutions alter the angular diameter distance, D_A(z_★), of the last scattering surface, whereas early-time solutions typically modify the comoving sound horizon, r_★ = r_s(z_★), providing the options listed below.

• Late times. Modifying H(z) in z ∈ [0, z_★] implies a different value of H₀, keeping the comoving sound horizon, r_★, and the angular diameter distance, D_A(z_★), fixed. This ensures that the angular scale of the sound horizon, θ_★, remains unchanged.

• Early times. Modifying r_★ prior to recombination by changing either z_★; the sound speed, c_s(z); or the Hubble parameter, H(z). Accordingly, H₀, as embedded in the expression of D_A, varies in order to keep θ_★ fixed.

The late-time approach appears quite unlikely, since H(z) is well constrained in that regime, and it turns out to be difficult to expect physical modifications. Conversely, typical early-time modifications tend to reduce r_★, reconciling, de facto, the locally measured value of H₀ with CMB observations by keeping θ_★ fixed (Schöneberg et al. 2022; Kamionkowski & Riess 2023; Poulin et al. 2023).

Motivated by the above considerations, we focused on early times and investigated an alternative approach to early dark-energy (EDE) models, without involving scalar fields. To do so, we addressed the H₀ tension predicting the existence of an additional barotropic fluid with a positive equation of state (EoS) that satisfies the Zel’dovic limit –namely matter with pressure – which is different from radiation and EDE models. We thus directly modified the comoving sound horizon, r_s, and H(z) due to the presence of this additional fluid before recombination. This procedure led to an effective increase in the Hubble constant, slightly modifying the ΛCDM model into a fast and direct extension. In particular, the additional barotropic fluid turns out to be subdominant with respect to dust, neutrinos, radiation, and so on. To show this, we computed exclusion plots, where the range of values for the fluid magnitude and EoS was carefully limited; this turned out to be substantially smaller than radiation, i.e., it did not substantially modify the CMB structure as well. Once the suitable priors were fixed, we performed a cosmological fit including neutrinos and radiation, in agreement with existing Planck results. The fit was worked out by analyzing our constructed paradigm, hereafter dubbed the Λ_ωs cold dark matter (Λ_ωsCDM) model, by means of a Markov chain Monte Carlo (MCMC) method, modifying the freely available CLASS code⁵ (Lesgourgues 2011); with this approach, we employed the ranges of values that cannot exclude the presence of a further fluid. Suitable findings, in a spatially flat universe, were thus preliminarily obtained, indicating that our slight extension of the ΛCDM model is not a priori disfavored with respect to the usual case. We find inconsistency with a stiff matter fluid (i.e., EoS equal to 1), but rather with an EoS smaller than that of radiation, i.e., ω_s ≲ 1/3, with a normalized density significantly smaller than that of radiation. Physical consequences and corresponding explanations of this matter-like fluid, inspired by pure stiff matter were thus explored in detail. Further, its fundamental origins, in terms of extra fields, are summarized in this paper.

The paper is structured as follows. In Sect. 2, we discuss the Hubble tension, summarizing the physical approaches and explaining our new treatment. In Sect. 3, we detail our analysis based on the use of a MCMC simulation and discuss the corresponding outcomes. Accordingly, we physically justify our fluid and the physics behind it. In Sect. 4, we conclude our work, emphasizing consequences and perspectives of our treatment.

2. Tackling the Hubble tension

Alleviating the H₀ tension should keep the angular scale of θ_★ fixed while modifying the comoving sound horizon, r_s, or the angular diameter distance, D_A. In particular, considering the standard approach, H₀ is effectively embedded within the calculation of r_★ = r_s(z_★), which is given by

$$\begin{array}{c}\hfill r_{★}={\displaystyle {\int }_{z_{★}}^{\infty }\frac{c_s(z)}{H(z)}\mathrm{d}z=\frac{c}{\sqrt{3}{H}_{★}}{\int }_{z_{★}}^{\infty }\sqrt{\frac{{\rho }_{★}}{\rho (z)}}\frac{\mathrm{d}z}{\sqrt{1+R}},}\end{array}$$(1)

where H_★ = H(z_★) and ρ_★ = ρ(z_★) are the Hubble parameter and the total density, ρ(z), at the recombination, respectively, and R = (3/4)(Ω_b/Ω_γ)(1 + z)⁻¹. The quantities Ω_b and Ω_γ are the present-day (normalized to the critical density) baryon and photon densities. Hence, when maintaining θ_★ fixed, if r_★ decreases, D_A(z_★) might also decrease. The latter is given by

$$\begin{array}{c}\hfill D_{\mathrm{A}}(z_{★})=\frac{c}{1+z}{\displaystyle {\int }_0^{z_{★}}\frac{\mathrm{d}z}{H(z)},}\end{array}$$(2)

leading to a “decrease” of H₀, inferred from local measurements. As late-time constraints disfavor changes in H(z) after recombination, we focused on early-time solutions, which are outlined in the following three classes of approaches, plus a fourth scenario that consists of our new fluid, suggesting a minimal corrected Λ_ωsCDM paradigm, inspired by matter with pressure.

2.1. Modifying the sound speed

The first solution changes the sound speed. Thus, before recombination, photons and baryons are tightly coupled through Thomson scattering and behave effectively as a single fluid. The corresponding adiabatic sound speed (i.e., at constant entropy S) definition is

$$\begin{array}{c}\hfill c_s^2=c^2{\left(\frac{\partial P}{\partial \rho }\right)}_S.\end{array}$$(3)

As the baryons are assumed to be nonrelativistic this photon-baryon fluid exhibits a total pressure due to the photon component only, P = ρ_γ/3, and a total energy density given by the sum of the two components, yielding ρ = ρ_γ + ρ_b. Thus, the sound speed is given by

$$\begin{array}{c}\hfill c_s^2=c^2{[3(1+R)]}^{-1}.\end{array}$$(4)

It is clear that altering this value leads to a modification of the comoving sound horizon. Specifically, if we decrease it by introducing new additional species into the plasma, the comoving sound horizon becomes smaller, increasing the inferred value of H₀ (Evslin & Sen 2018).

2.2. Varying the recombination redshift

The second way to heal the H₀ tension at early times is to vary z_★. In particular, changes in the recombination history directly impact the size of the sound horizon in the photon–baryon plasma. However, since the corresponding angular scale at decoupling, θ_★, is measured precisely, any variation in its physical scale must be balanced by a corresponding change in the angular diameter distance to the surface of the last scattering. The CMB’s sensitivity to the ionization history depends on the photon-scattering rate, which is governed by the ionization fraction

$$\begin{array}{c}\hfill X_e(z)=\frac{n_e(z)}{n_H(z)},\end{array}$$(5)

where n_e and n_H denote the number densities of free electrons and hydrogen nuclei, respectively. The recombination redshift is typically defined by the condition X_e(z_★) = 0.5. A modified recombination history leads to a different ionization fraction, thereby shifting z_★ and altering r_★ (Lynch et al. 2024a,b; Mirpoorian et al. 2025). In particular, a lower z_★ reduces r_s, which increases the value of H₀.

2.3. Resorting early dark energy

A third and very promising approach to alleviating the H₀ tension involves the presence of EDE. This term refers to any form of dark energy that becomes dynamically relevant before recombination, thus modifying the Hubble expansion rate and reducing the comoving sound horizon. The presence of a dark energy component with a time-dependent EoS, non-negligible at early times, has been extensively studied (Doran & Robbers 2006; Linder & Robbers 2008; Grossi & Springel 2009; Maggiore et al. 2011; Pettorino et al. 2013; Poulin et al. 2019; Niedermann & Sloth 2021; Braglia et al. 2020; Smith et al. 2021; Gómez-Valent et al. 2021; Sabla & Caldwell 2022), originally to investigate its impact on structure formation (Doran et al. 2001; Doran & Wetterich 2003; Wetterich 2004), and its effect on increasing the value of H₀ depends on the specific model considered. The impact of EDE is determined by the fractional contribution of the field to the total energy density, expressed as f_EDE(z) = ρ_EDE(z)/ρ(z), which increases over time. In particular, the scalar field initially remains frozen in the potential, resulting in a nearly constant energy density. Later, it is released from this configuration by a physical mechanism, such as a drop in Hubble friction below a critical threshold or a phase transition that alters the shape of the potential. Once the field becomes dynamical, its energy density dilutes faster than that of matter, leading to a rapid decay of its contribution to the total energy budget.

Therefore, the impact of EDE on the expansion rate is localized in redshift, occurring prior to recombination, and becomes negligible at late times. This allows it to reduce the comoving sound horizon, leaving the angular diameter distance unaltered. A successful alleviation of H₀ tension generally requires a fractional contribution of f_EDE ∼ 10%, peaking at redshifts of z ∼ 10³–10⁴ and subsequently diluting at a rate at least as fast as that of radiation (Kamionkowski & Riess 2023; Poulin et al. 2023).

2.4. Introducing matter with pressure

As fourth attempt, we propose a novel early-time resolution to the Hubble tension whereby a new barotropic fluid is introduced, i.e., a pseudo-relativistic matter characterized by an arbitrary EoS, which reduces the comoving sound horizon, simultaneously changing H(z) and increasing the inferred value of the Hubble constant. Our new setup is thus determined by a tightly coupled fluid with three main components:

• Photons, with P_γ = ρ_γ/3.

• Baryons, with negligible pressure, P_b ≈ 0.

• Relativistic matter, with P_s = ω_sρ_s, where ω_s > 0, and not necessarily ω_s = 1 (stiff matter).

Here, the total energy density and pressure are

$$\begin{array}{cc}& \rho ={\rho }_{\gamma }+{\rho }_b+{\rho }_s,\hfill \end{array}$$(6)

$$\begin{array}{cc}& P=P_{\gamma }+P_s=\frac{1}{3}{\rho }_{\gamma }+{\omega }_s{\rho }_s,\hfill \end{array}$$(7)

and the adiabatic sound speed becomes

$$\begin{array}{c}\hfill c_s^2=c^2{\left(\frac{\partial P}{\partial \rho }\right)}_S=\frac{c^2(\partial {\rho }_{\gamma }+3{\omega }_s\partial {\rho }_s)}{3(\partial {\rho }_{\gamma }+\partial {\rho }_b+\partial {\rho }_s)}.\end{array}$$(8)

When introducing a new dimensionless ratio involving the new fluid and the photon densities, W ≡ 3(1 + ω_s)ρ_s/(4ρ_γ), the modified sound speed before the recombination is

$$\begin{array}{c}\hfill c_s^2=\frac{c^2(1+3{\omega }_{s}W)}{3(1+R+W)},\end{array}$$(9)

indicating that this new species alters not only the Hubble parameter, but also the primordial plasma.

Barotropic fluids differ from a scalar field as the kinetic energy of the latter may be either positive or negative (see, e.g., Dymnikova & Khlopov 1998, 2000, 2001; Linder & Scherrer 2009; Steinhardt et al. 1999; Chiba 2006; Scherrer & Sen 2008b,a; Shahalam et al. 2015; Wolf & Ferreira 2023; Carloni & Luongo 2025). Ensuring the Zel’dovic limit, namely ω_s > 0, we ended up with the requirements listed below.

• The fluid has positive pressure and may contribute positively to the stress–energy tensor, unlike any pure dark energy (Mukhopadhyay et al. 2011; Luongo & Muccino 2018; D’Agostino et al. 2022).

• The fluid dilutes faster than cold dark matter and baryons as the Universe expands⁶, providing an energy density of ρ_s(z)∝(1 + z)^{3(1 + ω_s)}.

• Physically, it would contribute non-negligibly at early times to alter H₀ and remain subdominant at all stages of the Universe’s evolution.

Our proposed class of barotropic fluids represents a naive and classical example in theoretical cosmology, since it simply guarantees the Zel’dovich limit. Besides pure radiation, corresponding to ultrarelativistic particles, such as photons, neutrinos, etc., our Λ_ωSCDM model is inspired by the possible existence of genuine stiff matter, characterized by ω_sm = 1 (Chavanis 2015a; Grøn 2024). However, we explored the case of 0 < ω_s < 1/3 as it is more relevant to our study. Even though this scenario has not been directly observed yet, in principle it appears possible, as it corresponds to a matter that has a nonzero pressure and is pseudo-relativistic, but subdominant to radiation (see Fig. 1). Thus, if this fluid exists, one may wonder what differences there may be between this approach and a pure EDE contribution. To this end, the key differences between a positive-ω barotropic fluid and EDE are summarized in Table 1. In summary, even though both a barotropic fluid, with ω_s > 0, and EDE can impact the early Universe, they are fundamentally different in their origin and behavior. Mainly, a barotropic fluid – subdominant but present throughout the Universe’s evolution – follows a constant EoS and scales predictably with the expansion. Conversely, EDE models are constructed with specific dynamical mechanisms that make the energy density significant only during precise epochs, after which it disappears as it is jeopardized by unpleasant fine-tuning issues. Moreover, the EDE cosmology extends the six-parameter ΛCDM model, {ω_b, ω_cdm, H₀, A_s, n_s, τ_reio}, with three additional parameters: {Ω_φ, log(a_c),ϕ_i} (Poulin et al. 2019). In contrast, the Λ_ωsCDM model only introduces two extra parameters, {ω_s, Ω_s}, thereby reducing the overall model complexity.

Fig. 1. ΛωsCDM density parameters: Λ (dotted gray), pressureless matter (dashed blue), radiation (solid red), and matter with pressure (dotted-dashed green); vertical black lines (from left to right) mark the following redshifts: recombination, pressureless matter–radiation equivalence, and pressureless matter–matter with pressure equivalence.

3. Experimental setup

To prove the existence of our additional fluid with density, Ω_s(1 + z)^{3(1 + ω_s)}, we define the underlying background cosmology as an extension of the standard ΛCDM framework, in which the comoving distance is defined as

$$\begin{array}{c}\hfill D_{\mathrm{M}}(z)=c{\displaystyle {\int }_0^z\frac{\mathrm{d}{z}^{\prime }}{H(z^{\prime })},}\end{array}$$(10)

where the parameters describing the expansion history enter in the definition of the Hubble parameter:

$$\begin{array}{c}\hfill H(z)=H_0{[{\mathrm{\Omega }}_{\nu }(z)+{\displaystyle {\sum }_i}{\mathrm{\Omega }}_i{(1+z)}^{3(1+{\omega }_i)}]}^{1/2};\end{array}$$(11)

here, the set, i = {Λ, m, γ, s}, labels the cosmological constant, dust, radiation, and matter with pressure, respectively, exhibiting barotropic indexes as ω_i = { − 1, 0, 1/3, ω_s}. The energy density of neutrinos Ω_ν(z) scales as radiation (with EoS ω_ν = 1/3) at early times and as matter (with ω_ν = 0) at late times. Ensuring a spatially flat Universe, Ω_k = 0, naturally implies Ω_Λ = 1 − Ω_m − Ω_s − Ω_γ − Ω_ν. Thus, from Eq. (10) we can define the luminosity distance D_L(z) = (1 + z)D_M(z) and the diameter angular distance D_A(z) = D_M(z)/(1 + z).

3.1. Numerical results

We utilized the following datasets: the Pantheon+ sample of SNe Ia, the second data release (DR2) of baryonic acoustic oscillations (BAO) from the Dark Energy Spectroscopic Instrument (DESI), hereafter dubbed DR2 DESI-BAO, and the CMB data, indicated in detail below and in Appendix A.

3.1.1. Pantheon+ and SH0ES.

The Pantheon+ is a catalog of 1701 SNe Ia with redshifts of z ∈ [0, 2.3]. It comprises 18 different samples (Scolnic et al. 2022) and includes the released SH0ES Cepheid host distance anchors that provide the value of H₀ that is in tension with CMB determination (Brout et al. 2022).

When using SH0ES Cepheid host distances, the 1701 SN distance modulus residuals are given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mu }_i=\{\begin{array}{cc}{\mu }_i-{\mu }_{\mathrm{C},i}\hfill & i\in \text{Cepheid hosts}\hfill \\ {\mu }_i-{\mu }_{\mathrm{th}}(z_i)\hfill & \text{otherwise}\hfill \end{array},\end{array}$$(12)

in which μ_i are the inferred SN distance moduli and μ_C, i are the Cepheid-calibrated host-galaxy distance moduli by SH0ES; the model’s distance moduli are

$$\begin{array}{c}\hfill {\mu }_{\mathrm{th}}(z_i)=m_i-M=5log(D_{\mathrm{L}}(z_i)/10\mathrm{pc}),\end{array}$$(13)

where m_i are the rest-frame B-band apparent magnitudes of each SN Ia, and the rest-frame B-band absolute magnitude, M, can be viewed as a nuisance parameter and marginalized (Conley et al. 2011). Thus, the cosmological parameters can be constrained by maximizing the log-likelihood

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{SN}}=-\frac{1}{2}\mathrm{\Delta }{\mathit{\mu }}^T{C}^{-1}\mathrm{\Delta }\mathit{\mu },\end{array}$$(14)

where C = C_SN + C_C is the total covariance matrix of SNe and Cepheid host distances (accounting for both statistical and systematic errors)⁷.

3.1.2. DR2 DESI-BAO.

The DR2 DESI-BAO dataset is shown in Table A.1. We utilized a total of N_D = 13 measurements: six transverse distances, D_M; six Hubble rate distances, D_H (six data points); and one angle-averaged distance, D_V, all normalized over comoving sound horizon at the baryon drag redshift, z_d, i.e., r_d = r_s(z_d) (Abdul Karim et al. 2025). Breaking the r_d − H₀ degeneracy is mandatory to provide the explicit expression of r_d as a function of the model parameters given in Sect. 2. Here, we only fixed z_d = 1059.94 ± 0.30 (Planck Collaboration VI 2020). Provided that D_M is given in Eq. (10), the other two distance measurements are given by

$$\begin{array}{c}\hfill \frac{D_{\mathrm{H}}(z)}{r_d}=\frac{c}{r_{d}H(z)},\frac{D_V(z)}{r_d}=\frac{D_{\mathrm{H}}^{1/3}(z)D_{\mathrm{M}}^{2/3}(z)}{z^{-1/3}{r}_d}.\end{array}$$(15)

Thus, the total BAO log-likelihood is given by

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{BAO}}\propto -\frac{1}{2}{\displaystyle \sum _{i=1}^{N_D}}{[Y_i-Y(z_i)]}^T{C}_{\mathrm{BAO}}^{-1}[Y_i-Y(z_i)],\end{array}$$(16)

with Y_i denoting the DR2 measurements and Y(z_i) = {D_M(z_i)/r_d, D_H(z_i)/r_d, D_V(z_i)/r_d} given in Eqs. (10) and (15). Here, C_BAO indicates the BAO covariance matrix⁸.

3.1.3. CMB data.

We adopted the Planck 2018 CMB high-ℓ TTTEEE, low-ℓ TTEE, and lensing likelihoods. In particular, we used the Plik (Planck Collaboration V 2020) high-ℓ likelihood for TT over 30 ≤ ℓ ≤ 2508 and for TE and EE over 30 ≤ ℓ ≤ 1996. The low-ℓ TT and EE likelihoods cover 2 ≤ ℓ ≤ 29 and are based on the Commander algorithm and SimAll likelihood, respectively (Aghanim et al. 2020; Planck Collaboration IV 2020). The CMB Plik lensing likelihood is included following Planck Collaboration VIII (2020).

By accounting for all the above described datasets, we can deduce the model cosmological parameters by maximizing the total log-likelihood function defined as

$$\begin{array}{c}\hfill ln\mathcal{L}=ln{\mathcal{L}}_{\mathrm{SN}}+ln{\mathcal{L}}_{\mathrm{BAO}}+ln{\mathcal{L}}_{\mathrm{CMB}}\end{array}$$(17)

using the priors given by Table 2.

MCMC fitting was performed using our modified version of CLASS⁹. This code solves the Boltzmann equations for perturbations in the early Universe and accounts for the following components: photons, baryons, cold dark matter, neutrinos, dark energy (in the form of a cosmological constant Λ), and our extra matter fluid with pressure¹⁰. Specifically, in agreement with the Planck analysis, we assumed two massless neutrino species and one massive species with m_ν = 0.06 eV, leading to N_eff = 3.046 (Ade et al. 2019). In our analysis, chain convergence is assumed to be achieved when the Gelman–Rubin statistic satisfies R − 1 < 0.02 (Gelman & Rubin 1992).

The results are summarized in Table 3. The Λ_ωsCDM model yields a better fit to the data, with ${\chi }_{\mathrm{s}}^2=4093.06$, compared to χ_Λ² = 4115.48 for the ΛCDM model; this corresponds to $\mathrm{\Delta }{\chi }^2={\chi }_{\mathrm{s}}^2-{\chi }_{\mathrm{\Lambda }}^2=-22.42$. Then, to estimate the statistical preference, we computed the Bayesian evidence using the public code MCEvidence (Heavens et al. 2017)¹¹. A model comparison was done by computing Δlog B = log B_s − log B_Λ = +6.2, which yields decisive evidence in favor of the Λ_ωsCDM scenario based on the modified Jeffreys’ scale (Trotta 2008).

Using the results of the Λ_ωsCDM model from Table 3, we were able to compute the equivalence epochs among fluids. Figure 1 displays the density parameters, Ω_i(z) = Ω_i(1 + z)^{3(1 + ω_i)}, of all fluids¹². In particular, we mark the following epochs with vertical black lines (from left to right): (i) the recombination redshift, z_★ = 1089.92; (ii) the equivalence redshift between the pressureless matter and radiation (photons and neutrinos) equivalence redshift (z_eq ≈ 3600); and (iii) the equivalence redshift between matter with and without pressure (z_ms ≈ 68400). From the plot it also evident that our matter fluid with pressure is (a) subdominant with respect to pressureless matter and radiation (photons and neutrinos) at z_★, and (b) subdominant with respect to photons at late times and photons+neutrinos at early times. For the above reasons, all the other model parameters were not significantly modified –with respect to the ΛCDM case– by the presence of such a further contribution. The corresponding contour plots, obtained using GetDist, are portrayed in Fig. B.1.

3.2. Physical consequences

In our picture, forcing a positive EoS and a fluid in the Zel’dovic limit implies a result of

$$\begin{array}{c}\hfill {\omega }_s=\frac{1}{3}-ϵ,{\mathrm{\Omega }}_s(z)={\mathrm{\Omega }}_s{(1+z)}^{4-3ϵ},\end{array}$$(18)

with $ϵ=0.{040}_{-0.009}^{+0.003}$. This kind of fluid, with an EoS similar to radiation and here reinterpreted as matter with pressure, can mainly be attributed to specific cases of interest.

3.2.1. Unlikely physical scenarios

Below, we report two direct possibilities that appear to be unlikely.

3.2.1.1. Non-gauge fields such as the Proca Lagrangian.

In such a case, the Hubble tension can be solved including the presence of massive electromagnetic fields, where the action of photon mass, m_P, reduces the EoS. For example, (a) ensuring the spatial component of the vector potential, A, and the magnetic field, B, vanish, say A = B = 0, and (b) that for the scalar potential, A₀, and the electric field, E, the relation m_P²A₀² ≪ E² holds, we end up with

$$\begin{array}{c}\hfill {\omega }_{\mathrm{P}}\simeq \frac{1}{3}-\frac{4A_0^2{m}_{\mathrm{P}}^2}{3E^2},ϵ=\frac{4A_0^2{m}_{\mathrm{P}}^2}{3E^2}.\end{array}$$(19)

In a thermal background at temperature T, we have the following scalings: A₀ ≈ T and E² ≈ 2T⁴. Thus, using the above value of ϵ and the temperature at recombination T_★ = 0.26 eV, we obtain $m_{\mathrm{P}}=0.{064}_{-0.007}^{+0.002}$ eV, which is, however, excluded by the most conservative existing upper limit on the photon mass m_γ < 10⁻¹⁸ eV obtained from the solar-wind magnetic-field structure (Ryutov 2007).

3.2.1.2. Dark photons.

Highly debated and recently put under examination, dark photons are hypothetical massive gauge bosons associated with a dark-sector hidden U(1)_D symmetry. They may have a small kinetic mixing with massless photons and very weak or no direct coupling to standard model matter (Cline 2024). Massive dark photons behave as Proca fields:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{D}}\simeq \frac{1}{3}-\frac{4A_0^2{m}_{\mathrm{D}}^2}{3E^2},ϵ=\frac{4A_0^2{m}_{\mathrm{D}}^2}{3E^2},\end{array}$$(20)

but with much greater freedom, since dark photons are not subject to the same tight experimental bounds as standard photons. Therefore, for these constituents, we simply took the previous bound, $m_{\mathrm{D}}=0.{064}_{-0.007}^{+0.002}$ eV.

3.2.2. Most plausible physical interpretation: A thermal, semi-relativistic scalar component

In view of the above, our matter with pressure fluid may be interpreted either as a thermalized scalar field, whose nature can be associated with the quasi-quintessence, or as vector fields, in which the Proca field could represent an immediate and plausible explanation. We did not investigate alternatives making use of extended theories of gravity and/or modifications as they are being beyond the scope of this work (see, e.g., van der Westhuizen et al. 2025 and Paliathanasis 2025a,b).

The MCMC constraints indicate that the additional component behaves as a quasi-radiation fluid, with a constant ϵ > 0 and a small but non-vanishing energy density, which remains subdominant with respect to the standard radiation sector at recombination, roughly reading Ω_s/Ω_γ ∼ 1/3. In this regime, the fluid is naturally interpreted as an effectively relativistic species, whose pressure is mildly reduced by a finite mass scale.

Among the various possibilities discussed above, a thermal, semi-relativistic scalar component may be the most likely counterpart to associate with our fluid. In this case, the desired EoS is obtained in a standard kinetic-theory approximation. The simplest approach is to consider, for example, an additional purely relativistic species assumed to be a sterile neutrino. Sterile neutrinos imply severe bounds on the effective number of relativistic species, whereas scalar fields come as an extension of the generalized K essence described in Appendix C.

In the semi-relativistic approximation, we have

$$\begin{array}{c}\hfill {\omega }_{\mathrm{S}}\simeq \frac{1}{3}-\frac{m_{\mathrm{S}}^2/3}{{\kappa }^2{T}^2+m_{\mathrm{S}}^2},ϵ=\frac{m_{\mathrm{S}}^2/3}{{\kappa }^2{T}^2+m_{\mathrm{S}}^2},\end{array}$$(21)

where m_S is the mass and κ ≈ 2.7 for neutrinos (from Fermi-Dirac statistic) and κ ≈ 3.2 for scalars (from Bose-Einstein statistic). Both massive components behave as radiation with ω_S ≈ 1/3 when m_S ≪ T and as pressureless matter with ω_S ≈ 0 when m_S ≫ T. Limits on the mass, for example at the recombination, can be obtained as follows:

$$\begin{array}{c}\hfill m_{\mathrm{S}}\simeq \kappa T_{★}\sqrt{\frac{3ϵ}{1-3ϵ}}.\end{array}$$(22)

From the above discussion, it immediately appears evident that sterile neutrinos, or Majorana fermions more broadly, shall be accounted for in N_eff; moreover, since in our MCMC fits it is fixed to the standard value of N_eff = 3.046, any mass estimates would be physically incorrect. In addition, recent developments seem to indicate that sterile neutrinos appear quite disfavored, as they have not been detected by independent experiments so far (see, e.g., Grimm et al. 2025).

Conversely, for thermal scalar fields (namely bosons) there is no such limitation. We thus obtain

$$\begin{array}{c}\hfill m_{\mathrm{S}}=0.{31}_{-0.04}^{+0.01}eV.\end{array}$$(23)

This interpretation has two attractive features, which are summarized below.

• It naturally explains why the extra component is very close to radiation. Indeed, the fluid is simply a thermal species that is still semi-relativistic at z_★. Accordingly, it contributes to early-time expansion, and precisely to the sound horizon integral, without requiring EDE; but it is also constantly subdominant throughout the evolution of the Universe.

• It clarifies why the component does not behave as an additional clustering agent at late times. For a barotropic fluid with constant ω_s and adiabatic perturbations, the rest-frame sound speed is c_s, s² = ω_s; hence, pressure support prevents gravitational growth once modes enter the horizon, in close analogy with the photon sector. We faced this point later while dealing with Newtonian perturbations.

A more accurate study that could exclude dark photons and exotic species could be pursued in future developments. We would check whether, statistically speaking, our fluid may be favored to exist even at late times.

3.3. Comparison with Big Bang nucleosynthesis constraints

The existence of a new barotropic fluid not only alters the expansion rate of the Universe at any epoch, but could also affect the primordial abundances of light elements produced during Big Bang nucleosynthesis (BBN). In this work, we investigated the consistency with the observed ⁴He mass fraction value Y_p = 0.2450 ± 0.0050 (with statistical and systematic errors) and the CMB-inferred ratio between the number densities of photons and baryons η = (6.104 ± 0.058)×10⁻¹⁰ (Navas et al. 2024). We use simplified analytical calculations, excluding the possible identification of the extra fluid with an additional neutrino species as stated above, and we verified the consistency with the standard BBN framework when Ω_s = 0.

At temperatures of T ≫ 1 MeV, proton and neutron number densities are the same because the weak interaction rate Γ(T) = qT⁵, with q = 9.6 × 10⁻²² MeV⁻⁴ (Matei et al. 2025), is much larger than the expansion rate:

$$\begin{array}{c}\hfill H(T)=\frac{T^2}{M_p}{\left[\frac{8{\pi }^{3}g(T)}{90}\right]}^{1/2}{[1+{\delta }_s(T)]}^{1/2},\end{array}$$(24)

where g(T) is the number of relativistic degrees of freedom, which depends on the primordial plasma temperature, and

$$\begin{array}{c}\hfill {\delta }_s(T)=\frac{{\mathrm{\Omega }}_s}{{\mathrm{\Omega }}_r}{\left[\frac{g_S(T)}{g_S(T_0)}\right]}^{{\omega }_s-1/3}{\left(\frac{T}{T_0}\right)}^{3{\omega }_s-1}.\end{array}$$(25)

The ratio between the entropy degrees of freedom, g_S, computed at T, and the CMB temperature, T₀, appears as a consequence of the entropy conservation, g_S(T)T³a³ = const, because from MCMC fits we have the current extra-fluid density parameter at T₀, but in Eq. (24) it has to be evaluated at T. In fact, the temperature of the plasma is not an exact proxy of the scale factor, a = (1 + z)⁻¹, which is altered by changes in the number of relativistic species across the evolution of the Universe.

As the Universe cools down, the above competing rates become comparable –Γ(T_f) = H(T_f)– at the “freeze-out” temperature –T_f ≃ 0.8 MeV– for which protons and neutrons are locked in a ratio of 5 : 1. This is the onset of the BBN.

The first nuclear reaction involves the synthesis of deuterium nuclei. However, because of the extremely high abundances of high-energy photons, compared to baryons, deuterium and ⁴He syntheses are delayed until the plasma temperature drops well below the binding energy of a deuterium nucleus, B_D = 2.224 MeV. This is referred as to as the deuterium bottleneck. Thus, the production of deuterium nuclei starts when their creation and destruction rates become equal. This occurs at a temperature of T_d ≃ 0.06 MeV, which can be evaluated as

$$\begin{array}{c}\hfill T_{\mathrm{d}}\simeq -B_D{\{ln\eta +\frac{3}{2}ln\left(\frac{4\pi T_{\mathrm{d}}}{m_{\mathrm{n}}}\right)+ln\left[\frac{3\zeta (3)}{{\pi }^2}\right]\}}^{-1}\end{array}$$(26)

(Mukhanov 2005), where m_n is the neutron mass. Thus, the mass fraction of ⁴He can be evaluated using the neutron-to-proton number density ratio (n/p)_d evaluated across freeze-out and deuterium production phases; i.e.,

$$\begin{array}{c}\hfill Y_p=\frac{2{(n/p)}_{\mathrm{d}}}{1+{(n/p)}_{\mathrm{d}}},{\left(\frac{n}{p}\right)}_{\mathrm{d}}={\mathrm{e}}^{-\mathrm{\Delta }m/T_{\mathrm{f}}}{\mathrm{e}}^{-(t_{\mathrm{d}}-t_{\mathrm{f}})/{\tau }_{\mathrm{n}}},\end{array}$$(27)

where Δm = 1.293 MeV is the neutron–proton mass difference and τ_n ≈ 880 s is the neutron mean lifetime. The cosmic times at freeze-out and deuterium production, t_f and t_d, respectively, were computed from the relation t(T_x) = ∫_Tx^+∞dT/[T H(T)]. Moreover, in evaluating the cosmic times and δ_s(T), one has to take into account that at T < T_f neutrinos decouple from the other particles, while keeping the same temperature; afterwards, at T ≲ m_e, the e^± pairs (with masses of m_e) in the thermal background begin to annihilate, heating up the photon background relatively to the decoupled neutrinos. Here, we considered the e^± annihilation to be completely settled below an effective temperature of T_e ≃ 0.15 MeV. In view of these considerations, our calculations consider

$$\begin{array}{c}\hfill \begin{array}{ccc}T\ge T_{\mathrm{e}}:\hfill & \hfill g(T)=g_S(T)& =2+(7/4)(2+3)=10.75,\hfill \\ T<T_{\mathrm{e}}:\hfill & \hfill g(T)& =2+(7/4)N_{\mathrm{eff}}{(4/11)}^{4/3}=3.38,\hfill \\ \hfill & \hfill g_S(T)& =2+(7/4)N_{\mathrm{eff}}(4/11)=3.94.\hfill \end{array}\end{array}$$

The results are shown in Fig. 2, where in the standard BBN we recovered Y_p = 0.2463, whereas in our scenario we obtained the value Y_p = 0.02478 ± 0.0021, which is fully consistent with the observational bounds.

Fig. 2. BBN constraints in η–Yp plane for the ΛωsCDM model (blue curve with the light blue 1σ error bands), compared with the standard result (dashed black curve) and the observational constraints on Yp (green box) and η (gray box).

3.4. Linear perturbations

Our extra component was modeled as a barotropic perfect fluid with a constant and positive EoS, P_s = ω_sρ_s, with an index ω_s –according to our fits– close to, but slightly below, the photon value ω_γ = 1/3. Accordingly, the new species behaves as a “quasi-radiation” component: it contributes to the early-time background expansion and to the primordial plasma properties, while remaining subdominant with respect to the standard radiation sector.

It is important to clarify what we mean by “no clustering”. This does not imply the absence of perturbations. Rather, it means that, once a Fourier mode enters the horizon, the density contrast of the ω_s–fluid does not develop a gravitationally growing mode that could seed, or substantially contribute to, late-time structure formation. This behavior is analogous to radiation; i.e., photon perturbations exist and are observable (e.g., through acoustic physics), but they do not collapse into bound structures.

For a barotropic fluid with constant ω_s and adiabatic perturbations (no intrinsic entropy mode), the rest-frame sound speed equals the adiabatic one,

$$\begin{array}{c}\hfill c_{s,s}^2\equiv {\omega }_s;\end{array}$$(28)

so, for ω_s = 𝒪(10⁻¹), the pressure support is relativistically large. On subhorizon scales, the competition between pressure support and gravity leads to a damped acoustic behavior for δ_s ≡ δρ_s/ρ_s, with an effective frequency ∼c_s, sk/a. This is conveniently summarized by the Jeans scale estimate,

$$\begin{array}{c}\hfill k_{J,s}^2\sim \frac{4\pi G{\rho }_{\mathrm{tot}}}{{(c_{s,s}/a)}^2}\sim \frac{a^2{H}^2}{c_{s,s}^2}\Rightarrow {\lambda }_{J,s}\sim \frac{2\pi a}{k_{J,s}}\sim \frac{c_{s,s}}{H},\end{array}$$(29)

which enters the density-contrast equation,

$$\begin{array}{c}\hfill {\ddot{\delta }}_s+2H{\dot{\delta }}_s+(\frac{c_{s,s}^2{k}^2}{a^2}-4\pi G{\rho }_{\mathrm{tot}}){\delta }_s\simeq 0,\end{array}$$(30)

which was written in a Newtonian context; i.e., analogous to the more complicated fully relativistic version.

Our findings show that, using Eq. (28), the Jeans length is on order the Hubble radius. Hence, for subhorizon modes (k ≫ aH) we have

$$\begin{array}{c}\hfill \frac{c_{s,s}^2{k}^2}{a^2}\gg H^2,\end{array}$$(31)

and δ_s undergoes acoustic oscillations with Hubble damping rather than a growing gravitational instability. Consequently, our new fluid does not act as an additional clustering component at current times. Its crucial contribution lies in the early-time expansion and horizon-scale dynamics, consistently with our Boltzmann treatment of perturbations.

Future works will explore the fully relativistic versions of Eq. (30). However, the main outcome cannot depart from the one above reported, indicating that the Λ_ωsCDM paradigm does not depart significantly from the expectations of the standard cosmological model.

4. Final remarks

In this work, we focused on the early-time Universe to heal the observational H₀ tension and proposed an alternative view to EDE models, namely invoking the existence of an additional barotropic fluid component that is physically inspired by stiff matter-like fluids exhibiting relativistic properties differently from dust. In so doing, we ensured this fluid obeyed the Zel’dovic condition on the EoS; i.e., it had to possess ω_s > 0. This matter-like constituent was constructed with nonnegligible pressure, which is quite distinct from both radiation and standard EDE candidates, as well as pure dust-like counterparts or neutrinos. Thus, by introducing such a component into the early cosmological dynamics, we directly altered the expansion rate, H(z), and the comoving sound horizon, r_s, prior to recombination, modifying, moreover, the sound speed, c_s. This enabled an effective increase in H₀ through a minimal and direct extension of the ΛCDM paradigm, which we call the Λ_ωsCDM model. To prevent any significant modification of cosmic measurements, we required the additional barotropic component to remain subdominant with respect to radiation, matter, and dark energy throughout the entire cosmic-expansion history. Thus, we carried out a full cosmological parameter estimation using a MCMC simulation, involving Pantheon+ with SH0ES SNe, DR2 DESI-BAO, and the Planck 2018 CMB high-ℓ TTTEEE, low-ℓ TTEE and lensing likelihoods. To this end, we modified the CLASS code (Lesgourgues 2011) to include the dynamics of the extra fluid and performed the analysis in accordance with the latest Planck observational constraints. In this respect, we included radiation and neutrino contributions, employing the assumption of spatial flatness.

We showed that the proposed extension is favored with respect to the ΛCDM benchmark model, while our findings ruled out the case of a purely stiff fluid with EoS ω = 1, instead favoring a scenario in which the fluid exhibits an EoS slightly below that of photons, ω_s ≲ ω_γ, and a density parameter smaller than the photons’ background, which is Ω_s/Ω_γ ≃ 0.33. Indeed, we demonstrated that the additional fluid contribution remained significantly smaller than that of radiation and thus did not produce any substantial alteration in the CMB anisotropy power spectrum.

Afterwards, we discussed the physical implications of such a fluid, emphasizing its phenomenological consistency and theoretical plausibility. We thus propose that the fluid may be interpreted in view of a possible Proca-like field, where the effect of nonzero mass for the vector field can alter the EoS, departing from the genuine case given by radiation. We thus conclude that at least two vector fields may be present just before the last scattering surface, without, however, excluding more exotic scenarios involving dark radiation, axions, and so on. Based on these considerations, we evaluated the ⁴He mass fraction, Y_p, produced during the BBN with the inclusion of our extra fluid, and we verified the consistency with the observational bounds.

After this crucial point, we explored the linear Newtonian perturbations, showing analogy between our new fluid and pure radiation. Accordingly, we conclude that it is unlikely that our new fluid would influence the clustering without significantly altering observations at current times.

We also inferred constraints on equivalence and transition between epochs between pressureless matter and radiation (occurring at z_eq ≈ 3600) and between matter with and without pressure (occurring at z_ms ≈ 68400). The matter fluid with pressure is subdominant with respect to pressureless matter and radiation at z_★, and always subdominant with respect to radiation. These features show that all other constituents are not significantly modified by the presence of this kind of further contribution.

Future works will compare this matter with pressure fluid through other approaches that predict matter with pressure, by complicating the barotropic factors, or simply involving thermodynamic alternatives other than a constant and positive ω_s. Moreover, we will see the influence on structure formation of this fluid and, above all, its fundamental nature in more detail; thus, we will further explore the role of Proca-like fields in early-time cosmology. A direct comparison with additional constituents, such as neutrinos, relativistic matter, and so on. will also be the object of future speculations. Last but not least, further experimental analyses will be pursued with the aim of refining our bounds, taking a deeper look at perturbations, and adopting a fully relativistic first- and second-order analysis.

Acknowledgments

YC acknowledges Francesco Pace for valuable discussions and is grateful to Vivian Poulin for support during his stay at LUPM, where this work was finalized. OL acknowledges Maryam Azizinia for discussion and support by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19 )” – Next Generation EU (NGEU).

References

Appendix A: Datasets used in the analysis

Pantheon+ SN Ia and Cepheid host distance moduli, together with the covariance matrixes C_SN and C_C can be found at https://github.com/PantheonPlusSH0ES/DataRelease.

DR2 DESI-BAO with N_D = 13 measured values of D_M(z)/r_d, D_H(z)/r_d, and D_V(z)/r_d are taken from (Abdul Karim et al. 2025) and shown in Table A.1. The covariance matrix is provided by https://github.com/LauraHerold/MontePython_desilike.

CMB data, consisting of the Planck 2018 high-ℓ TTTEEE, low-ℓ TTEE, and lensing likelihoods, can be accessed via the Planck Legacy Archive https://pla.esac.esa.int/pla/#home.

Appendix B: Experimental results

The contour plots of the parameters for Λ_ωsCDM and ΛCDM models, with the corresponding mean values listed in Table 3, are shown in Fig. B.1.

Fig. B.1. Contour plots of ΛωsCDM and ΛCDM models with respect to the data points employed in this work. The contours have been obtained by fixing cosmic radiation, in a universe with a further fluid, constituted by matter, dark energy under the form of a cosmological constant, Λ, plus our barotropic fluid, whose density is given by Ωs(z) = Ωs(1 + z)3(1 + ωs), ensuring ωs > 0.

Appendix C: Alternative models of beyond-dust matter

It is possible to extend the EoS for matter, including pressure, by considering alternatives to the standard dust hypothesis. For example, barotropic fluids can be considered using effective unified dark energy-dark matter models (Dunsby et al. 2016), or working extensions of the Chaplygin gas that appear suitable at the level of background and perturbations (Dunsby et al. 2024a). Alternatives, quite more complicated, would imply the presence of barotropic fluids that act differently from standard ones, namely profoundly modifying the EoS, including complicated terms inspired by solid state physics, such as logotropic fluids (Boshkayev et al. 2021; Capozziello et al. 2023; Boshkayev et al. 2019), Anton-Schmidt approaches (Capozziello et al. 2018, 2019a), Murnagham EoS (Dunsby et al. 2024b) and so on.

C.1. Generalized K-essence models

Matter with pressure can be accounted involving a barotropic fluid with constant EoS (Luongo & Muccino 2022; Boshkayev et al. 2020), albeit a fundamental representation of such a prerogative can be obtained ensuring c_s = 0 within a scalar field Lagrangian, whose effective form can be described through a Lagrange multiplier as (Luongo & Muccino 2018), ℒ = λ Y(X,φ) + K₀ − V(φ), with the generalized kinetic term of the scalar field X, kinetic functions K₀ = const and Y, the Lagrangian multiplier λ, and the potential V, yielding both spontaneous symmetry breaking and chaotic inflation pictures (D’Agostino et al. 2022). The above paradigm defines density and pressure

$$\begin{array}{c}\hfill \rho =2X{\partial }_X\mathcal{L}-(K_0-V),P=K_0-V.\end{array}$$(C.1)

Distinguishing between baryonic (b) and cold dark matter (c) yields K₀ = K_b, 0 + K_c, 0 and ${\partial }_X\mathcal{L}=\lambda {\partial }_X(Y_b+Y_c)$.

After the transition induced by the symmetry breaking, the field φ reaches the minimum of the potential V₀ = −0.22 M_P⁴, where M_P is the Planck mass, and the magnitude of dark matter pressure K_c, 0 ≫ K_b, 0 can be adjusted to cancel out the vacuum energy density, healing the fine tuning issue associated with the cosmological constant Λ (Belfiglio et al. 2023, 2024). In this picture, Eqs. (C.1) become

$$\begin{array}{c}\hfill \rho \approx 2X\lambda {\partial }_X(Y_b+Y_c)-K_{b,0},P\approx K_{b,0}.\end{array}$$(C.2)

The above fluid is isentropic, leading to the identification 2Xλ ∂_X(Y_b+Y_c) ≡ (ρ_b+ρ_c)(1 + z)³ that describes the densities of pressureless dark matter and baryonic matter endowed with negative pressure, K_b, 0. For additional details, see Luongo & Muccino (2018).

C.2. Anton-Schmidt and logotropic fluids

The issue of negative pressure poses a problem to obtain direct experimental evidence of negative pressure in laboratory. However, in condensed matter physics, it is possible to argue a pressure that turns out to be locally negative. A notable example is given by the Anton–Schmidt EoS (Anton & Schmidt 1997), which empirically describes the pressure of a crystalline solid under isotropic deformation. For cosmological applications, it becomes

$$\begin{array}{c}\hfill P=A{\stackrel{\sim }{\rho }}^{{\gamma }_G+\frac{1}{6}}ln\stackrel{\sim }{\rho },\end{array}$$(C.3)

where γ_G denotes the so-called Grüneisen parameter, whose macroscopic definition is related to the thermodynamic properties of the material. Here, $\stackrel{\sim }{\rho }\equiv \rho /\overline{\rho }$ is the normalized density with respect to a reference density $\overline{\rho }$, and A represents the amplitude of the fluid. Remarkably, as γ_G → −1/6, the fluid becomes purely logotropic, i.e., described by

$$\begin{array}{c}\hfill P_l=Aln{\stackrel{\sim }{\rho }}_l,\end{array}$$(C.4)

initially proposed in Refs. (Chavanis 2015b, 2016), where the subscript l indicates its logotropic nature. The model has been motivated by the hypothesis that the same fluid acting as dark energy could also account for dark matter in galaxies, thus alleviating the cusp problem observed at the centers of spiral galaxies (Chavanis 2015b).

These two examples of matter with pressure, represented by Eqs. (C.3) and (C.4), also constitute possible unified dark energy models, albeit their validities have been severely criticized (see, e.g., Carloni et al. 2025; Boshkayev et al. 2021).

C.3. The Chaplygin and Murnaghan fluids

The Chaplygin gas was originally introduced as a simple prototype of unified dark energy models. This is achieved through an EoS that scales inversely with the density itself. However, such models are ruled out by observational data (Fabris et al. 2011; Carloni et al. 2025).

More recently, motivated by a two-fluid interacting Lagrangian framework, where one fluid transports vacuum energy, it has been proposed that the pressure can be expressed by the Murnaghan EoS,

$$\begin{array}{c}\hfill P=-A({\stackrel{\sim }{\rho }}^{-\alpha }-1),\end{array}$$(C.5)

which corresponds to a Chaplygin gas supplemented by a cosmological constant term, where A is the amplitude of the fluid with density normalized as $\stackrel{\sim }{\rho }\equiv \rho /\overline{\rho }$, and α is a constant. This paradigm was first introduced in (Dunsby et al. 2024a) and subsequently reviewed in (Dunsby et al. 2024b), demonstrating promising fits to cosmological observations.

Within all the above alternatives, it becomes clear that matter with pressure can be straightforwardly modeled by a constant EoS parameter, as realized in the two simplest cases, i.e., the generalized K-essence and stiff-like fluids. Hence, rather than complicating the CMB consequences adopting logotropic, Anton-Schmidt and/or Chaplygin gas, we focused above on our Λ_ωsCDM framework that suggests matter (1) with enhanced elastic properties, (2) not fully cold, and (3) distinct from neutrinos, while closely approximating radiation.

All Tables

All Figures

	Fig. 1. ΛωsCDM density parameters: Λ (dotted gray), pressureless matter (dashed blue), radiation (solid red), and matter with pressure (dotted-dashed green); vertical black lines (from left to right) mark the following redshifts: recombination, pressureless matter–radiation equivalence, and pressureless matter–matter with pressure equivalence.
In the text

	Fig. 2. BBN constraints in η–Yp plane for the ΛωsCDM model (blue curve with the light blue 1σ error bands), compared with the standard result (dashed black curve) and the observational constraints on Yp (green box) and η (gray box).
In the text

	Fig. B.1. Contour plots of ΛωsCDM and ΛCDM models with respect to the data points employed in this work. The contours have been obtained by fixing cosmic radiation, in a universe with a further fluid, constituted by matter, dark energy under the form of a cosmological constant, Λ, plus our barotropic fluid, whose density is given by Ωs(z) = Ωs(1 + z)3(1 + ωs), ensuring ωs > 0.
In the text