Cosmological Constraints on Ghost Dark Energy in the Brans-Dicke Theory by Using MCMC Approach

By using a Markov Chain Monte Carlo simulation, we investigate cosmological constraints on the ghost dark energy (GDE) model in the framework of the Brans-Dicke (BD) theory. A combination of the latest observational data of the cosmic microwave background radiation data from seven-year WMAP, the baryon acoustic oscillation data form the SDSS, the supernovae type Ia data from the Union2 and the X-ray gas mass fraction data from the Chandra X-ray observations of the largest relaxed galaxy clusters are used to perform constraints on GDE in the BD cosmology. In this paper, we consider both flat and non-flat universes together with interaction between dark matter and dark energy. The main cosmological parameters are obtained as: $\Omega_{\rm b}h^2= 0.0223^{+0.0016}_{-0.0013}$, $\Omega_{\rm c}h^2=0.1149^{+0.0088}_{-0.0104}$ and $\Omega_{\rm k}=0.0005^{+0.0025}_{-0.0073}$. In addition, the Brans-Dicke parameter $\omega$ is estimated as $1/\omega\simeq 0.002$.


I. INTRODUCTION
Accelerating expansion of the Universe [1,2] can be explained either by a missing energy component usually called "dark energy" (DE) with an exotic equation of state, or by modifying the underlying theory of gravity on large scales.
Among various models of DE, the so called ghost dark energy (GDE) has attracted a lot of interests in recent years. The origin of DE in this model comes from Veneziano ghosts in QCD theory [27][28][29][30]. Indeed, the contribution of the ghosts field to the vacuum energy in curved space or time-dependent background can be regarded as a possible candidate for DE [31,32]. The magnitude of this vacuum energy is of order Λ 3 QCD H, where H is the Hubble parameter and Λ QCD is the QCD mass scale. With Λ QCD ∼ 100M eV and H ∼ 10 −33 eV , Λ 3 QCD H gives the right order of magnitude ∼ (3 × 10 −3 eV ) 4 for the observed dark energy density [31]. The advantages of GDE model compared to other DE models is that it is totally embedded in standard model and general relativity, therefore one needs not to introduce any new parameter, new degree of freedom or to modify gravity. The dynamical behavior of GDE model in flat universe have been studied [33]. The study was also generalized to the universe with spacial curvature [34]. The instability of the GDE model against perturbations was studied in [35]. In [36,37] the correspondence between GDE and scalar field models of DE were established. In the presence of bulk viscosity and varying gravitational constant, the GDE model was investigated in [38]. Other features of the GDE model have been investigated in Refs. [39][40][41][42][43]. The cosmological constraints on this model have been considered by some authors [33,43,44]. Scalar tensor theories have been reconsidered extensively, recently. One important example of the scalar tensor theories is the BD theory of gravity which was introduced by Brans and Dicke in 1961 to incorporate the Mach's principle in the Einstein's theory of gravity [45]. This theory also passes the observational tests in the solar system domain [46]. In addition, BD theory can be tested by the cosmological observations such as the cosmic microwave background (CMB) and large scale structure (LSS) [47][48][49][50][51]. Since the GDE model have a dynamic behavior, it is more reasonable to consider this model in a dynamical framework such as BD theory. It was shown that some features of GDE in BD cosmology differ from Einstein's gravity [52]. For example, while the original DE is instable in all range of the parameter spaces in standard cosmology [35], it leads to a stable phase in BD theory [53]. In the framework of BD cosmology, the ghost model of DE has been studied [52]. It is also of great interest to see whether the GDE model in the framework of the BD theory is compatible with observational data or not.
In this paper, cosmological constraints on GDE in the BD theory (GDEBD) [52] theory is performed by using the Marko Chain Monte Carlo (MCMC) simulation. The used observational datasets are as follows: cosmic microwave background radiation (CMB) from WMAP7 [54], 557 Union2 dataset of type Ia supernova [55], baryon acoustic oscillation (BAO) from SDSS DR7 [56], and the cluster X-ray gas mass fraction from the Chandra X-ray observations [57]. To put the constraints, the modified CosmoMC [58] code is used.
The organization of this paper is as follows. In section II, we review the formalism of the GDE in the framework of Brans-Dicke cosmology. In section III the methods which are used in this paper to analyze the data are introduced. Section IV contains the results of the MCMC simulation and we conclude our paper in section V.

