Observational constraints on a phenomenological fR,∂R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\left( R,\partial R\right) $$\end{document}-model

This paper analyses the cosmological consequences of a modified theory of gravity whose action integral is built from a linear combination of the Ricci scalar R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document} and a quadratic term in the covariant derivative of R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document}. The resulting Friedmann equations are of the fifth-order in the Hubble function. These equations are solved numerically for a flat space section geometry and pressureless matter. The cosmological parameters of the higher-order model are fit using SN Ia data and X-ray gas mass fraction in galaxy clusters. The best-fit present-day t0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{0}$$\end{document} values for the deceleration parameter, jerk and snap are given. The coupling constant β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} of the model is not univocally determined by the data fit, but partially constrained by it. Density parameter Ωm0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{m_0}$$\end{document} is also determined and shows weak correlation with the other parameters. The model allows for two possible future scenarios: there may be either an eternal expansion or a Rebouncing event depending on the set of values in the space of parameters. The analysis towards the past performed with the best-fit parameters shows that the model is not able to accommodate a matter-dominated stage required to the formation of structure.


Introduction
The cosmological constant is the simplest solution to the problem of the presentday acceleration of the universe. It is consistent with all the available observational tests, from the constraints imposed by CMB maps [1][2][3][4] to BAO picks [5][6][7][8] to SN Ia constraints [9][10][11][12][13][14] and the estimates of fraction of baryons in galaxy clusters through X-ray detection [15,16]. Nevertheless, there is an uncomfortable feeling from the fact that is not satisfactory explained as a vacuum state of some field theory: the energy scale observed for and predicted by quantum field theory badly disagree [17]. This is one of the reasons why alternative explanations for the present-day cosmic acceleration have been put out since 1998, when it was first detected [18,19].
These proposals can be classified into two broad categories, each one related to a modification of a different side of the Einstein field equations of General Relativity [20].
The modified matter approach introduces a negative pressure (that counteracts gravity) through the energy momentum tensor T μν on the right-hand side of Einstein equations. Three types of models constitute this approach, namely quintessence models [21][22][23][24][25], k-essence solutions [26,27] and perfect fluid models (such as the Chaplygin gas [28,29]). The quintessence and k-essence models are based on scalar fields whose potential (in the former case) or the kinetic energy (in the later one) drive the accelerated expansion. The perfect fluid models use exotic equations of state (for the pressure as a function of the energy density) to generate negative pressure. The weakness of k-essence and quintessence proposals lies in the fact that they require a scalar field to be dominant today, when the energy scale is low. This requirement demands the mass of such a field φ to be extremely small. One would then expect long range interaction of this field with ordinary matter, which is highly constrained by local gravity experiments. The difficulty with the perfect fluid models is the lack of an interpretation to the exotic equations of state based on ordinary physics.
The second approach, called modified gravity, affects the left hand side of Einstein equations [30,31]. The subtypes of models in this category include braneworld scenarios [32,33], scalar-tensor theories [34,35], non-riemannian geometries [36][37][38], f (R) gravity [39][40][41][42][43][44] and f (T ) theories [45][46][47]. In the braneworld models there exists an additional dimension through which the gravitational interaction occurs and by which the possibility of a cosmic acceleration is accommodated. Scalar-tensor theories couple the scalar curvature R with a scalar field φ in the action S, allowing enough degrees of freedom for an accelerated solution. Non-riemannian geometries are another alter-native way to explore the geometrization of the gravitational interaction, but allowing new degrees of freedom describing different geometric structures than curvature. The f (R) gravity corresponds to models resulting from non-linear Lagrangian densities dependent on R. f (T ) theories which can be interpreted as an extension of teleparallel gravity [48], build cosmological models on manifolds whose connection is equipped with torsion instead of curvature.
All these mechanisms of generating the present-day cosmic dynamics are generally called dark energy models.
The two categories mentioned above concern the Cosmology obeying the principle of large-scale homogeneity for the universe. However, on the scales less than (or of the order of) 100 Mpc, the universe is clearly non-homogeneous, exhibiting voids and filaments of aggregate of matter [49]. As gravity is a non-linear theory, there is the possibility of these inhomogeneities in the energy content of the cosmos to engender a feedback, ultimately leading to a dynamics interpreted as a global cosmic accelerated expansion. These ideas work as an alternative to dark energy and appear in the context of inhomogeneous cosmologies [50,51], such as Buchert and Räsänen's back-reaction model [52][53][54], and Wiltshires timescape cosmology [55][56][57].
Our aim with the present paper is to contribute with the study of dark energy in the branch of modified gravity (therefore, we work in the context of the standard homogeneous cosmology). We explore the cosmological consequences of a particular model in the context of the f (R, ∂ R) gravity testing the model against some of the observational data available. f (R, ∂ R) gravity is a theory based on an action integral whose kernel is a function of the Ricci scalar R and terms involving derivatives of R. Ref. [58] is perhaps an early example of what we call f (R, ∂ R) gravity. A more recent example is the (non-local) higher derivative cosmology presented in Refs. [59][60][61]. These works intend to eliminate the Big Bang singularity through a bounce produced by a theory based on an action containing the regular Einstein-Hilbert term plus a String-inspired term. This last ingredient contains an infinite number of derivatives of R, a requirement needed to keep this modified gravity theory ghostfree. References [62,63] are f (R, ∂ R)-gravity prototypes concerned not with the past-incomplete inflationary space-time but with the cosmological constant problem. In this regard, they try to address a similar issue that interests us in the present work: the cosmological acceleration nowadays.
Few years ago, we proposed a gauge theory based on the second order derivative of the gauge potential A following Utiyama's procedure [64]. This led to a second field strength G = ∂ F + f AF (depending on the field strength F) and allowed us to obtain massive gauge fields among other applications, such as the study of gauge invariance of the theory under U (1) and SU (N ) groups when L = L A, ∂ A, ∂ 2 A . In a subsequent work, gravity was interpreted as a second order gauge theory under the Lorentz group, and all possible Lagrangians with linear and quadratic combinations of F and G− or the Riemann tensor, its derivatives and their possible contractions, from the geometrical point of view-were built and classified within that context [65]. In a third paper we selected one of those quadratic Lagrangians for investigating the cosmological consequences of a theory that takes into account R-and (∂ R) 2 -terms [66]. That Lagrangian appears in the action S below-Eq. (1)-and can be understood as a particular example of an f (R, ∂ R) gravity: where χ = 8π G (once we take c = 1), β is the coupling constant, ∂ μ is the derivative, and L M is the Lagrangian of the matter fields. As usual, g = det g μν . 1 An important feature of S is the presence of ghosts (see [61] and references therein). However, here this is not a serious issue because (1) should be considered as an effective action of a fundamental theory of gravity, this last one a ghost-free theory. There are some consistent proposals of extended theories of gravity in the literature. A good example is the String-inspired non-local higher derivative gravitation model presented in Ref. [60]. The ghost-free action of that model is: where M P is the Planck mass (M 2 P = 1/ (8π G), G is the gravitational constant), M * is the mass scale where the higher-derivative terms become important and In this context, our higher-order model must be treated as a low energy limiting case of (2) where the only relevant contribution of Eq. (3) for the present-day universe is the first term of the series, i.e. f 1 f n , ∀n 2. In this approximation, Eq. (2) is equivalent to the action integral (1) of our f (R, ∂ R)-cosmology up to a surface term. The equivalence of both models leads to a relation between the coupling constant β and the parameter f 1 ∼ M −2 * , namely Taking the variation of (1) with respect to the metric tensor g μν (cf. the metric formalism) gives the field equations with where = ∇ μ ∇ μ , ∇ μ is the covariant derivative and R μν is the Ricci tensor.
A four dimensional spacetime with homogeneous and isotropic space section is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) line element where the curvature parameter κ assumes the values +1, 0, −1 and a (t) is the scale parameter. Substituting the metric tensor from Eq. (7) along with the energymomentum tensor of a perfect fluid in co-moving coordinates, (ρ is the density and P is the pressure of the components fulling the universe), into the gravity field equations Eq. (5) leads to the two higher order Friedmann equations: where (12) which are written in terms of the Hubble function H = H (t), and where we have set κ = 0 as is largely favored by the observational data. The dot (such as inȧ) means derivative with respect to the cosmic time t and for any arbitrary i.e. H (5) is the fifth-order derivative of the Hubble function with respect to the cosmic time.
Equations (9) and (10) reduce to the usual Friedmann equations if β → 0. The action (1) is invariant under translations implying in covariant conservation of the energymomentum tensor, as demonstrated in Ref. [65] by means of Bianchi identities. Since we use an homogeneous and isotropic model described by a perfect fluid with (8), Such covariant conservation is important since it implies that light will follow null geodesics (for electromagnetic field minimally coupled). This way, the luminosity distance has the usual formula given by Eq. (33) shown subsequently. In Ref. [66], we considered simultaneously the set formed by Eq. (14) and a combination of Eqs. (9) and (10), namelẏ The system of Eqs. (14)- (15) was solved perturbatively with respect to the coupling constant β for pressureless matter. Using the observational data available, we have determined the free parameters of our modified theory of gravity and the scale parameter as a function of the time. The plot of a = a (t) indicates that the decelerated regime peculiar to the universe dominated by matter evolves naturally to an accelerated scenario at recent times. This makes it promising to perform a more careful analyzes of our f (R, ∂ R) model making a complete maximum likelihood fit to all parameters. As far as we know, this statistical treatment has only been done for models dealing with derivatives up to fourth-order and this is the first time this is proposed for a theory of gravity with higher derivatives. This way, besides checking consistency/inconsistency of the theory with experimental data, we establish an algorithm for doing this. We understand this is an important result of this work, since other theories/models with higher derivatives can be tested according to the guidelines established here.
In this paper we will solve the system of Eqs. (14)-(9) exactly, without using any approximation technique. Section 2 presents the script to obtain the complete solution to the higher order cosmological equations. Another improvement we shall do in the present work (when comparing with Ref. [66]) is to use the observational data available to fully test our model. In Sect. 3 we construct the likelihood function to be maximized in the process of fitting the model parameters using data from SN Ia observations [13] and from the measurements of the mass fraction of gas in clusters of galaxies [15]. The results of this fit are given in Sect. 4, where we show the present-day numerical values for the deceleration parameter q 0 = q (t 0 ) [67], and also of the snap s 0 = s (t 0 ), These last two functions are given in terms of higher-order time derivatives of the scale factor and it is only natural that they show up here. In Sect. 4.2 we perform an investigation of the dynamical evolution of the model, comparing its behaviour towards the past with the standard solutions for single component universes. Moreover, we show that our f (R, ∂ R)-model accommodates an eternal expansion or a Rebouncing scenario depending on the values of the coupling constant β. The Rebouncing event is the point where the universe attains its maximum size of expansion and from which it begins to shrink. Our final comments about the features of our model and its validity are given in Sect. 5.

