f(R) gravity modifications: from the action to the data

It is a very well established matter nowadays that many modified gravity models can offer a sound alternative to General Relativity for the description of the accelerated expansion of the universe. But it is also equally well known that no clear and sharp discrimination between any alternative theory and the classical one has been found so far. In this work, we attempt at formulating a different approach starting from the general class of f(R) theories as test probes: we try to reformulate f(R) Lagrangian terms as explicit functions of the redshift, i.e., as f(z). In this context, the f(R) setting to the consensus cosmological model, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda $$\end{document}ΛCDM model, can be written as a polynomial including just a constant and a third-order term. Starting from this result, we propose various different polynomial parameterizations f(z), including new terms which would allow for deviations from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda $$\end{document}ΛCDM, and we thoroughly compare them with observational data. While on the one hand we have found no statistically preference for our proposals (even if some of them are as good as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda $$\end{document}ΛCDM by using Bayesian Evidence comparison), we think that our novel approach could provide a different perspective for the development of new and observationally reliable alternative models of gravity.


Introduction
Almost 20 years have passed since distant Ia supernovae showed for the first time that the universe is expanding in an accelerated way [1][2][3]. Since then, this evidence has been supported by many other observations as cosmic microwave background anisotropies (CMB) [4] and large scale structure [5][6][7][8][9][10] but the question of the origin of this acceleration is still waiting for an answer. Despite the favorable experimental results which support the standard paradigm [4,11], the a e-mail: ruth.lazkoz@ehu.es b e-mail: maria.ortiz@ehu.eus c e-mail: vincenzo.salzano@usz.edu.pl ΛCDM model − a cosmological constant [12][13][14] plus cold dark matter in the framework of general relativity (GR) − there are some important theoretical shortcomings in it as well as tensions between data and theory [15], the latest one being the discrepancy [16] between the values of the Hubble constant (H 0 ) as measured by Planck [4] (i.e. assuming ΛCDM as the background cosmology) and the local measurements from Cepheids [17].
All these facts raise the need for the formulation of many other approaches and theories. Many authors do not leave the realm of GR, and explore dynamical alternatives to the cosmological constant, the so-called dark energy models, where a new component of the mass-energy tensor is added. These alternative settings can be accomplished in the most varied ways, from the theoretical side [18,19] to the phenomenological one, where the most common approach is a generalization of the cosmological constant itself by making the equation of state have additional parameters and/or geometry dependence [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. Many others try to give up on GR itself and propose entirely new alternative theories of gravity, also known as modified gravity theories. Even in this case, the number of ways in which GR can be extended and/or modified is extremely large; for a non-exhaustive list see, for example, [36,37], but the border between the dark energy and the modified gravity formulation is not clearly paved [38][39][40][41].
Among the plethora of alternative models proposed, in this work we will focus on the so-called f (R) framework, also known as fourth-order theories or extended theories of gravity [42][43][44][45][46]. The easiest and most basic approach to this kind of theories consists in replacing the Ricci scalar, R, appearing in the Hilbert-Einstein action of GR with a general function f (R). The f (R) proposals, of course, cannot be arbitrary, but they have to be able to both fit the cosmological data and to satisfy the Solar System constraints, given that on such scales GR has been experimentally tested and confirmed. Some examples of such viable models are in [47][48][49][50][51][52]. These theories have also been intensively analyzed, by comparing them with cosmological probes [53][54][55][56][57][58][59][60], with the internal dynamics of many different-scale gravitational structures [61][62][63][64][65][66][67], and with cosmological simulations [68][69][70][71][72][73][74], eventually putting weather too weak or too severe (i.e. very consistent with the ΛCDM limits) constraints on their viability.
The most common approach when one tries to accomplish this task is to first propose an f (R) expression at the level of the Lagrangian, with the requirement that it satisfies some priors (e.g. Solar System constraints), then to solve the corresponding dynamical equations, and finally to test them against the data. The main obstacle in this procedure is that the fourth order differential equations which come out from a very general f (R) Lagrangian are not analytically solvable. For cosmological applications, consider that the Friedmann equations can result to be quite complicated third order differential equation in H [75], the expansion rate function which is generally needed in order to calculate cosmological distances, and we might miss an easy analytical expression for it (see [76] for an application to the famous model from [48]). In order to overcome such difficulties, thus, weather one needs to fix a functional form for the f (R) function in order to obtain analytically manageable equations [53], or some phenomenological ansätze for an H compatible with the given f (R) have to be proposed [77]. Alternatively, one can go the other way around, that is, to recover analytically a specific f (R) model from any requested standard cosmology [78] or given any H (t) proposal [75,76], or to reconstruct numerically the f (R) from the data [79]. It has not to be forgotten, on the other hand, that for f (R) theories a scalartensor equivalence holds [80], meaning that a scalar field can be introduced, coupled to gravity (geometrically described by a standard f (R) = R) and following a given potential, which is completely equivalent to the given f (R) model.
In all the previous cases, anyway, there are some flaws, or limitations: in one way or another, we need to make assumptions at some step, so that the results will be always somewhat biased by these choices and will not be as general as they should or as we would like; in some cases, the relation between f (R) and H is generally approximate and some information might be lost (especially if advocating the scalar-tensor equivalence). One usual approach is to take into account that the f (R) cosmology should mimic the ΛCDM model as much as possible given that, at the end of the day, this model gives the best description to most of the cosmological data available nowadays. Thus, the f (R) cosmology should recover a matter dominated stage at high redshift, and should show accelerated expansion at low redshift, ideally, without a true cosmological constant, but through purely geometrical terms. Apart from setting some general limits based on these considerations, there really is not much more information one can provide a priori in order to yield sensible f (R) theories.
Our present study is precisely related to this topic. We want to explore this problem but using a different approach, a sort of practical poor cosmologist approach. The vast amount of probes in cosmology provides us with a big variety of measurements usually corresponding to specific values of the redshift; thus, we think it would be interesting to have gravity models which are described in terms of this observable, so that we are provided with a more straightforward way to test the models against the observations. Studies exist in the literature where the redshift formulation appears [77], however they are not focused on constructing a sensible f (z) model but rather on proposing ansätze both for the Hubble function H and for the f (R) function, and then fitting the parameters in order to do a cosmographic reconstruction. Our present analysis will seek a way to provide reasonable f (z) models since the beginning, i.e., from the action, so that at a later stage one can numerically solve the dynamics of the Universe, to test their validity and study their deviation with respect to the ΛCDM scenario.
We begin with Sect. 2 by introducing briefly a very general approach to f (R) gravity, with a description of the modified Friedmann equations and how we build f (R) as a function of the redshift z. The expression of the derivatives of f (R) with respect to R, and of R with respect to time are provided in terms of derivatives with respect to the redshift z. We then consider various different dark energy and modified gravity scenarios and calculate the corresponding high and low redshift limits to validate the choice of our proposals. In Sect. 3 we describe the observational data used in our analysis; in Sect. 4 we discuss the obtained results and in Sect. 5 we give some conclusions.

From f (R) to f (z)
We consider the most general f (R)-type modification to the Einstein-Hilbert action described by the action [44] where g is the determinant of the metric g μν and L m is the Lagrangian of any considered energy-matter component. In order to obtain the field equations one has to vary the action with respect to the metric field, g μν , ending up with where R μν is the Ricci tensor, and ∇ are respectively the d'Alembertian and Laplacian operators, T m μν is the stress energy tensor, and we define f R ≡ d f/d R. In the remainder of the paper we assume a background Friedmann-Lemaître-Robertson-Walker (FLRW) metric in spherical coordinates with c, the speed of light in vacuum and a(t), the scale factor, and we restrict our considerations to spatially flat spaces (k = 0), with matter and radiation as the only contribution to the stress-energy tensor. Thus we will consider f (R) as a purely geometrical contribution even if, as stated before, we can always find a scalar field, entering the stress-energy tensor, and completely equivalent to the original proposed f (R). Finally, we obtain the generalized Friedmann equations which govern the dynamics of the universe at large scales, namely [75] where · ≡ d/dt and we have defined together with the energy conservation equation for standard energy-matter perfect fluids: where the suffix X stands for matter, radiation or any fluid in the stress-energy tensor. For our goals, the next step is to perform a change of variables in order to have all the derivatives appearing in Eq. (4) in terms of the redshift. First, we note that by combining the Friedmann equations we obtain which is the usual definition for the Ricci scalar for homogeneous and isotropic flat FLRW spacetimes. From this, we calculate where the subindex z means derivative w.r.t. the redshift. Finally, one gets: In this work we will provide f (R) as f (z), we will solve Eq. (4) numerically for H , and we will compare it to observational data. In terms of the redshift, the first Friedmann equation Eq. (4) now reads: with Clearly, here we have a third order differential equation, for which we need to set initial conditions for H , H z and H 2z ; we will discuss in the next section how we choose such initial conditions.

Requirements for f (z) proposals
Once we have the equation for H , the next step is to provide a "good" f (z) model to test against data. We know that a spatially flat ΛCDM universe can be expressed as the f (R) theory: If we assume the universe filled with matter and radiation, then the first Friedmann equation can be cast into the following form: where we have defined the dimensionless density parameters as with i = m, r indicating, respectively, matter and radiation, and 1−Ω m −Ω r corresponding to the cosmological constant.
To get f Λ in terms of the redshift we can use Eq. (7), thus obtaining Here we have normalized by H 2 0 , the Hubble constant H 0 ≡ H (z = 0), just to simplify notation, and we will keep this notation for the whole analysis.
Taking this into account we would like to choose an f (z) which is somehow the simplest generalization of the latter expression, that is, a polynomial with more terms and not only a constant and a third-order one. In order to set some restrictions when choosing a specific model, a useful analysis is to study the high and low redshift limit of R and f (R) for the case of ΛCDM and other different models of dark energy or modified gravity theories. At the end of the day, even if ΛCDM may not be the model which really underlies our universe, it is the one which, so far, best describes it. Thus, any generalization, should have a behaviour which should follow its same trends. If we were able to detect some special feature, we could propose a reasonable f (z) according to this.
Actually, Eqs. (15)-(16) exactly provide us with such information; it is straightforward and somehow trivial to verify that the high and low redshift limits of ΛCDM for both the Ricci scalar and the f Λ (R) function (which in this case are identical, because of GR) read Let us generalize a little bit such scenario, considering the most popular and common dark energy model used in the usual references, the so-called Chevallier-Polarski-Linder [24,25] (CPL) parametrization. In this model the dark energy fluid has a dynamical equation of state which is linear in the scale factor: The well-known expression for the first Friedmann equation in this case is: from this, we can compute the Ricci scalar: As we are still working in the context of general relativity, f C P L = R and the CPL dark energy fluid is included in the stress-energy tensor. Let us notice that w(a = 1) = w 0 has to be negative today, in order to lead to the observed accelerated expansion of the universe; and we can further assume the strong prior w 0 + w a < 0, as confirmed by observations [4], which implies that the universe was matter dominated at early times. Accordingly, while for the low redshift limit one finds: = const.
Setting w 0 = −1 we would get the limit for the ΛCDM model. Another possibility is given by the "early dark energy" models [81,82]. We will focus on the model discussed in [83]. It is a phenomenological scenario which behaves as a cosmological constant at early times and decays rapidly at late times: where Ω ee is the fractional energy density of the early dark energy today, and a c = 1/(1 + z c ) is the critical value of the scale at which it shifts from the early-time behaviour to the late-time behaviour. Computing the Ricci scalar in terms of the redshift, one gets Taking the limits one finds: Thus, in these two cases, we qualitatively recover the same limits as ΛCDM. While this is somehow expected, because both the CPL parametrization and the early dark energy models are generalizations of the cosmological constant and contain it as a special case, we are also interested in exploring the limits of more radically different approaches. For example, we will consider the Ricci dark energy model [84] which belongs to the so-called Holographic Dark Energy models. The basic idea behind this class of models [85] is that our universe is in a sense finite and can be described by a two dimensional spherical holographic screen, thus there must be finite size effects. At least on a theoretical background, there is a big conceptual difference with respect to a "simple" cosmological constant; moreover, another interesting feature here is that ΛCDM is not explicitly included in this model. This difference, anyway, as it happens in many cases, does not necessarily translate into a difference in the quantitative description of the observational data. In this model the cosmological evolution is governed by from which the Ricci scalar reads with One can notice that the first term in R will dominate as far as γ < 2. And, actually, this parameter has been constrained with observational data [86] to be γ = 0.325 +0.009 −0.010 , so that one can write The low redshift limit exhibits, like in the previous cases, a constant behaviour: We also consider another well-known example of modified gravity, the Dvali-Gabadadze-Porrati (DGP) model [87], in which gravity leaks off the four dimensional Minkowski brane into the five dimensional bulk Minkowski space-time and such setting should yield a self-acceleration of the universe without introducing dark energy. The dynamics for this model is given by the modified Friedmann equation [88] As in the previous case, the ΛCDM model is not a sub-case within this theory. Plugging Eq. (33) into Eq. (7) we get: with and where r c is the cross-over scale that governs the transition between four-dimensional behavior and five-dimensional behavior. One can easily see that as z → ∞, the biggest contribution corresponds to while for the low redshift limit we have: with It is important to note that the four-dimensional part of the total DGP action has f DG P (R) = R. A generalization of this model can be found in [89].

Final proposals
Given the previous examples and reminding that for all of them we have f (R) = R, we can see that in many relevant cases in the literature these hold: then we can conclude that it is a reasonable choice to propose polynomial expressions with a maximum third-order term as extensions of the f Λ (R). Anyway, one could ask what would happen with more general f (R) models alternative to those described above. For example, in [54], the authors consider the model Following [54] it is easy to deduce that from which we have This result is not in contrast with previous high redshift limits. Again, in [54], another f (R) model is considered which given by Again, we obtain with A = Ω m H 2 0 α −1 . In this case, which implies that Thus, again, we have the same high redshift limit. Note also that in both the f (R) models we have just considered, the low redshift limit is not a proper constant, which means, they do not include a cosmological constant. Of course, we cannot check here all possible models, as this would be out of the purpose of this work, and would be a gigantic, yet useless, exercise, also because even the more commonly used, as the one in [48], can not have an analytical solution for H [76]. But it is clear, that the limits we have described so far, clearly depict a possible trend which we use as guideline for our proposals. Thus, we are going to choose different models to see how much each model deviates with respect to ΛCDM. Let us notice that some of the models will contain a constant term (so they might resemble a ΛCDM model) and others will not. The latter may have more interest if we want to provide fully alternative theories to the standard model of cosmology. Finally and all in all, the models we have decided to focus on are:

Data
We use the combination of various current observational data to constrain the f (R) = f (z) models described previously. In this section, we describe the cosmological observations used in this work. We will only consider the observational data related to the expansion history of the universe, i.e., those describing the distance-redshift relations. Specifically, we use the Type Ia Supernovae (SNeIa), the cosmic microwave background (CMB) distance priors, the Baryon Acoustic Oscillations (BAO) data, the expansion rate data from early-type galaxies (ETG), plus a prior on H 0 .

Hubble data
We use a compilation of Hubble parameter measurements estimated by the differential evolution of passively evolving early-type galaxies used as cosmic chronometers, in the redshift range 0 < z < 1.97, and recently updated in [90]. The corresponding χ 2 H estimator is defined as with σ H (z i ) the observational errors on the measured values H obs (z i ), θ the vector of the cosmological background parameters. Moreover, we will add a gaussian prior, derived from the Hubble constant value given in [91], H 0 = 69.6 ± 0.7.

Type Ia supernovae data
We used the SNeIa data from the JLA (Joint-Light-curve Analysis) compilation [92]. This set is made of 740 SNeIa obtained by the SDSS-II (Sloan Digital Sky Survey) and SNLS (Supernovae Legacy Survey) collaborations, covering the redshift range 0.01 < z < 1.39. The χ 2 in this case is defined as with ΔF = F theo − F obs , the difference between the observed and the theoretical value of the observable quantity for SNeIa, the distance modulus; and C SN the total covariance matrix (for a discussion about all the terms involved in its derivation, see [92]). The predicted distance modulus of the SNeIa, μ, given the cosmological model and two other quantities, the stretch X 1 (a measure of the shape of the SNeIa light-curve) and the color C, is defined as where D L is the luminosity distance given by with E(z) = H (z)/H 0 (following [92], only for SNeIa analysis we assume H 0 = 70 km/s Mpc −1 ) and c the speed of light measured here and now. In this case the vector θ b will include cosmologically-related parameters and three other fitting parameters: α and β, which characterize the stretchluminosity and color-luminosity relationships; and the nuisance parameter M B , expressed as a step function of two more parameters, M 1 B and Δ m : where M stellar is the mass of the host galaxy. Further details are given in [92].

Baryonic acoustic oscillations
The χ 2 B AO for Baryon Acoustic Oscillations (BAO) is defined as where the quantity F B AO can be different depending on the considered survey. We used data from the WiggleZ Dark Energy Survey, evaluated at redshifts 0.44, 0.6 and 0.73, and given in Table 1 of [6]; in this case the quantities to be considered are the acoustic parameter and the Alcock-Paczynski distortion parameter where D A is the angular diameter distance and D V is the geometric mean of the physical angular diameter distance D A and of the Hubble function H (z), and defined as We have also considered the data from the SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) DR12, described in [93] and expressed as where r s (z d ) is the sound horizon evaluated at the dragging redshift z d ; and r f id s (z d ) is the same sound horizon but calculated for a given fiducial cosmological model used, being equal to 147.78 Mpc [93]. The redshift of the drag epoch is well approximated by [94] The sound horizon is defined as: with the sound speed and with T C M B = 2.726 K. Finally, we have also added data points from Quasar-Lyman α Forest from SDSS-III BOSS DR11 [95]:

