\section{Introduction}
\label{Introduccion}

  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~\citep[see e.g.][]{Will}, where the gravitational field is
moderately weak \citep{Eddington, Pound, Reasenberg, Anderson, Chandler,Ciufolini,Sitter, Eubanks, Nordtvedt, Nordtvedt2, Nordtvedt3, Schiff, Shapiro}.
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 \lesssim 10^{-10} \textrm{m}/\textrm{s}^2 \), or equivalently when the mass-to-length
ratios are much smaller than those of the solar system~\citep{mendoza15},
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 \citep{staro,Nojiri11,Nojiri08,shamir,Odin20,Nor-Odin20,Nor-Odin-04,Nor-Odin-03},
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~\citep{frlm,Harko1,Harko3,Harko4,harko-lobo-book, Bertolami,barrientos18} 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 \citet{bernal11,mendoza13,barrientos16,barrientos18,Barrientos20}, 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~\citet{bernal11,mendoza13,Barrientos20} 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~\citet{Mash1,Mash2,Chicone1,Chicone2,Chicone3, Blome}. 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~\citep{Maggiore1,HehlM, Foffa}.
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 \( \sqrt{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
\citep{Podlubny, Handbook1,Handbook2,Handbook3,Handbook4,Handbook6,Handbook7}.  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~\ref{apen1} 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~\citep{shch1,shch2,shch3,shch4,Roberts,Vacaru,Rami}. 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~\citep{Frederico,Balenu07,Nabulsi,Herzallah,Baleanu05,Nabusi13}.

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~\ref{general-relativity}
we describe a few key results of standard cosmology and cosmographic
parameters.  In Section~\ref{fractional-friedmann} 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~\ref{cosmography}.  In Section~\ref{sub_resul} we
present our statistical results and finally in Section~\ref{discussion} we
state final remarks of the toy model developed in this article.

\section{Standard cosmology}
\label{general-relativity}

  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~\citep[see e.g.][and references therein]{Peacock,Misner,Longair,Peebles}.

The Friedmann--Lema$\hat{\textrm{\i}}$tre--Robertson--Walker (FRLW) metric describes an isotropic and homogeneous Universe,  and is given by :

\begin{equation}
    \mathrm{d}s^2=c^2
     \mathrm{d}t^2-a^2(t)\left(\frac{\mathrm{d}r^2}{1-k
     r^2}+r^2 \mathrm{d} \theta^2+r^2 \sin^2\theta \mathrm{d}\varphi^2
     \right),
    \label{FRLW}
\end{equation}

\noindent 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 \( \theta \) and \( \varphi \) are the polar and azimuthal angles respectively.
Substitution of the FLRW metric into the field equations of general relativity with a cosmological constant  \( \Lambda \)
yields two independent expressions for the time \( 00 \) and radial
\( 11 \) components respectively:

\begin{eqnarray}
\frac{\dot{a}^2+kc^2}{a^2}=\frac{8\pi G\rho +\Lambda c^2}{3},
\label{F1}\\
\frac{\ddot{a}}{a}=-\frac{4\pi G \rho-\Lambda c^2}{3},
\label{F2}
\end{eqnarray}

\noindent for a dust, i.e.~pressure-less Universe. In the previous
equations $\rho$ represents the matter density.
Using the  standard definition for the Hubble parameter:

\begin{equation}
H(t) := \frac{\dot{a}}{a},
\label{Hubble}
\end{equation}

\noindent equation~\eqref{F1} turns into:

\begin{equation}
    H^2=\frac{8\pi G \rho+\Lambda c^2}{3}-\frac{kc^2}{a^2}.
    \label{FriedH}
\end{equation}

\noindent The right hand side of this equation contains all the information for a late-time Universe's constituents: the curvature \( k \), the matter density \( \rho \) and the
cosmological constant \( \Lambda \).
The previous equation means that:

\begin{equation}
    1=\Omega_\text{M}+\Omega_\Lambda+\Omega_k,
\label{suma1}
\end{equation}

\noindent where:

\begin{equation}
    \Omega_\text{M} :=\frac{8\pi G \rho}{3H^2}, \quad \Omega_\Lambda :=\frac{c^2\Lambda}{3H^2},
    \quad \Omega_k := -\frac{kc^2}{H^2a^2},
    \label{densities}
\end{equation}

\noindent represent the matter, dark energy and curvature density parameters
respectively. Equation~\eqref{suma1} 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:

\begin{equation}
    q := -\frac{1}{H^2}\frac{\ddot{a}}{a},
\label{def_q}
\end{equation}

\noindent can be rewritten as:

\begin{equation}
    q=\frac{\Omega_\text{M}}{2}-\Omega_\Lambda,
    \label{deceleration}
\end{equation}

\noindent by means of equation~\eqref{F2}.

  Derivating  with respect to cosmic time \( t \) the second Friedmann
equation~\eqref{F2} yields:

\begin{equation}
    \frac{\dddot{a}}{a}-\frac{\ddot{a}\dot{a}}{a^2}=-\frac{4\pi G \dot{\rho}}{3}.
    \label{derivative1}
\end{equation}

\noindent We now introduce another cosmographic parameter, namely the jerk:

\begin{equation}
  j := \frac{1}{H^3}\frac{\dddot{a}}{a}.
\label{def_j}
\end{equation}

\noindent 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.~$\nabla_\mu T^{\mu\nu}=0$,
For a matter-dominated dust Universe, it follows that this relation yields:

\begin{equation}
    \dot{\rho}=-3H\rho.
    \label{timerho}
    \end{equation}

\noindent and so $\rho\propto a^{-3}$. Substitution of this last relation into
equation~\eqref{derivative1} yields:

\begin{equation}
    j=\frac{3}{2}\Omega_\text{M}-q,
    \label{j_q}
\end{equation}

\noindent which can be expressed in terms of energy densities only using relation~\eqref{deceleration}:

\begin{equation}
    j=\Omega_\text{M}+\Omega_\Lambda
    \label{jerk}
\end{equation}

  An additional derivative with respect to cosmic time in
equation~\eqref{derivative1} gives the following expression:

\begin{equation}
    \frac{\ddddot{a}}{a}-2\frac{\dddot{a}\dot{a}}{a^2}+2\frac{\ddot{a}\dot{a}^2}{a^3}-\left(\frac{\ddot{a}}{a}\right)^2
    =-\frac{4\pi G \ddot{\rho}}{3}.
    \label{derivative2}
\end{equation}

  The snap cosmographic parameter $s$ is defined as:

\begin{equation}
    s := \frac{1}{H^4}\frac{\ddddot{a}}{a}.
    \label{def_s}
\end{equation}

\noindent  The time derivative of equation~\eqref{timerho} yields
$\ddot{\rho}=-3(\dot{\rho}H+\rho\dot{H})=-3(-3\rho H^2+\rho\dot{H})$.
So, using equation~\eqref{derivative2}, together with the definitions of $H$, $q$, $j$
and $\dot{H}=-H^2(1+q)$, the snap parameter can be written as:

\begin{equation}
    s=-\frac{3}{2}\Omega_\text{M} (4+q)+2j+2q+q^2,
    \label{s_q_j}
\end{equation}

\noindent which in terms of the density parameters is:

\begin{equation}
s=-3\Omega_\text{M}-\frac{1}{2}\Omega_\text{M}^2+\frac{1}{2}\Omega_\text{M}\Omega_\Lambda+\Omega_\Lambda^2.
\label{snap}
\end{equation}

Equations \eqref{deceleration}, \eqref{jerk} and \eqref{snap} can be expressed in terms of the curvature
density parameter $\Omega_k$ \citep{Kun} instead of the dark energy density parameter $\Omega_\Lambda$
using equation~\eqref{suma1}. Since the current observations from Planck \citep{planck} strongly suggest that
we are living in a flat universe $k=0$, we will work with this value in what follows.  Therefore, equation~\eqref{suma1} simplifies to:
$1=\Omega_\text{M}+\Omega_\Lambda$ and so, the value for the jerk parameter
in general  relativity has a constant unitary value, i.e.~$j=1$.

\section{Fractional Friedmann equation}
\label{fractional-friedmann}

  In what follows, we explore the cosmological consequences of replacing
the integer time derivatives in the Friedmann equations~\eqref{F1} and~\eqref{F2}
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:

\begin{align}
\left(\frac{D^\gamma a}{Dt^\gamma}\right)^2 &=\kappa a^2 \left(\frac{8\pi G\rho +\Lambda c^2}{3}\right),
\label{fraction1}
\\
\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma a}{Dt^\gamma}\right) &=\kappa a \left(-\frac{4\pi G \rho-\Lambda c^2}{3}\right),
\label{fraction2}
\end{align}

