text:
# Is there Supernova Evidence for Dark Energy Metamorphosis ?

## 1 Introduction

One of the most tantalizing observational discoveries of the past decade has been that the expansion of the universe is speeding up rather than slowing down. An accelerating universe is strongly suggested by observations of type Ia high redshift supernovae provided these behave as standard candles. The case for an accelerating universe is further strengthened by the discovery of Cosmic Microwave Background (CMB) anisotropies on degree scales (which indicate $`\mathrm{\Omega }_{\mathrm{total}}1`$) combined with a low value for the density in clustered matter $`\mathrm{\Omega }_m1/3`$ deduced from galaxy redshift surveys. All three sets of observations strongly suggest that the universe is permeated by a relatively smooth distribution of ‘dark energy’ (DE) which dominates the density of the universe ($`\mathrm{\Omega }_{\mathrm{DE}}2\mathrm{\Omega }_m2/3`$) and whose energy momentum tensor violates the strong energy condition ($`\rho +3p0`$) so that $`w_{\mathrm{DE}}=p/\rho <1/3`$.

Although a cosmological constant ($`w=1`$) provides a plausible answer to the conundrum posed by dark energy, it is well known that the unevolving cosmological constant faces serious ‘fine tuning’ problems since the ratio between $`\rho _\mathrm{\Lambda }=\mathrm{\Lambda }/8\pi G`$ and the radiation density, $`\rho _r`$, is already a miniscule $`\rho _\mathrm{\Lambda }/\rho _r10^{54}`$ at the electroweak scale ($`T100`$ GeV) and even smaller, $`\rho _\mathrm{\Lambda }/\rho _r10^{123}`$, at the Planck scale ($`T10^{19}`$ GeV). This issue is further exacerbated by the ‘cosmological constant problem’ which arises because the $`\mathrm{\Lambda }`$-term generated by quantum effects is enormously large $`\rho _\mathrm{\Lambda }\text{ }\stackrel{>}{}m_{\mathrm{Pl}}^4`$, where $`m_{\mathrm{Pl}}1.2\times 10^{19}`$ GeV is the Planck mass (Zeldovich, 1968; Weinberg, 1989).

Although the cosmological constant problem remains unresolved, the issue of fine tuning which plagues $`\mathrm{\Lambda }`$ has led theorists to explore alternative avenues for DE model building in which either DE or its equation of state are functions of time. (Following Sahni et al.  (2003) we shall refer to the former as Quiessence and to the latter as Kinessence.) Inspired by inflation, the first dark energy models were constructed around a minimally coupled scalar field (quintessence) whose equation of state was a function of time and whose density dropped from a large initial value to the small values which are observed today (Peebles & Ratra, 1988; Wetterich, 1988). (‘Tracker’ quintessence models had the advantage of allowing the current accelerating epoch to be reached from a large family of initial conditions (Caldwell, Dave & Steinhardt, 1998).)

Half a decade after SNe-based observations pointed to the possibility that we may be living in an accelerating universe, the theoretical landscape concerning dark energy has evolved considerably (see the reviews Sahni & Starobinsky, 2000; Carroll, 2001; Peebles & Ratra, 2002; Sahni, 2002; Padmanabhan, 2003). In addition to the cosmological constant and quintessence, the current paradigm for DE includes the following interesting possibilities:

* Dark energy with $`w1`$ (Chiba, Okabe, & Yamaguchi, 2000; Caldwell, 2002; McInnes, 2002; Sahni & Shtanov, 2003; Alam & Sahni, 2002; Caldwell, Kamionkowski & Weinberg, 2003; Carroll, Hoffman & Trodden, 2003; Frampton, 2003; Frampton & Takahashi, 2003; Singh, Sami & Dadhich, 2003; Johri, 2003)
* The Chaplygin Gas whose equation of state drops from $`w=0`$ at high redshifts to $`w1`$ today (Kamenshchik, Moschella & Pasquier, 2001)
* Braneworld models in which the source for cosmic acceleration rests in the gravity sector rather than in the matter sector of the theory (Deffayet, Dvali & Gabadadze, 2002; Sahni & Shtanov, 2003; Maeda, Mizuno & Torii, 2003)
* Dark energy models with negative potentials (Felder et al. , 2002; Kallosh et al. , 2002; Alam, Sahni & Starobinsky, 2003)
* Interacting models of dark matter and dark energy (Amendola, 2000; Chimento et al. , 2003; Hoffman, 2003)
* Modified gravity and scalar-tensor theories (Boisseau et al. , 2000; Bertolami & Martins, 2000; Bartolo & Pietroni, 2000; Damour, Kogan & Papazoglou, 2002)
* Dark energy driven by quantum effects (Sahni & Habib, 1998; Parker & Raval, 1999)
* Dark energy with a late-time transition in the equation of state (Bassett et al. , 2002; Corasaniti et al. , 2003)
* Unified models of dark energy and inflation (Peebles & Vilenkin, 1999; Copeland, Liddle & Lidsey, 2001; Sahni, Sami & Souradeep, 2002) etc.

Faced with the current plethora of dark energy scenarios the concerned cosmologist is faced with two options:

(i) She can test every single model against observations,

(ii) She can take a more flexible approach and determine the properties of dark energy in a model independent manner.

In this paper we proceed along route (ii) and demonstrate that model independent reconstruction brings us face to face with exciting new properties of dark energy.

Applying the techniques developed in Saini et al.  (2000); Sahni et al.  (2003) to a new data set consisting of $`172`$ Supernovae from Tonry et al.  (2003) and an additional 22 Supernovae from Barris et al.  (2003) we show that the DE equation of state which best fits the data evolves from $`w0`$ at $`z1`$ to $`1.2\text{ }\stackrel{<}{}w\text{ }\stackrel{<}{}1`$ today. An evolving equation of state of DE is favoured by the data over a cosmological constant for a large region in parameter space.

## 2 Model independent reconstruction of dark energy

Supernova observations during the previous decade have been pioneered by two teams: The High-z Supernova Search Team (HZT) (Riess et al. , 1998) and the Supernova Cosmology Project (SCP) (Perlmutter et al. , 1999). The enormous efforts made by these two teams have changed the way cosmologists view their universe. A recent analysis (Tonry et al. , 2003) of 172 type Ia supernovae by HZT gives the following bounds on the cosmic equation of state (at $`95\%`$ CL)

$$1.48<w<0.72,$$
(1)
when the 2dFGRS prior $`\mathrm{\Omega }_mh=0.2\pm .03`$ is assumed (Percival et al. , 2001). A similar bound

$$1.61<w<0.78,$$
(2)
is obtained for a new sample of high-z supernovae by SCP (Knop et al. , 2003). 111It is interesting that, when no priors are set on $`\mathrm{\Omega }_m`$, the dark energy equation of state becomes virtually unbounded from below and has a $`99\%`$ confidence limit of being $`<1`$ ! (Knop et al. , 2003)

These results clearly rule out several DE contenders including a tangled network of cosmic strings ($`w1/3`$) and domain walls ($`w2/3`$). However a note of caution must be added before we apply (1) or (2) to the wider class of DE models discussed in the introduction. Impressive as the bounds in (1) & (2) are, they strictly apply only to dark energy having a constant equation of state since this prior was assumed both in the analysis of the supernova data set as well as in the 2dFGRS study (Tonry et al. , 2003; Knop et al. , 2003). Aside from the cosmological constant ($`w=1`$), the topological defect models alluded to earlier and the sine-hyperbolic scalar field potential (Sahni & Starobinsky, 2000; Urena-Lopez & Matos, 2000; Sahni et al. , 2003) no viable DE models exist with the property $`w=\mathrm{constant}`$. Indeed, most models of dark energy (Quintessence, Chaplygin gas, Braneworlds, etc.) can show significant evolution in $`w(z)`$ over sufficiently large look back times.

In this paper we shall reconstruct the properties of dark energy without assuming any priors on the cosmic equation of state. (The dangers of imposing priors on $`w(z)`$ have been highlighted in Maor et al.  (2002) and several of our subsequent results will lend support to the conclusions reached in this paper.)

### 2.1 Cosmological reconstruction of $`w(z)`$

Cosmological reconstruction is based on the observation that, in a spatially flat universe, the luminosity distance and the Hubble parameter are related through the equation (Starobinsky, 1998; Huterer & Turner, 1999; Nakamura & Chiba, 1999):

$$H(z)=\left[\frac{d}{dz}\left(\frac{d_L(z)}{1+z}\right)\right]^1.$$
(3)
Thus knowing $`d_L`$ we can unambiguously determine the Hubble parameter as a function of the cosmological redshift. Next, the Einstein equations

$`H^2`$ $``$ $`{\displaystyle \frac{8\pi G}{3}}[\rho _m+\rho _{\mathrm{DE}}],`$
$`q`$ $`=`$ $`{\displaystyle \frac{\ddot{a}}{aH^2}}={\displaystyle \frac{4\pi G}{3H^2}}{\displaystyle \underset{i}{}}(\rho _i+3p_i),`$ (4)
are used to determine the energy density and pressure of dark energy:

$`\rho _{\mathrm{DE}}`$ $`=`$ $`\rho _{\mathrm{critical}}\rho _\mathrm{m}={\displaystyle \frac{3H^2}{8\pi G}}(1\mathrm{\Omega }_\mathrm{m}(x)),`$
$`p_{\mathrm{DE}}`$ $`=`$ $`{\displaystyle \frac{H^2}{4\pi G}}(q{\displaystyle \frac{1}{2}}),`$ (5)
where $`\rho _{\mathrm{critical}}=3H^2/8\pi G`$ is the critical density of a FRW universe. The equation of state of DE $`w_{\mathrm{eff}}=p_{\mathrm{DE}}/\rho _{\mathrm{DE}}`$ follows immediately (Saini et al. , 2000)

$$w_{\mathrm{eff}}(x)=\frac{2q(x)1}{3\left(1\mathrm{\Omega }_\mathrm{m}(x)\right)}\frac{(2x/3)d(\mathrm{ln}H)/dx1}{1(H_0/H)^2\mathrm{\Omega }_{0\mathrm{m}}x^3},$$
(6)
where $`\mathrm{\Omega }_{0\mathrm{m}}=8\pi G\rho _{0\mathrm{m}}/3H_0^2`$, $`x=1+z`$. In quintessence models and in $`\mathrm{\Lambda }`$CDM, the equation (6) determines the true ‘physical’ equation of state of dark energy. However the subscript ‘eff’ in $`w_{\mathrm{eff}}`$ stresses the fact that this quantity should be interpreted as an ‘effective’ equation of state in DE models in which gravity is non-Einsteinian or in models in which dark energy and dark matter interact. Examples of the former include Braneworld models and scalar-tensor theories. It is well known that in a large class of Braneworld models the Hubble parameter does not adhere to the Einsteinian prescription (4) since it includes explicit interaction terms between dark matter and dark energy (Deffayet, Dvali & Gabadadze, 2002; Sahni & Shtanov, 2003). In this case the equation of state determined using (6) can still be used to characterize DE, but physical interpretations of $`w_{\mathrm{eff}}`$ need to be treated with caution. 222One way around this difficulty is to define observables solely in terms of $`H`$ and its derivatives (called ‘Statefinders’ in Sahni et al. , 2003). A detailed discussion of these issues can be found in Alam et al.  (2003).

One route towards the meaningful reconstruction of $`w(z)`$ lies in inventing a sufficiently versatile fitting function for either $`d_L(z)`$ or $`H(z)`$. The parameters of this fitting function are determined by matching to Supernova observations and $`w(z)`$ is determined from (3) and (6).333Alternatively one could apply an ansatz to $`w(z)`$ itself (Chiba & Nakamura, 2000; Weller & Albrecht, 2002; Corasaniti & Copeland, 2003; Gerke & Efstathiou, 2002; Maor et al. , 2002; Linder, 2003). See Alam et al.  (2003) for a summary of different approaches to cosmological reconstruction. Non-parametric approaches are discussed in Wang & Lovelace (2001); Huterer & Starkman (2002); Saini (2003); see also Daly & Djorgovsky (2003); Nunes & Lidsey (2003). Our reconstruction exercise will be based upon the following flexible and model independent ansatz for the Hubble parameter (Sahni et al. , 2003)