Cosmic microwave background data
The χ 2 C M B for Cosmic Microwave Background (CMB) is defined as where F C M B is a vector of quantities taken from [96], where Planck 2015 data release is analyzed in order to give the socalled shift parameters defined in [97]. They are related to the positions of the CMB acoustic peaks which depends on the geometry of the model considered and, as such, can be used to discriminate between dark energy models of the different nature. They are defined as: where we introduce the baryonic density parameter, Ω b . As before, r s is the comoving sound horizon, evaluated at the photon-decoupling redshift z * , given by the fitting formula [98]: with while r is the comoving distance defined as:

Monte Carlo Markov chain (MCMC)
In order to test the predictions of our theory with the available data, we implement an MCMC code in order to minimize the total χ 2 defined as As described in Sect. 2, we will need to solve a third-order differential equation in H , in order to recover this quantity and calculate all the required observational signatures. This means we need to specify the initial condition for H , H z and H 2z . Contrarily to what is done in some literature, where similar equations are solved assuming that the model behaves as ΛCDM from z → ∞ up to some finite high redshift value in order to assure a deviation from ΛCDM only on the narrow redshift range covered by the data, we have decided to leave more freedom to our models to adjust to observations. In particular, we will fix the "initial conditions" for our H (z) only at one single redshift point, with no restrictions to its behaviour for both larger and smaller redshift values. The parameters we need to specify are: z min : this redshift corresponds to the early universe regime, i.e. very high z. While in theory one should fix z = ∞, due to numerical issues we will take it as finite, and we set z = 10 10 ; z max : this value corresponds to the present time and it is set to z = 0. Both z min and z max define the redshift interval where the differential equation is solved; z pivot : this point is used to set some "initial conditions" of our differential equations system. In our case we impose a ΛCDM model on this single point. Note that in order to study possible influences of the initial conditions on the final results, we have performed our calculations for the values z piv = 10, 100.
Moreover, analysing carefully the range of values of the parameters in many simulations (by setting different initial conditions) we have decided to apply the following priors on the free parameters in our approach, i.e.: We have verified a posteriori that these choices are licit and do not induce any further strong restriction onto the parameters. Finally, in order to set up the reliability of one model against the other, we use the Bayesian Evidence, which is generally recognized as the most reliable statistical comparison tool even if it is not completely immune to problems, like its dependence on the choice of priors [99]. We calculate it using the algorithm described in [100]; as this algorithm is stochastic, in order to take into account possible statistical noise, we run it ∼ 100 times obtaining a distribution of values from which we extract the best value of the evidence as the median of the distribution. The Evidence, E, is defined as the probability of the data D given the model M with a set of parameters θ , E(M) = dθ L(D|θ , M) π(θ|M), where π(θ|M) is the prior on the set of parameters, normalized to unity, and L(D|θ , M) is the likelihood function.
Once the Bayesian Evidence is calculated, one can obtain the Bayes Factor, defined as the ratio of evidences of two models, M i and M j , Even if the Bayes Factor B i j > 1, one is not able yet to state how much better is model M i with respect to model M j . For this, we choose the widely-used Jeffreys' Scale [101]. In general, Jeffreys' Scale states that: if ln B i j < 1, the evidence in favor of model M i is not significant; if 1 < ln B i j < 2.5, the evidence is substantial; if 2.5 < ln B i j < 5, is strong; if ln B i j > 5, is decisive. Negative values of ln B i j can be easily interpreted as evidence against model M i (or in favor of model M j ). In [99], it is shown that the Jeffreys' scale is not a fully-reliable tool for model comparison, but at the same time the statistical validity of the Bayes factor as an efficient model-comparison tool is not questioned: a Bayes factor B i j > 1 unequivocally states that the model i is more likely than model j. We present results in both contexts for readers' interpretation.

