Constraints and cosmography of $\Lambda$CDM in presence of viscosity

In this work, we study two scenarios of the Universe filled by a perfect fluid following the traditional dark energy and a viscous fluid as dark matter. In this sense, we explore the simplest case for the viscosity in the Eckart formalism, a constant, and then, a polynomial function of the redshift. We constrain the phase-space of the model parameters by performing a Bayesian analysis based on Markov Chain Monte Carlo method and using the latest data of the Hubble parameter (OHD), Type Ia Supernovae (SNIa) and Strong Lensing Systems. The first two samples cover the region $0.01<z<2.36$. Based on AIC, we find equally support of these viscous models over Lambda-Cold Dark Matter (LCDM) taking into account OHD or SNIa. On the other hand, we reconstruct the cosmographic parameters ($q,j,s,l$) and find good agreement to LCDM within up to $3\sigma$ CL. Additionally, we find that the cosmographic parameters and the acceleration-deceleration transition are sensible to the parameters related to the viscosity coefficient, making of the viscosity an interesting physical mechanism to modified them.


I. INTRODUCTION
Nowadays, one of the challenges in cosmology is the understanding and description of the accelerated expansion phase of the Universe. Several cosmological observations give support to this phenomena; firstly it was confirmed by Supernovaes of the Type Ia (SNIa) [1] then, by the acoustic peaks of cosmic microwave background radiation (CMB) [2] and supported by direct measurements of the Hubble parameters (OHD) [3], Baryon acoustic oscillations (BAO) [4], and strong lensing systems (SLS) [5]. The simplest model, called Λ-Cold Dark Matter (LCDM), describes very well such cosmological observations and proposes a cosmological constant (Λ or CC) characterized by an equation of state (EoS) w = −1 to model the accelerated expansion of the Universe and dust matter (w = 0) to simulate the dark matter evolution at the background level. These two extra components correspond about ∼ 95% of the total [2], being the rest of components associated with baryons and relativistic species like photons and neutrinos. Also, in the literature these first two ingredients are known as dark energy (DE) and dark matter (DM).
Besides its successful at large scale, LCDM presents several problems at local scales, for instance, the well known missing satellite problem that refers to the dis-crepancy of about 10 times more dwarf galaxies obtained by the numerical simulations based on LCDM model and the observed ones in cluster of galaxies [6,7]. Also, the well-known core-cusp problem [8]. Furthermore, there are open questions concerning to the origin of the cosmological constant. In this vein, it is the concordance problem that consists in a disagreement of about 120 orders of magnitude in the CC value measured from the Quantum Field Theory point of view and the one obtained from cosmological measurements [9][10][11]. Additionally, the degeneracy problem which afflicts also the LCDM model, refers to the inability of measuring the energy-momentum of each component, instead the total one. In other words, this implies the inability to know if the dark sector is composed by one or several components 1 .
Several models have emerged in order to propose alternatives to the LCDM paradigm, for instance, braneworld [13], Chaplygin gases [14], Unimodular gravity [15,16], among others [17], have entered into the scene as a greater contenders, resolving conundrums that the LCDM cannot, moreover, scalar fields as DM [18][19][20][21], axion [22,23], etc, are an important approaches to resolve the problem of DM. In this vein, fluids with viscosity are great candidates not only to aboard the DM problem, but also, the DE problem from an unifying approach. Indeed, the viscous models could affront problems like the H 0 tension, and the problem associated with the matter fluctuation amplitude. This kind of models also have the characteristic of predicting an earlier transition to an accelerated phase in comparison with the standard cosmological model (see for example [24][25][26]). This is because the causative of the Universe acceleration could be related with a dynamical DE and not a cosmological constant as states the consensus model. Moreover, we cannot ignore the recent results of the Experiment to Detect the Global EoR Signature (EDGES) which detect an excess of radiation (not predicted by the standard model) in reionization epoch, specifically in z ≈ 17 [27], which can be boarded through the viscous model scenario.
In the non perfect fluid context, there are two kind of viscosity known as shear and bulk. The shear viscosity does not play an important role at late epochs of the Universe because it does not satisfy the cosmological principle as the bulk viscosity does. Then, at late times, it is of great interest to study the bulk viscosity in any of its formalisms and there are, typically, two ways to aboard it, using the Eckart theory [28] or the Israel-Steward (IS) theory [29]. Although, in the Eckart approach the propagation of the perturbations on the non perfect fluid occurs at infinite speed, it is a simpler theory than the IS formalism, allowing to study more complex forms of the viscosity. For instance, bulk viscosity coefficient has been modelled as a constant [30][31][32][33], polynomials [34][35][36], and hyperbolic [36,37] functions. Moreover, it allows easily to explore the presence of interacting terms in the viscous fluid [38]. Because the non perfect fluid should satisfy the near equilibrium condition of thermodynamics, the pressure of the fluid must be greater than the one generated by the viscosity. Then, to alleviate this condition, it is convenient to add an extra perfect fluid such as CC or any other. It is worth mention that this assumption has a price to pay in the viscous fluids because they loss the power to describe the dynamics of the Universe by unifying the DM and the DE in an unique viscous fluid.
Regarding to IS formalism, there is a small quantity of viscosity models that has been studied. For instance, some authors [39][40][41][42][43][44] consider the viscosity coefficient as ξ ∼ ρ s , where ρ is the energy density of the fluid and some solutions have been studied for the case s = 1/2. Nevertheless, to solve the mentioned near equilibrium condition of thermodynamics, [45] studies the evolution of the Universe by adding a fluid as CC to a viscous fluid (considering ξ ∼ ρ s ).
In this work, we study a Universe filled by two fluids under the Eckart formalism, a perfect fluid as DE mimicking the dynamics of the CC, while a non-perfect fluid as DM which is added to its EoS a viscosity term. We will assume two cases for the viscosity coefficient: a constant and a polynomial function of the redshift. It is interesting to remark that our model is reduced to LCDM model (when the relativistic species are negligibly) by turning off the viscosity. In this sense, we aim mainly to compare the cosmographic parameters of our viscous models and LCDM to understand the viscosity effects in LCDM. Additionally, we analyze correlations between the cosmography and the viscosity parameters. Hence, we perform a Bayesian Markov Chain Monte Carlo (MCMC) analysis to constrain the free model parameters using the largest samples of the observational Hubble parameter distance measurements [3], Supernovae Pantheon sample [46] and Strong Lensing Systems [5].
The manuscript is structured as follow: In Sec. II, it is presented the mathematical details of the viscous dark fluid model. Section III describes the OHD, SNIa and SLS samples, together with the joint analysis of the three previously mentioned samples, in Sec. IV we present the constraints and results, finally in Sec. V we give some discussions and conclusions.