\noindent for a flat Universe.  The constant $\kappa$, with dimensions of
$\textrm{time}^{2(1-\gamma)}$ has been introduced
into the fractional Friedmann equations in order to have dimensional coherence.
The left-hand side of equation~\eqref{fraction2} is written as such since in
general $D^\gamma D^\gamma \neq D^{2\gamma}$.

  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-$\Pi$ theorem for the dimensional  analysis \citep{sedov}, it follows
that:

\begin{equation}
  \kappa= A \left(\frac{a_0}{c}\right)^{2(\gamma-1)}.
\label{kappa}
\end{equation}

 \noindent 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:

\begin{equation}
   H^\star := \frac{1}{a}\frac{D^\gamma a}{Dt^\gamma}.
    \label{fracH}
\end{equation}

Since our intention is to express the Friedmann equation in terms of the density parameter,
we follow the procedures of Section~\ref{general-relativity} and so, dividing
equation~\eqref{fraction1}  by \( H^2 \) and using the previous definition
it follows that\footnote{\label{foot01} It is important to mention that the use
of the standard definitions for the density parameters in equation~\eqref{densities} implies
that in this extended fractional Friedmann cosmological model, the sum~\eqref{suma1} is no longer
valid.  This essentially occurs because the critical density to ``close'' a pure matter dominated
universe is no longer given by \( 3 H_0^2 / 8 \pi G \).  One can of course redefine the density
parameters in such a way that the sum~\eqref{suma1} holds.  Indeed, if:

\begin{align*}
    \Omega_\text{M}^\star &=   A \left(\frac{a_0}{c}\right)^{2(\gamma-1)}\Omega_\text{M}\left(\frac{H}{H^\star}\right)^2 \\
\intertext{and}
    \Omega_\Lambda^\star &=  A \left(\frac{a_0}{c}\right)^{2(\gamma-1)}\Omega_\Lambda\left(\frac{H}{H^\star}\right)^2.
\end{align*}

\noindent then:

\begin{displaymath}
  \Omega_\text{M}^\star + \Omega_\Lambda^\star = 1.
\end{displaymath}

\noindent However, in order to avoid more
confusion with new extended definitions we decided to keep the standard cosmological definitions
of equation~\eqref{densities} at the cost of breaking up the validity of relation~\eqref{suma1}.
}:

\begin{equation}
    \left(\frac{H^\star}{H}\right)^2 = A \left(\frac{a_0}{c}\right)^{2(\gamma-1)}\left(\Omega_\text{M}+\Omega_\Lambda\right).
    \label{intermedio1}
\end{equation}

\subsection{Matter dominated Universe}
In order to compute the term \( H^\star / 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 ($\Omega_\Lambda=0$). For this kind of
Universe, the following ansatz is proposed~\citep[see e.g.][]{Longair,liddle,dodelson}:

\begin{equation}
  a=a_1 t^n,
  \label{ansatz}
\end{equation}

\noindent where
\( a_1 \) is a constant, and
using the rules of fractional derivative for a power law given in Appendix~\ref{apen1}, the fractional Hubble
parameter \(  H^\star \)
is given by:

\begin{equation}
    H^\star=\frac{\Gamma(n+1)}{\Gamma(n+1-\gamma)}t^{-\gamma},
    \label{Hstar}
\end{equation}

\noindent where the exponent $\gamma$ that appears in
$t^{-\gamma}$ in the previous equation follows standard algebraic rules.
Also, the standard Hubble parameter for this scale factor
is: $H=nt^{-1}$. Thus, $H^\star/H$ can be written as:

\begin{equation}
    \frac{H^\star}{H}=\frac{\Gamma(n)}{\Gamma(n+1-\gamma)}t^{1-\gamma}
    =\frac{\Gamma(n)}{\Gamma(n+1-\gamma)}\left(\frac{H}{n}\right)^{\gamma-1}.
    \label{intermedio2}
\end{equation}

\noindent Substitution of this last result into
equation~\eqref{intermedio1}  yields:

\begin{equation}
    H=B\frac{a_0}{c}\left(\Omega_\text{M}\right)^{1/2(\gamma-1)},
    \label{intermedio3}
\end{equation}

\noindent where the constant $B$ is defined as:

\begin{equation}
    B:=\left[\frac{\Gamma(n+1-\gamma)}{\Gamma(n)}\right]^{1/(\gamma-1)} n
    A^{1/2(\gamma-1)}.
    \label{Bconstant}
\end{equation}

 At this point, we have complete freedom over $A$ and so,
in order to simplify the Hubble parameter
equation~\eqref{intermedio3}, the following choice is made:

\begin{equation}
    A=\left[\frac{\Gamma(n)}{\Gamma(n+1-\gamma)}\right]^2 n^{2(1-\gamma)}.
    \label{A}
\end{equation}

\noindent Therefore, the equation for the Hubble parameter is\footnote{
Equation~\eqref{H_final} can of course be rewritten in such a way so that
\( \Omega_\text{M}^\star = 1 \), where \( \Omega_\text{M}^\star \) is the
right-hand side of equation~\eqref{H_final} divided by \( H \), but as we
mentioned before we preferred to stay with the definitions used in
standard cosmology.  Furthermore, the use of equation~\eqref{H_final} is
very convenient since in it, the Hubble parameter \( H \) is  represented by
a simple  function of \( \Omega_\text{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~\ref{CDM}}:

\begin{equation}
    H=\frac{a_0}{c}\left(\Omega_\text{M}\right)^{1/2(\gamma-1)}.
    \label{H_final}
\end{equation}

\noindent In the previous equation, the definitions of the density parameters given
in~\eqref{densities}  were used.  Nonetheless it is possible to define
suitable new adequate density parameters in order to obtain an equation
similar to~\eqref{suma1}, as mentioned on Footnote~\ref{foot01}.
Also, the use of equuation~\eqref{H_final} is quite useful since
the definitions of the fractional
cosmographic parameters in equations~\eqref{frac_q}, \eqref{fraction_j}
and~\eqref{fraction_s},  involve the Hubble parameter $H$, making necessary
an explicit equation for such parameter. Note that the previous
equation is not valid for \( \gamma = 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 \( \gamma \neq 1 \).

  By defining a fractional deceleration parameter as:

\begin{equation}
    q^\star := -\frac{1}{aH^{2\gamma}}\frac{D^\gamma}{Dt^\gamma}
      \left(\frac{D^\gamma a}{Dt^\gamma}\right),
\label{frac_q}
\end{equation}

\noindent where $\gamma$ in $H^{2\gamma}$ is an exponent,
the second  fractional Friedmann equation~\eqref{fraction2}
for a matter dominated Universe can be
written as:

\begin{equation}
  q^\star H^{2\gamma}=\kappa\frac{4\pi G \rho}{3},
\label{fried_q_frac}
\end{equation}

\noindent which after dividing by \( H^2 \) yields:

\begin{equation}
  q^\star H^{2(\gamma-1)}=\kappa\frac{1}{2}
    \Omega_\text{M},
\label{fried_q_frac2}
\end{equation}

\noindent and so, using equation~\eqref{H_final} and \eqref{kappa} we find:

\begin{equation}
    q^\star = \frac{A}{2}.
\label{q_density}
\end{equation}

  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~\ref{general-relativity}. This procedure will involve
the cumbersome application of two consecutive fractional cosmic derivatives
(one for the jerk and another for the snap)
in equation~\eqref{fraction2}.
To avoid that, we apply a Last Step Modification procedure in
equations~\eqref{derivative1} and~\eqref{derivative2} 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~\eqref{derivative1} is given by:

\begin{eqnarray}
    \frac{1}{a}\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma a}{Dt^\gamma}\right)\right)
    -\frac{1}{a^2}\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma a}{Dt^\gamma}\right)\frac{D^\gamma a}{Dt^\gamma} \nonumber \\
    = -\kappa \frac{4\pi G }{3}\frac{D^\gamma \rho}{Dt^\gamma}.