Higher order cosmological model
The solution of our higher-order model shall be obtained from the integration of the pair of Eqs. (9) and (14). The definitions and enables us to express Eq. (9) as and the conservation of energy-momentum Eq. (14) aṡ which is readly integrated: It is standard to take a 0 = 1. Moreover, substituting the solution (23) into Eq. (21), leads to a fifth-order differential equation in a: with In order to solve Eq. (24), we define the non-dimensional variables which are related to the Hubble function, deceleration parameter, jerk and snap. In fact, when these variables are calculated at the present time, t = t 0 , they give precisely the well known cosmological parameters: In the set of variables (26)-(29), one variable is essentially the derivative of the previous one, therefore they can be used to turn the fifth-order differential equation Eq. (24) into five first order differential equations to be solved simultaneously. Following the Method of Characteristics, we choose the scale factor a as the independent variable of the model (instead of the cosmic time t), obtaining the system with where we have used the relationship between the redshift and the scale factor, The system of Eq. (30) is solved numerically, which is indeed a far more simpler task than dealing directly with Eq. (24). This leads to a function H (z) (essentially equal toȧ) for each set [M] = {q 0 , j 0 , s 0 , B, m 0 } of cosmological parameters. This function is used to obtain an interpolated function for the luminosity distance d L = d L (z, [M]) of our model-defined below-which is to be fitted to the data from SN Ia observations and from f gas measurements.

