A fractional matter sector for general relativity

In this work, we construct a fractional matter sector for general relativity. In particular, we propose a suitable fractional anisotropy function relating both the tangential and radial pressure of a spherically symmetric fluid based on the Grünwald–Letnikov fractional derivative. The system is closed by implementing the polytropic equation of state for the radial pressure. We solve the system of integro-differential equations by Euler’s method and explore the behavior of the physical quantities, namely, the normalized density energy, the normalized mass function, and the compactness.


Introduction
The construction of real models for compact stars satisfying Einstein field equations remains a great challenge and, although there exist a lot of conditions we can consider, the simplest one is assuming static and spherically symmetric anisotropic fluids.[1][2][3][4][5][6][7][8][9].In this case, the problem reduces to solving three equations with five unknowns, namely two metric potentials and the three thermodynamic quantities (the energy density, the radial, and the tangential pressure of the fluid).In general, the strategy to close the system is to consider an equation of state relating both the energy density and the radial pressure of the fluid, and certain anisotropy functions.Of course, any relation we could consider encodes our ignorance about the complex microscopic interaction occurring in the interior of the star.In other words, the equation of state and the anisotropy function are macroscopic effective quantities capturing what is really happening in the microscopic realm.Although, in principle, we could be able to construct microscopic a e-mail: econtreras@usfq.edu.ecmodels leading to a particular matter sector of the Einstein field equations, in this work we trace back "quantum effects" with the fractional calculus approach [10] (see also [11,12]).
In Ref. [10] the author presents how the fractional calculus modifies the Einstein equations and how difficult is to obtain solutions given that: i) arriving at the Einstein field equation given a particular parameterization of the metric requires the use of a fractional Leibnitz rule which differs from the classical one by the additions of a combinatory infinite series and ii) the non-locality leads to a set of integro-differential equations which are far from trivial to solve.
The first attempts to extract physics from the fractional Einstein equations can be found in the works by Munkhammar [13] and Vacaru [14] in astrophysics and Roberts [15] in cosmology (for applications in other contexts see [16][17][18][19], for example).It is clear that these works are focused on space-time as a fractional geometric structure.However, our intention here is completely different in the sense that we will consider a spacetime described by the classical Einstein equations and the fractional calculus enters only in the anisotropy of the matter sector.At this point, some comments are in order.First, it is worth mentioning that this strategy avoids the issue of interpretation of a "fractional manifold".Instead, we solve the classical Einstein's equations (which we know clearly how to interpret) sourced by an effective matter sector involving fractional calculus.Second, as claimed in [20], in the last years the interest in the so-called fractional quantum mechanics has increased considerably.In particular, in Ref. [20] the authors studied the implications of the fractional calculus in black hole thermodynamics by modifying the Wheeler-DeWitt equation introducing a fractional derivative in the quantization process in the coordinate representation.In this regard, our strategy of tracing back quantum effects through the fractional calculus approach is well-set (at least formally).
This work is organized as follows.In the next section, we briefly review the Einstein field equations for compact objects in static and spherically symmetric space-times.In section 3, we explain the numerical setup we developed to solve the system of integro-differential equations obtained, we propose the anisotropy function, and show the results.The last section is devoted to the final remarks and conclusions.

Einstein field equations: a brief introduction
Let us consider a spherically symmetric space-time with a line element given in Schwarzschild-like coordinates by, where ν and λ are functions of the radial coordinate only.The metric (2.1) satisfies the Einstein field equations given by, where encodes the matter content with, u µ = (e −ν/2 , 0, 0, 0), (2.4) the four-velocity of the fluid and s µ is defined as with the properties s µ u µ = 0, s µ s µ = −1 (we are assuming geometric units c = G = 1).The metric (2.1), has to satisfy the Einstein field equations (2.2), which are given by ) where primes denote derivative with respect to r.Furthermore, we shall consider that the fluid distribution is surrounded by the Schwarzschild vacuum, namely where M represents the total energy of the system.In order to match the two metrics smoothly on the boundary surface r = r Σ = constant, we require continuity of the first and second fundamental forms across that surface.As a result of this matching we obtain the wellknown result, where the subscript Σ indicates that the quantity is evaluated at the boundary surface.From the radial component of the conservation law, one obtains the generalized Tolman-Oppenheimer-Volkoff equation for anisotropic matter which reads, At this point is clear that the problem has reduced to solving three differential equations, namely, (2.6)-(2.8),for five unknowns, {ν, λ, ρ, P r , P ⊥ }.In this regard, we must implement two conditions to close the system.Alternatively, the system can be reduced to the structure equations in the following way.
Let us parameterize the g rr component of the metric as where m is the mass function.Now, after replacing (2.15) in (2.6) and (2.7) we arrive at from where with measures the anisotropy of the system.Now, the structure equations are (2.16) and (2.18) that must be solved with the following boundary conditions so the problem has been reduced to solve two differential equations with four unknowns, namely, {P r , ρ, m, ∆}.To close the system, in this work, we will provide a polytropic equation of state relating the radial pressure with the energy density (see [21][22][23][24][25][26] and references therein), and the anisotropy function ∆.The next step is to rewrite (2.16) and (2.18) by replacing (2.21) and a new set of dimensionless variables which will led us to the generalized Lane-Emden equations for anisotropic matter.To this end, let us first consider defining the variable ψ by where ρ c denotes the energy density at the center (from now on the subscript c indicates that the variable is evaluated at the center), so that (2.21) now reads where α = P rc /ρ c .Let us now introduce the following dimensionless variables from where (2.24) reads ξ 2 dψ dξ Now, in terms of the new variables, the other structure equation (2.16) reads, Equations (2.26) and (2.27) corresponds to the generalized Lane-Emden equations which can be integrated numerically with the conditions (2.28) after providing some anisotropy function ∆.In this work, we will represent the anisotropy function in terms of the Grünwald-Letnikov derivative that will be defined in the next section.