Results of the observational tests
The complete set of free parameters in our analysis is Ω m , Ω b , h, f i , α, β, M B , Δ m , where f i are the parameters corresponding to the various polynomial orders we have considered for each of the proposed parametrization of f (z). We will focus in our comments only on the matter density parameter Ω m and on the f i parameters, given that the other parameters are fully limited by imposed priors (e.g., Ω b and h), or are independent of the cosmological background (e.g., SNeIa parameters, α, β, M B , Δ m ). In Tables 1 and 2 we report the results obtained for such parameters in terms of the corresponding confidence levels, for the two different choices of Table 1 Results for z piv = 10 Model  Table 2 Results for z piv = 100 Model −0.03 for z piv = 10 and 3Ω m = 0.89 +0.02 −0.02 for z piv = 100, which perfectly agrees with the corresponding free estimations of f 3 . We can also calculate 6(1 − Ω m − Ω r ) = 4.14 +0.06 −0.09 for z piv = 10 and 6(1 − Ω m − Ω r ) = 4.21 +0.04 −0.04 for z piv = 100, which also agree with our free estimations of f 0 , mainly because of the larger errors on this parameter than on Ω m .
Before discussing generalizations and/or deviations from ΛCDM, let us note that the value of Ω m does not really change from one case to another; there is tiny trend toward smaller values for z piv = 100 than for z piv = 10, but all the values are perfectly consistent at 1σ level.