\label{frac_deriv}
\end{eqnarray}

  We define the fractional jerk parameter as:

\begin{equation}
    j^\star := \frac{1}{aH^{3\gamma}}\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma a}{Dt^\gamma}\right)\right),
\label{fraction_j}
\end{equation}

\noindent and substitute this together with  the definitions of the fractional deceleration~\eqref{frac_q}
and Hubble parameter~\eqref{fracH} into equation~\eqref{frac_deriv} to obtain:

\begin{equation}
    H^{3\gamma}j^\star + q^\star H^{2\gamma}H^\star= -\kappa \frac{4\pi G }{3}\frac{D^\gamma \rho}{Dt^\gamma}.
    \label{aux1.1}
\end{equation}

  In order to obtain an expression for $D^\gamma \rho/Dt^\gamma$
an equation for $\rho$ must be found. From equation \eqref{H_final} the following
relation is given:

\begin{equation}
    \rho=\frac{3}{8\pi G}\left(\frac{c}{a_0}\right)^{2(\gamma-1)}H^{2\gamma}.
    \label{rho}
\end{equation}

  With the ansatz $a=a_1t^n$, the fractional derivative of order $\gamma$
of the previous equation is:

\begin{equation}
    \frac{D^\gamma \rho}{Dt^\rho}=\frac{3}{8\pi G}\left(\frac{c}{a_0}\right)^{2(\gamma-1)}
    n^{2\gamma}\frac{\Gamma(1-2\gamma)}{\Gamma(1-3\gamma)}t^{-3\gamma},
    \label{rho_gamma}