Grünwald-Letnikov derivative and anisotropy function
The introduction of a fractional derivative to represent the anisotropy function will lead to a set of integrodifferential equations which are far from trivial to be solved.Among all the possibilities, in this work, we will use the Grünwald-Letnikov derivative which is advantageous to the numerical computation given that it is defined in a discrete way.The Grünwald-Letnikov derivative is defined as [27]: ) This definition involves a sum of values of the function f (x) at different points.At this point, a couple of comments are in order.First, it is worth mentioning that the Grünwald-Letnikov derivative can be obtained as a limiting case of the Riemann-Liouville derivative.More precisely, (3.1) converges to (see [28][29][30]) where s = [β] + 1, x > a and Γ is the gamma function.In this regard, the Grünwald-Letnikov derivative can be seen as a discrete approximation of the Riemann-Liouville derivative, and in the limit, as the step size goes to zero, the Grünwald-Letnikov derivative becomes a continuous fractional derivative.Second, for 0 < β < 1 the fractional derivative D β a f (x) reduces to the function when β → 0, namely D 0 a f (x) = f (x) and is proportional to the classical first derivative when β → 1, namely, D 1 a f (x) ∝ df (x)/dx.In this work, we propose the following anisotropy function where C and γ are constants, P rc is the pressure at the center of the fluid distribution and D β a is the Grünwald-Letnikov derivative.It is worth noticing that there are many ways to choose the anisotropy function.However, Eq. (3.3) makes the work as we will explain in what follows.First, note that to avoid divergences in the TOV equation the anisotropy must vanish at ξ = 0. Nevertheless, if we assume that the anisotropy function is proportional to the fractional derivative, this requirement is violated because the fractional derivative coincides with the function Ψ for β = 0 but ψ(ξ = 0) ̸ = 0.In this regard, the factor ξ in (3.3) ensures ∆ = 0 at the center of the star.Second, note that as the fractional derivative could take negative values, the extra term in parentheses in (3.3) ensures the positiveness of ∆ in whole the domain.At this point, it is worth emphasizing that Einstein field equations were obtained by assuming regular partial derivative and the fractional operator enters in the relation between the components of the matter sector.To be more precise, we are considering a regular geometric sector and a fractional matter sector which resembles theories with multi-fractional derivatives introduced in [10] to some extent.Even more, the fractional sector is only valid inside the star as the exterior geometry is described by the Schwarzschild vacuum solution.The construction we are doing here assumes that standard general relativity is valid outside the star (a neutron star observed at a distance where the weak field limit holds) but inside the star, the strong curvature and high energy density require to take into account non-local effects encoded in the fractional derivative entering in the anisotropy function.Probably, a more complete description showing the smooth transition between fractional and regular gravity is possible but it is not the main goal of this work.Now, replacing (3.3) in (2.26), the fractional Lane-Emden equation reads, To compute the approximate solutions for the differential equation, we implemented Euler's method, which can only be used because we are able to find an analytical expression for the derivative of the desired function.A general form of the algorithm can be described by the following recurrence relation where δ is the step size used for the numerical approximation.In our case, the function dψ dt (t n , ψ(t n )) is calculated numerically using all values previously computed.Note that this algorithm stops once the function ψ takes values less than zero, which is desirable in the present work because the function clearly has a root and we are only interested in the values before that root.
In figures 1 and 2, we show the behavior of ψ and η as a function of ξ for the parameters shown in the legend.Note that ψ decreases monotonously as expected.The root indicates the size of the star which increases with n.The numerical data reveals that the radius of the star is similar for the extreme cases (β = 0.1 and β = 0.9) and grows around the middle of the interval of the fractional index, namely β = 0.5.
Another quantity of interest that deserves to be discussed is the compactness of the star which measures how much mass can be packaged in a fixed radius and is defined as the ratio of the total mass of the star and the radius of the star, namely y = M/r Σ . (3.6) In figure 3 we show the compactness parameter as a function of n for different values of the fractional parameter β.First, we note that the solutions have high compactness with respect to neutron stars.Indeed, for the different values of the parameters, we find y ∈ (0.35, 0.39) which is over the average of y ≈ 0.2 for neutron stars.Second, we observe that the compactness of the star decreases as β increases, namely, the compactness decreases as we move from the function to its first derivative.In this regard, the fractional parameter allows controlling the degree of compactness of the star.Finally, as n increases, the compactness decreases except for β = 0. Using the fractional Grünwald-Letnikov derivative we have been able to define a new function form for the anisotropy (3.3), which can be interesting for building compact object models, so in figure 4 we report the behavior of ∆ as a function of the parameters shown in the legend of the same figure.In the interior of the fluid distribution, the local anisotropy is an increasing function for some values of β (as expected) but it modifies its behavior as we approach the surface as β grows.More precisely, for small values of this parameter (0 ≤ β < 0.5) the anisotropy function is monotonically increasing, but for larger values (0.5 ≤ β ≤ 1) modifies its behavior towards the outer regions of the compact relativistic object with the appearance of a local minimum.It is worth mentioning that the appearance of a local minimum in the anisotropy could be related to instabilities in the system.However, a complete study on instabilities requires a complete treatment involving the evolution equations of the system which is out of the scope of this work.Finally, we notice that by increasing the index associated with the radial polytrope (n), in general, the anisotropy increases although its behavior is similar in form.
We would like to conclude this section by estimating the numerical error in the numerical computations.Note that, the Euler method is to first order, which means that the global truncation error is proportional to the step size, δ.In this regard, we can effectively make the error arbitrarily small, with a small enough step size.To show that we are using a small enough δ, we will use the L 1 metric for continuous functions in the interval [0, 1]: In Table 1, we have computed the metrics between the interpolated functions obtained from different values of δ, where the next delta is always half of the previous one, and it can be seen that, at δ = 0.000488, the change obtained by using the next delta is less than 10 −3 .