II. INTERACTING GHOST DARK ENERGY IN THE BRANS-DICKE THEORY IN A NON-FLAT UNIVERSE
Let us first review the formalism of the interacting GDE in the framework of BD theory in a non-flat universe [52]. The action of the BD theory in the canonical form may be written [59] where R is the Ricci scalar and φ is the BD scalar field. Varying the action with respect to the metric g µν and the BD scalar field φ, yields where T M µν stands for the energy-momentum tensor of the matter fields. The line element of the Friedmann-Robertson-Walker (FRW) universe is where a(t) is the scale factor, and k is the curvature parameter with k = −1, 0, 1 corresponding to open, flat, and closed universes, respectively. Nowadays, there are some evidences in favor of closed universe with a small positive curvature (Ω k ≃ 0.01) [61]. Using metric (4), the field equations (2) and (3) reduce to where H =ȧ/a is the Hubble parameter, ρ D and p D are, respectively, the energy density and pressure of DE, and ρ m is the pressureless matter density which contains both dark matter (DM) and baryonic matter (BM) densities i.e. ρ m = ρ c + ρ b where ρ c and ρ b are the energy densities of dark matter and baryonic matter respectively. To be more general and because of some observational evidences [62,63], here we propose the case where there is an interaction between GDE and DM. In this case the semi-conservation equations reaḋ where Q represents the interaction term between dark matter and dark energy and here we assume that the baryonic matter is conserved separately. We assume Q = 3ξH(ρ m + ρ D ) with ξ being a constant. Such a choice for interacting term implies the the DE and DM component do not conserve separately while the total density is still conserved throughρ where ρ = ρ D + ρ m and P = P D . The ghost energy density is proportional to the Hubble parameter [31] ρ D = αH.
where α > 0 is roughly of order Λ 3 QCD and Λ QCD is QCD mass scale. Taking into account the fact that Λ QCD ∼ 100M eV and H ∼ 10 −33 eV for the present time, this gives the right order of magnitude ρ D ∼ (3 × 10 −3 eV) 4 for the ghost energy density [31].
Since the system of equations (5-7) is not closed, we still have another degree of freedom in analyzing the set of equations. As usual we assume the BD scalar field φ has a power law relation versus the scale factor, An interesting case is when ε is small whereas ω is high so that the product εω results of order unity [64,65]. In section IV we will consider the ωε = 1 condition for constraining the model by observational data. This is interesting because local astronomical experiments set a very high lower bound on ω [66]; in particular, the Cassini experiment implies that ω > 10 4 [46,48]. Now we take the time derivative of relation (13). We arrive aṫ Combining Eqs. (13) and (14) with the first Friedmann equation (5), we get As usual the fractional energy densities are defined as where Using (12) we can rewrite equation (19) as Based on these definitions, equation (15) can be rewritten as where Ω m = Ω c + Ω b and we have defined Next we take the time derivative of (15), after using (22), we finḋ Combining the above equation with Eqs. (8) and (12), we obtain the EoS parameter as The first and second derivatives of the distance can be combined to obtain the acceleration parameter q. It was shown that the zero redshift value of q 0 , is independent of space curvature, and can be obtained from the first and second derivatives of the coordinate distance [67]. It was argued that q 0 , which indicates whether the universe is accelerating at the current epoch, can be obtained directly from the supernova and radio galaxy data [67]. The acceleration parameter is given by Using (24)the acceleration parameter (26) is obtained as Finally, we obtain the equation of motions of GDE in BD theory. For this purpose, we first take the time derivative of relation (21). We findΩ Substituting q from (27) into equation (28) and using relation Ω ′ In the remaining part of this paper we will constrain the GDEBD model by using the most recent observational date in the three different physical models: model I which is the GDEBD model in a flat universe (ξ = 0 and Ω k = 0), model II is the interacting GDEBD model in a flat universe (ξ = 0 and Ω k = 0) and finally model III is the interacting GDEBD model in a non-flat universe (ξ = 0 and Ω k = 0).

