{'fraction_non_alphanumeric': 0.031053945267702547, 'fraction_numerical': 0.012997234152593825, 'mean_word_length': 4.334932821497121, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 8, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Toward unbiased estimations of the statefinder parameters

With the use of simulated supernova catalogs, we show that the statefinder parameters turn out to be poorly and biased estimated by standard cosmography. To this end, we compute their standard deviations and several bias statistics on cosmologies near the concordance model, demonstrating that these are very large, making standard cosmography unsuitable for future and wider compilations of data. To overcome this issue, we propose a new method that consists in introducing the series of the Hubble function into the luminosity distance, instead of considering the usual direct Taylor expansions of the luminosity distance. Moreover, in order to speed up the numerical computations, we estimate the coefficients of our expansions in a hierarchical manner, in which the order of the expansion depends on the redshift of every single piece of data. In addition, we propose two hybrids methods that incorporates standard cosmography at low redshifts. The methods presented here perform better than the standard approach of cosmography both in the errors and bias of the estimated statefinders. We further propose a one-parameter diagnostic to reject non-viable methods in cosmography.

With the use of simulated supernova catalogs, we show that the statefinder parameters turn out to be poorly and biased estimated by standard cosmography. To this end, we compute their standard deviations and several bias statistics on cosmologies near the concordance model, demonstrating that these are very large, making standard cosmography unsuitable for future and wider compilations of data. To overcome this issue, we propose a new method that consists in introducing the series of the Hubble function into the luminosity distance, instead of considering the usual direct Taylor expansions of the luminosity distance. Moreover, in order to speed up the numerical computations, we estimate the coefficients of our expansions in a hierarchical manner, in which the order of the expansion depends on the redshift of every single piece of data. In addition, we propose two hybrids methods that incorporates standard cosmography at low redshifts. The methods presented here perform better than the standard approach of cosmography both in the errors and bias of the estimated statefinders. We further propose a one-parameter diagnostic to reject non-viable methods in cosmography.

I. INTRODUCTION
The universe is currently undergoing a late-time speeding up phase [1,2]. The component responsible for the acceleration is commonly named dark energy and manifests a negative equation of state providing gravitational repulsive effects. Nonetheless, a totally satisfactory fundamental description does not exist and dark energy still continues challenging the concordance paradigm [3,4]. In addition, to characterize the dark energy evolution, one needs to postulate a model a priori. This leads to numerical constraints on the free coefficients of the model which are determined in a model dependent way. Hence, treatments which encapsulate different aspects of cosmology without calling any specific model become useful to understand whether dark energy evolves or not in time.
In the ΛCDM model, the accelerated expansion of the Universe is driven by a cosmological constant, while the clustering of matter at large scales is a consequence of the gravitational self-attraction of a stress-free and dustlike collection of particles dubbed cold dark matter. Although the ΛCDM model suffers from the cosmological constant and cosmic coincidence problems [5,6], the cosmological constant is physically well-motivated as the minimal modification to General Relativity consistent with general covariance. Perhaps more intriguing is the nature of the dark matter sector: for example, today observations are not able to tell if it is really cold, or it is warm.
The success of the ΛCDM model is much more appreciable in the early universe (see, e.g., figures 1 and 3 in [7]), but it is still not very well tested at late times, where there is large room for evolving dark energy and several behaviors for the dark matter. 1 Hence, approaches beyond the ΛCDM model are widely investigated by the community.
Among general model independent treatments, cosmography-on the background (hereafter cosmography) attempts to reconstruct the expansion history of the universe in a model independent way; see [13] for a recent review. To do so, it usually takes into account a set of parameters, derivatives of the scale factor, related to the statefinders [14,15] and generically named cosmo-graphic series. 2 Such a procedure enables us to consider the fewest number of assumptions as possible. Indeed, for the ΛCDM model at the background level, one needs three quantities only to address the late evolution: the Hubble rate H, the deceleration parameter q and its variation, related to the jerk parameter j. Furthermore, it is usually assumed a flat universe, which reduce the parameters to H and q. The need to go beyond these quantities has increased recently in order to constrain additional parameters with higher accuracy. This has led cosmologists to propose new diagnostics to investigate the corresponding phase spaces of coefficients [14][15][16].
The case of fitting theoretical models by using the statefinders has always attracted attention. More recently, it has been argued that the impossibility for the statefinders to constrain general models poses questions about its usefulness [28][29][30]. In this regard, our point of view is mostly different. Let us suppose to have a theoretical model M(Ψ; α i ) with Ψ denoting the fields of the model and α i its free parameters. For M to be well defined, each realization of the parameters, subjected to initial conditions, should give a unique Hubble diagram 3 ; and given that, it is as simple as taking derivatives to find the statefinders in that model. Thereafter, one can make comparisons to the measured statefinders, and accept or reject the realization α i ; clearly, one can restart the procedure with a new realization of the parameters, although at this point it could be a better idea to fit directly to the data. That is, to fit a general class of models, say e.g. the whole class of f (R) theories, using the statefinders is not always possible at least some extra assumptions are considered [31].
Although doable, the main objective of cosmography is not to constrain theoretical models -the statefinders are not data. Instead, it is to reconstruct the Hubble diagram as model independent as possible. Cosmography is a "top-down" approach to cosmology [25], that attempts to deduce its kinematics directly from observations; contrary to a "bottom-up" approach, that assumes the dynamics of a given model from the very beginning.
In the class of cosmographic methods lying on series expansions, the most spinous difficulty remains likely the 2 In our paper we refer to the cosmographic series with the name statefinders, which slightly differs from the introduced first in [14]. For details see Appendix A. 3 It is not always possible to obtain the Hubble function for arbitrary free parameters. Particularly in higher derivatives theories.
convergence problem [21]. It is intimately related to truncating series choosing the particular Taylor expansion to use. Indeed, all finite Taylor series diverge when z 1 leading to possible misleading outcomes. Some authors have attempted to overcome this obstacle by using auxiliary variables built up in terms of the redshift z [21,22,32].
Another major problem of the available cosmographical methods is that the estimation of the statefinder parameters is in general biased. For example, the authors of Ref. [28] conclude that estimations in cosmography are biased when expanding the luminosity distance as powers of the function of redshift y = z/(1 + z).
SC also suffers from this bias problem. To neatly observe it, we consider the following simple example. In this approach, the luminosity distance d L is given bỹ whered L ≡ H 0 d L /c with c = 1 the speed of light and a subindex "0" means that the statefinders are evaluated at z = 0. We fit the first two cosmographic parameters in SC to our "exact" simulated catalog with Ω m = 0.3 and w = −1 (see Sec. IV C below). From this fit we can derive the parameters of a wCDM model. In Fig. 1 we show a plot of the joint 2-dimensional posterior region Ω m -w. It can be observed that the true cosmology is out of the estimations of SC at 2σ, showing the large bias in SC.
Another issue in cosmography is due to correlations. In general, the statefinders are degenerated, for example in the flat ΛCDM model there is the linear degeneracy s 0, ΛCDM = −2 − 3q 0, ΛCDM . Since the true cosmological model should be close to ΛCDM, we expect that a successful cosmography method should follow this trend in a small neighborhood of its best fit. 4 In Sec. IV, we show that this is not the case for SC. We assume the reason behind it is that the third-degree polynomial form of the z 4 coefficient in Eq. (1) leads to fictitious degeneracies.
In this work, we extensively use simulated supernova catalogs to address the bias problem and to study the posterior distributions in SC. We show that the posteriors not only have wide dispersions, but also a large bias, non-tolerable for future observations (for example, [cite]). To overcome these problems, and motivated for the unwanted degeneracies in SC exposed above, we propose a new route that lies on expanding the Hubble function Here we fit only the first two cosmographic parameters to our simulated exact data built upon a wCDM model with Ωm = 0.3 and w = −1. We note that standard cosmography fails to obtain the true cosmology at 2σ. and use that expansion as an input for observables as the modulus distance. By using several tests, we indeed show that the dispersions and bias turn out to be smaller. We consider our approach by means of a hierarchy in which we split the data into redshift domains, speeding up the numerical computations. Due to the success of SC at low redshifts we further propose two hybrid methods. The paper is organized as follows: in Sec. II we highlight the treatment of cosmography, whereas in Sec. III, we present our proposed cosmographic approach. In Sec. IV, we go toward it and we test our model and SC with simulated data, while in Sec. V we estimate the statefinders parameters using the JLA and Union2.1 compilations. Finally, in Sec. VI we present our conclusions and perspectives.

II. COSMOGRAPHY OF THE UNIVERSE
Cosmography, or better the cosmographic method, stands for a model-independent treatment able to fix limits on universe's kinematics, without imposing a cosmological model a priori [21,22]. Particularly, cosmography makes possible to gather the universe expansion history at small redshift regimes by simply invoking the cosmological principle, without any additional requirements on Einstein's energy momentum tensor. In other words, postulating the homogeneity and isotropy, with spatial curvature somehow fixed, cosmography can account for the evolving dark energy contributions in Einstein's equations. To do so, one can understand if dark energy is composed by a pure cosmological constant or by some particular exotic fluid.
The strategy is based on expanding all cosmological observables of interest around present time. Even though this procedure is the most used, it is even possible to expand only the scale factor a(t) and then to rewrite all quantities in terms of it. In any cases, each expansions may be matched with cosmic data in order to give bounds on the evolution of each variable under exam. Hence, one gets numerical outcomes which do not depend on any particular requirements since only Taylor expansions are compared with data. Further, relating observations to theoretical predictions means to heal the degeneracy problem among cosmological models. Indeed, cosmography is therefore able to distinguish different classes of models which turn out to be compatible with cosmographic predictions and those ones that have to be discarded.
Bearing this in mind leads to consider cosmography as a powerful method that allows to study present-time cosmology and to describe the dynamical evolution of the universe. Following standard recipes lies on postulating a scale factor a ≡ 1/(1 + z) expanded today as: The cosmographic parameters are the terms entering the expansion (2) 5 , which represent model-independent quantities related to the Taylor expansion of Hubble's rate by: This gives: where we used the convention in which the subscripts represent the derivative orders with respect to the redshift z. In analogy, we can explicitly write Eq. (2) to have: where we have normalized the scale factor to a(t 0 ) = 1. 5 They have been explicitly reported in Appendix A.
Each term displays a specific meaning associated to dynamical properties of the universe. For example, q and j fix kinematic properties at lower redshift domains, since the value of q at a given time specifies whether the universe is accelerating or decelerating. Further j is intimately related to its variation in the past. In particular, 1. q 0 > 0 shows that the universe is expanding universe although it undergoes a deceleration phase, as for example in the case of pressureless matter dominated phase.
2. −1 < q 0 < 0 represents an expanding universe which speeds up as expected by current observations.
3. q 0 = −1 indicates that all the whole cosmological energy budget is dominated by a de Sitter fluid.
4. j 0 < 0 means that dark energy influences early time dynamics in the same way of late time evolution.
5. j 0 = 0 indicates that the acceleration parameter smoothly tends to a precise value, without changing its behavior as z → ∞.
6. j 0 > 0 implies that the universe acceleration started at a precise time during the evolution, associated to the transition redshift. In such a way, it provides the acceleration changes sign during time.
All those properties are valid since the variation of q and j are intertwined by the following relation: Today, we have j 0 = dq dz | 0 + 2q 2 0 + q 0 . Since we measure that −1 ≤ q 0 < −1/2, we clearly obtain: 2q 2 0 +q 0 > 0 and so, when q 0 < −1/2 the term j 0 is linked to the sign of the variation of q. Although highly predictive, cosmography is jeopardized by some drawbacks which limit its use. In particular: Degeneracy between coefficients: The whole list of independent parameters is q 0 , j 0 , s 0 , . . ., but measuring H 0 leads to degeneracy since it is unfortunately impossible to estimate H 0 alone by using measurements of distance modulus. In other words, H 0 degenerates with the rest of the parameters.
Degeneracy with spatial curvature: Spatial curvature must be fixed somehow, otherwise j and s parameters may be strongly influenced by its value, since spatial curvature degenerates with them. However, small deviations do not influence the simplest case Ω k = 0 and we extensively use the former case.
Systematics due to truncated series and convergence: Slower convergence in the best fit algorithm may be induced by choosing truncated series at a precise order, while systematics in measurements occur, on the contrary, if series are expanded up to a certain order. On the other side, the problem of truncated series lies also on determining the particular Taylor expansion which is better to use with precise data survey. Taylor series, however, are always evaluated around present time 6 , or better defined when z = 0. So all data sets exceed the bound z 0 giving, in principle, that all Taylor series do not converge when z 1. Combining different data would alleviate the systematics which is produced with the aforementioned theoretical problems, but cannot be considered exhaustive to handle high redshift data sets.
While degeneracies can be alleviated by refined analyses and combined techniques of cosmographic reconstructions, systematics and convergence are difficult to treat. In particular, to alleviate the convergence problem one can employ parameterizations of the redshift z, by means of auxiliary variables (Z new ). This technique seems to enlarge the convergence radius of Taylor's expansions to a wider sphere having radius Z new < 1 [21]. In such a picture, data lying within z ∈ [0, ∞) are rewritten into shorter (non-divergent) ranges. For example, . Another technique, more recent and likely suitable, is to consider rational approximations, such as the Padé approximants, in which the expansions are taken by rational functions which do not diverge as the redshift increases [18]. In order to reduce and alleviate systematics over cosmographic measurements, we show below our new scenarios which refine cosmographic results. We demonstrate that using our method the numerics can be clearly improved in a concise and suitable way.

III. NEW STRATEGIES TOWARD COSMOGRAPHY WITH HUBBLE EXPANSIONS
We start by expanding the Hubble function in Taylor series about redshift z = 0, with The first four coefficients, here named the eis, can be written in terms of the statefinders parameters (q 0 , j 0 , s 0 ) as The comoving distance is defined as the (comoving) distance a photon travels from a source at a redshift z to us, at z = 0, In our new approach to cosmography we use directly the Taylor expansion of H(z) in the comoving distance expression and integrate numerically to obtain the luminosity distance, that is From Eq. (7) we can note that at low redshifts not all of the powers in the E(z) expansion are important. Thus, in order to speed up the numerical computations, we estimate the eis parameters in a hierarchical manner. To this end we use Eq. (10) in redshift bins as For z > z high we expand in Taylor series the integrand ofd L (z k ) up to z 3 and analytically perform the integration. This last step is necessary for numerical stability, otherwise the tails of the posterior distributions become very noisy.
That is, for a supernova at redshift z k in a given simulated catalog, we used We performed preliminar numerics and found that a good choice for these redshift cuts is The numerical outcomes of this particular binning does not differ significatively from those obtained by a direct application of Eq. (10). We name the method of Eq. (11), the Eis method, or simply Eis. We further propose two hybrid methods: hybrid 1 which consists in the use of SC up to z mid , and beyond that redshift in using the Eis method with the integrand expanded. The other method is hybrid 2, which is a modification of hybrid 1 with z high = z mid = 0.4. The physical reasons behind the two hybrid approaches are essentially based on the fact that at low redshifts, SC constrains with small dispersions the deceleration and jerk parameters. The here involved equations for the hybrid methods are given in Appendix A.
In the following sections we use the code CosmoMC [33] to draw the likelihood distributions for all the methods studied here with a Metropolis-Hasting MCMC algorithm [34]. A module for CosmoMC is available at https://github.com/alejandroaviles/EisCosmography; also, all the simulated catalogs and further statistics can be found there.
Despite our numerics perform fits of the eis parameters (even for SC; see Eq. (A4)), we are primarily interested in the statefinders because its physical interpretation is more familiar. We are imposing flat priors on the eis parameters, which is not equivalent to flat priors on the statefinders. Nevertheless, this does not have a significative impact, as can be observed by comparing the triangular plot in Fig. 3 (bottom panel) to previous works on cosmography; see, e.g. [22,28].

IV. NUMERICAL ANALYSIS WITH SIMULATED DATA
In this section we test the performance of our Eis method by using simulated catalogs of supernovae Ia. To this end, we construct two kinds of simulated data based on fiducial ΛCDM and wCDM models. On each of these simulations we take 740 data distributed with the same redshifts and error bars of the observed peak magnitudes (∼ 0.12) as those of the JLA compilation [35]. We choose this catalog because it has a large amount of low redshift data providing a good inference of E 1 that acts as a leverage for a better estimation of the rest of eis parameters. We better discuss this point in Sect. (IV D).

A. The dispersed simulated data
The dispersed simulated data set consists in 100 simulations in which each supernovae in every catalog is obtained from fiducial ΛCDM models with Ω m drawn from a Gaussian distribution with mean 0.30 and standard de- viation of 0.034. 7 That is, in a single catalog, for each supernova at redshift z k we give the modulus distance a value µ k = µ(z k ; Ω m k ∈ N (0.30, 0.034)). Taking the mean value and the standard deviation of this Ω m distribution and using ΛCDM-statefinder relations, we expect that our results intersect the intervals q 0 = −0.55±0.051, j 0 = 1 and s 0 = −0.35 ± 0.153. This is exactly the case that we got for all the models involved into our analyses. One may be tempted to assume that the underlying, true, cosmology of the dispersed simulated data sets is the same for all of them, and given by Ω m = 0.3; see, e.g. [28]. The drawback of this approach is that the true cosmology is actually unknown for each catalog. Thus, in order to analyze how good the fits are, we must compare against ΛCDM fittings to the same simulations.
In Fig. 2 we show a histogram of the estimatedΩ m 's by performing these fits. The average value is Ω m = 0.307, the dispersion is σΩ m = 0.036, and the average of the standard deviations is σ Ωm = 0.012.
In Fig. 3 we show a triangular plot for the derived statefinders of one of our simulated catalogs, named SD 1, both for the SC and the Eis methods, revealing that the dispersions are smaller than in SC, most notably for the parameter E 3 (or s 0 ). For the E 1 parameter there is also an improvement over SC, nevertheless this is not as remarkable as for the other parameters, the reason behind this is that q 0 is constrained mainly by low redshift supernovae, as it is shown in Sec.IV D below, particularly in Fig. 8, and all the expansions considered in this work converge at low redshift. In Fig. 4 we display a zoom of the q 0 -s 0 2D joint posterior. The shadows show the confidence intervals (c.i.) for q 0 and s 0 derived from fitting the ΛCDM model. The solid black line is s 0 = −2 − 3q 0 , which corresponds to the allowed region in ΛCDM [this can be derived by setting w = −1 in Eqs. (A5)]. We notice that SC does not follow this degeneracy trend, while Eis can do it inside its 0.68 confidence region. 7 We fix the Hubble constant to H 0 = 70 km/s/Mpc. Nevertheless, it is irrelevant because supernovae data cannot estimate the Hubble constant. Therefore, in the numerical analysis we internally marginalize the combination 5 ln(c/H 0 ) + M b as described in [36].
Complementing Fig. 3, in Table I we show the 1dimensional marginalized posterior intervals for the simulated data SD 1. The average statistics of the standard deviations and mean posterior values for all the methods are shown in Table II.
In Fig. 5 we show the triangular plot derived from the fitting to the simulated catalog SD 1. We make note that the hybrid 1 method provides the smallest dispersions among the four methods studied here.
From the results of this subsection, we conclude that the approaches proposed in this work improve the standard deviations of the statefinders estimations from those obtained by SC.

B. Bias on the estimators
We use the bias of an estimatorθ defined as b θ = bias(θ) ≡θ − θ true (13) where θ true is the true value of the parameter θ. For the estimatedθ we use the mean value of the posterior distribution. We may use the maximum likelihood estimator instead, but our preliminary numerics show only slight differences despite the distributions are in general skewed.
As explained above, we do not know the true values of the parameters θ; thus, we will assume that the ΛCDM provides unbiased estimations for them. That is, we use θ true =θ ΛCDM .
By itself, the bias does not provide a complete information of how wellθ estimates θ, specially if we know only approximately the true value. For this reason, it is convenient to use complementary statistics. First, we implement the risk statistics [37,38] which penalize the bias with the standard deviation. Furthermore, for the whole cosmology, following [37,39] (see also [40,41] for applications in cosmology), we compute the bias statistics which roughly quantifies the slip from the χ 2 -statistics due to bias. Here F is the reduced Fisher matrix for the estimated parameters and b = (b E1 , b E2 , b E3 ) is the bias vector for eis parameters. 8 We note that this statistics requires the maximum likelihood estimator instead of the posterior mean. A smaller ∆χ 2 does not imply a smaller bias, this can be noted for the case of one single parameter, since ∆χ 2 = b 2 θ /σ 2 θ for this case. Therefore, we additionally compute the figure of merit (FoM), that we define as 8 ∆χ 2 is not an invariant between statefinders (sf ) and eis parameters because they are not linearly related and therefore sf = sf ( eis). Nevertheless, the differences in values are small. On the other hand, the determinant of the Fisher matrix, and hence the FoM, is an invariant because the determinant of the transformation matrix between eis and sf parameters is equal to −1.
The numerical factor 4π/3 is not common in the literature. We consider it so that the FoM coincides with the volume of the 3-dimensional ellipsoid defined by the covariance matrix.
Strictly, the FoM and ∆χ 2 bias statistics work only for multi-variate Gaussian distributions. By using them, we are implicitly approximating the confidence regions by ellipsoids. In Fig. 6 we show the projection for two of these on the q 0 -s 0 subspace. For comparison we also show the Markov chains obtained in the MCMC analysis, as well as the region allowed in the ΛCDM model.
We perform the four statistics for each simulated catalog and for each one of our models.
In Table II we show the average values of the bias and risk for both the eis and statefinder parameters. It can be noted that in these 1-parameter bias tests, Eis, as well as hybrid 1 and hybrid 2, perform better than SC.
For the whole 3-dimensional bias statistics, the average values over the 100 simulated catalogs for ∆χ 2 and FoM are: For ν = 3 parameters, the 1σ contour is ∆χ 2 ≤ 3.53 [39], thus for all the methods, on the average, the true value is well inside this region. We note SC is able to do it because the volume of its error ellipsoid, or FoM,  is very large in comparison with the other methods; cf. Fig. 6. We have seen in Fig. 1 that, when using 2 parameters, the true value lies outside the SC 2σ confidence region. This is because when the number of parameters is reduced, the bias is propagated from higher to lower order cosmographic parameters. In the upcoming section we will observe this effect more clearly by reducing the number of parameters to just one. On the overall analysis it seems the three new methods presented in this work have a similar performance. Thus, by virtue of the simplicity of Eis and the more appealing form of its contour plots, in the following we focus mainly to this method. tribute the exact value µ k = µ(z k ; Ω m , w) to the modulus distance. The aim of constructing unphysical exact simulated data is to compare the performance of our new method with the SC approach. Indeed, knowing the exact cosmology permits to compare the bias of the estimators for the statefinders in both approaches with higher precision.
Our results indicate that the width of the posterior distributions are similar to those of dispersed simulated data, with Eis providing smaller standard deviations than SC. Keeping this in mind, here we focus our attention to the statefinders' bias.
In Table III we report the bias in j 0 and s 0 , we do not show the bias on q 0 , which is similar for both methods and less than 2%. In Table IV we show the bias statistics ∆χ 2 , noting again that the values for both methods are comparable, although Eis has a smaller FoM.
We note from Table IV that there is a clear tendency to have larger bias for larger matter abundances for the Eis method. This could have been foreseen due to the fact that modulus distance curves µ(z; Ω m ) are more densely distributed for larger values of Ω m , and consequently the cosmographic parameters are more sensitive to small steps on larger matter abundances. The tendency is not clear for the case of SC.
To end up this section, with the Ω m = 0.3, w = −1 exact simulated catalog, we perform a test that lets us observe the consequences of the bias in a simple way and by one single parameter. The test consists in fixing the eis to the ΛCDM model. Hence we set j 0 = 1 and s 0 = −2 − 3q 0 , that can be translated to the eis parameters as E 2 = −E 2 1 + 2E 1 and E 3 = 3E 3 1 − 6E 2 1 + 2E 1 . We impose these conditions internally in the code, and let E 1 to be the only variable. In Fig. 7, we show the posteriors of Ω m for this test. The obtained means and standard deviations are The dispersions are similar for the five models. Albeit the best fit is quite different for SC. This is a consequence of the large bias of j 0 and s 0 that now has been absorbed by q 0 . Given that we are fitting to exact simulated data, this bias is intrinsic to the cosmographic method; a reduction of the error bars or, equivalently, adding more data over the same redshift domain, reduces the standard deviations of the estimations, but it has a small impact on the bias. Therefore, this one-parameter diagnostic provides a robust criterion to discard non-viable cosmography approaches. The methods proposed in this work are not free of this bias propagation effect, although it is much smaller.

D. Redshift distributions
We now explore different redshift distributions. This analysis shows the importance of having a large amount of low redshift supernovae to reduce the bias. We split the redshift range covered by the JLA catalog (z ∈ (0.01, 1.3)) in four bins divided by the redshift cuts z low , z mid z high as in Eq. (12). We explore three different distributions: zdist 1 has the same redshifts as the JLA compilation, zdist 2 is skewed to low z, and zdist 3 is skewed to the high z. The number of supernovae per bin is denoted by (N 1 , N 2 , N 3 , N 4 ) with N 1 the number of supernovae in the bin 0 < z < z low , N 2 the number of supernovae in the bin z low < z < z mid , N 3 the number of supernovae in the bin z mid < z < z high , and N 4 the number of supernovae in the bin z > z high . For zdist 2  With these distributions we construct exact simulated data based on a ΛCDM model with Ω m = 0.3 and errors σ µ = 0.118, which is the average of errors in the JLA catalog.
In Table V we show the marginalized 1-dimensional estimations, along with the 3-dimensional bias statistics ∆χ 2 and FoM. It is expected that the standard deviations of the higher statefinders are smaller when we add more supernovae at high redshifts, because the influence of variations in Ω m on d L is more evident at higher redshifts. Nevertheless it is clear that the bias is enhanced for the zdist 3 case, having it higher FoM and ∆χ 2 . Contrary, when considering larger amounts of low redshifts supernovae, the bias is considerably reduced, this is the case of zdist 1.
This is more evident if we use the one-parameter diagnostic of Sec. IV C. In Fig. 8 we show the plots for Ω m obtained by performing this test and in Table V we show the best fits and 0.68 c.i.. It is clear that the inclusion of a large amount of low redshift supernovae in zdist 1 reduces the bias, while for zdist 3 it becomes quite large, such that the true value is not attainable. Instead, the standard deviation is reduced as explained above.  In this section we fit the Eis and SC methods to real data. To this end we consider the two most used supernova type Ia compilations: the JLA [35] and the Union2.1 [42].
The major complication on real supernova catalogs is the presence of systematic errors, mainly due to photometric calibration and selection bias, but also to physical effects as photon absortion by intervening dust and gravitational lensing. Supernova systematics specially limits the fitting procedure when considering several parameters, as it is our case. Since these errors do not follow any specific physical model (e.g. ΛCDM) we expect departures on the bias obtained in the previous sections, specially for the parameter j 0 , which serves as a nulldiagnostic for flat ΛCDM -j 0 will not be necessarily underestimated (overestimated) for the Eis (SC) method, as it is the case for the simulated data. To partially alleviate these problems we could additionally use other types of data, as BAO or direct observations of the Hubble flow. Although straightforward, we delegate this endeavor for future investigations, since in this work we have only analyzed the response of the statefinders to supernovae data.
The numerical fit to JLA is quite slow due to the incorporation of nuisance parameters that are present also on the covariance error matrix. Thus, for every step on the Markov chains one has to invert the full covariance matrix. In the presence of nuisance parameters, the modulus distance to the k-th supernova on the compilation follows the linear model µ k = m B, k − (M B, k − α JLA X 1, k + β JLA C k ) where m B is the observed peak magnitude, M B the absolute magnitude, X 1 the time stretching of the light curve, and C the color at maximum brightness. In principle the absolute distance is the sum of two nuisance parameters M B = M 1 B + ∆ M , but since the errors do not depend on them we can internally marginalize, as it is already incorporated in the JLA likelihood module to CosmoMC. Therefore, we only estimate α JLA and β JLA nuisance parameters.
We were unable to sample adequately the tails of the posterior distributions for E 3 (or equivalently the snap) in the SC case using the JLA compilation. For that reason we warn the reader to take the JLA results for SC in Table VI with precaution. We also report the χ 2 statistics maximum values; the number of degrees of freedom (d.o.f.) for this compilation considering three free parameters and two nuisances is d.o.f. = 740 − 3 − 2 = 735.
The nuisance parameters for both SC and Eis have almost the same 0.68 c.i.: α JLA = 0.14 ± 0.01, β JLA = 3.11 ± 0.08, differing only in the third significative figure.
We perform a second fit to the Eis method, but with the nuisance fixed to their best fit values. We name this fit as Eis*.
The Union2.1 compilation is lacking of nuisance parameters, thus the fitting procedure is straightforward. Considering systematic errors the statefinders constraints at 0.68 c.i. are shown in Table VI. Additionally, we report the χ 2 statistics maximum value, which for this Again, for the case of SC we obtained non-conclusive posterior distributions for E 3 . We conclude that the available supernova data compilations alone are not able to fit SC when using more than two statefinder parameters.
In Table VI the central values are the means of the marginalized posterior distributions, and the errors denote the departures from the means of the lower and higher limits of the 0.68 c.i. For the Eis method fitting to the Union 2.1 compilation, due to the skewness of s 0 , its mean value is located out of its 0.68 c.i

VI. CONCLUSIONS
In this work, by accounting for Hubble expansions inside the luminosity distance, we tailored a method that leads to less biased estimations and smaller standard deviations than those in SC. We baptized our new method Eis, since it estimates the coefficients (named eis) of Taylor's expansion of the normalized Hubble rate E(z) = H(z)/H 0 . We focused on the first three eis parameters by using data spanning over the redshift interval 0 < z ≤ 1.2. From them, the deceleration q 0 , jerk j 0 , and snap s 0 statefinders parameters can be derived.
In order to speed up the computations, we further split all the numerical analyses in redshift bins. In so doing, the order of Hubble's expansion depends on the redshift of a single supernova data. We have checked that this binning procedure does not alter significatively the parameters' estimations. This speeding up is in average a factor of 2.5, and it is mainly due to the convergence of MCMC chains, which is attained within a less number of steps.
We further constructed two hybrid models that made use of SC for redshifts z < z mid = 0.4, whereas over the rest of the redshift domain they employed expansions of our Eis method.
By extensively considering simulated catalogs of supernovae built up through the JLA compilation, we also compute several bias statistics for cosmological models near the concordance model. Specifically, for each single parameter, we compute the bias and the risk statistics; and, in order to account for the correlations of the statefinders, we further use the ∆χ 2 bias statistics and the FoM.
We concluded that our methods provide less biased estimations than in the SC case. Moreover, the standard deviations of the posterior distributions are considerably smaller than in SC.
The other issue that we faced out was the one due to degeneracies. We showed that SC is not able to follow the ΛCDM degeneracies on the statefinders, even when the fitting is performed against an "exact" ΛCDM model simulated catalog. Meanwhile our method Eis is capable to do it. We assumed that the third-degree polynomial form of the z 4 coefficient in the SC luminosity distance is responsible for these fictitious degeneracies. Actually, this intuition motivated us in the first place to construct more viable methods in cosmography.
We further proposed a new test to reject non-viable methods in cosmography. It consists in building up an exact simulated catalog of supernovae by means of a given fiducial ΛCDM cosmology. Thereafter, all except the E 1 (or equivalently q 0 ) of the cosmographic parameters should be fixed into the code by using the degeneracies of the ΛCDM model. Finally, the estimations are performed. This is a simple 1-parameter diagnostic that addresses the propagation of bias from higher orders statefinders to the deceleration parameter. We showed how SC is not able to pass this test since the estimation of Ω m for it turns out to be more than 2 σ away from the true value, reflecting its highly biased estimations. We highlighted that the methods presented here are not completely free of this bias propagation, albeit the effects are much smaller than in SC.
Finally, we applied our method to real data, which provides the most stringent to date constraints of the statefinders by using supernovae data only; for either JLA or Union2.1 compilations.
We make publicly-available a module to the code CosmoMC that perform the MCMC numerical analysis for the cosmographic methods of this work at https://github.com/alejandroaviles/EisCosmography. There, we also uploaded all the simulated data as well as further tables and statistical files.
Future works will focus on investigating the same procedures, discussed in our work, for investigating bias and dispersions of cosmographic estimations for other kinds of surveys such as BAO and/or H(z) data. For H(z) data, in particular, we expect the results to be the same for both E is and SC methods. This consideration comes from the definitions of both the methods. In fact, when assuming H(z) data, SC needs the direct expansions of the Hubble rate. This fact turns out to give analogous outcomes with respect to our approach.
The luminosity distance for hybrid 2 method is the same as that for hybrid 1 but with z high = z mid = 0.4.
Other useful formulae are the relations between flat wCDM and the cosmographic parameters. These are given by and Appendix B: Binning Analysis in the Eis method Throughout this work we have assumed that a suitable choice for the redshift cuts is given by Eq. (12). We now check the goodness of this approach. In particular, we do not perform the analysis for the hybrid methods because the Eis paradigm works better beyond the 1-parameter test and turns out to be more appropriate.
First, we note that z low and z mid cuts are physically different than z high . Indeed, the three aforementioned bins make use of Eqs. (11) for the estimation of the eis parameters, whereas for z > z high we require the integrand Taylor expansion.
We highly emphasize that the choice of z low and z mid is tailored only for speeding up the numerical outcomes, and due to this fact we decide to test three models: the first with the same binning already discussed throughout the work, the second model without bounds over z low and z mid , implying the use ofd L (z) =d (3) L (z) for z < z high while finally the third model with nuisance parameters for the redshift cuts with uniform priors over the intervals 0 < z low,N < 0.2 and 0.3 < z mid,N < 0.7 (black curves). The three methods are compared with the ΛCDM model with Ω m = 0.30 by using an exact simulated catalog; we even employ the redshift distribution found in the JLA catalog, with errors σ µ = 0.118, as in zdist 1 of Sec. IV D.
We perform MCMC estimations and we show the 1-dimensional marginalized posterior distributions in Fig. 9. We note no differences among the estimations of each methods. Moreover, the nuisances do not show a preferred value. From this analysis we conclude that the cuts z low and z mid are not very relevant for the estimations.
However, the situation is quite different for the fourth bin, defined by z > z high . When it is not considered, the best fit does not differ significatively, but the tails of the posterior distributions become very noisy, specially for E 3 (or equivalently s 0 ). This behavior has been portrayed in the q 0 -s 0 contour plot of Fig. 10 and are traced back to the samples of E 3 and E 2 that can make the denominator ofd L in Eq. (11) very close to zero or even negative.
One may wonder that if considering only SNe redshifts which are smaller than the value of z high turns out to be more convenient for better estimations. This does not seem to be the case if considering three parameters (although considering only two, instead of three, becomes a viable strategy). In Fig. 10 we also plot the former behavior, revealing that high redshift SNe are necessary for the estimation of the third statefinder.
Finally, we test methods in which there is no z low and z mid cuts, and we choose z high to take different values from z = 0.1 to z = 1.3, comprising 120 different z high redshifts separated by intervals of size ∆z = 0.1. We perform the 1-dimensional parameter test of Sec. IV C against the same simulated catalog of Figs. 9 and 10. In Fig. 11 we show the 1-dimensional mean (solid line) and 0.68 (blue) and 0.95 (light blue) c.i. of the Ω m distributions for each redshift z high ; for comparison we also show the true Ω m value (Dashed line). Below z high ≈ 0.65 the plots show a constant trend both in the best fit around Ω m ≈ 0.312 and in the standard deviation σ Ωm ≈ 0.015, FIG. 9. Here we focus on the impacts of z low and z mid on our estimations. Red curves correspond to the choice of cuts used throughout this work. Blue curves are for estimation with z low = z mid = 0. Black curves correspond to estimations using z low and z mid as nuisance parameters. We note the differences between the three models are very small.
for intermediate redshifts 0.7 z high 0.9 there is a transition, and above z high ≈ 0.95 there is also a constant trend with Ω m ≈ 0.294 and σ Ωm ≈ 0.011. At our choice, used in this work, z high = 0.9, we get Ω m = 0.297 ± 0.012 in agreement with Eq. (18).
We also note that the noise present in the red curves of Fig.10 is not displayed in the 1-dimensional test; we assume that this is consequence of the fact that E 3 and E 2 were fixed to E 1 through their ΛCDM constraints. No SNe z > 0. 9 No z high cut z high = 0.9 FIG. 10. q0-s0 contour plots for three models with different behaviors at high z. Blue contours correspond to the standard E is method, the red ignores the z high cut, and the gray does not use SNe at redshifts higher than z = 0.9. This analysis forecasts the high importance of getting limits over z high .