II. VISCOUS DARK FLUID MODEL
In what follow, we summarize briefly the two components models following the mathematical formalism described in [33,47]. Thus, we consider a flat Universe (k = 0) under the FriedmannLemaîtreRobertsonWalker (FLRW) metric, where a(t) is the scale factor as function of the cosmic time and dΩ 2 = dθ 2 + sin 2 θdϕ 2 is the solid angle. In addition, the bulk viscosity term, Π, is introduced through the energy-momentum tensor as an effective pressurep, i.e., being h µν = g µν + u µ u ν and u µ = (1, 0, 0, 0) is the cuadrivelocity in the co-moving coordinate system,p = p + Π with p is the total barotropic pressure of the fluids presented in the Universe. In this vein, we consider a viscous dust-like matter, with EoS w = p/ρ = 0 coupled to a perfect fluid behaving as the cosmological constant, i.e., w = −1. Then, the Friedmann equations are where κ 2 = 8πG and H =ȧ/a is the Hubble parameter, ρ m and ρ de are the energy density of the viscous matter and dark energy, respectively. The effective pressure is p = −ρ de − 3ξH where we have inspired the form of Π to be proportional to the H through the bulk viscosity coefficient ξ = ξ(t) as analogy to the fluid dynamics which the viscous effects are proportional to the velocity. Based on the procedure presented in [33,47], the homogeneous solution (when Π = 0) of Eqs. (5)-(6) with respect to the redshift z is where i = m, de, and ρ i0 means the energy density at current epochs of the i-component. Notice that this solution corresponds to an Universe filled by two perfect fluids. Additionally, it is interesting to see that it is an approximation of the LCDM model which we have despised the relativistic species component at the background level. The viscosity effects will appear as correction terms in the general solution which is built as the sum of the homogeneous solution and the particular one.
In other words, we can express the general solution as It is important to remark that this form of the general solution is a simplification coming from the idea that the viscous term is the same for all the fluids, implying that the function u(z) is the same for any fluid. For more details see [47]. Then, it can be possible to find an expression for u(z) in function also with the bulk viscosity coefficient, through the following differential equation, as By defining the dimensionless viscosity coefficient λ(z) = ξ(z)H 0 /ρ cr , being ρ cr = 3H 2 0 /κ 2 the critical density, the above expression reads as where I 0 is an integration constant and we have defined , Ω i0 = ρ i0 /ρ cr , being i = m, de. Now, in order to solve the above integral, we need to propose a form of λ(z). By considering the form 9λ(z) = λ 0 + λ 1 (1 + z) n , [48], we have where λ 0 , λ 1 ans n are free parameters. Before to solve the latter expression, it is convenient to write the dimensionless Hubble parameter as function of the redshift, z, defined as where i = m, de and H 0 = 100 h km s −1 Mpc −1 and h is the dimensionless Hubble constant measured at current epochs. Then, by integrating and using (8) and (12), we obtain, where we have used E(0) = Ω(0) = 1, Ω(z) = Ω m0 (1 + z) 3 + Ω de0 , and 2 F 1 is the hypergeometric function. In this work, based on the results obtained on [48] we will set n = −2. For λ 1 = 0, we have the case when the bulk viscosity coefficient is constant; in this case, we obtain where sinh −1 is the inverse of hyperbolic sin function.
In the latter we can observe that the standard cosmology (LCDM) is recovery when λ 0 = 0. It is worth to notice that for these two cases, a future singularity appears for z → −1, crossing to the Phantom DE region. An evidence of this behavior it is discussed later in Sec. V, which coincide with previous studies provide by [49] about a Little Rip for this kind of models.