III. DATA FITTING METHOD
In this section we discuss the data fitting method in the Markov Chain Monte Carlo (MCMC) simulation to estimate the parameters of the model in section II using cosmological data.
To get the best fit values of the relevant parameters, the maximum likelihood method is used. The total likelihood function L total = e −χ 2 tot /2 is defined as the product of the separate likelihood functions of uncorrelated observational data with where SNIa stands for type Ia supernovae, CMB for cosmic microwave background radiation, BAO for baryon acoustic oscillation and gas stands for X-ray gas mass fraction data. Best fit values of parameters are obtained by minimizing χ 2 tot . In this paper we use the cosmic microwave background radiation data from seven-year WMAP [54], type Ia supernovae data from 557 Union2 [55], baryon acoustic oscillation data from SDSS DR7 [56], and the cluster X-ray gas mass fraction data from the Chandra X-ray observations [57]. In the rest of this section we discuss each χ 2 i in detail.
To obtain χ 2 CMB , we use seven-year WMAP data [54] with the CMB data point (R, l A , z * ). The shift parameter R, which parametrize the changes in the amplitude of the acoustic peaks is given by [68] where z * is the redshift of recombination (see (36)), c is the speed of light in vacuum, Ω m0 is the present value of the matter density parameter (a "0" subscript shows the present value of the related quantity), and E(z) ≡ H(z)/H 0 .
In addition, the acoustic scale l A , which characterizes the changes of the peaks of CMB via the angular diameter distance out to recombination is defined as [68] l The comoving distance r(z) is defined and the comoving sound horizon distance at recombination r s (z * ) is given by in terms of the sound speed c s (a), defined by The seven-year WMAP observations gives Ω γ0 = 2.469 × 10 −5 h −2 and Ω b0 = 0.02258 +0.00057 −0.00056 [54]. The redshift of recombination z * is obtained by using the fitting function proposed by Hu and Sugiyama [69] where Then one can define χ 2 CMB as χ 2 where C −1 CMB is the inverse covariant matrix. To obtain χ 2 SNIa , the SNIa Union2 data [55] is used which includes 577 type Ia supernovae. The expansion history of the universe H(z) can be given by a specific cosmological model. To test this model, we can use the observational data for some predictable cosmological parameter such as luminosity distance d L . Assume that the Hubble parameter H(z; α 1 , ..., α n ) is used to describe the Universe, where parameters (α 1 , ...α n ) are predicted by a theoretical cosmological model. For such a theoretical model we can predict the theoretical 'Hubble-constant free' luminosity distance as where E ≡ H/H 0 , z is the redshift parameter, and Then one can write the theoretical modulus distance where µ 0 = 5 log 10 (cH −1 0 /M pc) + 25. On the other hand, the observational modulus distance of SNIa, µ obs (z i ), at redshift z i is given by where m obs and M are apparent and absolute magnitudes of SNIa respectively. Then the parameters of the theoretical model, α i s, can be determined by a likelihood analysis by definingχ 2 SNIa (α i , M ′ ) in (30) as where the nuisance parameter, M ′ = µ 0 + M , can be marginalized over as Here we should mention an important point about using supernovae data to constrain the Brans-Dicke theories which have a varying gravitational coupling constant. Variations of gravitational coupling constant and apparent magnitude of supernovae are correlated as follows. The luminosity L of a supernova is powered by Nickel-35 mass which is proportional to the Chandrasekhar mass Moreover, the luminosity distance d L is the integral over the inverse Hubble parameter, which is proportional to G −1/2 . Therefore, the apparent magnitude m obs m obs = −2.5LogL + 5Logd L varies with a change in the gravitational coupling constant ∆G as On the other hand, in the Brans-Dicke theory we haveĠ where φ is the Brans-Dicke scalar field. In the slow roll approximation we can writė where q is the decelaration parameter. The average value of q between today and z ∼ 1.2 (the redshift when the SN measurement are probing the dark energy) is of order unity, and by using H = d ln a/dt, we can write By using Eqs. (46) and (49) we obtain In Union2 data set, the redshift interval is between 0.51 and 1.12, i.e. ∆ ln a ∼ 1, with the systematic error of order 0.03 in the measurement of apparent magnitude. Therefore, the systematic error can induce a bias roughly of order 0.3 on parameter ǫ, which is three order of magnitudes larger than the statistical errors, as we will discuss in the next section. Therefore, in order to constrain ǫ with a higher precision, we combine the supernovae data with other cosmological data sets as follows. For more detailed discussion on possible evolution of the gravitational constant from cosmological type Ia supernovae see [60]. The baryon acoustic oscillation data from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) [56] is used here for constraining model parameters. The data constrain parameter d z ≡ r s (z d )/D V (z), where r s (z d ) is the comoving sound horizon distance (see (34)) at the drag epoch (where baryons were released from photons) and D V is given by [70] The drag redshift is given by the fitting formula [71] z d = 1291(Ω m0 h 2 ) 0.251 1 + 0.659(Ω m0 h 2 ) 0. 828 where Then we can obtain χ 2 BAO by and its covariance matrix is given by [56] The ratio of X-ray gas mass to the total mass of a cluster is defined as the X-ray gas mass fraction [57]. The model fitted to the ΛCDM model is [57] f ΛCDM The elements in equation (60) The angular correction factor A is caused by changes in angle for the alternative theoretical model θ 2500 compared to θ ΛCDM 2500 , where η = 0.214 ± 0.022 [57] is the slope of the f gas (r/r 2500 ) data within the radius r 2500 (r 2500 is the radius of the gas core in Mpc/h units).
The bias factor b(z) in equation (60) contains information about the uncertainties in the cluster depletion factor b(z) = b 0 (1 + α b z) and the parameter γ accounts for departures from the hydrostatic equilibrium. The function s(z) = s 0 (1 + α s z) denotes the uncertainties of the baryonic mass fraction in stars with a Gaussian prior for s 0 , with s 0 = (0.16 ± 0.05)h 0. 5 70 [57]. The factor K describes the combined effects of the residual uncertainties, such as the instrumental calibration. A Gaussian prior for the 'calibration' factor is considered as K = 1.0 ± 0.1 [57].
Then χ 2 gas is defined as [57] with the statistical uncertainties σ fgas (z i ) and At the end of this section we should assert that the data points parameters of the CMB and BAO data sets which we use in this paper are the best fit values for ΛCDM and the error estimates are also based on the ΛCDM model. Therefore they are not completely accurate in this application. However they are the only parameters which we have to constrain our model.