\end{equation}

\noindent or in terms of the Hubble parameter and the matter density:

\begin{equation}
   \frac{D^\gamma \rho}{Dt^\gamma}=\frac{\Gamma(1-2\gamma)}{n^\gamma\Gamma(1-3\gamma)}H^\gamma\rho.
   \label{Drho}
\end{equation}

\noindent With this result, equation~\eqref{aux1.1} takes the following form:

\begin{equation}
    H^{2\gamma}j^\star + q^\star H^\gamma H^\star=
    -A \left(\frac{a_0}{c}\right)^{2(\gamma-1)} \frac{\Gamma(1-2\gamma)}{n^\gamma\Gamma(1-3\gamma)}\frac{4 \pi G}{3}\rho.
    \label{aux1}
\end{equation}

\noindent or:

\begin{eqnarray}
    H^{2(\gamma-1)}\left[j^\star + q^\star \frac{\Gamma(n)}{\Gamma(n+1-\gamma)n^{\gamma-1}}\right] &=& \nonumber \\
    -A \left(\frac{a_0}{c}\right)^{2(\gamma-1)}\frac{\Gamma(1-2\gamma)}{n^\gamma\Gamma(1-3\gamma)}\frac{\Omega_{\text{M}}}{2}.
    \label{aux2}
\end{eqnarray}

\noindent where equation \eqref{intermedio2} was used. Direct substitution of
equations \eqref{H_final} and \eqref{q_density} into the previous relation yields:

\begin{equation}
   j^\star=-\frac{A}{2 n^\gamma}\left[ \frac{\Gamma(1-2\gamma)}{\Gamma(1-3\gamma)}+
   \frac{\Gamma(n)}{\Gamma(n+1-\gamma)n^{-1}} \right].
\label{fraction_jone}
\end{equation}

  A fractional analogous of equation \eqref{derivative2} involves the term $D^\gamma_t D^\gamma_t \rho$.
Since the Leibniz's rule is a complicated relation (see Appendix~\ref{apen1}), it is easier to derive
equation \eqref{rho_gamma} instead of \eqref{Drho}.

  The fractional snap parameter is defined by:

\begin{equation}
    s^\star := \frac{1}{aH^{4\gamma}}\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma}{Dt^\gamma}
    \left(\frac{D^\gamma}{Dt^\gamma}\left(\frac{D^\gamma a}{Dt^\gamma}\right)\right)\right).
\label{fraction_ss}
\end{equation}

\noindent By performing a similar procedure  as the one to obtain the fractional jerk parameter it follows that:

\begin{eqnarray}
    s^\star&=& -A\left[\frac{\Gamma(1-2\gamma)}{2\Gamma(1-4\gamma)n^{2\gamma}}\right. \nonumber\\
    & &\left. + \frac{\Gamma(n)\Gamma(1-2\gamma)}{\Gamma(1-3\gamma)\Gamma(n+1-\gamma)n^{2\gamma-1}}
    -\frac{A}{2}\right].
\label{fraction_s}
\end{eqnarray}

\section{SN~Ia fits}
\label{cosmography}

  The accelerated expansion of the Universe was first inferred through
cosmological observations of SN~Ia as standard candles at all
redshifts. \citet{Perlmutter99} used 42 high-redshift supernovae to construct
an apparent magnitude-redshift diagram
in order to obtain the best values for the matter density parameter $\Omega_\text{M}$
with the introduction of a  dark energy density parameter $\Omega_\Lambda$ component.

  The model independent relation between the luminosity distance
\( d_\text{L}(z) \) as a function of redshift \( z  \) and the apparent magnitude
\( \mu \) is given by \citep{Peebles,Barrientos20}:

\begin{equation}
    \mu(z) = 5 \log \left[\frac{H_0 d_L(z)}{c}\right] - 5 \log h(z) + 42.3856,
    \label{modulus}
\end{equation}

\noindent where $H_0$ is the Hubble constant evaluated at the present epoch
\( t_0 \) and $h := H_0 / 100 \textrm{km/s} / \mathrm{Mpc} $ is the normalised Hubble constant, with the luminosity distance for a flat universe given by~\citep{Visser}:

\begin{equation}
  \begin{split}
    d_L(z) = \frac{c}{H_0} \left[z+\frac{1}{2}(1-q_0)z^2
    -\frac{1}{6} \left(1-q_0-3q_0^2+j_0 \right)z^3\right.
    \\
    \left. +\frac{1}{24} \left(2-2q_0-15q_0^2-15q_0^3+5j_0+10q_0 j_0+s_0\right)z^4 + ...\right],
  \end{split}
\label{distance}
\end{equation}