III. COSMOLOGICAL SAMPLES
In order to analyse the viscous models, we use three data samples provided by direct measurements of the Hubble parameters, Supernovae observations and strong lensing system. This section is devoted to describe them and report the constraints of the model parameters considering each sample and also by performing a joint analysis.

A. Observational Hubble parameter measurements
Currently, the direct way of measuring the accelerated expansion of the Universe is through the Hubble parameter using the differential age tools and BAO measurements. The largest sample that include these observations is compiled by [3] covering a range 0.07 < z < 2.36 with 51 points. We will refer to this sample as the Observational Hubble distance (OHD). Then, to constrain the parameter phase space, Θ = (h, Ω m , λ 0 , λ 1 ) and Θ = (h, Ω m , λ 0 ) for the model when λ is polynomial (setting n = −2) and constant, respectively, we build the χ 2 -function to be minimize as In the above expression, H i obs represents the observational Hubble parameter with its uncertainty σ i obs at the redshift z i . On the other hand, H th represents the theoretical expression related to the Eqs. (13) and (14) for the polynomial and constant form of the bulk viscosity coefficient (λ), respectively.

B. Type Ia Supernovae data
Apart of the OHD sample, it is useful to include the luminosity distance measurements obtained through Type Ia Supernovae (SNIa). The largest sample, collected by Pantheon [46], covers the redshift region 0.01 < z < 2.3 with 1048 measurements of the bolometric apparent magnitude. In order to compare it with our models, we compute the theoretical one as where, M is a nuisance parameter. The quantity d L (z), known as the dimensionless luminosity distance, is given by where c is the speed of light and E(z) the dimensionless Hubble parameter presented in Eqs (13) and (14). Then, we build the χ 2 -function as where Cov −1 refers to the inverse of the covariance matrix and m obs to the observed quantity of m.