IV. RESULTS
Finally we apply a Markov chain Monte Carlo simulation on the parameters of the GDEBD model by using the publicly available CosmoMC code [72]. From table IV one can see that the main cosmological parameters Ω b h 2 , Ω c h 2 , Ω DE and Ω k in all three models are compatible with the results of the ΛCDM model [61]. In addition, in the presence of interaction between DM and DE, the parameter ǫ decreases and so the Brans-Dicke parameter ω increases. The best fit value of the parameter ǫ = 1/ω in all three models is compatible with the results of other cosmological constraining works (however see the discussion following Eq. (43)). For example in [47] the authors by using the CMB temperature and polarization anisotropy data, found 1/ω ≃ 0.001. Wu and Chen in [50,51] by using the five-year WMAP and SDSS data obtained ω > 97.8. For other cosmological constraints on the BD theory see [48,49,[73][74][75][76]. This estimated value is also compatible with the results of the solar system tests of the scalar-tensor theories such as the Cassini experiment where it has been obtained ω > 10 4 [46,48]. The positive best fit value of parameter ξ describe a conversion of dark matter to dark energy although both in flat and non-flat universes, in 1-σ CL, an inverse conversion is possible as well. The interacting DE and DM models have been constrained by observational data by many authors with different parametrization of the interacting parameter Q [77][78][79][80][81][82][83]. He et al. in [82] have parametrized the interaction parameter as in the present paper although they have chosen the prior on parameter ξ as ξ = [0, 0.02]. They obtained the best fit value of parameter ξ as ξ = 0.0006 ± 0.0006.

V. CONCLUSION
In this paper we considered the cosmological constraints on the parameters of the GDE in the framework of BD theory by using a Markov Chain Monte Carlo simulation. We used the SNIa+ CMB+ BAO+X-ray gas mass fraction  data for the model fitting. The best fit values of the cosmological parameters in this model are compatible with the results of the ΛCDM model. In addition, we obtained the best fit values of parameters ǫ = 1/ω and ξ where ω is the BD parameter and ξ is the interacting parameter. The best fit values of these parameters are also compatible with the results of previous constraining works. However as we mentioned in section III, due to large systematic error in ǫ by using the supernovae data only, to constrain this parameter with a higher precision, we should combine the supernovae data with other cosmological data sets as the CMB and BAO data sets. In addition we should assert one more time that the data points parameters of the CMB and BAO data which we have used in this paper are the best fit values for ΛCDM and the error estimates are also based on the ΛCDM model. Therefore they are not completely accurate in this application. The numerical results can be improved in the future works by using more recent data such as nine-yaer WMAP [61] or the Planck [84] projects.  To produce these plots, SNIa+CMB+BAO+X-ray gas mass fraction data together with the BBN constraints have been used.