Extending Friedmann Equations Using Fractional Derivatives Using a Last Step Modiﬁcation Technique: The Case of a Matter Dominated Accelerated Expanding Universe

: We present a toy model for extending the Friedmann equations of relativistic cosmology using fractional derivatives. We do this by replacing the integer derivatives, in a few well-known cosmological results with fractional derivatives leaving their order as a free parameter. All this with the intention to explain the current observed acceleration of the Universe. We apply the Last Step Modiﬁcation technique of fractional calculus to construct some useful fractional equations of cosmology. The ﬁts of the unknown fractional derivative order and the fractional cosmographic parameters to SN Ia data shows that this simple construction can explain the current accelerated expansion of the Universe without the use of a dark energy component with a MOND-like behaviour using Milgrom’s acceleration constant which sheds light into to the non-necessity of a dark matter component as well.


Introduction
The Hilbert-Einstein field equations, a success of general relativity, work tremendously well at mass-to-length scales similar to those of the solar system see (e.g. [1]), where the gravitational field is moderately weak [2][3][4][5][6][7][8][9][10][11][12][13][14]. At very large mass-to-length ratios, where the gravitational field is very strong, the detection of gravitational waves produced by the interactions of compact objects such as black holes and neutron stars has shown a remarkable good agreement with the predictions of general relativity.
When the Newtonian acceleration of a test particle reaches values a 0 10 −10 m/s 2 , or equivalently when the mass-to-length ratios are much smaller than those of the solar system [15], the Hilbert-Einstein field equations cannot fit an enormous amount of astrophysical and cosmological data unless: (a) extra dark matter and/or dark energy components are added or (b) the curvature caused by the matter and energy requires extensions or modifications. In other words, the Hilbert-Einstein field equations do not correctly predict the observed results at those acceleration scales unless (a) or (b) are adopted. In what follows, we will only deal with case (b) for the accelerated expansion of the Universe at the present epoch.
If the gravitational field equations are to be extended, the first intuitive attempt consists on assuming a general f (R) function of the Ricci curvature scalar R as the Lagrangian in the gravitational action. Despite some f (R) theories have interesting results [16][17][18][19][20][21][22][23], there is currently not a full f (R) Lagrangian which solves all the shortcomings between observations and the gravitational theory. Through the years, more general actions have been proposed, using for example functions of several scalars built with the Riemann tensor and even ones in which couplings between R and the matter Lagrangian or the trace of the energy-momentum tensor [24][25][26][27][28][29][30] are adopted. These proposals are still being developed and investigated in full, and although interesting in principle, they have a small general inconvenient: the motion of free particles is not necessarily geodesic and as such, a fifth force appears naturally.
In Barrientos and Mendoza [30], Bernal et al. [31], Mendoza et al. [32], Barrientos and Mendoza [33], Barrientos et al. [34], the authors proposed general gravitational actions with curvature-matter couplings in order to obtain a feasible explanation for the MOdified Newtonian Dynamics (MONDian) behavior of gravitational phenomena. The general conclusions reached by the works of Bernal et al. [31], Mendoza et al. [32], Barrientos et al. [34] is that it is possible to recover the MONDian expression for the acceleration in the regime a < a 0 for a pure metric f (R) theory provided a non-local action construction. Non-locality in these attempts is introduced as an extra scalar field with dimensions of mass that can be conveniently thought of as the causally-connected mass to each point in space-time.
The idea of a non-local gravity and its implications at solar system and cosmological scales have been recently revisited by the works of Mashhoon [35,36], Chicone and Mashhoon [37][38][39], Blome et al. [40]. These theories claim that locality, which allows to treat an accelerated observer as an inertial one at each instant on his proper world-line in special relativity, has its limitations. Since the Einstein Equivalence Principle relates an observer in a local gravitational field with another accelerated one with no gravitation, those limitations are thus extended to gravitational theory. A general characteristic of non-locality is that fields are no longer given by their instant (or local) value but have a contribution attached to the history of the observer (in general terms the Lagrangian density is not localised and as such is not defined as a simple function of the space-time coordinates). Recent investigations have shown that these kind of proposals can simulate dark matter behaviour [41][42][43]. For Tully-Fisher scalings, MOND introduces a gravitational potential for a point mass source proportional to ln(r), which flattens rotation curves. This is also included in these non-local proposals, but they cannot fully explain the proportionality of the MOND force to √ G, where G is Newton's gravitational constant. For the non-local constructions, the relation between the acceleration and the gravitational constant turns out to be linear.
In the present article, we discuss another possible non-local approach in which fractional calculus is used. Fractional calculus, although still a curiosity for many, has proved to describe appropriately non-local space-time effects in a wide range of applications [44][45][46][47][48][49][50]. In general terms, fractional calculus extends the order of differentiation and integration from the natural to the real numbers (in fact to the complex numbers, but we are only going to work with real values in this work). In order to simplify the understanding of fractional calculus to the non-expert, we have included in Appendix A a simple introduction to fractional calculus.
The introduction of fractional derivatives in gravity is not a trivial task since there are several proposals in order to perform the required generalisations. From a pure mathematical point of view, fractional derivatives must induce a somehow fractional geometry in such a way that all the geometric entities involved in general relativity e.g., the connection, the covariant derivative, the Riemann tensor and the metric should be defined in terms of a fractional derivative order. Such an approach is known as the First Step Modification (FSM). A more practical and simple way consist in modifying the Hilbert-Einstein field equations for a given geometry, replacing the covariant derivative order by its analogous fractional derivative without going deeper in how such equation can be obtained. This last approach is usually called a Last Step Modification (LSM). In recent years the applicability of both approaches has been studied in gravity at a classical and and cosmological level [51][52][53][54][55][56][57]. An intermediate approach is to formulate a variational principle for a fractional order action. This latter approach is of particular interest for the scientific community, not exclusively of the gravitational area, since a fractional variational theory is general enough for applications in different scientific fields [58][59][60][61][62][63].
The introduction of fractional derivatives (with unknown derivative orders, which are to be fitted to observations) in the Friedmann equations means that either the gravitational Lagrangian contains fractional derivatives, or the order of variation has a non-integer value (or both). In general terms, if fractional derivatives are to be allowed in physical equations, then the order of differentiation should also be a parameter in the action. In other words, the principle of least action for field equations should somehow contain fractional derivatives. If this is not the case, another possibility is that the variation of the action includes a fractional operator. In a more complicated scenario, both the variation and the Lagrangian could contain fractional derivatives. By itself, this constitutes a very profound subject outside the scope of this work. In this article we force the introduction of fractional derivatives into the standard Friedmann equation of cosmology and see whether we can fit the order of derivatives using SN Ia observations. The goal of this article is to find out whether it would be possible to explain the accelerated expansion of the universe without the introduction of a dark energy component into the cosmological density budget in the Friedmann equation for a dust flat universe as observed today. We deal with the dark matter component by introducing MOND's acceleration constant a 0 into the dynamics of the universe and in principle one would expect that this will force the end result to mean that a non-baryonic dark matter component is not required. The end result is that the fractional derivative order used to generalise the Friedmann equations turns out to be coherent with the reported value of a full MONDian construction of gravity using fractional derivatives. Nevertheless, as it will be discussed in the article, in general terms with the studies performed we can not completely claim the non-necessity of a non-baryonic dark matter, but since MONDian effects are introduced into the problem we expect that component to have a null value.
The article is organised as follows. In section 2 we describe a few key results of standard cosmology and cosmographic parameters. In Section 3 we extend the Friedmann equations and other relations -including cosmographic parameters-to their fractional derivative counterparts. Since the derivative order is a free parameter of the proposal, using these fractional extensions we calibrate all free parameters using SN Ia observations for the accelerated expansion of the Universe in Section 4. In Section 5 we present our statistical results and finally in Section 6 we state final remarks of the toy model developed in this article.