SN Ia data
In order to relate the proposed model to the observational data, one has to build the luminosity distance d L [20], Our task is to solve the system Eq. (30) for obtaining H, then to integrate numerically the non-dimensional luminosity distance, for each set of cosmological parameters Function d h appears in the definition of the distance modulus, where d 0 is the reference distance of 10 pc. Observational data give μ r of the supernovae and their corresponding redshift z.
According to the likelihood method the probability distribution of the data (N SNIa is a normalization constant) is the one that minimizes where σ μ i is the uncertainty related to the distance modulus μ r i of the i-th supernova. The data set used in our calculations is composed with the 557 supernovae of the Union 2 compilation [13]. Uncertainty σ ext is given by i.e. it is the combination of the error σ v in the magnitude due to host galaxy peculiar velocities (of order of 300 km s −1 ), and the uncertainty from gravitational lensing σ lens = 0.093z (40) which is progressively important at increasing SN redshifts. Further detail can be found in Refs. [68,69], for example.

Galaxy cluster X-ray data
Measurements of the mass fraction of gas in galaxy clusters f gas can be used to constrain cosmological parameters. 2 The connection between them occurs because f gas measurements depend on the angular distance to the cluster as f gas ∼ d 3/2 . This behavior is obtained under hydrostatic equilibrium condition [70] which is satisfied in view of the conservation of the energy-momentum tensor. Usually, this analysis is done on Newtonian limit. For distances involved in f gas measurements [15,70] the Newtonian limit is valid for our model. Besides, the quantity f gas constrains efficiently the ratio b0 / m 0 and exhibit a lower dependence on the values of other cosmological parameters [16]. 3 According to Ref. [15], the mass fraction of gas relates to b0 / m 0 as where d A is the angular distance for κ = 0, .
The term d CDM The other factors in Eq. (41) are defined below: • K is a constant due to the calibration of the instruments observing in X-ray. Following Ref. [15], we take K = 1.0 ± 0.1. • Constant A measures the angular correction on the path of an light ray propagating in the cluster due to the difference between CDM model and the model that we want to fit to the data. Ref. [15] argues that for relatively low red-shifts (as those involved here) this factor can be neglected, i.e. taken as A = 1. • γ accounts for the influence of the non-thermal pressure within the clusters. According to Ref. [15], one may set 1.0 < γ < 1.1. Our choice is γ = 1.050 ± 0.029. • b (z) is the bias factor related to the thermodynamical history of the inter-cluster gas. Once again due to arguments given in Ref. [15], this factor can be modeled describes the fraction of baryonic mass in the stars. 4 Ref. [15] uses s 0 = (0. 16 and we will set α s = 0.000 ± 0.115. Analogously to Eq. (44), we define that will be required below. • The baryon parameter b0 is set to b0 h 2 = 0.0214 ± 0.002 (46) following the constraints imposed by the primordial nucleosynthesis on the ratio D/H [71].
The uncertainties in the set of parameters {b 0 , s 0 , K , b0 h 2 , α s , α b , γ } are propagated to build the uncertainty of f CDM gas . It is worth mentioning that the quantities b 0 , s 0 , K , b0 h 2 have Gaussian uncertainties. On the other hand, the parameters {α s , α b , γ } present a rectangular distribution, e.g. γ has equal probability of assuming values in the interval [1.0, 1.1]. The uncertainty of the parameters obeying the rectangular distribution is estimated as the half of their variation interval divided by √ 3.
Maybe it is useful to remark that in the case of s 0 , we adopted h Ref. [15] provides data for the fraction of gas f gas and its uncertainty σ f of 42 galaxy clusters with redshift z. These data can be used to determine the free parameter of our model by minimizing which should be maximized.

Combining SN Ia and f gas data sets
The majority of papers that present cosmological models tested through SN Ia data fit-see e.g. Refs. [25,68]-perform an analytical marginalization with respect to H 0 . On the other hand, constraint analyses using f gas data adopt a fixed value for H 0 . Our basic reference on f gas treatment, namely Ref. [15], choose H 0 = 70 km s −1 Mpc −1 . As we will consider SN Ia and f gas data sets in the same footing, it is convenient to adopt H 0 = 70 km s −1 Mpc −1 as a fixed quantity in our calculations (and not the most up-to-date value of (73.8 ± 2.4) km s −1 Mpc −1 found in [72]). We want to maximize the probability distributions given in Eqs. (36) and (48) simultaneously, which is done by maximixing or minimizing

Estimation of the free parameters of the f (R, ∂ R)-model
The parameters to be adjusted are of two types. The kinematic parameters are {q 0 , j 0 , s 0 } 5 ; the cosmological parameters are m 0 , B . Maximization of the like-   It is not surprising that small values of B (B < 10 −3 ) are disfavored by observations, because our model does not include a cosmological constant; it is therefore necessary to have a non-zero modification of GR to mimic its effect on the Hubble diagram.
The standard statistical analysis leads to the best-fit values of the parameters shown in Table 1.
In order to verify that the best-fits are good fits we plotted them against SN Ia and galaxy cluster data in Fig. 2. It is clear from Table 1 that our higher order cosmological model presents an acceleration at present time: the deceleration parameter q 0 is negative for all the values of the coupling constant B. This is a robust result indeed once we have got negative values for q 0 with 99.7 % confidence level.
Both SN Ia and f gas data sets give information about astrophysical objects at relatively low redshifts (z 1.4). That is the reason for the high values of the uncertainties on the parameter j 0 and s 0 . Positive values of j 0 mean that the tendency is the increasing of the rate of cosmic acceleration. Because the jerk parameter j 0 does not exhibit the same sign for all values of B, it is not possible to argue in favor of an overall behavior for the cosmic acceleration today. This statement is also applicable to the snap s 0 obtained with B = ±10 −1 . However, the uncertainties in the value of s 0 for B = −10 −3 do indicate that it is positive within 68 % CL. The negative value of s 0 for B = +10 −3 is also guaranteed within the same CL. This difference in the sign of s 0 for B = −10 −3 and B = +10 −3 leads to distinct dynamics for the scale factor towards the future, as we shall discuss in Sect. 4.2.
The best-fits for m 0 are statistically the same for all the chosen values of B. This is linked to the fact that the f gas data produce a sharp determination of this parameter. Notice that the value m 0 of our model is compatible with that of CDM model as reported in Ref. [13] within 95 % CL. Inspection of Table 1   The values of s 0 , although irrelevant for the evaluation of the parameters of interest, were included for completeness. The last four columns show the values of covariances σ q 0 j 0 , σ q 0 m 0 and correlation coefficients r q 0 j 0 , r q 0 m 0 The shape of the confidence regions for (q 0 , j 0 ) indicates that the deceleration parameter q 0 and the jerk j 0 are correlated parameters; this occurs for all the values of B. The four plots (q 0 , j 0 ) point to an accelerated universe at recent times within 95 % CL, but the change in the acceleration rate takes positive and negative values with equal probability.
Conversely, the shape of the plots of q 0 , m 0 in Figs. 3, 4, 5, 6 indicate that q 0 and m 0 are weakly correlated parameters. Notice that the difference between the maximum and the minimum values of m 0 in the range of the plots q 0 × m 0 is only 0.05 with 95 % CL: this confirms that the parameter m 0 is sharply determined by the best-fit procedure (something that is consistent with the fact that we have employed f gas data to build the likelihood function).
In order to quantify the correlation between the parameters (q 0 , j 0 ) and q 0 , m 0 , the covariances σ q 0 j 0 , σ q 0 m 0 and correlation coefficients r q 0 j 0 , r q 0 m 0 were evaluated according to Ref. [73]. The mean values and their respective standard deviations used to calculate the covariances and correlation coefficients were those given in Table 2. The same table shows the values obtained for the covariances and correlation functions for each value of B taken in this paper.
The values of the correlation coefficient presented in Table 2 confirm what was stated before: The parameters q 0 and j 0 are strongly (negatively) correlated while q 0 and m 0 are weakly correlated for all values of B considered here.
Comparison of the results presented in Tables 1 and 2 shows that there is no statistical difference between considering the mode and a confidence interval of 68 % (Table 1) from taking the mean and standard deviation. This emphasizes the quasi-gaussian behaviour of the probability function.

Behaviour of the f (R, ∂ R)-model with time
We can investigate the evolution of our higher order model in time by plotting the function H =ȧ/H 0 as a function of z. This is done setting the parameters at their best-fit results for each value of B.
There is agreement between our model-regardless the value of B-and the standard CDM model in the region 0 z 1. In particular, the curves of H (z) for B = −10 −1 and B = +10 −1 almost coincide: the difference between them is only noticeable at z 3. The transition from a decelerated phase to an accelerated expansion in the higher order model occurs at z 0.6 for all values of B.
Towards the future (t > t 0 ) the redshift takes negative values, −1 < z < 0. In this region our model present different dynamics depending on the value of the coupling constant. For B = −10 −1 , B = −10 −3 and B = +10 −1 , there is an onset of an ever growing acceleration eventually leading to a eternal accelerated expansion. On the other hand, the velocity in the higher order model with B = +10 −3 -black curve in Fig. 7-reaches a maximum value at z −0.15, the acceleration then changes the sign, leading to a positive deceleration parameter. This deceleration is sufficiently intense to produce a Rebouncing, where the universe reaches a maximum size and then begins to shrink. This happens at z −0.32.
This difference between eternal expansion and Rebouncing scenarios deserves a more careful statistical analysis. We shall consider this in a broader space of parameters, and not using only the best-fit values for the parameters of our model, as we have At last, it is interesting to explore the features of the model to the past. Our higherorder model is evolved from redshift z = 10 till z = 1, 000 for each value of B while varying parameters q 0 , j 0 , s 0 and 0 within the confidence interval of 2σ (see Table 2). The dynamics found through this analysis is equivalent to the behaviour of a standard Friedmann model with a fluid equation of state P = wρ constrained by 6 This shows that the proposed model does not reproduce a matter era in the interval 10 z 1000, if the domain of the cosmological parameters is the one established by our data fits. In fact, our model produces an age with characteristics close to a radiation era (w = 1/3).
The perturbation theory in the linear regime for the standard Friedmann cosmology shows that the matter inhomogeneities do not grow if the background dynamics is constrained by Eq. (51). The matter simply does not agglomerate [74]. Therefore, it is not possible to make the background observations (SNIa and f gas ) compatible with the process of structure formation in the context of our model (in the domain of parameters we have determined).

Final remarks
In this paper we have studied the phenomenological features of a particular modified gravity model built with an action integral dependent on the Ricci scalar R and also on the contraction ∂ μ R∂ μ R of its covariant derivatives. The Friedmann equations resulting from the field equation and the FLRW line element involve sixthorder time-derivatives of the scale factor a. We solved those differential equations numerically-but without any approximation-assuming a null cuvature parameter (κ = 0) and pressureless matter. Our higher order model depends on five parameters, namely [M] = q 0 , j 0 , s 0 ; m 0 , B . The present day value of the Hubble constant is taken as given-and not constrained by data fit-the reason for that being consistency with the f gas data analysis. Four of the five parameters [M] are determined by the best-fit to the data for the SN Ia d L (z) and fraction of gas in galaxy clusters d A (z) available from Refs. [13] and [15], respectively.
The shape of the likelihood function L in terms of B shows that it is not possible to determine this parameter by maximization. We were forced to select some values of B and proceed to the best-fit of the other parameters in the set [M]. It was clear that the m 0 is statistically independent of the initial conditions, regardless the values of B and the other best-fit parameters. In addition the sign of both j 0 and s 0 depend on B, so that the behaviour of the model towards the future is not unique.
The confidence regions for the pair of parameters (q 0 , j 0 ), q 0 , m 0 and ( j 0 , s 0 ) were built. They showed that q 0 and j 0 are strongly correlated parameters, whereas q 0 and m 0 are weakly correlated quantities. This was confirmed by the results in Table 2 . The plot of the confidence region for the pair ( j 0 , s 0 ) is relevant to understand the future evolution of our higher-order model, which can exhibit either an eternal expansion or a Rebouncing depending on the value of B and on the values of j 0 and s 0 .
If additional data for high z objects were used, as Gamma Ray Bursts, for instance, they might reduce the uncertainties in our estimations of the jerk and the snap, and help us to decide if the Rebouncing scenario is favored with respect to the eternal expansion one, or conversely.
A positive feature of the particular f (R, ∂ R)-model presented here is its consistency with the standard CDM dynamics at recent times. Maybe this is anticipated because the model has a great number of free parameters to be fit to the data. The high number of kinematic parameters (i.e. those involving acceleration and its first and second derivatives) is a consequence of the model being of the fifth-order in derivatives of H , a fact that could always make room for agreement with the observations. However, the freedom for bias is not unlimited. Even with all these free parameters there is a constraint to the dynamics of the model towards the past.
Besides this consistency with the standard CDM dynamics (at recent times) the model does not reproduce a large scale structure formation phase, within the domain established for the set q 0 , j 0 , s 0 , m 0 with our data fits. This shows that just fitting a model to actual data is not enough to falsify or prove the model. General essential characteristic must be present too, as the instabilities allowing the large scale structure formation. Therefore, the extra term presented in the action (1) should be ruled out as viable alternative for dark energy.
Although the model should be discarded, as a by product we have established an algorithm for fitting data for gravity theories with higher derivatives. This procedure can be used for testing other models of this type.