Testing the equation of state for viscous dark energy

Some cosmological scenarios with bulk viscosity for the dark energy fluid are considered. Based on some considerations related to hydrodynamics, two different equations of state for dark energy are assumed, leading to power-law and logarithmic effective corrections to the pressure. The models are tested with the latest astronomical data from Type Ia supernovae (Pantheon sample), measurements of the Hubble parameter $H(z)$, Baryon Acoustic Oscillations and Cosmic Microwave Background radiation. In comparison with $\Lambda$CDM model, some different results are obtained and their viability is discussed. The power-law model shows some modest results, achieved under negative values of bulk viscosity, while the logarithmic scenario provide good fits in comparison to $\Lambda$CDM model.


I. INTRODUCTION
Over the last years, the intriguing question behind the late-time acceleration over the cosmological expansion has drawn much attention over the scientific community, being one of the most challenge problems in theoretical physics. As the universe is observed to be approximately homogeneous and isotropic at large scales, the best description for the cosmological evolution within General Relativity (and also within other geometrical theories) is provided by the so-called Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime. Nevertheless, as can be easily shown trough the equations, an accelerating expansion in a FLRW universe requires generally an effective negative pressure fluid, which violates the energy conditions and has been called dark energy when talking on the late-time acceleration, although is also assumed for producing the so-called cosmic inflation at early times. Dark energy models have been deeply analysed, by considering each one's pros and cons, inherent of every theoretical model (see [1]). However, the main remaining problem strikes on the similar predictions provided by a very large number of dark energy models, including new fields or modifications of General Relativity. Hence, a great effort is being done to reduce the number of models by studying more complex features, going beyond in perturbation theory or getting more accurate constraints by the use of the great amount of incoming data.
Moreover, as the main property for the dark energy fluid lies on the negativity of its pressure in order to achieve an accelerating expansion, a plausible scenario is provided by a viscous fluid, since its pressure is affected by a bulk viscosity term, as is well described in hydrodynamics, such that the fluid may keep the energy conditions satisfied, but providing an effective negative pressure. Such possibility has been widely analysed in the literature as well (for a review see [2]) and some realistic scenarios have been proposed where the well established knowledge of hydrodynamics is applied to cosmology [3]. In general, most of the analysis consider bulk viscosity as a possibility of assuming a viscous fluid, as can keep the conditions on homogeneity and isotropy, widely contrasted by the observations. In this sense, some authors have dealt with the possibility of incorporating viscosity to dark matter [4], which can lead to unify dark matter and dark energy under the same fluid, as for instance in the case of the Chaplygin gas [5][6][7] or a logotropic fluid [8]. Some other models consider a proper dark energy fluid with viscosity [9][10][11]. Such possibility may lead to a fluid with a negative pressure that may even cross the phantom barrier [9,10,[12][13][14][15]. In addition, the viscosity terms may play an essential role during the early time inflationary stage [2,16].
In the present paper we consider two models for dark energy with bulk viscosity, in the same line as proposed in some works previously. Here our models are based on some phenomenological considerations or inspired by condensed matter physics. In that sense, we study a model whose viscosity depends on the powers of the energy density and the Hubble parameter, being considered as effective corrections to a perfect dark energy fluid, while the second model is inspired in Anton-Schmidt's equation of state for crystalline solids ( [17]). We test the model by using some recent observational data, and techniques developed in some previous papers [18][19][20][21]. The observational datasets include the latest Pantheon sample [22] of Type Ia supernovae (SNe Ia), estimations of the Hubble parameter H(z), observational manifestations of baryon acoustic oscillations (BAO) and cosmic microwave background radiation (CMB) [23][24][25][26]. By using the likelihood, we obtain the best fits for the free parameters of the models and compare to ΛCDM model. The paper is organized as follows: section II is devoted to a brief description of dynamical equations and the models for the bulk viscosity. In Section III, we provide the observational datasets used along the paper for testing the models, which correspond to SNe Ia, H(z), BAO and CMB under investigation. Section IV is devoted to the results of the analysis of the models. Finally, in Section V we summarise the results of this work.

II. BACKGROUND
Let us start by introducing the basis of the paper. Here we are assuming a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, which can be expressed in co-moving coordinates as follows: Here we are interested in the analysis of viscous fluids, such that the bulk viscosity ζ is introduced as an effective contribution to the pressure [2,9,10]: Note that ζ(H) depends on time and can be written in terms of the Hubble rate H =ȧ/a as we are assuming a FLRW spacetime. Hence, the effective pressure p includes the viscous term B(H), which would satisfy the transportation equation, and may be extended to more general forms that contain derivatives of the Hubble parameter B = B(a, H,Ḣ) (see Ref. [2,9,10]). Here, we are considering two known models for the bulk viscosity (2).
In the first scenario, we also assume a non constant equation of state but depending on powers of the energy density w = w(ρ) while the viscosity is given by ζ(H) ∼ H 2β−1 . Then, the effective EoS yields [9,10,15]: The second model is inspired in Anton-Schmidt's equation of state for crystalline solids ( [17]), which includes a logarithmic-corrected power-law fluid and the the same viscous term as above ζ(H) ∼ H 2β−1 , leading to the following EoS [14]: For both models, A, α, B, β, ρ * are constants, essentially the free parameters of the model, while α = γ G + 1 6 , where γ G is the Grüneisen parameter. The density ρ * shows the limit where standard pressure vanishes, and can be identified with the Planck density ρ P = c 5 /( G) (see Ref. [27,28]. For aesthetic, we assume the units such as the speed of light reduces to unit c = 1. In the sections below, we will test the viability of models (3) and (4) by confronting their predictions with recent observational data, coming from different sources. Besides the dark energy component ρ x , we will also take into account the other two components that play an important role along the cosmological evolution, dust (baryons and cold dark matter) ρ m and radiation ρ r , such that the total energy density can be expressed as follows: In addition, we also assume that there is no interaction among the three components, such that they satisfy the continuity equation independently,ρ Let us now consider the scenario with the three matter components (5) with energy densities ρ m , ρ x , ρ r , where ρ x is governed by the EoS given in (3) or (4). In particular, by the power-law EoS (3), the expression for the pressure leads to p x = −ρ x +Ãρ α x +BH 2β . As usual in a FLRW Universe, cold dark matter ρ m and radiation ρ r evolve according to their continuity equations respectively (6): Here the index 0 refers to magnitudes measured at the present time t 0 while the scale factor at the present time is set as the unity, a(t 0 ) = 1. In this notation, the redshift z of a luminous object is z = a −1 − 1. On the other hand, it is convenient to rewrite the equations of state (3) and (4) for the viscous component in the following form ("logarithmic" model); (9) where the dimensionless cosmological parameter Ω x is defined as usual by the ratio among the dark energy density and the critical density: where the Hubble constant and the critical density are given by H 0 = H(t 0 ) and ρ cr = 3H 2 0 κ 2 , respectively, while the constant κ 2 ≡ 8πG. In the logarithmic model (9) we consider only the case α = 1 because of too large number of model parameters.
Then, by using the Einstein field equations, together with the FLRW metric (1), the corresponding FLRW equations are obtained: where ρ = ρ i corresponds to (5). Recall that the continuity equation (6) can be retrieved by combining (11) and (12), such that is not an independent equation, as usual in covariant theories. We can rewrite the FLRW equations (11) and (12) in terms of the cosmological parameters (10) as follows: whereas the continuity equation (6) can be expressed as function of the scale factor instead of the cosmic time for the models (8) and (9): Here Hence, by solving the system of equations (13) and (14), the cosmological evolution is obtained in terms of the scale factor for some values of the free parameters, together with the corresponding expressions for radiation and dust in terms of the scale factor, given in Eq. (7).
In order to simplify the model and reduce the number of free parameters, we are considering a flat FLRW universe, such that the curvature is assumed to be zero and the corresponding cosmological parameter leads to Ω k = 0. Moreover, we will also fix the cosmological parameter for radiation density Ω 0 r through the ratio among baryons and radiation as provided by Planck [24]: Since this value is rather small, the radiation density ρ r is assumed as negligible when fitting our models with SNe Ia, H(z) and BAO observations in the range 0 < z ≤ 2.36, as usual in most of the analysis of this kind. The component ρ r becomes important just at high redshifts, at which radiation density turns out essential to deal with is the corresponding CMB observational data for redshifts z ≃ 1000. Hence, by setting the spatial curvature to be zero Ω k = 0 and by fixing the radiation-matter ratio (15), the free parameters for both models (8) and (9) turn out: , Note that the large number of parameters N p is a lack of strength for any model in comparison with other cosmological scenarios, as the ΛCDM model, since the increasing number of free parameter may lead to a loss of information and to weaker constraints on the free parameters, which may flag the corresponding theoretical model from the point of view of information criteria [29,30]. However, we consider the Hubble constant H 0 as a nuisance parameter and reduce the effective number to N p = 5. In addition, we will also show that the power-law model (8) with α = 1, provides similar fits and errors as the case when α is considered as a free parameter, which together with the low correlation among the parameters, gives reliable results and constraints on the models.

III. OBSERVATIONAL DATA
Let us introduce now the data that will be used to test and comapre the models (8) and (9). These sets of data include the largest recent catalogue of Type Ia supernovae (SNe Ia), the so-called Pantheon sample [22], baryon acoustic oscillations (BAO) data [23,31], estimations of the Hubble parameter H(z) [32] and parameters from the Cosmic Microwave Background radiation (CMB) [33,34].
Here we use the technique of minimising the likelihood, where we assume a Gaussian distribution for the free parameters: The Pantheon SNe Ia catalogue [22] includes n SN = 1048 data points with redshifts 0 < z i ≤ 2.26 and their corresponding distance moduli µ obs i . Then, the theoretical models are compared with the data by calculating the theoretical value of the distance modulus µ th (z; Ω 0 m , λ i ) for each set of the free parameters: where λ i are the free parameters of the theoretical model and D L (z; Ω 0 m , λ i ) is the free luminosity distance, which is given by: Hence, χ 2 function yields: Here C SN is the 1048 × 1048 covariance matrix [22]. For any set of the model parameters (16) we solve the system of equations provided in (13) and (14), obtaining the Hubble parameter H(z), and consequently the luminosity distances (19) and the distance moduli (18). We also marginalise the χ 2 SN function over the nuisance parameter H 0 [18][19][20][21].
Baryonic Acoustic Oscillations (BAO) are provided by the analysis of galaxy clustering and the following two magnitudes can be compared with the observational data [23]: where whereas r s (z d ) is the comoving sound horizon at the end of the baryon drag era z d , which corresponds to a peak in the correlation function of the galaxy distribution. As in previous works (see Refs. [18,19]), here we use 17 BAO data points for d z (z) and 7 data points for A(z) from Refs. [31] estimated for galaxy clusters with mean redshifts z = z i and represented in Table I.  , which shows a dependence on the Hubble parameter r s (z d ) ∼ H −1 0 and leaves d z (z) Hubble free. Then, the χ 2 function for the BAO fits (21) is where ∆d, ∆A are vector columns with elements ∆d i = d obs . . ), C d and C A are the covariance matrices for correlated BAO data [31] described in Ref. [20].
In addition, the Hubble parameter H(z) data are given by N H = 31 data points H obs (z i ) from Refs. [32] for redshifts 0 < z < 2, whose χ 2 function yields: χ 2 function Here we use only data [32] estimated by the method of differential ages (cosmic chronometers), where the values for the Hubble parameter at different redshifts are estimated at differential ages ∆t of galaxy clusters with certain differences ∆z of redshifts. These estimations are not correlated with the BAO data points [31] at the level (23).
Finally, we will also use the CMB parameters for testing our models. Unlike the datasets coming from SNe Ia, BAO and H(z) observations, measured for 0 < z ≤ 2.36, the CMB observational parameters are related with the photon-decoupling epoch z * = 1089.90 ± 0.25 [24,26]. Hence, as we are dealing with high redshifts here, the radiation density is not negligible and enters in the equations through the radiation-matter ratio X r = ρ 0 r /ρ 0 m in the form (15). We use the CMB parameters released by Planck [24,25] in the following form [33,34]: where the comoving sound horizon r s at z * is calculated as The current baryon fraction Ω 0 b is considered as the nuisance parameter and it is marginalized over ω b = Ω 0 b h 2 and H 0 in the χ 2 CMB function The following data are provided in [34] from Planck collaboration [25]: x P l = R P l , ℓ P l A , ω P l b = 1.7448 ± 0.0054, 301.46 ± 0.094, 0.0224 ± 0.00017 (27) which are given with free amplitude for the lensing power spectrum. The covariance matrix C CMB = C ij σ i σ j , C 12 = 0.53,C 13 = −0.73,C 23 = −0.42 and other details are described in [34] and also in Refs. [18,19].