$$H(x)=H_0\left[\mathrm{\Omega }_mx^3+A_0+A_1x+A_2x^2\right]^{\frac{1}{2}},$$
(7)
where $`x=1+z`$. This ansatz for $`H(z)`$ is exact for the cosmological constant $`w=1`$ ($`A_1=A_2=0`$) and for DE models with $`w=2/3`$ ($`A_0=A_2=0`$) and $`w=1/3`$ ($`A_0=A_1=0`$). It has also been found to give excellent results for DE models in which the equation of state varies with time including quintessence, Chaplygin gas, etc. (Sahni et al. , 2003; Alam et al. , 2003). The ansatz (7) is equivalent to the following expansion for DE

$$\rho _{\mathrm{DE}}=\rho _{0\mathrm{c}}(A_0+A_1x+A_2x^2+A_3x^3),$$
(8)
where $`\rho _{0\mathrm{c}}=3H_0^2/(8\pi G)`$ is the present day critical density. The condition $`A_30`$ allows $`\rho _{\mathrm{DE}}`$ to mimic the properties of dark matter at large redshifts ($`A_31`$ follows from large scale structure constraints). From (7) and (8) we find $`\mathrm{\Omega }_m=\mathrm{\Omega }_{0\mathrm{m}}+A_3`$, i.e. the value of $`\mathrm{\Omega }_m`$ in (7) can be slightly larger than $`\mathrm{\Omega }_{0\mathrm{m}}`$ in this case.

Substituting (7) into the expression for the luminosity distance we get

$$\frac{d_L(z)}{1+z}=\frac{c}{H_0}_1^{1+z}\frac{dx}{\sqrt{\mathrm{\Omega }_mx^3+A_0+A_1x+A_2x^2}}.$$
(9)
The parameters $`A_0,A_1,A_2`$ are determined by fitting (9) to supernova observations using a maximum likelihood technique. This ansatz has only three free parameters $`(\mathrm{\Omega }_m,A_1,A_2)`$ since $`A_0+A_1+A_2=1\mathrm{\Omega }_m`$ for a flat universe. A note of caution: since the ansatz (8) is a truncated Taylor expansion in $`x=1+z`$ its range of validity is $`z\text{ }\stackrel{<}{}\mathrm{few}`$, consequently the ansatz-derived $`H(z)`$ and $`d_L(z)`$ should not be used at higher redshifts.

Note that the weak energy condition for dark energy $`\rho _{\mathrm{DE}}0,\rho _{\mathrm{DE}}+p_{DE}0`$ has the following form for the ansatz (7) :

$$A_0+A_1x+A_2x^20,A_1+2A_2x0,$$
(10)
provided we assume that the $`\mathrm{\Omega }_mx^3`$ term in (7) is totally due to non-relativistic dark matter and does not include any contribution from dark energy. The demand that the WEC (10) be satisfied for all $`x0`$ (i.e. in the past as well as in the future) requires $`A_0,A_1,A_2`$ to be non-negative. However, the demand that the WEC (10) be satisfied in the past ($`x1`$) but not necessarily in the future, leads to the somewhat weaker constraint

$$A_1+2A_20,A_20.$$
(11)
(Models in which $`\rho _{\mathrm{DE}}(z)<0`$ for $`z<0`$ and which violate the WEC in the future, have been discussed in Felder et al.  (2002); Kallosh et al.  (2002); Alam, Sahni & Starobinsky (2003).)

The presence of the term $`\mathrm{\Omega }_mx^3`$ in (7) has two important consequences: (i) It ensures that the the universe transits to a matter dominated regime at early times ($`z1`$), (ii) It allows us to incorporate information (available from other data sets) regarding the current value of the matter density in the universe. This information can be used to perform a maximum likelihood analysis with the introduction of suitable priors on $`\mathrm{\Omega }_{0\mathrm{m}}`$. In further analysis we will assume that the $`\mathrm{\Omega }_mx^3`$ term in (7) does not include any contribution from dark energy.

We have also studied simple extensions of the ansatz (7) by adding new terms $`A_1x^1`$ and $`A_4x^4`$. The $`A_1x^1`$ term allows $`w(z)`$ to become substantially less than $`1`$, thereby providing greater leeway to phantom models. The $`A_4x^4`$ term allows DE to evolve towards equations of state which are more stiff than dust ($`w=0`$); its role is therefore complementary to that of $`A_1x^1`$. Despite the inclusion of these new terms, our best fit to the supernova data presented below does not change significantly(choosing $`A_1=0.0003`$ and $`A_4=0.008`$), which points to the robustness of the ansatz (7) for the given data set.

We should add that our reason for choosing an ansatz to fit $`H(z)`$ rather than some other cosmological quantity was motivated by the fact that the Hubble parameter is directly related to a fundamental physical quantity – the Ricci tensor, and is therefore likely to remain meaningful even when other quantities (such as the equation of state) become ‘effective’. (This happens for instance, in the case of the braneworld models of dark energy discussed in Deffayet, Dvali & Gabadadze (2002); Sahni & Shtanov (2003).)

The rationale for choosing a three parameter ansatz for $`H(z)`$ is the following. The observed luminosity distance determined using type Ia supernovae is rather noisy, therefore in order to determine the Hubble parameter from $`D_L(z)`$ and following that the equation of state, one must take two derivatives of a noisy quantity. This difficulty can be tackled in two possible ways: (i) either one smoothes the data over some interval $`\mathrm{\Delta }t`$ (binning is one possibility), or (ii) we may choose to smooth ‘implicitly’ by parameterizing $`H(z)`$ through an appropriate fitting function. The number of free parameters $`𝒩`$ in the fit to $`H(z)`$ will be related to the smoothing interval $`\mathrm{\Delta }z`$ through $`\mathrm{\Delta }z=z_{\mathrm{max}}/𝒩`$. Increasing $`𝒩`$ implies decreasing $`\mathrm{\Delta }z`$ which results in a rapid growth of errors through $`\mathrm{\Delta }H(z)(\mathrm{\Delta }z)^{3/2}`$, and $`\mathrm{\Delta }w(z)(\mathrm{\Delta }z)^{5/2}`$ (Tegmark, 2002), therefore in order not to loose too much accuracy in our reconstruction we considered 3 parameter fits for $`H(z)`$ in our paper (these correspond to 2 parameter fits for $`w(z)`$).

We now test the usefulness of the ansatz (7) in reconstructing different dark energy models. The ansatz returns exact values for $`\mathrm{\Lambda }`$CDM, and $`w=1/3`$, $`w=2/3`$ quiessence models. In figure 1 we show the accuracy of the ansatz (7) when applied to several other dark energy models such as tracker quintessence, the Chaplygin gas and super-gravity (SUGRA) models. We plot the deviation of $`\mathrm{log}(d_LH_0)`$ (which is the measured quantity for SNe) obtained with the ansatz (7) from the actual model values. Clearly the ansatz performs very well over a significant redshift range for $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$ (Also see appendix B). In fact, in the redshift range where SNe data is available, the ansatz recovers these models of dark energy with less than $`0.5\%`$ errors. However it would be appropriate to add a note of caution at this point. Although figure 1 clearly demonstrates the usefulness of the ansatz for some DE models, its performance vis-a-vis other models of DE is by no means guaranteed. By its very construction the ansatz (7) is expected to have limitations when describing models with a fast phase transition (Bassett et al. , 2002) as well as rapidly oscillating quintessence models (Sahni & Wang, 2000). (The ansatz (7) can give reasonable results even for these models provided the resulting DE behaviour is suitably smoothed.) For this reason, although the bulk of our analysis will be carried out using (7), we shall supplement it when necessary with other fitting functions, which will provide us with an independent means with which to test the robustness of our reconstruction exercise.

Methodology :