Models 2-4
Now, considering case by case, let us consider model (2), which generalizes ΛCDM by adding intermediate powers corresponding to the f 1 and f 2 parameters. Given the very good agreement of ΛCDM with cosmological background data, we expect small deviations from it, if any. Actually, we can see how f 1 and f 2 are O(10 −2 ) − (10 −3 ) at 1σ confidence level. What is even more important, is that they are also consistent with zero, at 1σ for z piv = 10 and maximum at 2σ for z piv = 100. Thus, we can conclude that any f (z) with this form is basically indistinguishable from ΛCDM at the present stage of observational data.
We also note that the simultaneous presence of both f 1 and f 2 seems to be tied to some degeneracy between the two parameters. In fact, for the models (3) and (4), where only one of the two parameters is present, the corresponding likelihoods are much more regular and far less noisy than in the case of model (2).
Anyway, we also stress that models (2)-(4) have in common the presence of a constant term f 0 which, in some way, can bias the possible detection of any departure from ΛCDM because it can always be considered as representing a cosmological constant (on the limit z → 0). Actually, note that also for models (2)-(4) the same relations we have described for the ΛCDM case still hold. Eventually, a departure or agreement for them could tell us about the effective reliability and weight of the terms f 1 and f 2 . In particular, we find that the relation 3Ω m = f 3 is always perfectly satisfied by all models for both the pivot redshift values we have considered. The same holds true for the condition 6(1 − Ω m − Ω r ) = f 0 , except for model (2), which exhibits a value of f 0 clearly different from all the other cases, even if with larger errors.
Thus, we are not in the position to establish if the role of f 1 and f 2 is really effective and/or needed, or not: such relations should not hold in models (2)-(4), because for them f (z) = f Λ (z), but we find a good agreement, so that one might be led to think we are not actually detecting any deviation from the ΛCDM model. Moreover, the new parameters f 1 and f 2 are very small so that their role is really less significant than those of f 0 and f 3 . But we cannot avoid to comment that these results strongly depend on the present observational status; future more precise data could help to discriminate between one model or another.

