Bayesian analysis of bulk viscous matter dominated universe

In our previous works, we have analyzed the evolution of bulk viscous matter dominated universe with a more general form for bulk viscous coefficient, $\zeta=\zeta_{0}+\zeta_{1}\frac{\dot{a}}{a}+\zeta_{2}\frac{\ddot{a}}{\dot{a}}$ and also carried out the dynamical system analysis. We found that the model reasonably describes the evolution of the universe if the viscous coefficient is a constant. In the present work we are contrasting this model with the standard $\Lambda$CDM model of the universe using the Bayesian method. We have shown that, even though the viscous model gives a reasonable back ground evolution of the universe, the Bayes factor of the model indicates that, it is not so superior over the $\Lambda$CDM model, but have a slight advantage over it.


Introduction
Many observations lead to the conclusion that the present universe is accelerating [1,2,3,4,5,6]. The reason for this acceleration was attributed to the dominant presence of a new cosmic component called dark energy. The ΛCDM model came out as the most successful one for explaining this late time acceleration of the universe. In this model the cosmological constant is being considered as the dark energy. But the model is plagued with severe drawbacks. The foremost among them is the cosmological constant problem and is about the discrepancy between the observed and predicted values of the dark energy density, which is of the order of 120. The other is the coincidence problem, the mysterious coincidence between the energy densities of the dark energy and dark matter component during the current epoch of the universe in spite of their completely different evolution history. This motivates a large class of models with varying dark energy density [7,8,9,10,11,12,13,14]. Perfect fluid models like Chaplygin gas model [15,16] would be an alternate suggestion, due to their ability to explain both the deceleration and late acceleration by a single cosmic component, which thus effectively leads to a unification of the dark matter and dark energy sectors. There were also attempts to study this phenomenon by modifying the geometry part of the gravity theories, like f (R) gravity [17,18,19], f (T ) gravity [20,21], Gauss-Bonnet theory [22], Lovelock gravity [23], Horava-Lifshitz gravity [24], scalar-tensor theories [25], braneworld models [26] etc.
As in the case of the Chaplygin gas model, another possibility of the unified description of both dark energy and dark matter arises in the dissipative fluid models. It has been shown that the early inflationary period of the universe can be due to the presence of an imperfect fluid with bulk viscosity [27,28,29,30,31]. This motivates the study of the disspiative cosmologies in the context of the late acceleration of the universe [32,33,34,35,36]. In [33], by considering a single cosmic component, which is the dark matter with bulk viscosity, ζ(ρ) = αρ m with α and m being constants, the authors have shown that, the universe can make a transition from a decelerating phase to a late accelerating phase and ultimately to a de Sitter epoch. Inspite of this good background evolution, the model have come across with some negative aspects while analysing the structure formation. For instance, in reference [33] with ζ = αρ −0.4 the authors have shown that the density perturbation would rapidly be damped out, which adversely affect mainly the CMBR. It may be due to the power factor of the density −0.4, which was obtained by constraining the model with old supernovae luminosity data by Riess and others [1]. At around the same time, in reference [35], the authors have considered a constant bulk viscous dark matter dominated universe, with 0 < ζ < 3 and predicts that the universe began with a Big-Bang, followed by a decelerated expansion epoch and later transition into an accelerated epoch. Later these authors [36] extended their model by taking varying bulk viscosity of the form ζ = ζ 0 +ζ 1 H and found that it shows a background evolution close to that of the standard ΛCDM model. In [37], the authors have shown that the data from Planck CMB observation and different LSS observations prefer small but non-zero amount of viscosity in cold dark matter fluid.
In [38], we have analyzed the evolution of bulk viscous matter dominated universe with a more general form for bulk viscous coefficient, ζ = ζ 0 + ζ 1ȧ a + ζ 2ä a , such that the viscosity depends on both the velocity and acceleration of the expansion, where a is the scale factor of expansion of the universe. We have studied the model by contrasting it with Union compilation of supernovae data and found that it predicts the transition to the late accelerating phase. A similar work has been done in reference [39], however, in constraining the parameters (ζ 0 , ζ 1 , ζ 2 ), the authors fixed either ζ 1 or ζ 2 as zero while evaluating the other. In a later work [40] we found that a general evolution of the universe with all the conventional phases (radiation epoch, matter dominated epoch and later acceleration) can be described without any causality violation only when the bulk viscous coefficient is a constant, ζ = ζ 0 .
In the present paper we concentrate on the Bayesian analysis of the bulk viscous model corresponding to the different forms of the bulk viscous coefficient. We intend to compare these models in general with the standard ΛCDM model and also among themselves. Bayesian model comparison method is commonly used in the context of cosmological model selection [41,42,43,44,45,46,47]. More details regarding Bayesian model comparison are discussed in the section 2.
The paper is organized as follows. In Section 2, we discuss the basic details regarding the Bayesian analysis. In section 3, we introduce the bulk viscous model of the universe. In section 4, we extract the bulk viscous parameters corresponding to different cases of our model and discuss the results of the Bayesian analysis of the viscous model and finally, in section 5, we present our conclusion.