Standard Cosmology
In this section we briefly summarise a few standard concepts of relativistic cosmology that will be extended to their fractional derivatives counterparts later on. Many of the results presented in this section can be found elsewhere (see e.g., [64][65][66][67] and references therein).
The Friedmann-Lemaître-Robertson-Walker (FRLW) metric describes an isotropic and homogeneous Universe, and is given by: where c is the speed of light, a(t) is is the cosmological scale factor as a function of the cosmic time t, r is the radial distance, k is the curvature and θ and ϕ are the polar and azimuthal angles respectively. Substitution of the FLRW metric into the field equations of general relativity with a cosmological constant Λ yields two independent expressions for the time 00 and radial 11 components respectively: for a dust, i.e., pressure-less Universe. In the previous equations ρ represents the matter density. Using the standard definition for the Hubble parameter: Equation (2) turns into: The right hand side of this equation contains all the information for a late-time Universe's constituents: the curvature k, the matter density ρ and the cosmological constant Λ. The previous equation means that: where: represent the matter, dark energy and curvature density parameters respectively. Equation (6) is a convenient normalisation for all the energy constituents of the universe, so that their sum equals one. We now define a few cosmographic parameters. To begin with, the deceleration parameter: can be rewritten as: by means of Equation (3). Derivating with respect to cosmic time t the second Friedmann Equation (3) yields: We now introduce another cosmographic parameter, namely the jerk: In order to express the jerk in terms of the density parameters, we use the fact that the covariant divergence of the energy-momentum tensor vanishes, i.e., ∇ µ T µν = 0, For a matter-dominated dust Universe, it follows that this relation yields: and so ρ ∝ a −3 . Substitution of this last relation into Equation (10) yields: which can be expressed in terms of energy densities only using relation (9): An additional derivative with respect to cosmic time in Equation (10) gives the following expression: The snap cosmographic parameter s is defined as: The time derivative of Equation (12) yieldsρ = −3(ρH + ρḢ) = −3(−3ρH 2 + ρḢ). So, using Equation (15), together with the definitions of H, q, j andḢ = −H 2 (1 + q), the snap parameter can be written as: which in terms of the density parameters is: Equations (9), (14) and (18) can be expressed in terms of the curvature density parameter Ω k [68] instead of the dark energy density parameter Ω Λ using Equation (6). Since the current observations from Planck [69] strongly suggest that we are living in a flat universe k = 0, we will work with this value in what follows. Therefore, Equation (6) simplifies to: 1 = Ω M + Ω Λ and so, the value for the jerk parameter in general relativity has a constant unitary value, i.e., j = 1.