\noindent where $q_0$, $j_0$ and $s_0$ are the cosmographic parameters evaluated at the present epoch.
Substitution of  equations~\eqref{deceleration}, \eqref{jerk} and~\eqref{snap} into the previous relation yields
an expression for the distance modulus $\mu$ as a function of
$z$ and $\Omega_\text{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~\eqref{distance} and the
cosmographic parameters defined in terms of the ordinary time derivatives
can be replaced by their fractional analogues, so instead of using equations~\eqref{deceleration},
\eqref{jerk} and \eqref{snap} for the cosmographic deceleration, jerk and snap parameters,
we use their fractional extension given in~\eqref{q_density}, \eqref{fraction_jone}
and~ \eqref{fraction_s}. This generalisation requires
the Hubble parameter $H_0$ to be  given  by equation~\eqref{H_final}.
The SN~Ia distance modulus-redshift data
was taken from the Supernova Cosmology Project (SCP)
Union 2.1~\citep{SCP-Union2012}\footnote{ The data can be
downloaded  at the following  location: \url{http://supernova.lbl.gov/Union/figures/SCPUnion2.1_mu_vs_z.txt},
or in the cosmology
tables section of \url{http://supernova.lbl.gov/Union/} }.

  In order to perform the statistical fits, we used the free software
\textit{gnuplot} (\url{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 \textit{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 $\chi$-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^\star \), \( j^\star \) and \( s^\star \) we used the Hubble parameter
\( H \) and not \( H^\star \).  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^\star \) 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~\eqref{modulus}.
That equation contains the Hubble parameter \( H(z) \) through the definition of \( h \).
As mentioned on Appendix~\ref{CDM}, for the case of the standard \( \Lambda \)-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^\star \) instead of \( H \) on equation~\eqref{modulus}, but then we would be required
to know by observations the value of \( H_0^\star \).  In order to avoid this it is then
best to leave equation~\eqref{modulus} as it is and take advantage of relation~\eqref{H_final},
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^\star \), \( j^\star \)  and \( s^\star \) in
equations~\eqref{frac_q},~\eqref{fraction_j} and~\eqref{fraction_ss}
could have been made in terms of \( H^\star \) instead of \( H \).
But this just means to rewrite
the standard cosmological equations~\eqref{suma1},~\eqref{deceleration},~\eqref{jerk}
and~\eqref{snap} with their star counterparts. This will not have any departure
from the standard $\Lambda$CDM
cosmological model.

\section{Results}
\label{sub_resul}

The fractional cosmographic parameters are functions
that depend on $n$ (via $A$), $\gamma$ and \( \Omega_M \).
Therefore, an extensive search for a coherent
fit with the simple \textit{gnuplot} routine with those
three parameters as free variables was made. Setting up the initial
values as follows: $n=2.6$, $\gamma=2.1$ and $\Omega_\text{M}=4.5$,
yields quite good convergence for the routine. The obtained results
are reported in Table~\ref{2_parameters}.

Although the error in \( n \) is \( \sim 90\% \), the mean plotted
best distance modulus-redshift function in Figure~\ref{fig_gamma_free}
seems to be in a good agreement with the observations.

\begin{table}
	\begin{center}
		\begin{tabular}{|c|c|}
		\hline
		$\gamma$ & $1.7254\pm 0.043$\\
		\hline
		$\Omega_\text{M}$  & $11.6273\pm 1.583$ \\
		\hline
		$n$ & $0.9306 \pm 0.846$\\
		\hline
		\end{tabular}
		\hspace{1.0cm}
        \begin{tabular}{|c|c|c|c|}
		 \hline
         & $\gamma$ &  $\Omega_\text{M}$ &n\\
         \hline
         $\gamma$ & 1.000 & &\\
         $\Omega_\text{M}$ & 0.999 & 1.000 &\\
         $n$ & 0.995 & 0.995 & 1.000\\
         \hline
		\end{tabular}
	\end{center}
\caption[Table01]{Left: Best fit results for the fractional derivative model with three
free parameters: the fractional  derivative order \( \gamma \), the
matter density parameter \(\Omega_\text{M} \) and the power for the scale factor $n$
presented with their corresponding errors for the initial values: $n=2.6$, $\gamma=2.1$ and $\Omega_\text{M}=4.5$.
The right	panel shows the correlation matrix for the best
fit values reported. The SSR and the $p$-value for this model are:
\(636.015 \) and  $0.04$
respectively}.
\label{2_parameters}
\end{table}

\begin{figure}
  \includegraphics[scale=0.7]{fit_ng_free_Hstar2.eps}
  \caption[Figure01]{Apparent magnitude \( \mu \) vs.~redshift \(
z \) Hubble diagram from the Union 2.1 SNe~Ia data (dots with their
corresponding error bars) and the best fit from our model. The solid
line represents the distance modulus $\mu(z)$ from the best fit to the
data of the model.}
	\label{fig_gamma_free}
\end{figure}

  The envelope curves that represent the statistical errors above and below the
solid mean curve in Figure~\ref{fig_gamma_free} 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~\ref{2_parameters},
the Hubble constant $H_0$ given by equation \eqref{H_final} has the following numerical value:

\begin{equation}
    H_0=66.95 \,\mbox{km/s/Mpc},
    \label{H0result}
\end{equation}

\noindent which is in a great agreement with the value reported by Planck \citep{planck}
(see Appendix~\ref{CDM}).

Due to limitations of the \textit{gnuplot}'s fit routine and the divergences in the
range $x<0$ of the $\Gamma(x)$ function, the convergence of the fit is highly dependent
on the initial values of the free parameters $n$, $\gamma$ and $\Omega_\text{M}$.
The results shown in Table~\ref{2_parameters} 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$, $\gamma=2.6$ and $\Omega_\text{M}=2.1$. The obtained values are reported
in Table \ref{table_error}. The distance modulus-redshift plot for the values
obtained in this fit is shown in Figure \ref{errors}. 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.

\begin{table}
	\begin{center}
		\begin{tabular}{|c|c|}
		\hline
		$\gamma$ & $1.4937\pm 0.0003$\\
		\hline
		$\Omega_\text{M}$  & $5.4220\pm 0.0243$ \\
		\hline
		$n$ & $0.5539 \pm 0.0046$\\
		\hline
		\end{tabular}
		\hspace{1.0cm}
        \begin{tabular}{|c|c|c|c|}
		 \hline
         & $\gamma$ &  $\Omega_\text{M}$ &n\\
         \hline
         $\gamma$ & 1.000 & &\\
         $\Omega_\text{M}$ & 0.273 & 1.000 &\\
         $n$ & -0.140 & 0.587 & 1.000\\
         \hline
		\end{tabular}
	\end{center}
\caption[Table02]{Left: Best fit results for the fractional derivative model with three
free parameters: the fractional  derivative order \( \gamma \), the
matter density parameter \(\Omega_\text{M} \) and the power for the scale factor $n$
presented with their corresponding errors for the initial values: $n=1.3$, $\gamma=2.6$ and $\Omega_\text{M}=2.1$.
The right	panel shows the correlation matrix for the best
fit values reported. The SSR and the $p$-value for this model are:  \(573.544 \) and  $0.50$
respectively}.
\label{table_error}
\end{table}

\begin{figure}
  \includegraphics[scale=0.7]{fit_error.eps}
  \caption[Figure02]{Apparent magnitude \( \mu \) vs.~redshift \(
z \) Hubble diagram from the Union 2.1 SNe~Ia data (dots with their
corresponding error bars) and the best fit from our model. The solid
blue line represents the distance modulus $\mu(z)$ from the best fit to the
data of the model. The other solid curves represents the envelope of the statistical
errors.}
	\label{errors}
\end{figure}

For this set of initial values the resulting asymptotic error of the free parameters
is considerably less than those reported in Table\ \ref{2_parameters}. The parameters $\gamma$ and
$\Omega_\text{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 \textit{gnuplot}. With these results, the Hubble constant has the following value:

\begin{equation}
    H_0=68.37 \,\mbox{km/s/Mpc}.
    \label{H0result2}
\end{equation}

\section{Final remarks}
\label{discussion}

  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~\citep{andrea}.    This result is in excellent agreements with the recent work
by~\citet{Barrientos20} since the Universe at the present epoch is in the
deep MOND regime\footnote{
  \citet{Milgrom1} noticed a coincidental relation between the acceleration
\( a_0 \) and \( H_0 \): \( a_0 \approx c H_0 / 2 \pi \).  Since the Hubble
radius \( r_H = c / H_0 \) and the Hubble mass is \( M_H \approx c^3 / 2 G H_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 \( \approx G M_H / r_H^2 \approx \pi a_0 \) \citep[see
e.g.][]{Milgrom:2008,Bernal:2011a}.  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 \citet{Barrientos20}, 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 \( \lesssim a_0 \).

  Strictly speaking, the fitting procedure is only giving information about the total
matter density \( \Omega_\textrm{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~\ref{fractional-friedmann}, the introduction of Milgrom's acceleration constant
\( a_0 \) in equation~\eqref{kappa} 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 \( \sim 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 \(\Lambda\)-dominated Universe,
equation~\eqref{ansatz} 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~\eqref{intermedio1} is a relation for $H^\star$.
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.

\section*{Acknowledgements}
The authors acknowlodge the excellent comments made by three reviewers for
useful suggestions that  improved the final version of the article.
This work was supported
by DGAPA-UNAM (IN112019) and CONACyT (CB-2014-01 No.~240512)
grants. EB thanks support from a Postdoctoral fellowship granted by CONACyT.
EB, SM and PP acknowledge economic support from
CONACyT (517586, 26344, 14558).

\bibliographystyle{aipnum4-1}