Bayesian model comparison
Various models have been proposed to interpret the cosmological observational data which eventually add to our understanding of the evolution of the universe. So there exist, in fact many models explaining the expected evolution of the universe. Contrasting these models among themselves to select the better ones is essential for understanding the true evolution of the universe. Bayesian statistical approach [47,48,49] is an effective tool to compare the new models with the standard ΛCDM model and also among themselves. The basic approach of this method is originated from the theory of random variables. In general, the relative merit of a random variable can be obtained by calculating the basic probability of it among the ensemble of values obtained theoretically or through repeated observations. But in cosmology repeated observations are virtually impossible. Here, what one can often do is to form hypothesis or a theory. For making the decision regarding the viability of such a proposed theory one have to assign certain probability to it in contrast to other theories existing for the same purpose. It is in this stage the Bayesian theory help us, so as to assign probability for a certain hypothesis by considering the observational data already available to us. Due to the acquisition of more data, one can in fact adjust the plausibility of the hypothesis using Bayesian theorem. This method have been adopted by many in the past, for instance, Jaffe [50] and Hobson et al. [51] have analysed the relative merits of certain cosmological models. Also John and Narlikar [47] have compared a simple cosmological model with scale factor a(t) ∝ t with standard and inflationary models of the universe. In many models one does not have a prior knowledge about the model parameters for assigning the corresponding probability and in such cases one often starts with a flat prior for the parameter.
According to Bayes's theorem [52], the posterior probability p(H i |D, I) of a hypothesis H i , given the data D and assuming any other background information I to be true, is given as, where p(H i |I) is the prior probability, i.e., the probability of H i given I is true and p(D|H i , I) is the likelihood for the hypothesis H i , which is the probability for obtaining the data D provided the hypothesis H i and I are true. The factor p(D|I) helps in normalization.
In Bayesian model comparison, we take the ratios between the posterior probabilities for different models. Let M i and M j be the two models which we need to compare, then using Bayes theorem the ratio between their posterior probability O ij can be written as, Since p(D|M i , I) for the data D is the likelihood for the model M i , we re-notate it with L(M i ), then the equation (2) becomes, If the background information I does not give any preference to a model over any other, then the prior probabilities becomes equal, so that, where B ij is called the Bayes factor and is thus the ratio of the likelihood of the two models. This factor helps to compare two models with reference to their power in predicting the given data, hence it can be taken as a summary of the evidence provided by the data in favor of one model over the other [53]. If 1 < B ij < 3, then the model M i is not worth more than a bare mention. If 3 < B ij < 20, the strength of evidence of the model M i is positive. If 20 < B ij < 150, the evidence is strong and if B ij > 150, it is very strong [48].
For a model having one or more free parameters such as α, β,...etc, it's likelihood L(M i ) can be evaluated as where p(α, β, ...|M i ) is the prior probability for the set of parameter values α, β, ... for the model M i to be true and L i (α, β, ...) is the likelihood for the combination of the parameters in the model and is usually taken as [48], where χ 2 i (α, β, ...) is the conventional χ 2 -function. If we assume that the model M i has two parameters α and β having flat prior probabilities in some range, [α, α + ∆α] and [β, β + ∆β], respectively, then the likelihood of the model is given as where is called the marginal likelihood for the parameter α. Similarly one can also find the marginal likelihood for the parameter β. This can be extended with the number of parameters. For instance, if we have three parameters, α, β, γ, then the marginal likelihood of a parameter, say, α, can be evaluated as,