C. Strong Lensing System
Finally, we also use the latest compilation of the strong lensing systems (SLS) provided by [5] which consider only systems where the lens are early type galaxies. With a total of 205 points, the sample covers the redshift region 0.0625 < z l < 0.9280 for the lens galaxy and 0.196 < z s < 3.595 for the source. To constrain cosmological parameters, it is useful to built the chi-square function given by [50] In the latter, D th is the theoretical angular diameter distance ratio defined by being the angular diameter distance to the source, Similarly, D ls means the angular diameter distance between the source and the lens galaxy, i.e., the same previous equation but now evaluated in the region z s < z < z l . The observational counterpart, D obs , is built as where θ E is known as the Einstein radius and σ is the velocity dispersion of the lens DM halo. Its uncertainty is estimated by where δθ E and δσ are the uncertainty of θ E and σ, respectively. Notice that δD obs does not consider correlation between θ E and σ. Following [5] we choice an absolute uncertainty on θ E of 0.05 for the data points without a reported uncertainty.

D. Joint analysis
In order to minimize the total χ 2 -function for each model given by we perform a Bayesian MCMC analysis based on emcee module [51]. After achieving a lower value than 1.1 in the Gelman-Rubin criteria [52] for each free model parameter to stop the n-burn phase, we obtain 5000 chains with 250 steps, each one to explore the confidence region taking into account a Gaussian prior on the Hubble constant h and on DE density Ω m0 according to the Planck results [2]. Additionally, we select flat priors for λ 0 and λ 1 based on [48]. Table I summarizes the priors used in the Bayesian analysis.