Models 5-8
More interesting are, from some point of view, models from (5) to (8) which explicitly avoid the introduction of the constant term f 0 so that they are not reducible to ΛCDM in any way. On the one hand, the first point to note is that, even in these cases, the parameters f 1 and f 2 , when included, must be very small in the light of observations, ranging from O(10 −2 ) to O(10 −7 ). On the other, we see that the parameters corresponding to the lowest order powers 1/2 and 1/4 are very well constrained and the likelihood profiles have a very regular Gaussian shape. Furthermore, they seem to provide an equally good fit for the data with respect to ΛCDM with the same number of theoretical parameters.
Thus, the point to be understood here is: are these low order terms a real possible alternative to the cosmological constant, or should we consider them as a "smeared" version of the standard case only due to the low (for such type of discrimination) accuracy of the data? In fact, we have to consider that the power of the redshift we are considering might be low enough so that, due to the observational constraints we are using, they are actually mimicking the effective behaviour of a constant.
A possibly obvious trend is that the higher the order, the smaller the value of the corresponding parameter f i is (seen by comparing f 0 to f 1/4 and f 1/2 ); but this trend might be expected, given that f 1/2 and f 1/4 are by construction coupled to redshift dependent terms which could compensate the magnitude value with time dependence.

Bayesian evidence
As we have mentioned before, apart from the information one can extract from the values of the parameters which best agree with the data, we have also computed the Bayes factor for each one of our proposed models against ΛCDM, the reference model. This statistical tool should provide us information about the reliability of each model. In general, we can see how | ln B i Λ | < 1 so that there is no statistical significant preference for one model with respect to another.