Final comments and conclusion
In this work, we constructed self-gravitating spheres with fractional anisotropy in the sense that it involves a fractional derivative.For convenience, we worked with the Grünwald-Letnikov derivative which is the discrete version of the Riemann-Liouville one.The anisotropy function was constructed based on two basic requirements it must fulfill, namely: (i) vanish at the center of the star and (ii) positive.In principle, we can construct any anisotropy function satisfying the above requirements but the one we used here made the work.Let us keep in mind that two additional conditions are required to solve Einstein's field equations and our anisotropy serves as one of these essential conditions.This definition for ∆ involves a sum of values for the energy density function ψ taken at different points, so that the anisotropy has a kind of "memory", in the sense that it reminds the previous point as each step h is carried out.This, can rather, be interpreted as a nonlocal equation of state.In this sense, many of the conditions, applied to the scope of General Relativity, that can perfectly be considered to provide the extra necessary information, are given by non-local equations of state.The relevant and well known concept of complexity, introduced by L. Herrera [31] is a recent example.This condition may be regarded as a non-local equation of state, that can be used to obtain non-trivial configurations with zero complexity, somehow similar to the one proposed some years ago in [32].The connections between this new fractional way to describe the anisotropy function and previous non-local studies may be a rich field of study that we reserve for future works.
We implemented Euler's method to integrate the Lane-Emden equation and obtain that the solution is physically acceptable for n ∈ (0.01, 0.25) in the whole range of the fractional parameter β.We observed that out of that interval for n the solution ψ fails to be a decreasing function.Another interesting point that deserves to be mentioned is the behavior of the compactness parameter that is sensitive to the values of the fractional parameter β.In particular, we found that the compactness decreases as β grows so the fractional parameter serves as a controller of the ratio mass/radius of the star.
Before concluding this work, it is worth commenting that, as an alternative to the Grünwald derivative we could consider the discrete version of a Riesz-like operator [30,33] that can be written as where A and B are reals and D β a + and D β a − stand for the right and the left Riemann-Liouville derivative.Another possibility is to take into account the so-called discretization of the bilateral derivative defined in [10].To be more precise, we could construct an extension of the Grünwald-Letnikov derivative containing the discrete version of either the Weyl or Liouville derivatives [34] so that, a combination of both leads to the bilateral derivative in [10] in the infinitesimal limit.In any case, it could be interesting to solve the system of integro-differential equations in this work based on those operators mentioned previously and compare the results with the ones obtained here.However, such a task is out of the scope of this work and we leave this for future developments.
We would like to conclude this work by saying that, in principle, we can construct any anisotropy function based on fractional derivatives such that the solution admits other values of n and/or an increase of the compactness of the star with β.However, these and other issues are left to future developments.

5
AcknowledgementsE.C is supporteb by Poligrant N • 17946.The authors acknowledge the referees for their valuable comments regarding the implementation of the fractional calculus in the context of General Relativity.

Table 1
Table of converging distances as δ decreases