For our primary reconstruction, we use a subset of 172 type Ia Supernovae, obtained by imposing constraints $`A_V<0.5`$ and $`z>0.01`$ on the 230 SNe sample, as in the primary fit of Tonry et al.  (2003). For the ansatz (9), we require to fit four parameters: ($`H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2`$). We may use prior information on $`H_0`$ ($`H_0=72\pm 8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, Freedman et al.  (2001)) and $`\mathrm{\Omega }_{0\mathrm{m}}`$ ($`\mathrm{\Omega }_{0\mathrm{m}}h=0.2\pm 0.03`$, Percival et al.  (2001)). 444One should note however that the prior on $`\mathrm{\Omega }_{0\mathrm{m}}`$ is not model independent since it relies on the $`\mathrm{\Lambda }`$CDM model to project from redshift space to real space. Results coming from the use of this prior should therefore not be taken too literally in the present context. See Kunz et al.  (2003) for an interesting discussion of related issues.

The measured quantity for this data is $`y=<\mathrm{log}(d_LH_0)>`$, therefore the likelihood function is given by

$``$ $`=`$ $`𝒩\mathrm{exp}\left({\displaystyle \frac{\chi ^2}{2}}\right),`$ (12)
$`\chi ^2`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{172}{}}}\left({\displaystyle \frac{y_iy_{\mathrm{fit}}(H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)}{\sigma _i}}\right)^2,`$ (13)
where $`𝒩`$ is a normalisation constant. Therefore, the probability distribution function in the four-space $`(H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ is

$$\mathrm{P}(H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)\mathrm{exp}\left(\frac{\chi ^2}{2}\right)\mathrm{Pr}(\mathrm{\Omega }_{0\mathrm{m}}h)\mathrm{Pr}(H_0).$$
(14)
where $`\mathrm{Pr}`$ refers to the priors applied on the parameters of the system.

Our goal is to reconstruct cosmological parameters such as the equation of state $`w(z)=w(z;\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$, therefore we marginalise over $`H_0`$ and obtain the probability distribution function in the $`(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ space:

$$\stackrel{~}{P}(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)=\mathrm{P}(H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)𝑑H_0.$$
(15)
In order to do this, we have to define the bounds of a four-dimensional volume in $`(H_0,\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$. The bounds of $`H_0`$ are taken at $`5\sigma `$ of the HST prior. For $`\mathrm{\Omega }_{0\mathrm{m}}`$, the natural choice is $`0\mathrm{\Omega }_{0\mathrm{m}}1`$. It is not immediately obvious what the bounds should be for $`A_1,A_2`$. We choose a sufficiently large rectangular grid for $`A_1,A_2`$ (roughly corresponding to $`6\text{ }\stackrel{<}{}w_0\text{ }\stackrel{<}{}5`$) which includes most known models of dark energy. This bound is merely a matter of convenience and does not affect our results in any way. After marginalisation, we have a grid in $`(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ space on which $`\stackrel{~}{P}(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ is specified at each point. We may now proceed in two ways. Firstly, we may choose to fix $`\mathrm{\Omega }_{0\mathrm{m}}`$ at a suitable constant value (e.g., $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$) thereby obtaining a grid in the $`(A_1,A_2)`$ plane with $`P`$ (the probability if $`\mathrm{\Omega }_{0\mathrm{m}}`$ is known to be an exact value) defined at each point. For a particular redshift, we may then calculate $`w(z;\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ at each point of the grid. This would yield results that would hold true if $`\mathrm{\Omega }_{0\mathrm{m}}`$ were known exactly. Instead of using the exact value of $`\mathrm{\Omega }_{0\mathrm{m}}`$, we may use the prior information about it available to us ($`\mathrm{\Omega }_{0\mathrm{m}}h=0.2\pm 0.03`$), and calculate $`w(z;\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ at each point of a three-dimensional grid, the probability $`\stackrel{~}{P}`$ at each point being known. Therefore, at any given redshift $`z`$, $`w(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$ can be tagged with a numerical value $`\stackrel{~}{P}(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$. Starting from the best-fit $`w(z)`$ (the value at the peak of the probability distribution), we may move down on either side till $`34\%`$ of the total area is enclosed under the curve, thus obtaining asymmetric $`1\sigma `$ bounds on $`w(z)`$. The $`2\sigma ,3\sigma `$ bounds can be similarly obtained.

Results :

We first show preliminary results for which the matter density is fixed at a constant value of $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$. A detailed look at the $`\chi ^2`$ surface in the ($`A_1,A_2`$) plane (figure 2) reveals the existence of two minima in $`\chi ^2`$, a shallower one close to $`\mathrm{\Lambda }`$CDM ($`A_1=0.177,A_2=0.119,w_0=1.03,\chi ^2=1.0402`$), and a deeper minimum at $`A_1=4.360,A_2=1.829,w_0=1.33,\chi ^2=1.0056`$. We would like to draw the readers attention to the fact that imposing the prior $`w(z)1(z0)`$ amounts to disallowing a significant region of parameter space (the unshaded region in figure 2). Consequently an analysis which assumes $`w(z)1`$ loses all information about the region $`2\sigma `$ around the deeper minimum ! Since we have no reason (observational or theoretical) to favour either minimum over the other, we shall always choose the deeper minimum as our best-fit in all the subsequent calculations.

In the figure 3, we plot the deviation of the squared Hubble parameter $`H^2/H_0^2`$ from $`\mathrm{\Lambda }`$CDM over redshift for the best-fit. We note that the quantity $`H^2/H_0^2`$ has a simple linear relationship with the parameters of the fit (Eq 7), therefore the errors in this quantity increase with redshift. Another quantity of interest is the energy density of dark energy. For this ansatz, $`\stackrel{~}{\rho }_{\mathrm{DE}}=\rho _{\mathrm{DE}}/\rho _{0\mathrm{c}}=A_0+A_1x+A_2x^2`$ (where $`\rho _{0\mathrm{c}}=3H_0^2/8\pi G`$ is the present day critical density). The figure 4 shows the logarithmic variation of $`\stackrel{~}{\rho }_{\mathrm{DE}}`$ with redshift. In this figure too the errors increase with redshift. An interesting point to note is that initially, dark energy density decreases with redshift, showing the phantom-like nature ($`w<1`$) of dark energy at lower redshifts of $`z\text{ }\stackrel{<}{}0.25`$, while at higher redshifts, the dark energy density begins to track the matter density. Before moving on to the second derivative of the luminosity distance (e.g., the equation of state) we may obtain more information from the dark energy density by considering a weighted average of the equation of state :

$$1+\overline{w}=\frac{1}{\mathrm{\Delta }\mathrm{ln}(1+z)}(1+w(z))\frac{dz}{1+z},$$
(16)
where $`\mathrm{\Delta }`$ denotes the total change of the variable between integration limits. This quantity can be elegantly expressed in terms of the difference in energy densities over a range of redshift as

$$1+\overline{w}=\frac{1}{3}\frac{\mathrm{\Delta }\mathrm{ln}\stackrel{~}{\rho }_{\mathrm{DE}}}{\mathrm{\Delta }\mathrm{ln}(1+z)}.$$
(17)
Thus the variation in the dark energy density depicted in figure 4 is very simply related to the weighted average equation of state !

In table 1 we show the values of $`\overline{w}`$ obtained using different ranges in redshift for our best-fit with corresponding $`1\sigma `$ and $`2\sigma `$ errors. We have taken the ranges of integration to be approximately equally spaced in ln$`(1+z)`$, with the upper limit set by the furthest supernova known at present. The values of $`\overline{w}`$ may be calculated using the equation (16) (which uses the second derivative of the luminosity distance), or they can simply be read off from figure 4 using equation (17). From this table, a “metamorphosis” in the properties of dark energy occurring somewhere between $`z0`$ and $`z1`$ can be clearly seen (note that, effectively, one needs to differentiate $`d_L(z)`$ only once to come to this conclusion).

We now reconstruct the equation of state of dark energy which, for the ansatz (7), has the form

$$w(z)=1+\frac{x}{3}\frac{A_1+2A_2x}{A_0+A_1x+A_2x^2},x=1+z.$$
(18)
(Note that, since $`w(z)`$ was derived from the ansatz (7), its domain of validity is $`z\text{ }\stackrel{<}{}\mathrm{few}`$.) In figures 5 (a), (b), (c), we show the evolution of the equation of state $`w(z)`$ with redshift for different values of $`\mathrm{\Omega }_{0\mathrm{m}}`$. The $`1\sigma `$ and $`2\sigma `$ limits are shown in each case. The $`\chi ^2`$ per degree of freedom for the best-fit for the different cases is given in Table 2. We find that for $`0.2<\mathrm{\Omega }_{0\mathrm{m}}<0.4`$, the behaviour of the equation of state does not change significantly with change in the matter density. However, for larger values of $`\mathrm{\Omega }_{0\mathrm{m}}`$, a *smaller* current value of $`w_0=w(z=0)`$ is preferred. In all three cases considered, the present value of the equation of state is $`w_0\text{ }\stackrel{<}{}1.2`$ for the best-fit, and the equation of state rises steeply from $`w\text{ }\stackrel{<}{}1.2`$ to $`w0`$ with redshift. In fact, the behaviour of $`w`$ appears to be extremely different from that in $`\mathrm{\Lambda }`$CDM ($`w=1`$). We note here that, for this analysis, the errors on $`w`$ appear to decrease with redshift. This may appear counter-intuitive, since there are fewer SNe at higher redshifts, but this is merely a construct of the fact that $`w`$ depends non-linearly on the parameters of the ansatz (see appendix A).

Quintessence models satisfy the weak energy condition (WEC) $`\rho +p=\dot{\varphi }^20`$ and it would be interesting to see how the imposition of the WEC as a prior on the equation of state will affect the results of our analysis. We therefore perform the same analysis as above with the added constraint $`w_01`$ (note that this implies $`w(z)1`$ for all $`z0`$ for our fitting function of $`H(z)`$ provided $`A_20`$). The results are shown in the figures 6 (a), (b), (c). We see that in this case the errors are larger and the evolution of the equation of state with redshift follows a much gentler slope. Such an equation of state would be largely consistent with the cosmological constant model. (These results are in broad agreement with an earlier analysis of Saini et al.  (2000) in which a smaller SNe data set was used and a different ansatz for the luminosity distance was applied.)

In Table 2, we show how the $`\chi ^2`$ for the best-fit evolving dark energy models compare with that for $`\mathrm{\Lambda }`$CDM. We find that $`\chi _{\mathrm{\Lambda }\mathrm{CDM}}^2>\chi _{\mathrm{best}\mathrm{fit}}^2`$ always. For $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$, the value of $`\chi _{\mathrm{\Lambda }\mathrm{CDM}}^2`$ is just within $`2\sigma `$ of the best fit $`\chi ^2`$, but for $`\mathrm{\Omega }_{0\mathrm{m}}=0.2`$, or $`\mathrm{\Omega }_{0\mathrm{m}}=0.4`$, $`\chi _{\mathrm{\Lambda }\mathrm{CDM}}^2`$ is outside the $`2\sigma `$ limits of the best-fit. It is also noteworthy that when the prior $`w(z)1(z0)`$ is used, the best-fit model has a slowly evolving equation of state with $`w_0=1`$ and the $`\chi ^2`$ for the best-fit becomes smaller for a smaller value of the matter density. When no priors are assumed on $`w`$, the trend reverses, and better fits are obtained for larger values of $`\mathrm{\Omega }_{0\mathrm{m}}`$. From this it appears that at least at $`1\sigma `$ the evolving dark energy model is favoured over $`\mathrm{\Lambda }`$CDM, and it does as well, if not better at the $`2\sigma `$ level, depending upon the value of the present-day matter density.

Using Priors on $`\mathrm{\Omega }_{0\mathrm{m}}`$ :

Instead of assuming an exact value for $`\mathrm{\Omega }_{0\mathrm{m}}`$, which is somewhat optimistic given the present observational scenario, we may use the 2dF prior on $`\mathrm{\Omega }_{0\mathrm{m}}`$ and calculate $`w(z)`$ as a function of $`(\mathrm{\Omega }_{0\mathrm{m}},A_1,A_2)`$. It should be noted here that the 2dF error bars on $`\mathrm{\Omega }_{0\mathrm{m}}h`$ have been calculated using 2dF data in conjunction with CMB, and this calculation assumes a $`\mathrm{\Lambda }`$CDM model, therefore this prior should be used more as a benchmark for the value of $`\mathrm{\Omega }_{0\mathrm{m}}`$ rather than as an absolute when considering evolving dark energy models. The resultant “marginalised” $`w`$ is shown as a function of the redshift in figure 7 (a). The nature of the equation of state for the analysis with the added prior $`w(z)1(z0)`$ is shown in figure 7 (b). We find that the general nature of evolution of the equation of state is not changed by adding this extra information on the matter density. If no priors are assumed on the equation of state to begin with, $`w(z)`$ still rises sharply from $`w_0\text{ }\stackrel{<}{}1`$ up to $`w0`$ at maximum redshift and the analysis appears to favour a fast-evolving equation of state of dark energy over the standard $`\mathrm{\Lambda }`$CDM model. If a prior $`w1`$ is assumed, then the marginalised equation of state is more consistent with the cosmological constant. From this we see that marginalisation over $`\mathrm{\Omega }_{0\mathrm{m}}`$ does not lead to any significant change in our results. In the subsequent sections, we will show our results for $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$.

From the above analysis, we find, therefore, that our results change significantly depending upon whether or not the prior $`w1`$ is imposed. We saw earlier that in the absence of any prior on $`w(z)`$, the best-fit equation of state rose from $`w\text{ }\stackrel{<}{}1`$ at $`z=0`$ to $`w0`$ at $`z1`$. By imposing a prior on the equation of state, we effectively screen off a sizeable part of the parameter space (see figure 2), and therefore the reconstruction is forced to choose its best-fit away from the true minima of the $`\chi ^2`$ surface. The effect of imposing a prior on $`w(z)`$ is therefore to make the best-fit $`w(z)`$ grow much more slowly with $`w=1`$ being preferred at $`z=0`$. Our results show that the reconstructed equation of state with the prior $`w1`$ is in good agreement with a cosmological constant at the $`68\%`$ CL. However, if no prior is imposed, then the steeply evolving dark energy models are favoured over the cosmological constant at $`1\sigma `$, and are at least as likely as the cosmological constant at the $`2\sigma `$ level.

Age and Deceleration Parameter of the Universe:

We may also use this ansatz to calculate other quantities of interest, such as the age of the universe, $`t(z)`$, and the deceleration parameter, $`q(z)`$:

$`t(z)`$ $`=`$ $`{\displaystyle _{1+z}^{\mathrm{}}}{\displaystyle \frac{dx}{xH(x)}}`$ (19)
$`q(z)`$ $`=`$ $`{\displaystyle \frac{\ddot{a}}{aH^2}}{\displaystyle \frac{H^{}}{H}}x1.`$ (20)
where $`x=1+z`$.

In figure 8 we plot the evolution of the age of the universe with redshift. We find that the best-fit age of the universe today is $`t_0=12.8`$ Gyrs if the Hubble parameter is taken to be $`H_0=72\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, which is slightly lower than the age of a $`\mathrm{\Lambda }`$CDM universe, $`t_0=13.4`$ Gyrs (both values are for $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$). At the $`2\sigma `$ level, the age of the universe today would vary between $`11.2t_013.6`$ Gyrs.

Figure 10 shows the evolution of the deceleration parameter with redshift. We find that the behaviour of the deceleration parameter for the best-fit universe is quite different from that in $`\mathrm{\Lambda }`$CDM cosmology. Thus, the current value of $`q_00.9`$ is significantly lower than $`q_00.55`$ for $`\mathrm{\Lambda }`$CDM (assuming $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$). Furthermore the rise of $`q(z)`$ with redshift is much steeper in the case of the best-fit model, with the result that the universe begins to accelerate at a comparatively lower redshift $`z0.3`$ (compared with $`z0.7`$ for $`\mathrm{\Lambda }`$CDM) and the matter dominated regime ($`q1/2`$) is reached by $`z1`$.

Using Priors on Age of the Universe :

Important consistency checks on our best-fit Universe may be provided by observations of the age of the Universe. Unfortunately, estimates of the age of the universe from different methods can produce widely varying results one reason for which is that estimates of the Hubble parameter itself can vary significantly. For instance, the HST key project yields $`H_0=72\pm 8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, while studies of the Sunyaev-Zeldovich effect in galaxy clusters give a significantly lower value $`H_0=60\pm 10\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (Krauss, 2001). Estimates of the ages of the oldest globular clusters suggest $`t_0=12.5\pm 2.5`$ Gyrs, at the $`95\%`$ confidence level (Krauss, 2001; Krauss & Chaboyer, 2001; Gnedin, Lahav & Rees, 2001; Hansen et al. , 2002; Gratton et al. , 2003; Marchi et al. , 2003) and this age estimate is consistent with several other measurements including observations of eclipsing spectroscopic binaries (Thompson et al. , 2001; Chaboyer & Krauss, 2002), results from radioactive dating of a metal-poor star (Cayrel et al. , 2001) and WMAP data (Spergel et al. , 2003) (see also Alcaniz, Jain and Dev (2002)). The results from the WMAP experiment suggest $`t_0=13.4\pm 0.3`$ Gyrs with a Hubble parameter $`H_0=72\pm 5\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, for $`\mathrm{\Lambda }`$CDM cosmology (which satisfies the WEC). Adding SDSS and SNe Ia data to WMAP, Tegmark et al.  (2003) find an age of $`t_0=14.1_{0.9}^{1.0}`$ Gyrs for a slightly closed $`\mathrm{\Lambda }`$CDM universe with $`H_0=66_{6.4}^{6.7}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Although these results cannot be carried over to evolving dark energy models including those implied by our best-fit reconstruction (which violate the WEC) they provide an indication of the range within which the age of the universe might vary. Keeping in mind these various results, we use two different priors on the Hubble parameter: $`H_0=72\pm 8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ ($`1\sigma `$ bound from HST; Freedman et al.  (2001)), and $`H_0=66\pm 6\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (approximate bound from WMAP, SDSS, SNe Ia; Tegmark et al.  (2003)). For each case, we choose three different Gaussian priors on the present age of the universe: $`t_0=13\pm 1\mathrm{Gyrs},14\pm 1\mathrm{Gyrs},\mathrm{and}15\pm 1\mathrm{Gyrs}`$ respectively, and perform the reconstruction for a $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$ universe. The results are shown in the figure 9. We find that, for a Hubble parameter of $`H_0=72\pm 8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and with an additional prior on the age of the universe $`t_0=13\pm 1`$ Gyrs, the best-fit remains nearly the same, showing a rapid evolution of the equation of state from $`w0`$ at $`z1`$ to $`w1.2`$ at $`z=0`$, and the errors become narrower. As the age is increased, the best-fit equation of state evolves more slowly, and the $`\chi _{\mathrm{dof}}^2`$ also increases (see Table 3). For the prior $`H_0=66\pm 6\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, we find that the lowest $`\chi _{\mathrm{dof}}^2`$ is obtained for the age prior of $`t_0=14\pm 1`$ Gyrs, which once again matches our best-fit. It should be noted that the errors are smaller in all cases, even though the $`\chi ^2`$ may be larger. We must remember that the addition of a new prior which is consistent with the underlying dataset would lead to a natural reduction in errors. However, the addition of a prior inconsistent with the dataset would lead to a shift of the likelihood maximum as well as a reduction in errors, and the results would then fail to reflect the actual information present in the dataset. That this is happening here for the higher values of age can be seen from the fact that although the errors are reduced, the $`\chi _{\mathrm{dof}}^2`$ is actually larger. Therefore priors from other observations should be added prudently to ensure that they do not lead to incorrect representation of the data. Since there is as yet no clear model independent consensus on the age of the universe, the results we obtain in this section should be interpreted with a degree of caution.

Figure 10 shows the evolution of the deceleration parameter with redshift. We find that the behaviour of the deceleration parameter for the best-fit universe is quite different from that in $`\mathrm{\Lambda }`$CDM cosmology. Thus, the current value of $`q_00.9`$ is significantly lower than $`q_00.55`$ for $`\mathrm{\Lambda }`$CDM (assuming $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$). Furthermore the rise of $`q(z)`$ with redshift is much steeper in the case of the best-fit model, with the result that the universe begins to accelerate at a comparatively lower redshift $`z0.3`$ (compared with $`z0.7`$ for $`\mathrm{\Lambda }`$CDM) and the matter dominated regime ($`q1/2`$) is reached by $`z1`$.

### 2.2 Robustness of Results :

Based on the above analysis, it is tempting to conclude that the dominant component of the universe today is dark energy with a steeply evolving equation of state which marginally violates the weak energy condition. (Of course, the less radical possibility of weakly time dependent dark energy satisfying the weak energy condition remains an alternative, too.) However, before any such dramatic claims are made, we need to check if our results are in any fashion a consequence of inherent bias in the statistical analysis itself, or in the sampling of the data. We therefore perform the following simple exercises to satisfy ourselves of the robustness of our results.

Using Different Subsets of Supernova Data :

In an attempt to understand how the nature of the reconstructed equation of state is dependent on the distribution of data, we perform the reconstruction exercise on different samples of data. We have confined ourselves to the case where $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$ for these exercises. Firstly, we may exclude the SCP data points from the 172 SNe primary fit, leading to a subsample of 130 SNe. We call this the HZT sample. Figures 11(a) and 12(a) show the result of performing the analysis on this subsample without any constraints. The $`\chi ^2`$ per degree of freedom for the best-fit is $`\chi _{\mathrm{HZT}}^2=0.9707`$, which is lower than $`\chi _{\mathrm{\Lambda }\mathrm{CDM}}^2=0.9939`$ for this sample. In this case we find that, though the error bars are slightly larger, overall the dark energy density behaves in the same way as before (compare figure 11(a) with figure 4), showing phantom like ($`w<1`$) behaviour at lower redshifts and tracking matter at higher redshifts. The equation of state of dark energy also evolves much in the same way as when the entire sample is used (compare figure 12(a) with figure 5(b)), starting at $`w_0\text{ }\stackrel{<}{}1`$ and evolving rapidly to $`w0`$. We may also use a sample complementary to this sample, where all the SCP data points published till date are considered, along with the low redshift Calan-Tololo sample. This leads to a sample of 58 SNe (Perlmutter et al. , 1999; Knop et al. , 2003), which we call the SCP sample. Using this sample, we obtain the figures 11(b) and 12(b). The best-fit has a chi-squared per degree of freedom: $`\chi _{\mathrm{SCP}}^2=1.0147`$, lower than $`\chi _{\mathrm{\Lambda }\mathrm{CDM}}^2=1.0369`$ for this sample. We find that here too, the dark energy density initially decreases and then starts tracking matter. The equation of state shows signs of rising steeply at low redshifts, but since the highest redshift in this sample is $`z=0.86`$, the behaviour of $`w`$ beyond this redshift cannot be predicted, therefore the apparent flattening out of the curve beyond a redshift of one cannot be seen in this case. For both these subsets of data, we may repeat the exercise using the prior $`w(z)1(z0)`$. The results obtained for the equation of state, as seen in figures 13(a), (b), are once again commensurate with the results obtained earlier for the full sample (figure 6(b)). We may therefore conclude from this exercise that subsampling the data does not significantly affect our results, and the steep evolution of the equation of state of dark energy is not a construct of the uneven sampling of the supernovae, but rather, reflects the actual nature of dark energy.

Testing our Ansatz against fiducial dark energy models :

The crucial question of course is whether the reconstructed equation of state of dark energy depends upon the ansatz which is used in the exercise, i.e. , whether the behaviour of the equation of state merely reflects a bias in the ansatz itself. In this section we show how the ansatz performs in recovering dark energy models whose equation of state is known, from simulated data. This ansatz was demonstrated to work extremely well when simulations of SNAP data were used (Alam et al. , 2003). However, simulation of SNAP-like data is an optimistic exercise, since data of this quality is unlikely to be available in the near future. We now demonstrate the accuracy with which the ansatz can recover the fiducial background cosmological model if data is simulated using present-day observational standards. In figures 14 (a), (b), (c), we show how well the ansatz recovers the equation of state for three fiducial models (assuming $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$):

(a) a quiessence dark energy model with a constant equation of state: $`w=0.5`$,

(b) a generalised Chaplygin gas model with $`p=A/\rho ^\alpha `$: with $`\alpha =0.5`$ and the present-day equation of state $`w_0=0.9`$, which would give rise to an effective equation of state

$$w(z)=\frac{|w_0|}{|w_0|+(1|w_0|)(1+z)^{3(1+\alpha )}},$$
(21)
and

(c) a model with a linearly evolving equation of state: $`w(z)=w_0+w_1z`$, with $`w_0=1.2,w_1=1`$.

(For DE models with $`w=1,2/3,1/3`$ the ansatz is exact therefore we don’t show the results for these cases.)

We find that in all three cases, the fiducial model lies within the $`68\%`$ confidence limits around the best-fit $`w(z)`$. Based on this result, we claim that within the $`1\sigma `$ error bars, the reconstructed equation of state represents the true properties of dark energy when we use real data.

Using other Ansatz :

It is also important to check whether the results of our reconstruction can be replicated using other ansatz such as fits to the luminosity distance or the equation of state. Many different fits have been suggested in the literature (see for example Huterer & Turner (1999), Saini et al.  (2000), Weller & Albrecht (2002), Gerke & Efstathiou (2002)). Here we choose the fit suggested in Linder (2003) in which the equation of state of dark energy is expanded as

$$w(z)=w_0+\frac{w_1z}{1+z}.$$
(22)
The luminosity distance can therefore be expressed as

$$\frac{d_L(z)}{1+z}=\frac{c}{H_0}_1^{1+z}\frac{dx}{\sqrt{\mathrm{\Omega }_{0\mathrm{m}}x^3+\mathrm{\Omega }_X}},$$
(23)
where $`\mathrm{\Omega }_X=(1\mathrm{\Omega }_{0\mathrm{m}})x^{3(1+w_0+w_1)}\mathrm{exp}[3w_1(\frac{1}{x}1)]`$.

We find that for this fit, the errors in the equation of state get larger with redshift, however this fit too demonstrates that the equation of state of dark energy increases rapidly with redshift (figure 15(a)) when no priors are assumed on the equation of state (EOS). The $`\chi ^2`$ per degree of freedom at the best-fit is $`\chi _{\mathrm{dof}}^2=1.0298`$. When the prior $`w(z)1(z0)`$ is invoked, the best-fit EOS remains very close to the $`\mathrm{\Lambda }`$CDM model (figure 15(b)). Therefore, from this ansatz, we may make the statement that at low redshifts, the equation of state of dark energy shows the same signs of rising steeply with redshift if no priors are assumed on the equation of state, thus supporting our earlier results. The large errors in the equation of state at redshifts of $`z\text{ }\stackrel{>}{}0.5`$ however make it difficult to make any definitive statements about the behaviour of dark energy at high redshifts.

A limitation of the fit (22) is that it is unable to describe very rapid variations in the equation of state. An ansatz which accommodates this possibility has been suggested in Bassett et al.  (2002)

$$w(z)=w_i+\frac{w_fw_i}{1+\mathrm{exp}(\frac{zz_t}{\mathrm{\Delta }})},$$
(24)
where $`w_i`$ is the initial equation of state at high redshifts, $`z_t`$ is a transition redshift at which the equation of state falls to $`w(z_t)=(w_i+w_f)/2`$ and $`\mathrm{\Delta }`$ describes the rate of change of $`w(z)`$.

The resulting luminosity distance has the form:

$$\frac{d_L(z)}{1+z}=\frac{c}{H_0}_1^{1+z}\frac{dx}{\sqrt{\mathrm{\Omega }_{0\mathrm{m}}x^3+\mathrm{\Omega }_X}},$$
(25)
where $`\mathrm{\Omega }_X=(1\mathrm{\Omega }_{0\mathrm{m}})\mathrm{exp}[3_0^{x1}(1+w(z))𝑑z/(1+z)]`$.

The results for the analysis using this fit to the equation of state are shown in fig. 16. We find that when the reconstruction is done without any priors on the equation of state (figure 16(a)), the best fit is remarkably close to the result for ansatz (7) (figure 5(b)). The $`\chi ^2`$ per degree of freedom at the minimum is $`\chi _{\mathrm{dof}}^2=1.0175`$ for this fit. The errors in this case are somewhat larger, especially at high redshift. If we constrain $`w(z)1`$, then as before, the evolution of the equation of state is much slower (figure 16(b)). So the reconstruction using this ansatz appears to confirm our earlier results.

The above exercises lead us to conclude that our results are neither dependent on the nature of the statistical analysis nor on the manner in which the SNe data is sampled. It therefore appears that dark energy with a steeply evolving equation of state provides a compelling alternative to a cosmological constant if data are analysed in a prior-free manner and the weak energy condition is not imposed by hand.

### 2.3 Reconstructing dark energy using a new Supernova sample

As this paper was nearing completion, a new dataset consisting of 23 type Ia SNe was released by the HZT team (Barris et al. , 2003). It is clearly important to check whether or not these new data points corroborate the findings reported in the previous sections. Accordingly, we use a subset of 200 type Ia SNe with $`A_V0.5`$ from the 230 SNe sample of Tonry et al.  (2003), and 22 SNe with $`A_V0.5`$ from the new sample to obtain a best-fit for our ansatz with $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$. We then plot the magnitude deviation of our best-fit universe from an empty universe with $`(\mathrm{\Omega }_{0\mathrm{m}},\mathrm{\Omega }_\mathrm{\Lambda })=(0.0,0.0)`$ in order to illustrate how well our model fits the data (figure 17). For clarity, we plot the median values of the data points. We obtain medians in redshift bins by requiring that each bin has a width of at least 0.25 in log$`z`$ and contain at least 20 SNe. For comparison, we also plot an $`\mathrm{\Lambda }`$CDM $`(\mathrm{\Omega }_{0\mathrm{m}},\mathrm{\Omega }_\mathrm{\Lambda })=(0.3,0.7)`$ model, as well as OCDM and SCDM models. From figure 17 we see that our dark energy reconstruction is a much better fit to SNe beyond $`z0.8`$ than $`\mathrm{\Lambda }`$CDM. At low redshifts ($`z0.1`$) the agreement between data and the two models is rather marginal. We now add 22 of the new supernovae (rejecting one with $`A_V>0.5`$) to our existing dataset of 172 supernovae and perform DE reconstruction on this new dataset of 194 SNe, assuming $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$ and no other priors. The resultant figure 18 is similar to the figure 5(b), with slightly smaller errors and has a best-fit $`\chi _{\mathrm{dof}}^2=1.015`$. The above exercises point to the robustness of results reported in previous sections, and indicate that evolving dark energy agrees well with the full data set containing 194 type Ia SNe.

## 3 The accelerating universe and the Energy Conditions

The energy conditions :

* Strong energy condition: $`\rho +3p0`$ (SEC),
* Weak energy condition: $`\rho 0`$, $`\rho +p0`$ (WEC)

play a vitally important role in our understanding of the accelerating universe, both in the context of inflation and dark energy. We therefore consider it worthwhile to review certain key developments which deepened our understanding of these issues.

In an expanding FRW universe the SEC implies that the universe decelerates while the WEC forbids the pressure from becoming too negative. Additionally, in the 1960’s and early 1970’s it was noted that energy conditions play a crucial role in the formulation of the singularity theorems in cosmology. Indeed, one of the necessary conditions for the existence of an initial/final singularity in big bang cosmology is that matter satisfies both the SEC and WEC (Hawking & Ellis, 1973).

By the late 1970’s it became clear that not all forms of matter satisfy the energy conditions. Perhaps the best example of a form of matter which satisfied the weak energy condition but violated the strong one is the cosmological constant, introduced into cosmology by Einstein in 1917. In addition, the vacuum expectation value of the energy momentum tensor, $`T_{ik}_{\mathrm{vac}}`$, which describes quantum effects (particle production and vacuum polarization) in an expanding universe could, in certain cases, violate both WEC and SEC (Birrell & Davies, 1982; Grib, Mamaev & Mostepanenko, 1980). (For certain space-times, such as de Sitter space, the vacuum energy momentum tensor generates a cosmological constant since $`T_{ik}_{\mathrm{vac}}=\mathrm{\Lambda }g_{ik}`$.) Thus by the late 1970’s it was well known that neither of the energy conditions could be held as being sacrosanct.555The importance of quantum effects to dark energy model building has been emphasised in Sahni & Habib (1998); Parker & Raval (1999).

The 1980’s, as we all know, led to great advances in the development of the inflationary paradigm. The inflaton field mimics the behaviour of a cosmological constant over sufficiently small intervals of time and therefore violates the SEC. Early dark energy models were based on inflaton-type scalars which coupled minimally to gravity (quintessence). Quintessence violates the SEC but respects the WEC. Precisely because of the latter property, not any experimentally obtained $`d_L(z)`$ is compatible with quintessence, as emphasized in Sahni & Starobinsky (2000). (The same observation holds for $`H(z)`$, since the latter can be derived from $`d_L(z)`$ using equation (3).) Clearly if observations do indicate that the WEC is violated by DE then more general (WEC-violating) models for DE should be seriously considered. One example of WEC-violating DE is provided by scalar-tensor gravity. Scalar-tensor models contain at least two functions of the scalar field (dilaton) describing dark energy. As shown in Boisseau et al.  (2000), these two functions, namely, the scalar field potential and its coupling to the Ricci scalar $`R`$, are sufficiently general to explain any $`H(z)`$ obtained from observations.

The WEC can also be effectively violated in DE models constructed in braneworld cosmology. It has recently been shown that such models, with $`w_{\mathrm{eff}}<1`$, are in excellent agreement with supernova data (Alam & Sahni, 2002). Since the field equations in these models are derived from a higher dimensional Lagrangian the unusually rapid acceleration of the four dimensional universe arises because of the full five dimensional theory and not because of matter which continues to satisfy the energy conditions and whose density remained finite and well behaved at all times (Sahni & Shtanov, 2003). This behaviour is in contrast to phantom, which assumes a conventional ‘perfect fluid’ form for the energy-momentum tensor and therefore contains pathological features such as an energy density which diverges in the future and a sound speed which is faster than that of light (Caldwell, 2002; Sahni & Shtanov, 2003).

The fact that the observed luminosity distance (derived from supernova observations) is better fit by dark energy violating the WEC than either quintessence or a cosmological constant was first noticed by Caldwell (2002). Caldwell called this ‘phantom energy’ and showed that larger values of $`\mathrm{\Omega }_{0\mathrm{m}}`$ ($`\mathrm{\Omega }_{0\mathrm{m}}\text{ }\stackrel{>}{}0.2`$) implied increasingly more negative values for the equation of state ($`w\text{ }\stackrel{<}{}1`$) of phantom. 666Note, however, that we will not consider the theoretical model of phantom matter based on a ghost scalar field proposed in this paper since, as is well known, it is unstable with respect to particle creation (particle + antiparticle of the ghost scalar field plus particle + antiparticle of all usual matter fields) and to the loss of spatial homogeneity at both quantum and non-linear classical levels. Caldwell’s results have since been confirmed by larger and better quality SNe data sets – for instance Knop et al.  (2003) find that, in the absence of priors being placed on $`\mathrm{\Omega }_{0\mathrm{m}}`$, the DE equation of state has a $`99\%`$ confidence limit of being $`<1`$ ! Both Caldwell (2002) and Knop et al.  (2003) however work under the assumption that the equation of state of dark energy is unevolving, so that $`w=`$ constant.

In this paper we have shown that, suspending the WEC prior and allowing the dark energy equation of state to evolve brings out dramatically new properties of dark energy. Thus the dark energy model which best fits the SNe observations has an equation of state which rapidly evolves from $`w(z)\text{ }\stackrel{<}{}1`$ at present ($`z=0`$) to $`w(z)0`$ at $`z1`$. Dark energy therefore appears to have properties which interpolate between those of dark matter (dust) at early times and those of a ‘phantom’ ($`w\text{ }\stackrel{<}{}1`$) at late times.

## 4 Conclusions and discussion

This paper reports the model independent reconstruction of the cosmic equation of state of dark energy in which no priors are imposed on $`w(z)`$. In the literature the imposition of various priors frequently precedes the analysis of observational data sets. Such a procedure is well founded and entirely justified when priors are dictated by complementary information such as orthogonal observations coming from different data sets. However, on occasion the use of priors is justified on ‘theoretical grounds’ and in this case one must be careful so as not to prejudge nature. (Compelling theoretical reasons might well reflect our own particular conditioning or set of prejudices !) In the case of the analysis of type Ia supernova data, the priors most frequently used have been $`w=`$ constant and $`w1`$. Both confine DE to within a narrow class of models. Moreover, as shown in Maor et al.  (2002), the imposition of such priors on the cosmic equation of state can, on occasion, lead to gross misrepresentations of reality.

In this paper we do not impose any priors on $`w(z)`$ and reconstruct the equation of state of dark energy in a model independent manner. In this case our best fit $`w(z)`$ evolves from $`w\text{ }\stackrel{<}{}1`$ at $`z=0`$, to $`w0`$ at $`0.8\text{ }\stackrel{<}{}z\text{ }\stackrel{<}{}1.75`$ (the upper limit is set by observations). This result is robust to changes in the value of $`\mathrm{\Omega }_{0\mathrm{m}}`$ and remains in place within the broad interval $`0.1\mathrm{\Omega }_{0\mathrm{m}}0.5`$. Our reconstruction clearly favours a model of DE whose equation of state metamorphoses from $`w=0`$ in the past to $`w1`$ today. An excellent example of a model which has this property is the Chaplygin gas (Kamenshchik, Moschella & Pasquier, 2001). However, in this model dark energy does not violate the weak energy condition (if it was not already violated initially). Our results also lend support to the dark energy models discussed in Bassett et al.  (2002); Corasaniti et al.  (2003) in which the DE equation of state shows a late-time phase transition. An interesting example of an evolving DE model in which $`w_0<1`$ at present whereas $`w(z)>1`$ at earlier times is provided by the braneworld models (called BRANE1) examined in Sahni & Shtanov (2003) which have been shown to agree very well with current supernova observations (Alam & Sahni, 2002).

It is also conceivable that the observed rapid growth in the EOS might characterise ‘unified’ models of dark matter (DM) and dark energy (DE). We end this paper with a small speculation on this last possibility. Since the nature of both DM and DE is currently unknown, it may be that a mechanism exists which converts DM (with $`w=0`$) into DE (with $`w1`$) in regions with sufficiently high density contrast $`\delta \rho /\rho 1`$. (This would happen if, for instance, the rate of conversion of DM into DE depended upon $`(\delta \rho /\rho )^x,x1`$, etc.) Since the conversion of DM to DE is confined to high peaks of the density field this process will not occur uniformly in the entire universe but will be restricted to regions occupying a small filling fraction ($`FF`$) ($`FF1`$ for regions with $`\delta \rho /\rho 1`$; see for instance Sheth et al.  (2003) and references therein). This process could commence as early as $`z1020`$ when the first peaks in a CDM model collapse. Since DE does not cluster and since $`\rho _{\mathrm{DE}}/\rho _{\mathrm{DM}}`$ grows rapidly as the universe expands, DE from high density regions ($`FF1`$) will spread at the speed of light, percolating through the entire universe ($`FF1`$) by $`z1`$. Since the creation of DE is tagged to the formation of structure, this model may not encounter the ‘coincidence problem’ which plagues other scenarios of DE including quintessence. (However this model might have problems in producing a sufficiently homogeneous and isotropic distribution of dark energy on the largest scales.) The concrete mathematical framework for a phenomenological model of this kind will be worked out in a companion paper.

In summary, evolving DE models have been shown to satisfy SNe observations just as well (if not better) than the cosmological constant. Our best fit equation of state, in the absence of any priors, evolves from $`w(z)\text{ }\stackrel{<}{}1`$ at $`z=0`$ to $`w(z)0`$ at $`z1`$. Indeed, figure 17 shows that our best fit EOS is better able to account for the relative brightness of supernovae at $`z\text{ }\stackrel{>}{}0.8`$ than $`\mathrm{\Lambda }`$CDM. However, the evolution in $`w(z)`$ is much weaker if the prior $`w(z)1`$ ($`z0`$) is imposed. Due to the paucity of SNe data beyond $`z=1.2`$ (till date, there is only a single data point beyond $`z=1.2`$, SN1999bf at $`z=1.75`$) it is not clear whether $`w(z)0`$ is a stable asymptotic value for the reconstructed DE equation of state at high redshifts.777An alternative explanation for the relative brightness of SNe at these redshifts, say, by gravitational lensing (Barber et al. , 2000) could clearly alter the high-z properties of our best-fit. New supernova data at $`z\text{ }\stackrel{>}{}1`$ from ongoing as well as planned surveys (SNAP) combined with data from other cosmology experiments (CMB, LSS, S-Z survey’s, lensing, etc.) are bound to provide important insights on the nature of dark energy at high redshifts. Our results clearly throw open exciting new possibilities for dark energy model building.

Acknowledgments:

We would like to thank John Tonry for several important clarifications and for help in preparing figure 17. We also thank Sarah Bridle, Pier-Stefano Corasaniti, Alessandro Melchiorri, Yuri Shtanov and Lesha Toporenskii for their useful comments on an earlier version of the paper. One of us (VS) acknowledges useful discussions with Salman Habib and Daniel Holz.

UA thanks the CSIR for providing support for this work. AS was partially supported by the Russian Foundation for Basic Research, grant 02-02-16817, and by the Research Program “Astronomy” of the Russian Academy of Sciences.

## Appendix A Propagation of errors

We have seen that the error bars on $`w(z)`$ for the analysis using ansatz (7) are non-monotonic with redshift. Low redshift behaviour of the equation of state affects the luminosity distance at all higher redshifts, while high redshift behaviour effects fewer such distances. This leads to an expectation that high-$`z`$ behaviour of the equation of state should be poorly constrained as opposed to the low-$`z`$ behaviour. This seems to contradict the behaviour seen in our figures. To investigate if this could be explained by our specific method of error analysis we describe the Fisher matrix error bars below and show that they are almost identical to what we obtain in our method.

In an analysis which uses an ansatz with $`n`$ parameters $`p_i`$, the Fisher information matrix is defined to be

$$F_{ij}\frac{^2𝐋}{p_ip_j},$$
(26)
where $`𝐋=\mathrm{log}`$, $``$ being the likelihood . For an unbiased estimator, the errors on the parameters will follow the Cramér-Rao inequality : $`\mathrm{\Delta }p_i1/\sqrt{F_{ii}}`$.

Since the likelihood function is approximately Gaussian near the maximum likelihood (ML) point, the covariance matrix for a maximum likelihood estimator is given by

$$(C^1)_{ij}\frac{^2𝐋}{p_ip_j}.$$
(27)
The Fisher information matrix is therefore simply the expectation value of the inverse of the covariance matrix at the ML-point.

Given the covariance matrix, the error on any cosmological quantity $`Q(p_i)`$ is given by :

$$\sigma _Q^2=\underset{i=1}{\overset{n}{}}\left(\frac{Q}{p_i}\right)^2C_{ii}+2\underset{i=1}{\overset{n}{}}\underset{j=i+1}{\overset{n}{}}\left(\frac{Q}{p_i}\right)\left(\frac{Q}{p_j}\right)C_{ij}.$$
(28)
Thus the nature of the errors on a quantity will depend essentially on the manner in which it is related to the parameters of the system.

We now consider how errors propagate for different cosmological quantities for the polynomial fit to dark energy which we have used for most of the results in this paper :

$$H^2/H_0^2=\mathrm{\Omega }_{0\mathrm{m}}x^3+A_0+A_1x+A_2x^2,x=1+z,$$
(29)
where $`A_0=1\mathrm{\Omega }_{0\mathrm{m}}A_1A_2`$. If $`\mathrm{\Omega }_{0\mathrm{m}}`$ is held constant then the parameters of the system are $`(A_1,A_2)`$.

We obtain the covariance matrix in $`(A_1,A_2)`$ from the ML analysis, and then using equation (28), calculate the errors on cosmological quantities of interest. For example, the errors on the quantity $`\mathrm{\Delta }H^2=(H^2H_{\mathrm{\Lambda }CDM}^2)/H_0^2`$ are given by :

$$\sigma _{\mathrm{\Delta }H^2}^2(x)=(x1)^2[C_{11}+2(x+1)C_{12}+(x+1)^2C_{22}].$$
(30)
Although the term $`C_{12}`$ is negative we find that $`\sigma _{\mathrm{\Delta }H^2}^2(x)`$ still increases with redshift. This is shown in the figure 19. The errors shown are approximately similar to those obtained in figure 3.

The corresponding errors on the equation of state can be calculated using equations (18) and (28), and has the somewhat more complicated expression :

$$\sigma _w^2(x)=\frac{x^2[f_1^2C_{11}+2f_1f_2C_{12}+f_2^2C_{22}]}{9[1\mathrm{\Omega }_{0\mathrm{m}}+A_1(x1)+A_2(x^21)]^4},$$
(31)
where

$`f_1`$ $`=`$ $`1\mathrm{\Omega }_{0\mathrm{m}}A_2(x1)^2,`$
$`f_2`$ $`=`$ $`2x(1\mathrm{\Omega }_{0\mathrm{m}})+A_1(x1)^2,`$
and $`A_1,A_2`$ are the mean values of the parameters. Although in this case it is difficult to predict the behaviour of error bars, after substituting the numerical values we obtain the error bars that are shown in figure 20. This figure can be compared to the figure 5(b), having almost identical errors.

This shows that the nature of our error bars is not an artifact of our specific method of error analysis. However, as shown in figure 15, a two parameter expansion in $`w(z)`$ shows monotonically deteriorating errors in $`w(z)`$ with the redshift, while the expansion in $`H_2(z)`$ shows errors that improve with redshift (figure 5(b)). This indicates that the nature of error bars might be affected by which quantity is being approximated. In the limit of infinite terms in the expansion of various quantities all the methods should produce identical result. The practical need for truncating these expansions make these approximations slightly different from each other. More specifically, we require setting of priors

$`f(z)`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}a_nz^n`$ (32)
$`a_n`$ $`=`$ $`0;(n>N_p)`$ (33)
where $`f(z)`$ could be $`H(z)`$, $`w(z)`$ or any other physical quantity and $`N_p`$ is the chosen number of parameters. The non-linear priors in the above equation make different finite expansions inequivalent. Since we do not know for certain if the underlying model for the accelerating expansion involves an energy component with negative pressure in a FRW setting we are forced to choose one of the alternatives for approximations. We hope that with increasingly high quality data the effect of such truncations will eventually disappear.

## Appendix B Reconstruction of Other Dark Energy Models

We have seen in the figure 1 that the ansatz (7) works well for several physically motivated models of quintessence, Chaplygin gas and SUGRA. In this section we take this exercise further and see how well it can reconstruct some of the other fits to dark energy known in literature. In figures 21 (a) and (b), we show results for simulations using $`\mathrm{\Omega }_{0\mathrm{m}}=0.3`$ and two different fits to the equation of state of dark energy :

(a) The fit suggested in Linder (2003) : $`w(z)=w_0+w_1z/(1+z)`$. For this we consider three sets of values in order of increasing evolution of $`w(z)`$ : (a) $`w_0=0.8,w_1=0.5`$, (b) $`w_0=1.0,w_1=1.0`$, and (c) $`w_0=1.2,w_1=2.0`$, and,

(b) The non-perturbative $`w(z)`$ suggested in Corasaniti & Copeland (2003) and Corasaniti et al.  (2003), which has the parameters $`w_Q^0`$ (the dark energy equation of state today), $`w_Q^m`$ (the dark energy equation of state at the matter dominated epoch), $`z_c`$ (the redshift where equation of state changes from $`w_Q^m`$ to $`w_Q^0`$), and $`\mathrm{\Delta }`$ (the width of transition). For the simulation we again use three sets of values in order of increasing growth rate of $`w(z)`$ : (a) $`w_Q^0=0.8,w_Q^m=0.5,z_c^m=1.5,\mathrm{\Delta }=0.1`$, (b) $`w_Q^0=1.0,w_Q^m=0.2,z_c^m=0.6,\mathrm{\Delta }=0.07`$, and (c) $`w_Q^0=1.2,w_Q^m=0.1,z_c^m=0.11,\mathrm{\Delta }=0.03`$.

We find that in both cases, the ansatz recovers the measured quantity to within $`0.5\%`$ accuracy in the redshift range important for SNe observations. Thus we find that even for fits for which the ansatz does not return exact values, it can recover cosmological quantities to a high degree of accuracy.


format: text

text:
# Superfluid current disruption in a chain of weakly coupled Bose-Einstein Condensates

## References


text:
# Mathematical models of haploinsufficiency

## I. Introduction

Complex multicellular organisms are in general diploids, i.e., each cell in an organism contains two copies of the full set of genes in contrast to haploids in which each cell contains a single copy of the genome. Genes provide the blueprint for the synthesis of proteins which perform essential functions in cells. If one copy of a gene is mutated, there is approximately a $`50\%`$ reduction in the level of proteins synthesized. In many cases this does not lead to observable changes and normalcy is retained. A common interpretation of haploinsufficiency (HI) is that it occurs when half normal levels of proteins are insufficient for completing particular tasks, leading to specific types of diseases. More generally, HI may occur when the level of proteins synthesized falls below a critical level for the onset of some desired activity.

There is presently an extensive literature on the genetic and biomedical aspects of HI  but mathematical models exploring the origin of HI and the issues related to it are practically non-existent. It is by now well-accepted that stochastic processes have considerable effect on patterns of gene expression in cells . Cook et al  have studied the role of stochastic gene expression in HI by constructing a minimal model of gene expression and using numerical techniques to simulate the model. Their major finding is that when one of the two genes in a diploid organism is inactivated due to mutations, there is an increased susceptibility to stochastic initiations and interruptions of gene expression. As a result, the number of proteins produced may transiently fall below the desired level giving rise to HI. Both increased gene copy number and faster gene expression kinetics reduce expression noise, thus enhancing the possibility of a stable outcome.

A large number of diseases are caused by mutations in genes encoding proteins called transcription factors (TF). More than $`30`$ different human maladies have been attributed to TF HI . TFs regulate gene expression by binding at the promoter region of the gene to be expressed. Cooperative interactions among the TFs favour the formation of bound TF complexes (oligomers). The TFs interact at only one site or at multiple sites of the promoter. A simple mathematical model has been proposed to explore HI in systems involving cooperative assembly of TFs . Such multimeric complexes are essential for initiation of gene expression in many eukaryotic systems. The model explores the relationship of fractional oligomerization $`Y`$ with the free ($`[S]`$) as well as total concentrations of TFs ($`[S_0]`$). The TFs oligomerise to form a bound complex. The curves $`Y`$ versus $`[S]`$ and $`[S_0]`$ have sigmoidal shapes. Due to the characteristic S shape of a sigmoid, a small change in the TF concentration around the inflection point (the point at which the tangent to the curve has the maximum slope) gives rise to a significant change in the magnitude of $`Y`$. Thus, if there are two TF-encoding genes and one of these becomes silent, the level of TFs produced may fall below the inflection point of the sigmoid and consequently the magnitude of $`Y`$, the fractional oligomerization, is considerably decreased. This results in reduced expression from the target gene, giving rise to TF HI if the amount of proteins synthesized falls below a threshold level.

In Section 2 of this paper, we extend the minimal model of Cook et al  to investigate the influence of bound complexes of TFs on the initiation of gene expression. In Section 3, we study the stochastic version of the minimal model and its extensions to elucidate the role of stochasticity in HI. We derive analytical expressions for the quantities determined numerically by Cook et al .

## 2. Deterministic Model

The model is an extension of the minimal model of gene expression studied by Cook et al . A brief description of the model is as follows. A gene can be in two possible states: inactive ($`G`$) and active ($`G^{}`$). Random transtions occur between the states $`G`$ and $`G^{}`$ according to the first order reaction kinetics

$$G\begin{array}{c}k_a\\ \\ k_d\end{array}G^{}\begin{array}{c}j_p\\ \end{array}p\begin{array}{c}k_p\\ \end{array}\mathrm{\Phi }$$
(1)
where $`k_a`$ and $`k_d`$ are the activation and deactivation rate constants. The corresponding half-times are $`T_a=\frac{log2}{k_a}`$ and $`T_d=\frac{log2}{k_d}`$ respectively. In the active state $`G^{}`$, the gene synthesizes a protein ($`p`$) with the rate constant $`j_p`$. The protein product degrades with a rate constant $`k_p`$ and the associated half-time is $`T_p`$. The protein degradation product is represented as $`\mathrm{\Phi }`$. We now assume that activation to the state $`G^{}`$ is brought about by an inducing stimulus $`S`$, e.g., TFs. The reaction scheme in the presence of the stimulus is given by

$$G+S\begin{array}{c}k_1\\ \\ k_2\end{array}GS\begin{array}{c}k_a\\ \\ k_d\end{array}G^{}\begin{array}{c}j_p\\ \end{array}p\begin{array}{c}k_p\\ \end{array}\mathrm{\Phi }$$
(2)
where $`GS`$ represents the bound complex of $`G`$ and $`S`$ from which transition to the active state $`G^{}`$ occurs. If $`n_G`$ is the total concentration of genes then

$$n_G=[G]+[GS]+[G^{}]$$
(3)
where $`[G]`$, $`[GS]`$ and $`[G^{}]`$ denote the concentrations of genes in the states $`G`$, $`GS`$ and $`G^{}`$ respectively. In the steady state, we have

$$\frac{d[G]}{dt}=k_2[GS]k_1[G][S]=0$$

so that

$$\frac{[G][S]}{K_s}=[GS]$$
(4)
where $`K_s=\frac{k_2}{k_1}`$ is the equilibrium dissociation constant. From (3) and (4), we get

$$[GS]=\frac{n_G[S]/K_s}{1+[S]/K_s}[G^{}]\frac{[S]/K_s}{1+[S]/K_s}$$
(5)
Also, in the steady state,

$$\frac{d[G^{}]}{dt}=k_a[GS]k_d[G^{}]=0$$
(6)
From (5) and (6), the expression for $`[G^{}]`$ in the steady state is given by

$$[G^{}]=\frac{n_Gk_a\frac{[S]/K_s}{1+[S]/K_s}}{k_a\frac{[S]/K_s}{1+[S]/K_s}+k_d}$$
(7)
The reaction scheme in (1) leads to the expression

$$[G^{}]=\frac{n_Gk_a}{k_a+k_d}$$
(8)
in the steady state. Expression (7) and (8) are equivalent on defining effective activation and deactivation rate constants:

$$k_a^{^{}}=k_a\frac{[S]/K_s}{1+[S]/K_s},$$

$$k_d^{^{}}=k_d$$
(9)
We now assume the inducing stimulus to be TFs. In the simplest approximation, $`n`$ individual TFs oligomerise to produce an active complex $`S_n`$ according to the reaction scheme

$$nS\begin{array}{c}K\\ \end{array}S_n$$
(10)
The $`n`$ TFs interact all at once to give rise to the bound complex $`[S_n]`$, i.e., we ignore the formation of dimers, tetramers, ….. etc. Let $`[S_0]`$ be the initial concentration of TFs. Then

$$[S_0]=[S]+n[S_n]$$
(11)
where $`[S]`$ and $`[S_n]`$ are the concentrations of free TFs and the bound TF-complex respectively. The global equilibrium constant $`K`$ is given by

$$K=\frac{[S_n]}{[S]^n}=\frac{[S_n]}{([S_0]n[S_n])^n}$$
(12)
The fractional oligomerization is defined as

$$Y=\frac{n[S_n]}{[S_0]}=\frac{n[S_n]}{[S]+n[S_n]}$$
(13)
Using (11) and (12), $`Y`$ can further be written as

$$Y=\frac{[S]^n}{\frac{[S]}{Kn}+[S]^n}$$
(14)
and

$$Y=\frac{([S_0]n[S_n])^n}{\frac{([S_0]n[S_n])}{Kn}+([S_0]n[S_n])^n}$$
(15)
Fig.1 shows the curve, fractional oligomerization $`Y`$ versus $`[S_0]`$ ( Eq. (15) ) for $`n=6,K=2`$ and $`[S_n]=0.2`$. The curve has the well-known sigmoidal shape.

We now replace $`[S]`$ by $`[S_n]`$ in the reaction scheme described by Eq. (2), i.e., we assume that the TF-oligomer $`S_n`$ binds to a gene in the inactive state $`G`$ to give rise to the bound complex $`GS_n`$. Transition to the active state $`G^{}`$ occurs from the intermediate state $`GS_n`$. The concentration $`[G^{}]`$ in the steady state is obtained from (7) by replacing $`[S]`$ by $`[S_n]`$ where $`[S_n]=K[S]^n`$ ( Eq. (12) ). One finally obtains

$$[G^{}]=\frac{n_Gk_a\frac{K[S]^n/K_s}{1+K[S]^n/K_s}}{k_a\frac{K[S]^n/K_s}{1+K[S]^n/K_s}+k_d}$$
(16)
The concentration of proteins in the steady state is given by

$$[p]=\frac{j_p}{k_p}[G^{}]$$
(17)
From (14), $`[S]^n`$ can be written as

$$[S]^n=\frac{Y}{1Y}\frac{[S]}{Kn}$$
(18)
From (16) and (17) and for $`k_a=k_d`$, we get

$$[p]=n_G\frac{j_p}{k_p}\frac{aY}{1+2aYY}$$
(19)
where $`a=\frac{[S]}{nK_s}.`$

Fig. 2 shows the protein concentration $`[p]`$ versus fractional oligomerization $`Y`$ for $`k_a=k_d`$ (Eq.(19)), $`n_G=2`$, $`j_p=0.5k_p`$, $`K_s=1.2`$ and $`n=6`$.

From Eqs. (11), (16), and (17) and for $`k_a=k_d`$, the concentration of protein $`[p]`$, as a function of the total concentration $`[S_0]`$ of TFs, can be written as

$$[p]=n_G\frac{j_p}{k_p}\frac{([S_0]n[S_n])^n}{2([S_0]n[S_n])^n+K_s/K}$$
(20)
Fig. 3 shows the plot $`[p]`$ versus $`[S_0]`$ (Eq. (20)) for $`n=6`$, $`K=2`$, $`K_s=1.2`$ and $`[S_n]=0.2`$. Figs. 1, 2 and 3 provide a possible explanation for the origin of HI. Suppose two TF-encoding genes produce TFs of total concentration $`[S_0]=4`$. If one of the genes is inactivated due to mutations, the total concentration $`[S_0]`$ falls to the value $`2`$. For the parameter values corresponding to Fig. 1, the fractional oligomerization $`Y`$ has the value $`0.797`$. The TFs form a bound complex $`S_n`$ ($`n=6`$) which then activates the gene synthesizing the protein $`P`$. The concentration of proteins $`[p]`$ corresponding to $`Y=0.797`$ and for $`[S]=0.8`$ is given by $`[p]=0.233`$ (Fig. 2). The same value $`[p]`$ is obtained from Fig. 3 with $`[S_0]=2.`$ If both the encoding genes are active, the values of $`[S_0]`$, $`Y`$ and $`P`$ are $`[S_0]=4,`$ $`Y=1.0`$ and $`[p]=0.5`$ respectively. Thus, for one gene, the protein level is reduced by more than half. If the level falls below a threshold, the amount of proteins synthesized is not sufficient for the execution of a particular task. This gives rise to HI, the manifestation of which is in the form of a disease.

## 3. Stochastic Approach

Let us first consider the model described in Section 2 in the absence of an inducing stimulus. Let $`n_{tot}`$ be the total no of genes and $`n_0`$ , $`n_1`$ ($`n_{tot}=n_0+n_1`$), the number of genes in the inactive ($`G`$) and active ($`G^{}`$) states respectively. In the stochastic model, a gene makes random transitions between the inactive and active states with $`k_a`$ and $`k_d`$ being the activation and deactivation rate constants. In the active state, protein production and degradation occur with the rate constants $`j_p`$ and $`k_p`$ respectively. Let $`p(n_1,n_2,t)`$ be the probability that at time $`t,`$ $`n_1`$ genes are in the active state $`G^{}`$ and the number of protein molecules is $`n_2.`$ The rate of change of the probability with respect to time is given by the Master Equation

$$\begin{array}{cc}\frac{p(n_1,n_2,t)}{t}=& k_a[(n_{tot}n_1+1)p(n_11,n_2,t)(n_{tot}n_1)p(n_1,n_2,t)]\\ & +k_d[(n_1+1)p(n_1+1,n_2,t)n_1p(n_1,n_2,t)]\\ & +j_p[n_1p(n_1,n_21,t)n_1p(n_1,n_2,t)]\\ & +k_p[(n_2+1)p(n_1,n_2+1,t)n_2p(n_1,n_2,t)]\end{array}$$
(21)
For each rate constant, there is a gain term which adds to the probability and a loss term which subtracts from the probability.

We now use the standard approach in the theory of stochastic processes  to determine the average number of activated genes $`<n_1>`$ and proteins $`<n_2>`$ in the steady state and the variances thereof. Define the generating function

$$F(z_1,z_2,t)=\underset{n_1,n_2}{}z_1^{n_1}z_2^{n_2}p(n_1,n_2,t)$$
(22)
In terms of the generating function, the Master equation (21) becomes

$$\frac{F}{t}=k_an_{tot}(z_11)Fk_a(z_11)z_1\frac{F}{z_1}k_d(z_11)\frac{F}{z_1}+j_p(z_21)z_1\frac{F}{z_1}k_p(z_21)\frac{F}{z_2}$$
(23)
In the steady state $`\frac{F}{t}=0.`$ The following properties of the generating function are used in subsequent calculations:

$$F_{z_1=1,z_2=1}=1$$
(24)
$$<n_1>=\frac{F}{z_1}_{z_1=1,z_2=1},<n_2>=\frac{F}{z_2}_{z_1=1,z_2=1}$$
(25)
where $`<n_1>`$ is the mean number of active genes, i.e., genes in the state $`G^{}`$ and $`<n_2>`$ is the same for proteins. Furthermore,

$$\begin{array}{c}\frac{^2F}{z_1^2}_{z_1=1,z_2=1}=<n_1^2><n_1>\\ \frac{^2F}{z_2^2}_{z_1=1,z_2=1}=<n_2^2><n_2>\end{array}$$
(26)
Hence the variances around the mean levels are given by

$$\begin{array}{c}Var_{n_1}=<n_1^2><n_1>^2=\frac{^2F}{z_1^2}_{z_1=1,z_2=1}+<n_1><n_1>^2\\ Var_{n_2}=<n_2^2><n_2>^2=\frac{^2F}{z_2^2}_{z_1=1,z_2=1}+<n_2><n_2>^2\end{array}$$
(27)
Successive differentiation of Eq. (23) ($`\frac{F}{t}=0.`$) with respect to $`z_1`$ and $`z_2`$ gives rise to linear equations for successively higher moments. The equations may be solved to obtain, in particular, the mean and the variance. For example, differentiating Eq. (23) with respect to $`z_1`$ and $`z_2`$ and then putting $`z_1,z_2=1`$, one obtains expressions for the mean.

The mean and variance are given by

$$<n_1>=\frac{n_{tot}k_a}{k_a+k_d}$$
(28)
$$Var_{n_1}=<n_1>\frac{k_d}{k_a+k_d}$$
(29)
$$<n_2>=<p>=<n_1>\frac{j_p}{k_p}=\frac{j_p}{k_p}\frac{n_{tot}k_a}{k_a+k_d}$$
(30)
$$Var_{n_2}=<n_1>\frac{j_p}{k_p}[1+\frac{j_pk_d}{(k_a+k_d)(k_a+k_d+k_p)}]$$
(31)
As in Ref. 3, temporal quantities are scaled relative to the product half-life $`T_p=\frac{log2}{k_p}.`$ Let $`T_a=\frac{log2}{k_a}`$ and $`T_d=\frac{log2}{k_d}`$ be the times for half-maximal gene activation and deactivation respectively. The times $`T_a`$ and $`T_d`$ are scaled relative to $`T_p`$. Some of the results obtained in Ref. 3, using numerical simulation techniques, can readily be derived from the analytical expressions in (28)-(31). Stochasticity introduces random fluctuations around the mean protein level and variance gives a measure of the spread. Let $`T_a=T_d=T_p/4`$, i.e., $`k_a=k_d=\alpha k_p`$ with $`\alpha >0`$. As $`\alpha `$ increases, one has faster expression kinetics and from (31) it is easy to verify that variance is reduced, i.e., the expression noise is less. The mean product level (Eq. (30)) is, however, independent of $`\alpha `$. With increase in $`j_p`$, i.e., the protein synthesis rate, the variance increases. Let us now consider the case when the net expression rate of $`n_{tot}`$ genes is distributed to one single gene so that the mean protein level remains the same. From (30)

$$<n_2>=\frac{j_p}{k_p}\frac{n_{tot}}{2}=\frac{j_p^{^{}}}{k_p}\frac{1}{2}$$
(32)
where $`j_p^{^{}}=j_pn_{tot}`$ is the expression rate when only one gene is considered. Since $`j_p^{^{}}>j_p`$, the gene copy number is reduced from $`n_{tot}`$ to $`1`$. Similarly, when the net expression rate of $`n_{tot}`$ genes is distributed to a larger number genes, say, from two to four, the variance is reduced.When one of two genes is inactivated due to mutations, the average protein level in the steady state is reduced by $`50\%`$. This may still be higher than the threshold level required for protein activity. Due to the variance around the mean level, the number of proteins may transiently fall below the threshold giving rise to HI. The occurrence of HI further becomes more probable for slower expression kinetics as then the variance is increased in magnitude. For stochastic gene expression in the presence of an inducing stimulus, say, TF’s, we use the effective model with the activation/deactivation rate constants given in Eq. (9). The expressions for the mean and the variance are the same as in Eqs. (28)-(31) but with $`k_a`$, $`k_d`$ replaced by $`k_a^{^{}}`$ and $`k_d^{^{}}`$ respectively.

We now derive expressions for the probability distributions of protein levels in the steady state. To do this, we consider a simpler stochastic model in which the only stochasticity arises from the random transitions of a gene between the inactive and active states. In each state of the gene, the concentration of proteins evolves deterministically according to the equation

$$\frac{dx}{dt}=\frac{j_p}{X_{max}}zk_px=f(x,z)$$
(33)
where $`z=1`$ $`(0)`$ when the gene is in the active (inactive) state and $`x=\frac{X}{X_{max}},`$ $`X`$ and $`X_{max}`$ being the protein concentration at time t and the maximum protein concentration respectively. We note that $`X_{max}=\frac{j_p}{k_p}.`$ Let $`p_j(x,t)`$ $`(j=0,1)`$ be the probability density function when $`z=j.`$ The total probability density function is

$$p(x,t)=p_0(x,t)+p_1(x,t)$$
(34)
The rate of change of probability density is given by

$$\frac{p_j(x,t)}{t}=\frac{}{x}[f(x,j)p_j(x,t)]+\underset{kj}{}[W_{kj}p_k(x,t)W_{jk}p_j(x,t)]$$
(35)
where $`W_{kj}`$ is the transition rate from the state $`k`$ to the state $`j`$ and $`W_{jk}`$ is the same for the reverse transition. The first term in Eq. (35) is the so called “transport” term representing the net flow of the probability density. The second term represents the gain/loss in the probability density due to random transitions between the state $`j`$ and other accessible states. In the present case, Eq. (35) gives rise to the following two equations:

$$\frac{p_0(x,t)}{t}=\frac{}{x}(k_pxp_0(x,t))+k_dp_1(x,t)k_ap_0(x,t)$$
(36)
$$\frac{p_1(x,t)}{t}=\frac{}{x}\{(\frac{j_p}{X_{max}}k_px)p_1(x,t)\}+k_ap_0(x,t)k_dp_1(x,t)$$
(37)
The Master equation (Eq. (21)) provides a full stochastic description of all the processes associated with gene expression, namely, gene activation and deactivation, protein synthesis and degradation. The only stochastic events considered by Cook et al  are those related to gene activation and deactivation. In their model, protein synthesis from the active gene and protein degradation occur in a deterministic manner. Eqs. (36) - (37) describe the scenario studied by Cook et al. The steady state distribution in this case is given by

$$p(x)=Cx^{(\frac{k_a}{k_p}1)}(1x)^{(\frac{k_d}{k_p}1)}$$
(38)
where $`C`$ the normalization constant is given by the inverse of a beta function.

$$C=\frac{1}{B(\frac{k_a}{k_p},\frac{k_d}{k_p})}$$
(39)
Since the probability density function is known, the mean protein level and its variance can be calculated in a straightforward manner. The mean protein level is identical to that obtained from the Master equation (Eq. (21)) whereas the variance is underestimated as stochasticity is taken into account only at the levels of gene activation and deactivation. Fig. 4 shows the plot of $`p(x)`$ versus $`x`$ for slow ($`T_a=T_d=T_p/4`$) and fast ($`T_a=T_d=T_p/40`$) gene expression kinetics. In the latter case, the distribution is significantly narrower, i.e., faster kinetics lead to a reduction in the variance. The same conclusion is reached from the Master equation approach. From the full width at half maximum of the broader distribution, one finds that $`x`$ ranges from 0.25 to 0.75, 0.5 being the mean value. Thus, it is probable that the protein level falls below the threshold for desired activity giving rise to HI. Let $`x_{thr}`$ ($`<1`$) be the threshold value of $`x`$. The probability that $`x`$ is greater than $`x_{thr}`$ is

$$p(x>x_{thr})=1\frac{_0^{x_{thr}}x^{(\frac{k_a}{k_p}1)}(1x)^{(\frac{k_d}{k_p}1)}𝑑x}{_0^1x^{(\frac{k_a}{k_p}1)}(1x)^{(\frac{k_d}{k_p}1)}𝑑x}$$
(40)
$$=1\frac{_{k_px_{thr}^{\frac{k_a}{k_p}}{}_{2}{}^{}F_{1}^{}[1\frac{k_d}{k_p},\frac{k_a}{k_p},1+\frac{k_a}{k_p},X_{thr}]}}{k_aB(\frac{k_a}{k_p},\frac{k_d}{k_p})}$$
(41)
where $`{}_{2}{}^{}F_{1}^{}(a,b,c;z)`$ is the hypergeometric function.

Let $`x_{thr}`$ be $`0.25`$. The probability $`p(x>x_{thr})`$ is computed using Mathematica for both slow and fast gene expression kinetics. The values of $`p(x>x_{thr})`$ in the slow and fast cases are $`0.9294`$ and $`0.9999`$ respectively. Since in the latter case, the probability that the protein level exceeds the threshold is higher, the chance of HI occurrence is correspondingly lower.

In the presence of an inducing stimulus, say, TFs, the probability of activation above a threshold is again given by (40) with $`k_a`$ and $`k_d`$ replaced by $`k_a^{^{}}`$ and $`k_d^{^{}}`$ (Eq. (9)) . Fig. 5 shows $`p(x>x_{thr}),`$ $`x_{thr}=0.25`$, versus $`S/K_s`$ in a semi-logarithm plot for both slow and fast kinetics. In the fast case, a substantially steeper curve is obtained leading to enhanced signal discrimination, i.e., a more predictable response in a gradient of inducing signal. As shown by Cook et al , the signal discrimination ability increases with gene copy number. One can thus speculate that diploid organisms utilise stochastic expression kinetics, preferably fast, for signal discrimination and are susceptible to degraded signal discrimination due to a reduction of gene copy number in the haploid state. Mutations in the subset of genes which generate a response to signaling gradients in diploid organisms may be the cause of some HI syndromes associated with these systems. Fig. 6 shows the same plot as in Fig. 5 but now the TF’s form bound complexes with $`n=6.`$ In Eq. (9), $`[S]`$ is replaced by $`[S_n]=K[S]^n`$ ( Eq. (12)). One now finds that the slopes of the curves for slow and fast kinetics are similar. Thus as $`n`$, the number of TF’s forming the bound complex, increases, the distinction between slow and fast gene expression kinetics, as regards their signal discrimination ability, becomes less pronounced.

## 4. Summary

In this paper, we have studied simple mathematical models to explore the possible origins of HI. In Section 2, we have considered a deterministic model in which a complex of $`n`$ TFs binds at the appropriate region of DNA to initiate gene expression in eukaryotes. The concentration of proteins synthesized versus the total concentration of TFs (Fig. 3) is a sigmoid. Due to the S-shape of the curve, the protein level may fall below a threshold when one of the two genes synthesizing the TFs is mutated, resulting in a $`50\%`$ reduction in the total concentration of TFs. The absence of required protein activity can give rise to HI. In Section 3, we have studied simple stochastic models of gene expression and shown that due to random fluctuations around the mean protein concentration in the steady state, the protein level may fall below the threshold even though it does not do so in the deterministic case. The variance, a measure of the spread around the mean protein level, is reduced with increasing gene copy number and faster expression kinetics. The variance increases if the rate constant $`j_p`$ associated with protein synthesis is increased. In the case of one gene, we have further calculated the probability that the concentration of proteins exceeds a threshold in the absence as well as the presence of an inducing stimulus. In the latter case, faster gene expression kinetics give rise to a sharper response to changing stimulus concentrations. As shown by Cook et al , this is also true when the gene copy number is increased. Thus the signal discrimination ability of diploid organisms may be impaired in the haploid state. When the inducing stimulus is a bound complex of $`n`$ TFs, the distinction between slow and fast gene expression kinetics becomes less with increasing $`n`$. To sum up, we have considered both deterministic as well as stochastic models of gene expression and indicated possible scenarios for the occurrence of HI.

ACKNOWLEDGEMENT

R. K. was supported by the Council of Scientific and Industrial Research, India under Sanction No. 9/15 (239) / 2002 - EMR - 1.