Summary and conclusions
In this work we have introduced a different approach to convert general f (R) theories in f (z) models, which should be more easily connected to observational data and thus could shed some light on the explanation of the accelerated expansion of the universe and provide a more straightforward formulation to perform tests against observations.
The most attractive scenario would be to be able to solve the Friedmann equations analytically to find H (z) solutions by proposing some f (z) or viceversa. However, as we have mentioned in the introduction, this formulation is not possible in general, neither when studying the problem in terms of the Ricci scalar nor in terms of the redshift.
Then, what we have done has been to start from some f (z) ansätze and perform the analysis numerically. By studying some of the most known dark energy models we have been able to find some general potentially interesting features in order to shed some light on the expression of our proposals and restrict their arbitrariness by imposing some physically well-sounded and expected trends. In particular, we have found that a simple algebraic expression of a polynomial including just a constant and a third-order term can describe the general behaviour whatever f (R) = f (z) model is expected to have in order to mimic ΛCDM at both high and low redshift. We need to stress that this polynomial does not need to include powers higher than order three, as we have seen by analysing different dark energy models. Even if such analysis cannot be exhaustive, we know that ΛCDM works fine with most of the data we have and we expect that, if any, deviations from it might be small.
Thus, we have selected eight proposals varying the number of free parameters in order to analyze their viability according to the data. In particular, we have included middle-order terms in between the constant and the third-order one, checking for their compatibility with observations and their statistical soundness. And we have also proposed models with no constant term at all, so that they cannot be reduced at a ΛCDM scenario in any way.
A general conclusion we can extract from this work is that there are some f (z) polynomial models which are competitive with ΛCDM at the background level. We were especially interested in models (5) and (7) because they explicitly do not include a Λ-like term and, in fact, we have evidence indicating that these models are as reliable as ΛCDM when it comes to analyse observational data. Even though the Bayesian Evidence does not claim any significant difference with respect to ΛCDM, we still think this is an interesting result, which could be a useful guide when formulating and studying manageable alternative models of gravity. Moreover, we want to stress and keep in mind that f (z) theories will probably not be the definitive answer to explain why the universe is expanding at an accelerated rate; but that was not really our objective here. Our goal was rather to provide a different perspective on how to relate the theoretical proposal of a modified gravity with observationally related quantities, and how to extract any kind of useful information which might contribute to the development of observationally reliable theories of gravity.