Bulk viscous FLRW Universe
We can now proceed towards the Bayesian comparison of the different bulk viscous models of the universe. First, we briefly discuss the bulk viscous model based on the reference [38]. The basic model we consider is a spatially flat universe, dominated with bulk viscous matter, described by the Friedmann equations, where we have taken 8πG = 1, ρ is the density of the matter component of the universe, a(t) is the scale factor and an overdot represents the derivative with respect to cosmic time t. We are neglecting the radiation components as it doesn't have any decisive role in the late evolution of the universe. The conservation equation of the cosmic component is then, According to Eckart theory [54], the effective pressure of the bulk viscous fluid is given as where P is the normal kinetic pressure and ζ is the coefficient of bulk viscosity. The matter component in the late universe is non-relativistic hence it is usually taken as, P = 0. Then the contribution to effective pressure is only due to the negative viscous pressure. The coefficient ζ is basically a transport coefficient, hence it would depend on the dynamics of the cosmic fluid. We consider the most general form for the bulk viscous coefficient ζ, which is a linear combination of the three terms [38,39,40,55,56], the first term ζ 0 , is a constant, the second is proportional to the velocity and the third is proportional to the acceleration in the expansion of the universe. Using the Friedmann equations and conservation equation, we can obtain the expression for the Hubble parameter H for the general form of ζ as [38], where we define the dimensionless bulk viscous parametersζ 0 = 3ζ 0 H 0 ,ζ 1 = 3ζ 1 andζ 2 = 3ζ 2 .
The above equation can be integrated to obtain the scale factor as where t 0 is the present cosmic time andζ 12 =ζ 1 +ζ 2 . For the best estimates of the parameters, it was shown in reference [38] that, in the limit t − t 0 → ∞ the scale factor tends to H 0 (t−t 0 ) , corresponding to the de Sitter phase, while in the limit t − t 0 → 0 the scale factor tends to a(t) → 1 + 3−ζ 12 3−ζ 12 , corresponding to an early decelerated expansion. Thus from the evolution of the scale factor, it is clear that the universe undergoes transition from the deceleration phase to accelerated phase and the transition redshift is found to be around 0.49. The model also predicts the present deceleration parameter around -0.68. However the age predicted was relatively low.
In reference [40], we have done the phase space analysis of the model with bulk viscous matter as the dominating cosmic component. Three choices for the bulk viscous coefficient, (i)ζ = ζ 0 + ζ 1ȧ a + ζ 2ä a , (ii)ζ = ζ 0 + ζ 1ȧ a and (iii)ζ = ζ 0 have been adopted. It was found that for all the three cases, it predicts a prior unstable decelerated epoch, and a later stable accelerating epoch, similar to the de Sitter phase. However, when the radiation component is also taken into account, the model support the conventional evolution of the universe, only for the case ζ = ζ 0 . The other two cases doesn't predicts a prior radiation dominated phase and conventional decelerated matter dominated phase of the universe respectively.

