PhysRevD.88.103528

Nonparametric reconstruction of dynamical dark energyvia observational Hubble parameter dataHao-Ran Yu,Shuo Yuan,and Tong-Jie Zhang*Department of Astronomy,Beijing Normal University,Beijing100875,China(Received10August2013;published26November2013)We study the power of current and future observational Hubble parameter data(OHD)on non-parametric estimations of the dark energy equation of state(EoS),wðzÞ.We propose a new method by conjunction of principal component analysis and the goodness offit criterion to reconstruct wðzÞ,ensuring the sensitivity and reliability of the extraction of features in the EoS.We also give a new error model to simulate future OHD data,to forecast the power of future OHD on the EoS reconstruction.The result shows that current OHD,despite their lesser quantity,give not only a similar power of reconstruction of dark energy compared to the result given by type Ia supernovae,but also extend the constraint on wðzÞup to redshift z’2.Additionally,a reasonable forecast of future data in greater quantity and better quality greatly enhances the reconstruction of dark energy.DOI:10.1103/PhysRevD.88.103528PACS numbers:98.80.Es,95.36.+xI.INTRODUCTIONThe accelerating expanding Universe is explained bydark energy[1,2]and one of the main challenges of themodern physical cosmology is to study the nature of darkenergy[3].Its evolution and dynamical properties are char-acterized by its equation of state(EoS)w p= ,the ratioof its pressure and energy density.A cosmological constantis usually the simplest model to explain dark energy.Its EoSdoes not depend on redshift z and remains constant:wðzÞ¼À1.Variations of dark energy models,such as quintessence, phantom,quintoms,etc.,are also used to explain darkenergy.Studying dark energy in a parametrized way,suchas Chevallier-Polarski-Linder parametrization[4,5],mayinduce misleading results due to prior assumptions of func-tion forms of the EoS;then we may do null-test diagnosticson dark energy[6]or use model-independent analyses toreconstruct the EoS nonparametrically[7–9].This has been achieved by astronomical data,such as luminosity distances d L of type Ia supernovae(Ia SN)[10,11].Although we have plenty of Ia SN data,they require the second derivative of d L with respect to redshift z and a complicated form of reconstruction equation.Meanwhile,observational Hubble parameter data(OHD)HðzÞreconstruct wðzÞvia a simpler way[Eq.(2)]and require only thefirst derivative of HðzÞ,so this method is more stable and less sensitive to the error of the data.During the last decade,OHD have proved to be very powerful in the constraint of cosmological parameters [12];they correspond to the cosmic expansion history directly and can be measured by various methods,e.g., cosmic chronometers[13],baryonic acoustic oscillation (BAO)peaks[14],and even other proposed Sandage-Loeb (SL)probes[15,16]and the standard siren by gravitational waves.On the way toward obtaining more data from observations,to simulate sets of future OHD appropriately is also instructive at the current stage in exploring the quality of EoS reconstruction in the near future.To do this,we need an error model giving the redshift distribution,the offsets from the theoretical value,and the sizes of their error bars, such as in[17].Here we construct a new and more accurate error model for next generation surveys,carefully estimat-ing the potential power of observation in different methods. For either observational or our simulated data,nonpara-metrically reconstructing HðzÞand wðzÞfrom data confronts the tradeoffs between extracting more features and avoiding the overfitting of errors—if we use fewer parameters,we are losing details of the underlying model,while using more parameters increases the detectability of such details,but it brings the risk of being polluted by the errors.Principle component analysis(PCA)is usually used to exclude error-induced oscillations in the reconstructing processes[10], and information criteria[18]are often used for choosing models[11].We alsofind that the goodness offit(GoF)is an effective probe of overfitting.The conjunction of the PCA and GoF criterion works naturally to deal with this tradeoff. The rest of the paper is organized as follows.We present our reconstruction method in Sec.II and briefly illustrate the current OHD and forecast future OHD in Sec.III.We show the results of our method in Sec.IV and conclude in Sec.V.II.METHODWe briefly show the mathematics of the reconstruction method in Sec.II.To test the method we use the new error model for simulated data,which is shown in Sec.II B,and the GoF criterion is discussed in Sec.II C.A.ReconstructionIn general relativity,the expansion of the Universe is affected by the densities of the components and the state of dark energy:*Corresponding author.tjzhang@PHYSICAL REVIEW D88,103528(2013)H ðz Þ2¼H 20½ M ð1þz Þ3þ K ð1þz Þ2þ Ãg ðz Þ ;(1)where H 0, M , K , Ãare the Hubble constant,current energy densities in pressureless nonrelativistic matter,cur-vature,and a smooth dark energy component,respectively,and M þ K þ ü1;and g ðz Þ¼exp ½3R z 0ð1þw ðz 0ÞÞ=ð1þz 0Þd z 0 .Here we neglect the current energy density of radiation R and assume the Universe to be dominated by nonrelativistic matter and dark energy.By expressing w ðz Þin terms of H ðz Þand H 0ðz Þ d H=d z ,we get the recon-struction equation:w ðz Þ¼3H 2À2ð1þz ÞHH 0À K H 20ð1þz Þ23H 20ð1þz Þ2½ M ð1þz Þþ KÀ3H 2:(2)Given a set of H ðz Þdata y ,we fit them with a smooth analytical function and then use it and its first derivative to reconstruct w ðz Þusing Eq.(2).Assume the data y ¼ðy 1;y 2;...;y n ÞT with covariance matrix C .By choos-ing a set of N (N <n )primary basis functions such as polynomials,rational functions,or wavelets,X ¼ðx T 1;x T 2;...;x TN Þ,where x i is a row vector representing the i th basis,we fit y by linear least squares.Simply minimize the weighted squared residual R ¼ðy ÀX ÞT C À1ðy ÀX Þby solving @R =@ ¼0,and we get the coefficient vector¼ðX T C À1X ÞÀ1X T C À1y :(3)With enough realizations from the error model,we diagonalize the inverse covariance matrix of by finding its N eigenvalue and eigenvectors:F ¼Cov À1ð Þ¼E ÃE T ,where E ¼ðe T 1;e T 2;...;e TN Þrearranges the N pri-mary basis functions to N new orthogonal eigenbasis U ¼X E with the corresponding eigenvalues diag ðÃÞ¼ð 1; 2;...; N ÞT .The first M (M N )eigenmodesU M ¼X E M ,where E M ¼ðe T 1;e T 2;...;e TM Þ,the principal components of the data,are usually of our interest,reflect-ing the main features of data.Conversely,the rest compo-nents are induced by the random offsets in the data and may induce overfittings.After the rearrangement,we solve Eq.(3)again and get new coefficients for the new basis:M ¼ðU T M C À1U M ÞÀ1U T M C À1y :(4)Because we have enough realizations to estimate Cov ð Þ,and any single new realization or observational data,hav-ing similar configurations,their true features can still be extracted and the nuisance features are excluded.With the reliable smooth H ðz Þ,getting w ðz Þis straightforward from Eq.(2),and its errors can be estimated from realizations.B.An optimal error modelUsually simulated offsets and errors are regarded as Gaussian distributed [17].More generally,the error,arandom variable ,is determined by several factors t ¼ðt 1;t 2;...;t k ÞT ,where these factors are assumed to be Gaussian distributed with covariance Æt .These factors are different sources of errors that jointly contribute to with a quadratic form 2¼t T D t ,and the matrix D de-notes the relation between the factor and the final error,up to second order accuracy.In such case, has a Nakagami m distribution f m ðx ;m; Þ(a scaling transform of generalized- distribution with parameters Æt and D )[19]and has a complicated form.We can never estimate the full contribution of these factors (weather,telescope,device,recording during observation,as well as data re-duction and systematic errors,etc.);however,we simplify the problem by assuming t i $N ð0;1Þ,Æt ¼D ¼I :they are independent,having the same importance,and contrib-ute to 2additively.In this case, $ k .More realisti-cally,there are few dominant factors of much greater importance (Æt ÞI )that jointly contribute to the final error.This effectively causes a reduction in the degree of freedom,k !k 0(k 0<k ),and a scaling transform ( k 0!A k 0).Practically,neglecting the mathematical details,we find that the m distribution f m ðx ;m; Þ(m ¼k 0=2, ¼A 2=k 0)well matches the relative error H =H ’s distribution and does not depend on redshift z .The redshift distribution is generated according to the configuration of each method of measurement.For a large enough sample,we find that uniformly distributed or evenly spaced samples are good approximations.For errors and offsets,we assume that the data points are indepen-dent,so the relative error is successively generated from the corresponding distribution =H ?ðz Þ$f i ,where f i is the relative error distribution for the i th kind of measure-ment.This means that we are ready for a measurement without knowing the true value y ?¼H ?ðz Þ,but the quality (error bar)is predetermined.For h f i i (1,the resulting measurement should have the distribution y $N ðy ?; Þ,and it gives the simulated data.C.Goodness of fitThe remaining problem is how to choose the number of primary basis functions N and the number of reserved principal components M appropriately.The choice of ðN;M Þdirectly relates to the upper limit of the ability of how completely we can detect the features of the under-lying model.If we suppose that the fitted model is just the underlying model,the offset y Ày ?,having the distribu-tion of N ð0; Þ,contributes to the residual R ?by ðN ð0;1ÞÞ2for each independent data point,and R ?has 2distribution R ?$ 2n .Thus,h GoF ?iR ?n ¼ðy Ày ?ÞT C À1ðy Ày ?Þn ¼1:(5)Here we use ?to denote that the fitted model is replaced bythe underlying model,and also define R ?=n as the good-ness of fit for underlying model GoF ?.Going back to theHAO-RAN YU,SHUO YUAN,AND TONG-JIE ZHANGPHYSICAL REVIEW D 88,103528(2013)resultingfitted model,we define the GoF in our circum-stance as1GoFR MnÀ1;(6)where R M is the least square R min given by thefirst M principal components U M:R M¼ðyÀU M MÞT CÀ1ðyÀU M MÞ:(7) Statistically,h GoF i should also be unity.If it is lessthan1,it indicates the overfitting of offsets due to the measurement error.If GoF is larger than1by a greatamount,it means that we have not yet captured the fullfeatures of the data;i.e.,the model is too simple tofit allthe features.We use GoF as the indicator to determine thecomplexity of ourfitted model.Here,the complexity of thefitted model is less than or equal to the true underlyingmodel,because errors and offsets lower the detectability ofthe underlying model.More specifically,for a very com-plex oscillating wðzÞ,as an example,it always allows itsgenerated data to have GoF?’1,and we are supposed to use higher N and M tofit its complex features.However,ifthe generated data are dominated by noise,such configu-ration(N and M)leads to illness infitting low quality data,because the errors and offsets are more dominant than thetrue features,and thus are amplified,and any subtler fea-tures are covered.As a result,the reconstructed model is no longer the underlying model and usually has much smaller GoF(GoF<GoF?’1).For a given large enough N,full usage of all eigenmodes(M¼N)leads to the overfitting,and we degrade thecomplexity of the model by reducing M consecutivelywhile examining h GoF i over realizations,until wefindthe last M to have h GoF i>1.It is still safe to use thisgrade of complexity.An insufficient number of primary basis functions N causes two consecutive M’s to have a great difference in h GoF i,skipping the range of appropriate fitting(h GoF i’1),since the eigenmodes are poorly deter-mined from a limited number of primary functions,or the primary functions X are chosen improperly.In these cases we should either increase N or choose more suitable X. Without clear knowledge of the essence of dark energy, the underlying complexity of wðzÞresults in innumerable forms of HðzÞ’s,so generally there is not an optimal choice of primary function basis X.However,we can still choose some popular,reasonable models to see how different X’s have an effect on thefiually polynomials are more effective than other rational functions for not-too-complexmodels,and they form stable eigenmodes that are invariantfor different values of N and various underlying HðzÞ,even if the underlying HðzÞis very oscillatory.III.A V AILABLE AND FUTURE HUBBLEPARAMETER DATACurrent OHD are obtained primarily by the method ofcosmic chronometers[20–23].Other methods to extractHðzÞare by the observations of BAO peaks[14,24]and Ly- forest of luminous red galaxies(LRGs)[25],the lastof which extended the current OHD deeply to z¼2:3. For the real data applied to our method,it is straightfor-ward to test their GoF?.Assuming the independence of28 measurements,to calculate the residual between the data and afiducial theoreticalÃCDM(cold dark matter with a constant cosmological constantÃ)model gives GoF?¼0:62<1,meaning that for independent data with current error level,their offsets would have been larger.It indicates that the assumption of independence for the current28 OHD is not proper.However,we do not have a good estimation of the off-diagonal elements of the OHD’s covariance matrix.To weigh the power of current OHD and compare with simulations,we rebin the data from LRGs that have close redshifts,and getfinally15mea-surements of OHD with smaller error bars(see Fig.1). Now for the binned data GoF?¼1:02,which is close to independent measurements.We have a larger residual contributing to R statistically from each binned datum, and it is reasonable from the error.Although this is only a rough estimation,it is still more accurate than the28data points with no knowledge of their covariance.The binned and unbinned error distributions for OHD well match the shape of m distribution(see the left panel of Fig.1).We also assume the relative error of future data to have such distribution.There are several ways to enhance the quality of OHD,ofwhich deeper-redshift,more-complete-sky-coverage LRGsurvey and spectroscopic observations of those identifiedLRGs give remarkable improvement.For example,the2SLAQ2has provided a LRG catalogue with the redshiftrange from0.3to0.9with180deg2coverage of the sky.Considering a future LRG survey with more than half skycoverage with redshift range z$2:3,it may give several millions of LRGs and effectively enlarges the number of OHD data points or lowers their ing a LRG sample binning strategy[20],we can extract up to100 OHD measurements with20%of the present error level. Other methods also have potential:[24]gives three OHD by processing the WiggleZ3survey andfitting the BAO peak parameter and the2D power spectrum at three red-shift slices.The raw WiggleZ samples are158741galaxies1Note that,although it is usual to scale RM by the degree offreedom ¼nÀMÀ1,here we instead use nÀ1to scale R M for the estimated underlying model from n data samples,because various choices of primary basis functions and inertial complex-ity of the underlying model lead to different numbers of pa-rameters needed.Even given afixed underlying model,resulting N and M depend on X.Thus,the expression of GoF should not be scaled by anything in terms of N or M in this circumstance.2/.3.au/site/.NONPARAMETRIC RECONSTRUCTION OF DYNAMICAL...PHYSICAL REVIEW D88,103528(2013)in the range z ¼ð0:2;1:0Þwith 800deg 2coverage of the sky.We can also forecast a more-than-half-sky survey like WiggleZ,which offers up to five million targets and thus lowers the error level to several percent of the present level.The more ambitious experiment for Sandage-Loeb signal (SLS),proposed by [15,16],may further deepen the red-shift range of OHD.It measures the quasar Ly- forest in a separation of a few decades with an extremely large telescope like the European Extremely Large Telescope (E-ELT).4The upcoming CODEX (COsmic Dynamics and EXo-Earth experiment)5is based on the E-ELT and offers a measurement of SL signal.If successful,it can help us to explore the redshift from 2to 5covering the ‘‘redshift desert’’and give useful data for the expanding history of the Universe [26].The CODEX group provided a full design for observing the SL signal and the prediction of the statistical error of SLS [27].We use the SL signal,with an error estimation by [28](a 15year observational inter-val is assumed),as optional simulated data to study their impact on the result.IV .RESULTSWe set cosmological parameters M , K , Ã,and H 0as in the last Planck data release [29].For the diagnostic EoS models,we choose ÃCDM with w Ãðz Þ¼À1and two other arbitrary models w 1ðz Þ¼À1=2þerf ðln ð2z=e ÞÞand w 2ðz Þ¼À1À0:31sin ð12ln ð1=ð1þz ÞÞÞ,where in model 1,w 1smoothly varies from À1to 0as z increases,while in model 2,w 2is very oscillatory.The theoretical H ðz Þcurves based on these two models are shown in dash-dotted lines in the top right panel of Fig.1.We also show the dashed blue line,representing the best-fit analytical H ðz Þcurve by using 15binned OHD.Its reconstruction of w ðz Þis shown in the bottom right panel of Fig.1.We do not show errors for these two dashed blue lines—usually Monte Carlo realizations are run on each data point,y MC $N ðy; Þ,to get the statistical properties.However,recall that y $N ðy ?; Þ,so y MC $N ðy ?;ffiffiffi2p Þ:its error is amplified,and also we have only one realization—real observation,y MC ,is biased by y Ày ?,due to the cosmic variance.For this reason we can only get an error estimation based on a supposed underlying model,e.g.,ÃCDM ,as-suming y is just one realization of y ?.We calculate the statistics from other realizations of y ?,y MC $N ðy ?; Þ,and see if the reconstruction from the real data is within the error region of y MC .We use our error model to simulate data and do this Monte Carlo realization.Here we use only cosmic chronometer data,with 100independent measure-ments,20%of the present error level,evenly distributed on 0<z <2:3.The results are shown with expectations (or-ange lines)and 1 of their error regions (translucent orange areas)in Fig.1.For H ðz Þ,it reasonably covers the H ðz Þby underlying ÃCDM .w ðz Þis confident when z &1:5,while beyond this range the reconstruction is biased due to the scarceness of data,and when z ’2it is hopeless because H 0is poorly paring the result (blue dashed line)with the error region,we still see an obvious $1:5 deviation from ÃCDM at 1:6&z &2,and it does not correlate with the bias of the central line of the error region.Although we note that this feature is also similar to theresultFIG.2(color online).Top left:Reconstruction of w ðz Þby OHD’s redshift and error configurations,but assuming model 1.Also shown as the blue dashed line is the reconstruction of the OHD set,which is obviously inconsistent with model 1.The remaining three panels show the reconstruction of w ðz Þfrom simulated future data,with the underlying models being ÃCDM ,model 1,and model2.FIG.1(color online).OHD.The purple and green of unbinned 28OHD and curves are their best fits with an m distribution.Their counterpart CDF and fitted CDF are shown in the inset.Right:Binned OHD and their reconstructions on H ðz Þ(top)and w ðz Þ(bottom).The shaded orange areas show the 1 error regions if we assume a ÃCDM model.Dashed blue lines are the results from binned OHD.4/public/teles-instr/e-elt.html .5http://www.iac.es/proyecto/codex/.HAO-RAN YU,SHUO PHYSICAL REVIEW D 88,103528(2013)from Ia SNe[11],it could not be verified by few data with low quality around1<z<2and a single datum(better quality)at z¼2:3.Here a single datum at z¼2:3[25], which is by BAO features in the redshift range2:1<z< 3:5,plays an important role in determining H and H0,which sensitively tune the resulting wðzÞ.Thus,further,more precise measurements at high redshift are needed.Our forecast of future data has proved to be able to verify or refute this deviation.The top left panel of Fig.2shows another Monte Carlo simulation,withÃCDM replaced by model1,and we can see that the reconstructed line by OHD is obviously excluded by this model.The remaining three panels in Fig.2show the results from simulated future OHD with the underlying models beingÃCDM,model1,and model2,respectively.Note that,for varying qualities of data and underlying model complexity,ðN;MÞare auto-matically adjusted,and are no longer suitable for current OHD,so we do not include current reconstructions.We can see that future OHD are able to reconstruct wðzÞvery accurately,even for very oscillatory models.V.DISCUSSION AND CONCLUSIONWe propose a new method by combining principal com-ponent analysis(PCA)and the goodness offit(GoF)crite-rion to reconstruct a dark energy equation of state wðzÞwith observational Hubble parameter data.We also used a new error model to simulate the error distributions of future OHD and get forecasted simulated data by estimating po-tential surveys and data acquisition methods.In the GoF criterion analysis,we calculate the residual between current 28OHD measurements and the concordanceÃCDM model andfind that the residual GoF? R?=n¼0:62(n¼28) is far from unity,which implies that these measurements should not be considered independent and thus it is im-proper to use a diagonal covariance matrix in the linear least squares.We bin the data by combining nearby redshift dataand get15newer measurements of OHD.They seem to beindependent—their GoF?¼1:02.We note that this is only a rough estimation,and generally we need to get the covari-ance between the measurements.However,there is noestimated covariance for OHD,so manually making themindependent is just for comparing the power of current andfuture data on the EoS reconstruction.We use our method to reconstruct EoS wðzÞby thecurrent OHD,and discover a feature of deviation fromÃCDM at z>1:5.The Ia SN data give similar results[11].While the quality of current data cannot verify it,future data from our simulations greatly enhance the result(Fig.2)and are able to confirm or deny this deviation.Inthe analysis we use only cosmic chronometer data.We alsouse simulations with the data from other sources(BAO,SLsignals),and the result does not improve much,because thequantity of three BAO data andfive SL data dominatesonly a small fraction of total100simulated data.However,with a small number of cosmic chronometer data(say,20),adding several SL signal data lowers the error of recon-struction at z’2,but beyond this redshift the result is stillnot confident.The reason is that adding a few high-redshiftdata helps to determine H0better at z’2but is still unableto determine H0well at2&z&5.In such a case we mayuse a derivative prior in Eq.(3)or use a Gaussian process[30]to help to determine H0.ACKNOWLEDGMENTSThis work was supported by the National ScienceFoundation of China(Grant No.11173006),the Ministryof Science and Technology National Basic Science pro-gram(project973)under Grant No.2012CB821804,and the Fundamental Research Funds for the CentralUniversities.[1] A.G.Riess et al.,Astron.J.116,1009(1998).[2]S.Perlmutter et al.,Astrophys.J.517,565(1999).[3]Y.Wang and M.Tegmark,Phys.Rev.Lett.92,241302(2004).[4]M.Chevallier and D.Polarski,Int.J.Mod.Phys.D10,213(2001).[5] E.V.Linder,Phys.Rev.Lett.90,091301(2003).[6]V.Sahni,A.Shafieloo,and A.A.Starobinsky,Phys.Rev.D78,103502(2008).[7]G.-B.Zhao,R.G.Crittenden,L.Pogosian,and X.Zhang,Phys.Rev.Lett.109,171301(2012).[8]R.A.Daly and S.G.Djorgovski,Astrophys.J.612,652(2004).[9]R.A.Daly,S.G.Djorgovski,K.A.Freeman,M.P.Mory,C.P.O’Dea,P.Kharb,and S.Baum,Astrophys.J.677,1(2008).[10] D.Huterer and G.Starkman,Phys.Rev.Lett.90,031301(2003).[11] C.Clarkson and C.Zunckel,Phys.Rev.Lett.104,211301(2010).[12]Z.-L.Yi and T.-J.Zhang,Mod.Phys.Lett.A22,41(2007).[13]R.Jimenez and A.Loeb,Astrophys.J.573,37(2002).[14] E.Gaztan˜aga,A.Cabre´,and L.Hui,Mon.Not.R.Astron.Soc.399,1663(2009).[15] A.Sandage,Astrophys.J.136,319(1962).[16] A.Loeb,Astrophys.J.Lett.499,L111(1998).[17] C.Ma and T.-J.Zhang,Astrophys.J.730,74(2011).[18] A.R.Liddle,Mon.Not.R.Astron.Soc.377,L74(2007).[19]J.Sheil and I.O’Muircheartaigh,J.Roy.Stat.Soc.C Appl.Stat.26,92(1977).NONPARAMETRIC RECONSTRUCTION OF DYNAMICAL...PHYSICAL REVIEW D88,103528(2013)[20]J.Simon,L.Verde,and R.Jimenez,Phys.Rev.D71,123001(2005).[21] D.Stern,R.Jimenez,L.Verde,S.A.Stanford,and M.Kamionkowski,Astrophys.J.Suppl.Ser.188,280(2010).[22]M.Moresco et al.,J.Cosmol.Astropart.Phys.8(2012)006.[23] C.Zhang,H.Zhang,S.Yuan,T.-J.Zhang,and Y.-C.Sun,arXiv:1207.4541.[24] C.Blake et al.,Mon.Not.R.Astron.Soc.425,405(2012).[25]N.G.Busca et al.,Astron.Astrophys.552,A96(2013).[26]P.-S.Corasaniti,D.Huterer,and A.Melchiorri,Phys.Rev.D75,062001(2007).[27]J.Liske et al.,Mon.Not.R.Astron.Soc.386,1192(2008).[28]M.Martinelli,S.Pandolfi,C.J.A.P.Martins,and P.E.Vielzeuf,Phys.Rev.D86,123001(2012).[29]Planck Collaboration,arXiv:1303.5076.[30]M.Seikel,S.Yahya,R.Maartens,and C.Clarkson,Phys.Rev.D86,083001(2012).HAO-RAN YU,SHUO YUAN,AND TONG-JIE ZHANG PHYSICAL REVIEW D88,103528(2013)。