Fractional Friedmann Equation
In what follows, we explore the cosmological consequences of replacing the integer time derivatives in the Friedmann Equations (2) and (3) with fractional derivatives. The idea is to fit the unknown derivative order to SN Ia observations. Following this path, we write down the fractional Friedmann equations as: for a flat Universe. The constant κ, with dimensions of time 2(1−γ) has been introduced into the fractional Friedmann equations in order to have dimensional coherence. The left-hand side of Equation (20) is written as such since in general D γ D γ = D 2γ . Since our target is to work with a Friedmann fractional model with no dark matter, we introduce Milgrom's acceleration constant a 0 as a fundamental physical quantity for the description of gravitational phenomena at cosmological scales. With this and since the velocity of light c and Newton's gravitational constant G are also fundamental, using Buckingham-Π theorem for the dimensional analysis [70], it follows that: where A is a dimensionless constant. Under the idea of fractional orders on the derivative, we can adapt the cosmographic parameters in a natural way as follows. To begin with, we define the fractional Hubble parameter as: Since our intention is to express the Friedmann equation in terms of the density parameter, we follow the procedures of Section 2 and so, dividing Equation (19) by H 2 and using the previous definition it follows that: At this point it is important to mention that the use of the standard definitions for the density parameters in Equation (7) implies that in this extended fractional Friedmann cosmological model, the sum (6) is no longer valid. This essentially occurs because the critical density to "close" a pure matter dominated universe is no longer given by 3H 2 0 /8πG. One can of course redefine the density parameters in such a way that the sum (6) holds. Indeed, if: then: However, in order to avoid more confusion with new extended definitions we decided to keep the standard cosmological definitions of Equation (7) at the cost of breaking up the validity of relation (6).