Bayesian Analysis of Bulk viscous models
A fairly detailed description of the method of Bayesian analysis was given in section 2. In this section we are going into the Bayesian analysis of the different bulk viscous models. For this we consider the following cases separately, where viscosity is depending on both the velocity and acceleration of the expansion of the universe.
where viscosity is depending only on velocity of the expansion of the universe apart from a constant additive part ζ 0 . This is equivalent to ζ = ζ 0 + ζ 1 ρ s , with s = 1/2.
3. ζ = ζ 0 , where viscosity is pure a constant where viscosity only has the velocity dependent term and is equivalent to ζ = ζ 1 ρ s , with s = 1/2.
where viscosity is depending on acceleration apart from an additive constant.
The best estimated values of the parameters (ζ 0 ,ζ 1 ,ζ 2 ) corresponding to the cases 1, 2 and 3 are extracted in the reference [38,40] and that of the cases 4 and 5 are extracted using the same set of data and are given in the table 1. The data used is the SCP "Union" Type Ia Supernova data [57] composed of 307 data points from 13 independent data sets and the method used is χ 2 minimization technique.  [38,40] Parameters Bulk viscous models ζ = ζ 1ȧ a ζ = ζ 0 + ζ 2ä ã The χ 2 function is constructed as, where µ k is the observational distance modulus for the k-th Supernova (obtained from the data) with red shift z k , σ 2 k is the variance of the measurement, n is the total number of data and µ t is the theoretical distance modulus for the k-th Supernova with the same redshift z k , which is given as where, m and M are the apparent and absolute magnitudes of the SNe respectively. d L is the luminosity distance and is defined as where c is the speed of light. After obtaining the χ 2 , we evaluate the marginal likelihood, using equation (9), and likelihood, using equation (7), for all the five cases of the model. We kept ΛCDM model as the reference model in order to compare the bulk viscous models and calculate the Bayes factor using equation (4). The marginal likelihood of the parametersζ corresponding to the five cases of bulk viscous models are shown in figures 1, 2, 3, 4 and 5, respectively.  We obtained the Bayes factor for three different priors for the parametersζ 0 ,ζ 1 ,ζ 2 . Since there is no prior information regarding any of the viscous parameters, we assume flat prior for all the three parameters.. Each prior corresponds to a range of values ofζ's for a fixed  likelihood. For choosing the prior forζ's we have plotted exp(−χ 2 /2) as a function of a giveñ ζ i.e.,ζ 0 ,ζ 1 , orζ 2 , by keeping the other two as constants (equal to that corresponding to the minimum of χ 2 ). For instance prior I, corresponds to the range ofζ's for likelihood of about 1 * 10 −70 , prior II corresponds to the range ofζ's for likelihood of about 1 * 10 −80 and prior III corresponds to the range ofζ's for likelihood of about 1 * 10 −90 . The results are tabulated in the Table 2. It can be seen from the table that for the first two cases the Bayes factor depends on the prior, while for the remaining three case such strong dependence on the priors are not evident. For case 2, i.e., when ζ = ζ 0 + ζ 1ȧ a , there is an increase in Bayes factor with prior and it exceeds 3, giving more evidence for its increasing strength. The important thing to be noted here is about the value of the Bayes factor for the respective model. For the cases ζ = ζ 0 and ζ = ζ 1ȧ a , the factor is much less than one. While for other cases the values are above one. As per the standard classification, it can be mentioned that, the models for which the Bayes factor is in between 1 and 3, can have only a very feeble advantage over the standard ΛCDM model, however is not worth more than a bare mention.   The SNe data that we have used contains the magnitude of supernovae in the red-shift range 0.015 < z < 1.55. The data predicts a transition from an early decelerated epoch to the late acceleration at a redshift of around z ∼ 0.5 [38]. Among the full data set, the low redshift data within the range 0.015 < z < 0.5 is often used for deducing the current value of the Hubble parameter. The high redshift data, which were obtained with small interference with the background and with high accuracy, are considered to be the best part of the data. For a critical analysis we have repeated our computation using the high redshift part of the data, corresponding to the red shift range, 0.5 < z < 1.55. For this range, we have extracted the values of the parameters,ζ's corresponding to the above 5 cases using the χ 2 minimization technique. The results are tabulated in table 3.
The Bayes factor of the bulk viscous models corresponding to the five cases for high redshift data are tabulated in the Table 4. Here prior I corresponds to the range ofζ's for likelihood of about 1 * 10 −38 and prior II corresponds to the range ofζ's for likelihood of about 1 * 10 −45 .
The models corresponding to the cases 1 and 5, where the viscous coefficient depends on the acceleration of the expansion, have Bayes factor less than one for both the priors. As a result these two cases are not worth of any mention against the standard model. This indicates that the dependence of viscosity on the acceleration is not so sensitive. Riess et. al. have found that  the magnitude of acceleration is small, since the distance of the high redshift supernovae were on average only 10% − 15% farther than expected in a universe with mass density parameter Ω m ∼ 0.3. Such a small acceleration would not have any observable effect on the transport coefficients like that of viscosity. For cases 2, 3 and 4, the Bayes factors are larger than one and it increases slightly with prior. Among these, the third case is of constant viscosity, while for second and fourth cases, the viscous coefficient depends on the velocity of expansion of the universe. As seen from the Table 4, the Bayes factors for the cases 2, 3 and 4 are all in the range 1 < B ij < 3 and it seems quite difficult to discriminate between them. All these three cases are thus qualified to have bare mention against the the ΛCDM model.
dynamical behavior that the model describes all the conventional phases of the universe and asymptotically tend towards a stable de Sitter epoch, provided the bulk viscosity of the cosmic fluid is a constant. In the present work we have contrasted the bulk viscous model of the universe with the standard ΛCDM model using the method of Bayesian analysis. We have first extracted the viscous parameters corresponding to the following five cases of bulk viscous model, (1) ζ = ζ 0 + ζ 1ȧ a + ζ 2ä a , (2) ζ = ζ 0 + ζ 1ȧ a , (3) ζ = ζ 0 , (4) ζ = ζ 1ȧ a , (5) ζ = ζ 0 + ζ 2ä a , using the "Union" data of Supernovae type Ia. We have obtained the Bayes factor for all the five cases, see table 2. For the full supernovae data set, the results indicate that the model corresponding to case 2, i.e., ζ = ζ 0 +ζ 1ȧ a have a Bayes factor just above 3, and thus have slight advantage over the ΛCDM model compared with other cases. For the model corresponding to cases 1 and 5 , the Bayes factor is just above one and can just have a bare mention, in contrast to the standard model. All other cases, especially the case 3, with constant viscosity, seems to fail in standing against the standard model. For more reliable result, we restrict to supernovae data with relatively high redshifts, z > 0.5, which were obtained with less background interference and hence are more reliable in making predictions regarding the evolution near the transition. The results consequent to this have a marked deviation from the previous one, such that the Bayes factor for the constant bulk viscosity ζ = ζ 0 (case 3) and models corresponding to cases 2 and 4, where the viscous coefficient depends on the velocity of expansion, are having a slight advantage over other cases when compared with the standard ΛCDM model. Since Bayes factors of the cases 2, 3 and 4, are all in the range 1 < B ij < 3, it is difficult to discriminate among themselves. However it was shown in reference [40] that only the case 3 will have asymptotically stable end de Sitter phase. Taking account of this, it can be concluded that, among the cases 2, 3 and 4, which are having almost same Bayes factor, the case 3, corresponding to the one with stable end de Sitter phase, can be preferred over the other cases.