IV. RESULTS AND DISCUSSION
Here the above SNe Ia, BAO, H(z) and CMB datasets are used to constrain the models (8) and (9) described in section II, through the analysis of the parameter space (16) in order to obtain the best fit and the confidence regions for each of the free parameters. The CMB observations (26) with small errors σ i (27) produce the most strict limitations in the parameter space in comparison with other data. For that reason, we analyse separately the χ 2 function obtained after fitting the free parameters with SNe Ia, BAO and H(z) data at redshifts 0 < z ≤ 2.36 [18,19]: Whereas we estimate the total χ 2 tot separately: where χ 2 CMB corresponds to redshifts near z * ≃ 1100.
Let us start by calculating the χ 2 Σ3 function (28) for the power-law model (8), with the free parameters as given in (16). The results are shown in Fig. 1, particularly the A − α contour plot is depicted in the top-left panel, where we have minimised the χ 2 function over the other parameters and have calculated the difference among the absolute minimum and its variation as a function of A − α: Here m abs Σ3 = min all χ 2 Σ3 is the absolute minimum of χ 2 Σ3 over all its parameters A, α, Ω 0 m , B, β which in this case is m abs Σ3 ≃ 1085.35. In the top-left panel of Fig. 1 we depict the two-parameter distribution ∆χ 2 Σ3 (A, α), where blue lines represents the values of ∆χ 2 Σ3 , as indicated in the figure. The dependence of min χ 2 Σ3 on A and α is very weak and the whole depicted area in the A − α plane lies in the 1σ confidence region. At the bottom panels, this weak dependence is also shown for the corresponding one-parameter distributions, where is minimised over all the other parameters.
As shown in Fig. 1, the minimum value for ∆χ 2 Σ3 in the A − α plane lies within the area with large negative values of A and α ∼ 1 (see the top-left panel), so for convenience, we can fix the value for α without loss of information and effectiveness when minimising χ 2 for this model, Indeed, the minimum for χ 2 Σ3 under the assumption (31)  By fixing the value for α as given in (31), the remaining free parameters are (recall we have marginalise over H 0 ): and its EoS (8) is reduced to The top-right panel of Fig. 1 shows the Ω 0 m − B plane of the two-parameter distribution min χ 2 Σ3 (Ω 0 m , B) (minimised over the other parameters) for the model (8) for a varying α (blue filled contours) and for the restricted case α = 1 (green contours). One can see that the contours of 1σ, 2σ and 3σ confidence regions for both cases do not differ much. This similarity is the most striking, when we compare the one-parameter distributions min other χ 2 Σ3 (Ω 0 m ) in the panel below: the corresponding green and blue dash-dotted lines practically coincide. The ΛCDM model is also depicted for comparison with our model (the black dashed line). Note that the power-law model (8) transforms into the ΛCDM model for Ω x = const = Ω Λ , which corresponds to the particular case A = B = 0. In addition, the difference among the cases α = 1 and α ∈ R becomes more remarkable in the bottom-right panel, where min χ 2 Σ3 depends on B. As shown in Fig. 1 and Table II, both parameters {A, B} are unbounded, since the function χ 2 Σ3 extends its 1σ region up to A → −∞ and B → +∞, respectively. Hence, the model (33) with α = 1 provides effectively the same results as the general model (8) with α ∈ R, since there is no correlation among α and the other parameters, as shown in Fig. 1. From here on, we assume α = 1. In addition, motivated by the behaviour of χ 2 Σ3 in the parameter space, we redefine the parameters as A * and B * , which are related to A and B by: In Fig. 2 we investigate in detail the EoS (33) with α = 1 in the Ω 0 m − A * , Ω 0 m − B * and β − B * planes including the CMB data (26), (27): the corresponding 1σ, 2σ, 3σ contour plots (top panels) are depicted for the function (29) χ 2 tot = χ 2 Σ3 + χ 2 CMB by red lines, the red diamonds show the local minimum points of χ 2 tot . Green filled contours and green dash-dotted lines in all panels correspond to the function χ 2 Σ3 (for the Ω 0 m − B * plane these contours were shown in Fig. 1).
As mentioned above, the CMB observational data, given in (26) and (27), impose the most severe constraints, specially for the parameter Ω 0 m , as obtained after analysing the corresponding χ 2 tot = χ 2 Σ3 + χ 2 CMB (see Fig. 2). This is connected to the observational data R = 1.7448 ± 0.0054, which is proportional to Ω 0 m . One can see that the function (29) has the additional local minimum at Ω 0 m ≃ 0.282, A * ≃ −0.55, B * ≃ 0.38, but is not the global minimum as shown in Fig. 2).
In the middle row of Fig. 2, the one-parameter distributions of the type χ 2 Σ3 (p j ) and χ 2 tot (p j ) are depicted for the 4 parameters p j = Ω 0 m , A * , B * and β. The corresponding likelihoods L Σ3 (p j ) and L tot (p j ) are shown in the bottom panels. These functions are obtained for L tot as follows: where the function is marginalised over all the other free parameters, being m abs tot the absolute minimum for χ 2 tot .
In the middle-left and bottom-left panels we compare these results with the corresponding distributions for the ΛCDM model for χ 2 Σ3 (Ω 0 m ), L Σ3 (Ω 0 m ) (black dashed lines) and for χ 2 tot (Ω 0 m ), L tot (Ω 0 m ) (brown lines). One can see that for the model (33) with α = 1, the absolute minimum for χ 2 Σ3 is essentially lower than the corresponding value for the ΛCDM model, but it is not true for χ 2 tot , when including the CMB data. In addition, we should also mention the the optimal values for A and B go to −∞ and +∞ respectively. Note that the domain B > 0 corresponds to negative viscosity ζ in Eq. (2). The values for the absolute minimum and the best fits (with 1σ errors) of the free model parameters for the model (33) are given in Table II. The results of the logarithmic model (9) are also included. The best fits and 1σ errors are calculated via the distributions χ 2 (p j ) or L(p j ).
On the other hand, the logarithmic model (9) behaves in another way and shows better results (see Table II and in Fig. 3). The minimums for χ 2 Σ3 and χ 2 tot are essentially smaller than in the ΛCDM model and the power-law case. Indeed, in Table II we can compare, for example, min χ 2 tot ≃ 1084.05 for the logarithmic model with the corresponding ΛCDM minimum 1089.03. Moreover, unlike the power-law scenario, these minimums in the logarithmic model (9) are achieved at finite values of A and B with B < 0, corresponding to positive viscosity ζ. FIG. 3: Contour plots for the logarithmic model (9). As above, we depict the confidence regions when considering χ 2 Σ3 (filled contours) and for χ 2 tot (red lines). In the bottom panels, the distributions for Ωm − B * , β * − B * and Ω * − A * are also depicted.
for the logarithmic model (9). The points of minimum are labeled as blue circles for χ 2 Σ3 and as red diamonds for χ 2 tot . These functions behave non-trivially in some domains of the parameter space as can be seen in Fig. 3. In particular, in the contour plots for β * − B * and for Ω * − A * (the top-center and top-right panels) the borders of the 1σ and 2σ confidence region are not regular, whereas the 3σ domain lies beyond. These unusual behaviour can be also seen in the corresponding one-parameter distributions min other χ 2 (p i ) of Fig. 3, particularly these functions decrease at large negative values of B * and A * and positive values of B * . However, here the local minimum coincide with the global minimum for the χ 2 function, unlike the power-law model.
Note that the narrow peak of χ 2 tot for the ΛCDM model in the bottom-left panels of Figs. 2 and 3 is connected with the CMB parameter R ∼ Ω 0 m in Eqs. (26) and (27), and with the number of free parameters for the flat ΛCDM model, that is Ω 0 m and the nuisance parameter H 0 , as shown in the FLRW equation: Nevertheless, the logarithmic model shows a good behaviour in comparison to ΛCDM, as shown in Fig. 3 and Table  II.   (33) with α = 1 and the logarithmic model (9), when considering χ 2 Σ3 = χ 2 SN + χ 2 BAO + χ 2 H (the Pantheon SNe Ia, BAO and H(z) data) and for χ 2 tot = χ 2 Σ3 + χ 2 CMB (including the CMB data) in comparison with the flat ΛCDM model. The table also shows the min χ 2 and the 1σ errors for the model parameters. Here β = sinh(β * ), similar to the relations (34).