Matter Dominated Universe
In order to compute the term H /H in the previous expressions, a further assumption about the scale factor a must be made. In standard cosmology, to get how the scale factor evolves as function of t, the different constituents of the Universe are treated separately. In this work, we are going to proceed that way. Thus, from now on we restrict our study to a matter dominated Universe (Ω Λ = 0). For this kind of Universe, the following ansatz is proposed (see e.g., [66,71,72] where a 1 is a constant, and using the rules of fractional derivative for a power law given in Appendix A, the fractional Hubble parameter H is given by: where the exponent γ that appears in t −γ in the previous equation follows standard algebraic rules. Also, the standard Hubble parameter for this scale factor is: Thus, H /H can be written as: Substitution of this last result into Equation (23) yields: where the constant B is defined as: At this point, we have complete freedom over A and so, in order to simplify the Hubble parameter Equation (30), the following choice is made: Therefore, the equation for the Hubble parameter is: In the previous equation, the definitions of the density parameters given in (7) were used. Equation (33) can of course be rewritten in such a way so that Ω M = 1, where Ω M is the right-hand side of Equation (33) divided by H, but as we mentioned before we preferred to stay with the definitions used in standard cosmology. Furthermore, the use of Equation (33) is very convenient since in it, the Hubble parameter H is represented by a simple function of Ω M . A similar relation is not found in standard cosmology and a value for the Hubble constant needs to be known somehow when fitting the standard cosmological model to SN Ia data as shown in the Appendix B.
Also, the use of Equation (33) is quite useful since the definitions of the fractional cosmographic parameters in Equations (34), (39) and (48), involve the Hubble parameter H, making necessary an explicit equation for such parameter. Note that the previous equation is not valid for γ = 1 and since we are going to use that result on many of the rest of the article, in here and in what follows we demand γ = 1.
By defining a fractional deceleration parameter as: where γ in H 2γ is an exponent, the second fractional Friedmann Equation (20) for a matter dominated Universe can be written as: which after dividing by H 2 yields: and so, using Equation (33) and (21) we find: In order to find expressions for the fractional jerk and snap cosmographic parameters as functions of the density parameters only, the natural way would be to follow an analogous procedure as the one described in Section 2. This procedure will involve the cumbersome application of two consecutive fractional cosmic derivatives (one for the jerk and another for the snap) in Equation (20). To avoid that, we apply a Last Step Modification procedure in Equations (10) and (15) with correct definitions for the fractional jerk and snap parameters. For simplicity and coherence, we will continue to use in what follows a Last Step Modification and so, the fractional equivalent of relation (10) is given by: We define the fractional jerk parameter as: and substitute this together with the definitions of the fractional deceleration (34) and Hubble parameter (22) into Equation (38) to obtain: In order to obtain an expression for D γ ρ/Dt γ an equation for ρ must be found. From Equation (33) the following relation is given: With the ansatz a = a 1 t n , the fractional derivative of order γ of the previous equation is: or in terms of the Hubble parameter and the matter density: With this result, Equation (40) takes the following form: or: where Equation (29) was used. Direct substitution of Equations (33) and (37) into the previous relation yields: A fractional analogous of Equation (15) involves the term D γ t D γ t ρ. Since the Leibniz's rule is a complicated relation (see Appendix A), it is easier to derive Equation (42) instead of (43).
The fractional snap parameter is defined by: By performing a similar procedure as the one to obtain the fractional jerk parameter it follows that:

SN Ia Fits
The accelerated expansion of the Universe was first inferred through cosmological observations of SN Ia as standard candles at all redshifts. Perlmutter et al. [73] used 42 high-redshift supernovae to construct an apparent magnitude-redshift diagram in order to obtain the best values for the matter density parameter Ω M with the introduction of a dark energy density parameter Ω Λ component.
The model independent relation between the luminosity distance d L (z) as a function of redshift z and the apparent magnitude µ is given by [34,67]: where H 0 is the Hubble constant evaluated at the present epoch t 0 and h := H 0 /100 km/s/Mpc is the normalised Hubble constant, with the luminosity distance for a flat universe given by [74]: where q 0 , j 0 and s 0 are the cosmographic parameters evaluated at the present epoch. Substitution of Equations (9), (14) and (18) into the previous relation yields an expression for the distance modulus µ as a function of z and Ω M . Since the values for the distance modulus and redshift are given by observations, a statistical fit will return the best estimated values for both density parameters. Since our intention is to see whether the fractional Friedmann model presented above can account for SN Ia observations, we made the assumption that the Last Step Approximation can be employed on Equation (50) and the cosmographic parameters defined in terms of the ordinary time derivatives can be replaced by their fractional analogues, so instead of using Equations (9), (14) and (18) for the cosmographic deceleration, jerk and snap parameters, we use their fractional extension given in (37), (46) and (48). This generalisation requires the Hubble parameter H 0 to be given by Equation (33). The SN Ia distance modulus-redshift data was taken from the Supernova Cosmology Project (SCP) Union 2.1 [75]. The data can be downloaded at the following location: http: //supernova.lbl.gov/Union/figures/SCPUnion2.1_mu_vs_z.txt, or in the cosmology tables section of http://supernova.lbl.gov/Union/.
In order to perform the statistical fits, we used the free software gnuplot (www.gnuplot. info) to obtain the best values for the free parameters of our model. The calibration can be directly performed since the constants c, a 0 and j 0 are known and so the corresponding equations for H 0 , q 0 and s 0 become functions of the density parameters and the fractional order. We used the fit function command in gnuplot for the calibration of the free parameters. This command uses non-linear and linear least squares methods and is able to fit the required function through the empirical data and provides the correlation matrix between the parameters, the number of iterations employed for the converged fit, the final sum of the squares of residuals (SSR), the best fit value for the parameters and its asymptotic standard error, the p-value for the χ-square distribution and the root mean squares of residuals.
At this point it is very important to note that in the definition of the cosmographic parameters q , j and s we used the Hubble parameter H and not H . With the use of the Last Step Modification technique we have no formal way to decide which one to use. We have used H for the simplicity it produces but of course at first sight it seems more natural that the choice H would be the correct choice.
In order to see that the choice H is the only viable one we proceed as follows. The final relation we would like to solve is the redshift-magnitude relation (49). That equation contains the Hubble parameter H(z) through the definition of h. As mentioned on Appendix B, for the case of the standard Λ-CDM cosmology a knowledge of H 0 is required to perform the fit to SN Ia data since there is no independent equation to relate H(z) with some of the parameters of the model. In the case of fractional calculus cosmology studied in this article we could have chosen to use H instead of H on Equation (49), but then we would be required to know by observations the value of H 0 . In order to avoid this it is then best to leave Equation (49) as it is and take advantage of relation (33), i.e., in our model we do not need to know a priori the value of H 0 to perform the fit to SN Ia data.
Also, it is important to note that the definitions of the csomographic parameters q , j and s in Equations (34), (39) and (47) could have been made in terms of H instead of H. But this just means to rewrite the standard cosmological Equations (6), (9), (14) and (18) with their star counterparts. This will not have any departure from the standard ΛCDM cosmological model.

Results
The fractional cosmographic parameters are functions that depend on n (via A), γ and Ω M . Therefore, an extensive search for a coherent fit with the simple gnuplot routine with those three parameters as free variables was made. Setting up the initial values as follows: n = 2.6, γ = 2.1 and Ω M = 4.5, yields quite good convergence for the routine. The obtained results are reported in Table 1. Although the error in n is ∼90%, the mean plotted best distance modulus-redshift function in Figure 1 seems to be in a good agreement with the observations. The envelope curves that represent the statistical errors above and below the solid mean curve in Figure 1 have not been drawn since they turn out to be quite wide due to the somewhat large error in n.
With the mean values reported in Table 1, the Hubble constant H 0 given by Equation (33) has the following numerical value: which is in a great agreement with the value reported by Planck [69] (see Appendix B). Due to limitations of the gnuplot's fit routine and the divergences in the range x < 0 of the Γ(x) function, the convergence of the fit is highly dependent on the initial values of the free parameters n, γ and Ω M . The results shown in Table 1 are obtained using a wide range of initial values. Another fit with a good convergence of the fit is given by the initial values: n = 1.3, γ = 2.6 and Ω M = 2.1. The obtained values are reported in Table 2. The distance modulus-redshift plot for the values obtained in this fit is shown in Figure 2. The envelope curves that represent the statistical errors are present, but due to the small asymptotic errors such envelope is cover by the mean curve.  For this set of initial values the resulting asymptotic error of the free parameters is considerably less than those reported in Table 1. The parameters γ and Ω M do not have a significant improvement, unlike the parameter n whose error drastically decreased. Another relevant difference resides in the p-value, this statistical indicator increased, but is acceptable to within the ranges of gnuplot. With these results, the Hubble constant has the following value: H 0 = 68.37 km/s/Mpc.

Final Remarks
In this article we have generalised the standard Friedmann equations of cosmology using fractional derivatives by performing a Last Step Modification approach. The obtained fractional Friedmann equations and the formulae of the corresponding fractional cosmographic parameters in terms of the cosmological density parameters were used to explain the current accelerated expansion of a dust and zero curvature Universe, without the introduction of a dark energy component. We introduced Milgrom's acceleration constant into the built expressions with the intention of not requiring any non-baryonic dark matter component. The statistical fitting of the free parameters in the model shows an excellent agreement with a fractional derivative order of 3/2, which has recently been shown to be the required fractional order to explain MOdified Newtonian Dynamics (MOND) phenomenology [76]. This result is in excellent agreements with the recent work by Barrientos et al. [34] since the Universe at the present epoch is in the deep MOND regime. Indeed, Milgrom [77] noticed a coincidental relation between the acceleration a 0 and H 0 : a 0 ≈ cH 0 /2π. Since the Hubble radius r H = c/H 0 and the Hubble mass is M H ≈ c 3 /2GH 0 it then follows that the Newtonian gravitational acceleration experienced by a test galaxy at a distance r H of a point mass source M H is ≈ GM H /r 2 H ≈ πa 0 (see e.g., [78,79]). In other words, the Newtonian gravitational acceleration at the present epoch in the Universe is approximately MONDian. As such, one would expect that if no dark matter components are introduced into the cosmological energy budget, then a MONDian description of the Universe at the present epoch is required. The introduction of fractional derivatives into physics usually means that non-locality is taking place at some level and it gets more noticeable at large scales. As noted by Barrientos et al. [34], there are an infinite number of non-local relativistic theories of gravity with curvature-matter couplings that have MOND on their weak field limit of approximation. All this is pointing into the direction that MOND is most probably a non-local phenomenon which occurs at sufficiently small mass to squared length scales, i.e., at acceleration scales a 0 .
Strictly speaking, the fitting procedure is only giving information about the total matter density Ω M which includes baryonic and non-baryonic components. In this sense the model presented shows that a simple introduction of fractional Friedmann equations in cosmology can account for a dark energy component. However, as mentioned in Section 3, the introduction of Milgrom's acceleration constant a 0 in Equation (21) was done with the intention that the dark matter component would be taken into account by a MOND-like prescription. The obtained value of the fractional derivative order of ∼3/2 in this article and the fact that we expect MONDian effects to be relevant in the present epoch dynamics of the Universe, further reinforce this result as mentioned in the previous paragraph.
The present work is restricted to a matter dominated Universe and puts aside the possibility of a dark-energy dominated Universe. The main reason lies in the complexity of the resultant equations. In fact, for a Λ-dominated Universe, Equation (27) is not longer valid. Thinking in analogy to the standard cosmological model, we can expect that the scale factor for a dark-energy dominated era has an exponential dependence on t. However, the fractional derivative of an exponential function is a cumbersome expression that makes the fit extremely difficult. The case of a Universe where both components are taken into account is equally complex since our definitions for the cosmographic parameters require an equation for the Hubble parameter and Equation (23) is a relation for H . In order to have a more solid proposal, a more profound model needs to be built, introducing fractional derivatives in the principle of least action. Such topic goes beyond to the toy model introduced in this work.
The toy model presented in this article about the current expansion of the Universe using fractional derivatives using the Last Step Modification is to be taken with care. It represents a first exploration onto whether there could be any interesting aspects of gravitational theory to be more deeply investigated since the field equations do not come from a variational principle. We intend to go beyond the present work to construct a formalism in terms of fractional derivatives for geometrical objects in general relativity and in formulating a fractional calculus of variations with applications to gravitational theory.

Acknowledgments:
The authors acknowlodge the excellent comments made by three reviewers for useful suggestions that improved the final version of the article. EB thanks support from a Postdoctoral fellowship granted by CONACyT. EB, SM and PP acknowledge economic support from CONACyT (517586, 26344, 14558).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Fractional Calculus: A Simple Introduction
The question of whether it is possible to take half the derivative of a function goes back to a letter written by Leibniz to L'Hôpital in 1665. Since then the subject of fractional calculus, as it is now known, has attracted the interest of mathematicians like Abel, Liouville and Riemann, who developed the basis of the theory. In recent years, fractional derivatives have proved to be useful in many applications ranging from anomalous diffusion, finance to modeling of viscoelastic materials. As expected, when considering the order of differentiation to be an integer, the fractional derivative coincides with the standard definition. However, for non-integer values it has an intrinsic non-local character, which explains its applicability in complex phenomena in which long range interactions or memory effects are present. Before providing the general definition of fractional derivative, we give a few examples motivating this notion. Take for instance Computing the first derivative we obtain The successive repetition of this procedure n times yields the n-th derivative: The above expression makes perfect sense for a real number k, if we can define its factorial in a meaningful way. The Γ function provides such a generalization, and the above expression can be written as where we have replaced n by α to emphasize the fact that now the order of differentiation can be taken to be any real number. A similar reasoning can be applied to other basic examples such as the exponential or trigonometric functions. In such a way, it is possible to define the fractional derivative or integral for functions expressed in their Taylor or Fourier series expansions. However, a more general and useful definition is based on the Riemann-Liouville representation formula of a function. More precisely, the integral operator of order α is defined as Again, this definition can be justified by applying the standard integral operator m times and integrating by parts and then making sense of the definition for a real order of integration. Less intuitive, but more appropriate in formulating initial value problems than other definitions of fractional derivatives is the Caputo proposal given by for m − 1 < α < m, m ∈ N and where f (m) denotes the standard derivative of f . For a systematic presentation of fractional calculus the reader is referred to Podlubny [44].
There are two properties of fractional derivatives that are important to mention since they depar from the common sense in ordinary calculus. The first one is the Leibniz's rule for the product of two functions. Oldham and Spanier [80] proved that the most general expression is given by: The second one is the composition rule: With the intention of showing the robustness of the fitting procedure we are using in the article with the fit routine of gnuplot, we show in this appendix that for the case γ = 1 we recover the values of Ω M and Ω Λ reported in the standard ΛCDM cosmological model. To do so, we cannot simply fix the value γ = 1 on all fractional equations since Equation (30) is not valid for such value. Therefore, the non-fractional Equations (6), (9), (14) and (18) must be employed. Unlike the fractional case, Equation (6) is a constriction and not an equation for the Hubble parameter H 0 . Since the distance modulus µ(z) relation (49) has an explicit dependence on this parameter, a value must be introduced. There are several values reported in the literature for H 0 that constitutes the so called "cosmological tension". We performed fits to SN Ia data using the values reported by the Planck Collaboration [69]: H Planck 0 = 67.36 km s −1 Mpc −1 , a value from SN Ia observations: H SN Ia 0 = 63 km s −1 Mpc −1 (see Figure 9 in [73] and references therein) and a local Hubble constant [81]: H loc 0 = 75.35 km s −1 Mpc −1 . Due to condition given by (6) the number of free parameters is reduced to one. Ω M was chosen to be this free parameter. The results obtained by the fit are shown in Table A1.   Figure A1 shows the Hubble diagram for the ΛCDM cosmological model using the Planck value for the Hubble parameter.