IV. CONSTRAINTS AND RESULTS
In this section we present and discuss our results obtained in the Bayesian analysis. Figure 1 shows the constrained phase space of the parameters for CVM (top panel) and PVM (bottom panel), respectively at 68% (1σ), 95% (2σ), and 99.7% (3σ) confidence level (CL).  both viscous models using OHD, SNIa, SLS and the joint analysis, respectively. The reported uncertainties correspond to 1σ CL. We find good agreement to data according to the chi-square value of the models, being the worse fit of them to SLS data with χ 2 = 602.3 (603.6) for the CVM (PVM). It is worth mentioning that CVM and PVM are better approximations to LCDM because consider a DE component with EoS w = −1 and a viscous matter component. By setting the λ = 0 (or ξ = 0), we recover the LCDM with DE and DM components, which the relativistic species are negligible. In other words, we are studying as a first approximation the consequences of the viscosity effects in the LCDM model under the Eckart formalism. In this sense, it is convenient to compare statistically our models with LCDM model. Instead of using χ 2 criteria, it is more convenient to use others such as the Akaike information criterion (AIC) [53,54] and Bayesian information criterion (BIC) [55] because they allow to compare models with different degree of freedom. They are defined as AIC = χ 2 + 2k and BIC = χ 2 + k log(N ) where k is the number of free parameters, and N is the size of the data sample.  lar support for both models in the AIC. Additionally, we find that the CVM is equally supported than ωCDM, the Chevallier-Polarski-Linder (CPL), and the Jassal-Bagla-Padmanabhan (JBP) models [5] and no evidence against as well. In contrast, we find a slightly less support for PVM than previous mentioned parameterizations, and a strong evidence against PVM over them.  We reconstruct the cosmographic parameters for both viscous models, CVM and PVM, in the redshift region 0 < z < 2.5. Figure 2 displays the Hubble and deceleration parameters for CVM (top panel) and PVM (bottom panel) respectively. We find a good agree-ment with LCDM (star markers) within 3σ CL along the mentioned range. Additionally, we calculate the corresponding deceleration-acceleration transition redshift values, obtaining z t = 0.683 +0.025 −0.022 and z t = 0.755 +0.063 −0.051 .
When we compare them with LCDM value (z LCDM t = 0.642 +0.014 −0.014 ), we obtain a deviation of 1.9σ and 2.2σ. Figure 4 shows the reconstruction of the high order cosmographic parameters (j, s, l) and their uncertainties at 1σ and 3σ. The jerk parameter gives us information about the dynamics of the DE EoS, corresponding to j = 1 to w = −1 for the DE component. In this sense, although our models are consistent up to 3σ with LCDM, our best fitting values indicate an effective dynamical DE EoS. In other words, we observe a deviation to LCDM in the jerk parameter due to the viscosity contributions.
On the other hand, it is well known that the snap and lerk parameters and higher cosmographic parameters does not have a well established physical meaning, however they are an important part of the Taylor series of the Hubble parameter in cosmography, giving us more precision in the preferred model by observations. Notoriously, the viscosity generates important differences (for snap and lerk) in the PVM case in comparison with LCDM, mainly at lower redshift, which in turn means that dark energy is dynamic and not constant in this kind of models. This is consistent with cosmographic studies where the LCDM is not the preferred model (see Refs. [56][57][58][59] for details).
Furthermore, we report the cosmographic parameters at current epochs, q 0 = −0.568 +0.018 −0.102 for a Universe filled by a single nonperfect fluid with polynomial viscosity, achieving a deviation of 5.3σ (3.7σ) relative to CVM (PVM). It is interesting to observe that PVM is a generalization of this polynomial single fluid model because PVM includes a DE component but has a better approximation to LCDM. In summary, Fig. 5 compares the q 0 and z t obtained for several models reported in the literature (see these Refs. [14,38,48,[60][61][62][63][64][65][66][67] for more details about the models and cosmological data used). The vertical band represents ±1σ around the central value of LCDM. It is interesting to observe that models with interactions and viscosity have a deceleration-acceleration transition earlier than LCDM and most of the models are in good agreement with the expected value q 0 for LCDM. Additionally, we estimate the effective EoS at current times, obtaining values of w ef f 0 = −0.712 +0.012 −0.014 and w ef f 0 = −0.648 +0.043 −0.038 for CVM and PVM respectively. These values have a deviation from one reported in [14] (Chaplygin-like model) of 3.95σ and 3.36σ respectively. Furthermore, we find a deviation of 5.66σ and 1.43σ to the value presented in [43] (dissipative dark fluid model), and good agreement (up to 2.1σ) with the one reported in [63].  Figure 3 displays the evolution of E(z) 2 as function of (1 + z) 3 for both viscous models. In this diagram, the black line represents the LCDM, and the region z > 0 over LCDM corresponds to a quintessence behavior and the region z > 0 below LCDM line to a phantom one. In this vein, it is interesting to observe that CVM behaves as phantom DE for any time, presenting a divergence in the future. Although, PVM also presents a singularity in the future (z = −1), PVM presents a transition from phantom to quintessence around z ≈ 2.2, conserving such behavior close to z = −1, before to finish in a Big Rip. It is worth mentioning that a Big Rip at the future is a typical final state presented in the viscous models (see for instance [47,49]).
On the other hand, we estimate the correlation between cosmographic parameters, z t , and those proper of the models through the formulae corr(x, y) = Cov(x, y)/σ x σ y where Cov(x, y) is the covariance coefficient between x and y, and σ x (σ y ) is the standard deviation of x (y). We find a strong correlation (abs(corr)> 0.7) between those related to the viscosity (λ 0 and λ 1 ) and q 0 , and between z t . In other words, the viscosity terms allow us to modify both physical quantities that characterize the Universe. Additionally, we also find a negative strong correlation between (λ 0 , λ 1 ) and the high order cosmographic parameters (j 0 , s 0 , l 0 ) at current epochs.
Finally, although it has been shown (see for instance [48,63]) that with an unique dissipative fluid it is possible to explain the accelerated expansion of the Universe, we split the dark sector into two fluids, dark matter and dark energy, to satisfy the near equilibrium condition required by thermodynamics, i.e., the viscous pressure Π must satisfy the condition |Π/p| 1 at least at current epochs, where p is the total equilibrium pressure of the fluids. On the other hand, by requiringä > 0 for late times, the condition −Π > p + ρ/3 must be satisfied, where ρ is the total energy density of the fluids. The latter is fulfilled because an acceleration phase occur as it is shown in the q(z) reconstruction. Figure 6 shows the evolution of |Π/p| over the redshift range −1 < z < 2.5. It is interesting to observe that the both viscous models are far equilibrium in the past, but have a trend to go to near equilibrium at current epochs. Additionally, due to the divergence at z = −1, CVM and PVM ends far from an equilibrium point.

V. SUMMARY
In this work, we study the Universe filled by one perfect fluid modelling the DE component and a non-perfect fluid describing the matter. While the DE is characterized by the traditional EoS w = −1, the matter follows an effective EoS given by w m = −3ξ(t)H(t)/ρ m (t), where ξ is the viscosity coefficient. As we mentioned, for ξ = 0 we recover the cosmology of LCDM when the relativistic species are negligible. In this vein, we constrain the free parameters assuming that ξ is constant and is a polynomial function of the redshift using the latest measurements and compilations of the Hubble parameter, type Ia supernovae and strong lensing. We find good agreement to data according to χ 2 -value presented in Table II. Additionally, we compare statistically our models with LCDM using AIC and BIC criteria, and we find that the OHD (SNIa) sample prefers equally our viscous model than LCDM. In contrast, the joint (OHD+SNIa+SLS) data give strong evidence against the viscous models over LCDM. Based on the joint analysis, we present the dynamics of the four main cosmographic parameters (q,j,s,l) in the region 0 < z < 2.5, in order to elucidate the differences and advantages regards to the LCDM model. In this context, the reconstruction of j for both models indicates an effective dynamical DE, besides the causative of the universe acceleration in this model, is proposed as a CC behavior. Therefore j parameter point us the effects of the viscosity, which characterize the model under study. Moreover, we find strong correlations between the viscosity parameters (λ 0 and λ 1 ) and those parameters that characterized the Universe dynamics (q 0 and z t ), allowing modify them according to the viscosity property.
Finally, we observe that both viscous models finish in a Big Rip state at z → −1 as it is typical for this kind of models. In addition, we observe that the PVM presents a slightly transition to a quintessence region at late times, i.e., a viscous coefficient as function of the redshift may change the behavior of the Universe from phantom to quintessence. Further studies should be developed to address this result, which will be presented elsewhere, and also to extend the model by including a radiation component, in order to study the effects of the viscosity at the CMB epochs and the consequences of these dissipative effects at the perturbative level as in those presented in [68][69][70][71]. In particular, the authors in [70] establish the constrictionξ < 0.24 (at 2σ) at z = 0 for a constant viscosity. In this sense, we estimate an upper bound of ξ = 9λ(0) < 0.149 andξ = 9λ(0) < 1.617 at 95% CL for CVM and PVM respectively. We can observe that the first is tighter than the one reported in [70], while our result for PVM is less restricted than the mentioned. Additionally, our results for PVM are also in agreement within 3σ with the one obtained by [69] (for ξ ∼ ρ s ), and we obtain a tighter constraint for CVM.