V. CONCLUSIONS
Along the paper we have considered two cosmological scenarios where dark energy is assumed to be described by a viscous fluid, through bulk viscosity, what leads to an effective pressure that can explain the late-time accelerating expansion. For that, and inspired on some hydrodynamics considerations, we have explored two different EoS for viscous dark energy: the power-law model (8), precisely, its variant (33) with α = 1 and the logarithmic model (9). By using data from Supernovae Ia, BAO, H(z) measurements and CMB, we have analysed the viability of these scenarios and compared to ΛCDM model.
Our analysis shows that the power-law model (33) behaves well, also in comparison with ΛCDM model, when considering the restricted set of observational data that excludes CMB data, as depicted in Figs. 1 and 2 and summarised in Table II). Actually, the model (33) provides a slightly lower minimum when considering χ 2 Σ3 than ΛCDM model, but higher errors for Ω m , and weak constraints on the free parameters {A, B}, since they show no bound from above/below, which may lead to a negative viscosity, as given in (2), particularly the model (33) achieves the best χ 2 values at the non-physical domain B → +∞ that corresponds to large negative viscosity ζ in Eq. (2). The other free parameter of the model α is very well constrained at α ∼ 1 and shows no correlations with the other free parameters, such that we have assumed α = 1 for a great part of our calculations, as depicted in Figs. 1 and 2. However, when assuming CMB data, the model (33) provides larger values for the minimum of χ 2 tot , but similar to the one given by ΛCDM model. This means that the viscous term as a power law behaves well at late-times, but shows some issues when increasing the covered region of the cosmological evolution. In addition, the same problems with the parameters {A, B} remain in this case (see Table II), and one can not obtain better constraints for both.
Unlike the power-law model (33), the logarithmic model (9) has no these drawbacks: it provides essentially lower values for min χ 2 Σ3 and min χ 2 tot , which are achieved at reasonable values of the free parameters, and the constraints on each parameter are well defined and limited (see Table II). The values of the minimums for χ 2 show that the model (9) fits better every set of observational data, in comparison to the power-law model and the ΛCDM model. However, despite the contour plots and statistical distributions in Fig. 3 show a well defined 1σ region, the errors increase much when one goes to confidence regions of upper σ, which may be interpreted as some lack of information on the free parameters. For instance, the analysis of χ 2 tot provides the best value for B as B = sinh B * = −0.465 +0.21 −0.275 , which corresponds to positive optimal values of viscosity ζ, strongly depending on H because for negative β, but if one increases the confidence region, B may take values that lead to a negative viscosity, and unconstrained model. In any case, the logarithmic model (9) seems to provide a very well description of the cosmological evolution at any redshifts, that is also when including CMB data.
Hence, we have explored the existence of a viscous dark energy fluid by using the last observational data coming from different sources and by considering some theoretical models for the viscosity terms that play a role in other areas of hydrodynamics. Our results show that the right viscosity term can provide better fits in comparison to other models, such that one should keep analysing this possibility by going beyond the analysis of the cosmological background evolution.