Self-gravitating systems in Extended Gravity

Starting from the weak field limit, we discuss astrophysical applications of Extended Theories of Gravity where higher order curvature invariants and scalar fields are considered by generalizing the Hilbert-Einstein action linear in the Ricci curvature scalar $R$. Results are compared to General Relativity in the hypothesis that Dark Matter contributions to the dynamics can be neglected thanks to modified gravity. In particular, we consider stellar hydrostatic equilibrium, galactic rotation curves, and gravitational lensing. Finally, we discuss the weak field limit in the Jordan and Einstein frames pointing out how effective quantities, as gravitational potentials, transform from one frame to the other and the interpretation of results can completely change accordingly.


Introduction
The today observed Universe appears spatially flat and undergoing an accelerated expansion. Several observational data probe this pictures [1][2][3][4][5][6] but two unrevealed ingredients are needed in order to achieve this dynamical scenario, namely Dark Matter (DM) at galactic and extragalactic scales and Dark Energy Other motivations to modify gravity come from the issue of a full recovering of the Mach principle which leads to assume a varying gravitational coupling. The principle states that the local inertial frame is determined by the averaging process of motion of distant astronomical objects [40]. This fact implies that the gravitational coupling can be scale-dependent and related to some scalar field. As a consequence, the concept of inertia and the Equivalence Principle have to be revised. For example, the Brans-Dicke theory [41] is a preliminary attempt to define an alternative theory of gravity: it takes into account a variable Newton gravitational coupling, whose dynamics is governed by a scalar field non-minimally coupled to the geometry. In such a way, Mach's principle is better implemented [41][42][43].
As already mentioned, corrections to the gravitational Lagrangian were already considered by several authors [13,16,17] soon after the GR was proposed. From a conceptual viewpoint, there are no reason a priori to restrict the gravitational Lagrangian to a linear function of the Ricci scalar, minimally coupled with matter [44]. Since all curvature invariants are at least second order differential, the corrective terms in the field equations will be always at least fourth order. That is why one calls them higher order terms (with respect to GR).
Due to the increasing complexity of the field equations, the main amount of works is related to achieve some formally equivalent theories, where the reduction of the order of field equations can be obtained by considering metric and connections as independent fields [44][45][46][47][48]. In addition, many authors exploited the formal relation to scalar-tensor theories to make some statements about the weak field regime [49]. Moreover other authors discussed a systematic analysis of such theories in the low energy limit [49][50][51][52][53][54][55][56][57][58][59]. In particular, Fourth Order Gravity has been studied in the Newtonian limit (weak field and small velocity) and in the Minkowskian limit (e.g. gravitational waves) [60]. In the former case, one finds modifications of gravitational potential, while in the latter, massive gravitational wave modes come out [61]. However, the weak field limit of such theories have to be tested against realistic self-gravitating structures. Galactic rotation curves, stellar hydrodynamics and gravitational lensing appear natural candidates as test-bed experiments [62][63][64][65].
In this paper, starting from Newtonian and post-Newtonian approximations of Extended Gravity, we match the outcomes with some typical astrophysical structures. The Newtonian limit of the general Fourth Order Gravity is considered by generalizing f (R) models with generic functions containing other curvature invariants, as Ricci square (R αβ R αβ ) and Riemann square (R αβγδ R αβγδ ). The spherically symmetric solutions of metric tensor yet present Yukawa-like behaviors, but, in general, one has two characteristic lengths. The Gauss-Bonnet invariant G = R 2 − 4R µν R µν + R µνλσ R µνλσ , and the Weyl invariant C αβγδ C αβγδ can be reduced to the above cases.
The plan of the review is the following. In Sec. 2 we summarize the basic ingredients of the Newtonian and the post-Newtonian limits [56]. Sec. 3 is devoted to the Newtonian limit of Fourth Order Gravity [57], while, in Sec. 3.1, we discuss the solutions generated by pointlike sources with and without gauge conditions. A brief analysis of the free parameters of the theory is performed.
In Sec. 4.1 we briefly review the classical hydrostatic problem for stellar structures and in Sec. 4.2 we derive the modified Poisson equation in the framework of the Newtonian limit of f (R)-gravity. Then the modified Lané-Emden equation is obtained in Sec. 4.3. Analytic and numerical solutions are discussed [66].
Galactic rotation curves are considered in Sec. 5. The main properties of galactic metric are discussed and, in particular, the violation of the Gauss theorem due to the corrections to the Newtonian potential. In Sec. 5.1, we discuss models for the mass distribution of galactic components. A qualitative and quantitative comparison between the experimental data and the theoretical predictions of rotation curve [65] is obtained for our Galaxy and for NGC 3198 in Sec. 5.2 by fixing numerically the free parameters of the theory. The motion of massive particles is discussed in Sec. 5.3: analogies and differences between general Fourth Order Gravity and R n -gravity are pointed out. Fourth Order Gravity is also considered with and without DM. The outcomes are compared with GR. The rotation curves results higher or lower than the same curves in GR, if the squared Ricci scalar contribution is dominant or not with respect to the contribution of the squared Ricci tensor. However by modifying gravity, also the spatial description of DM could undergoes modifications and the free parameters of model assume different values.
Gravitational lensing is analyzed in the case of pointlike sources (Sec. 6.1) and extended sources (Sec. 6.2). An important outcome is found in the case of thin lenses: the bending of photons in f (R)-gravity is the same as in GR [67,68]. If we add the squared Ricci tensor contribution, we obtain a shift of the image positions (Sec. 6.3).
Finally we take into account the weak field limit of scalar-tensor gravity [69] in the Jordan frame (sec. 7.1) and compare it with the analogous in the Einstein frame (sec. 7.2). Specifically, we consider how physical quantities, like gravitational potentials, derived in the Newtonian approximation for the same scalar-tensor theory, behave in the Jordan and in the Einstein frame. The approach allows to discriminate features that are invariant under conformal transformations and gives contributions in the debate of selecting the true physical frame.The particular case of f (R)-gravity is considered in Sec. 7.3.
Discussions and conclusions are drawn in Sec. 8.

Newtonian and post-Newtonian approximations
It is worth discussing some general issues on the Newtonian and post-Newtonian limits. Basically there are some general features one has to take into account when approaching these limits, whatever the underlying theory of gravitation is.
If one consider a system of gravitationally interacting particles of massM , the kinetic energy 1 2Mv 2 will be, roughly, of the same order of magnitude as the typical potential energy U = GM 2 /r, with M ,r, andv the typical average values of masses, separations, and velocities of these particles. As a consequence:v (for instance, a test particle in a circular orbit of radius r about a central mass M will have velocity v given in Newtonian mechanics by the exact formula v 2 = GM/r.) The post-Newtonian approximation can be described as a method for obtaining the motion of the system to higher than first order (approximation which coincides with the Newtonian mechanics) with respect to the quantities GM /r andv 2 assumed small with respect to the squared light speed c 2 . This approximation is sometimes referred to as an expansion in inverse powers of the light speed.
The typical values of the Newtonian gravitational potential U are nowhere larger than 10 −5 in the Solar System (in geometrized units, U/c 2 is dimensionless). On the other hand, planetary velocities satisfy the conditionv 2 U , while 1 the matter pressure p experienced inside the Sun and the planets is generally smaller than the matter gravitational energy density ρU , in other words 2 p/ρ U . Furthermore one must consider that even other forms of energy in the Solar System (compressional energy, radiation, thermal energy, etc.) have small intensities and the specific energy density Π (the ratio of the energy density to the rest-mass density) is related to U by Π U (Π is ∼ 10 −5 in the Sun and ∼ 10 −9 in the Earth [70]). As matter of fact, one can consider that these quantities, as function of the velocity, give second order contributions Therefore, the velocity v gives O(1) terms in the velocity expansions, U 2 is of order O(4), U v of O(3), U Π is of O(4), and so on. Considering these approximations, one has and Now, particles move along geodesics where ds = g αβ dx α dx β is the relativistic distance. The equation (5) can be written in details as where the set of coordinates adopted is x µ = (t, x 1 , x 2 , x 3 ). In the Newtonian approximation, that is vanishingly small velocities and only first-order terms in the difference between g µν and the Minkowski metric η µν , one obtains that the particle motion equations reduce to the standard result The quantity 1 − g tt is of order GM /r, so that the Newtonian approximation gives d 2 x i /dt 2 to the order GM /r 2 , that is, to the orderv 2 /r. As a consequence if we would like to search for the post-Newtonian approximation, we need to compute d 2 x i /dt 2 to the orderv 4 /r. Due to the Equivalence Principle and the differentiability of spacetime manifold, we expect that it should be possible to find out a coordinate system in which the metric tensor is nearly equal to the Minkowski one η µν , the correction being 1 We consider here on the velocity v in units of the light speed c. 2 Typical values of p/ρ are ∼ 10 −5 in the Sun and ∼ 10 −10 in the Earth [70]. expandable in powers of GM /r ∼v 2 . In other words one has to consider the metric developed as follows 3 g tt (t, x) 1 + g (2) tt (t, x) + g (4) tt (t, x) + O(6) where δ ij is the Kronecker delta, and for the controvariant form of g µν , one has In evaluating Γ µ αβ we must take into account that the scale of distance and time, in our systems, are respectively set byr andr/v, thus the space and time derivatives should be regarded as being of order Using the above approximations (8), (9), we have, from the definition of Christoffel symbols Γ µ αβ = 1 2 g µσ (g ασ,β + g βσ,α − g αβ,σ ), The Ricci tensor component are 3 The Greek index runs from 0 to 3; the Latin index runs from 1 to 3.

R
(2) tt,n − 1 2 g (2) mn g (2) tt,mn with and , respectively, the Laplacian and the gradient in flat space. The inverse of the metric tensor is defined by means of the equation (13) with δ α β the Kronecker delta. The relations among the higher than first order terms turn out to be Finally the Lagrangian of a particle in presence of a gravitational field can be expressed as proportional to the invariant distance ds 1/2 , thus we have : which, to the O(2) order, reduces to the classic Newtonian Lagrangian of a test particle L New = (1 + g As matter of fact, post-Newtonian physics has to involve higher than O(4) order terms in the Lagrangian.
An important remark concerns the odd-order perturbation terms O(1) or O(3). Since, these terms contain odd powers of velocity v or of time derivatives, they are related to the energy dissipation or absorption by the system. Nevertheless, the mass-energy conservations prevents the energy and mass losses and, as a consequence, prevents, in the Newtonian limit, terms of O(1) and O(3) orders in the Lagrangian. If one takes into account contributions higher than O(4) order, different theories give different predictions. GR, for example, due to the conservation of post-Newtonian energy, forbids terms of O(5) order; on the other hand, terms of O(7) order can appear and are related to the energy lost by means of the gravitational radiation.
On the matter side, i.e. right-hand side of the field equations we start with the general definition of the energy-momentum tensor of a perfect fluid with additional energy Π T µν = (ρ + Πρ + p)u µ u ν − p g µν (16) The explicit form of the energy-momentum tensor can be derived as follows where ρ is the density of energy.
3. The Newtonian limit of f (X, Y, Z)-gravity Let us start with a general class of fourth order theories given by the action where f is an unspecified function of curvature invariants X ≡ R, Y ≡ R αβ R αβ and Z ≡ R αβγδ R αβγδ . The term L m is the minimally coupled ordinary matter contribution. In the metric approach, the field equations are obtained by varying (18) with respect to g µν . We get The convention for Ricci's tensor is R µν = R σ µσν , while for the Riemann tensor is R α βµν = Γ α βν,µ + .... The affinities are the usual Christoffel symbols of the metric. The adopted signature is (+ − −−) (for details, see [71]). The trace of field equations (19) is the following Some authors considered a linear Lagrangian containing not only X, Y and Z but also the first power of curvature invariants R and R αβ ;αβ . Such a choice is justified because all curvature invariants have the same dimension (L −2 ) [72]. Furthermore, this dependence on these two last invariants is only formal, since from the contracted Bianchi identity (2R αβ ;αβ − R = 0) we have only one independent invariant. In any linear theory of gravity (the function f is linear) the terms R and R αβ ;αβ give us no contribution to the field equations, because they are four-divergences. However if we consider a function of R or R αβ ;αβ by varying the action we still could have the four-divergences but we would have the contributions of sixth order differential. For our aim (we want to stay in framework of fourth order differential equations) the most general theory of fourth order is the action (18).
The Newtonian limit analysis starts from the development (8) of the metric tensor. Then we set the metric as where we note that the presence of function g (2) ij (t, x) is superfluous here for our aim but it will be fundamental about the gravitational lensing. The curvature invariants X, Y , Z become The function f can be developed as and analogous relations for partial derivatives of f are obtained. From the interpretation of stress-energy tensor components (17) as energy density, momentum density and momentum flux, we have T tt , T ti and T ij at the various orders where T

(N )
µν denotes the term in T µν of orderM /r 3vN . In particular T tt is the density of rest-mass, while T (2) tt is the non-relativistic part of the energy density. From the lowest order of field equations (19) we have Not only in f (R)-gravity [55,56,73] but also in f (X, Y, Z)-theory a missing cosmological component in the action (18) implies that the space-time is asymptotically Minkowskian. The equations (19) and (20) at O(2) -order become 4 By introducing the quantities we get three differential equations for curvature invariant X (2) , ttand ij-component of Ricci tensor R The solution for curvature invariant X (2) in third line of (28) is where G 1 (x, x ) is the Green function of the field operator − m 1 2 . The solution for g (2) tt , by remembering R tt , is 4 We used the properties: 2R αβ ;αβ − R = 0 and R µ αβ where G 2 (x, x ) is the Green function of the field operator − m 2 2 .
The expression (30) is the "modified" gravitational potential (here we have a factor 2) for f (X, Y, Z)-gravity. The solution for the gravitational potential Φ = g (2) tt /2 has a Yukawa-like behaviors ( [55,57]) depending by a characteristic lengths on whose it evolves.
The ij-component of Ricci tensor in terms of metric tensor (12) can be expressed by using the harmonic gauge condition (g αβ Γ µ αβ = 0) as R ij . Then the general solution for g (2) ij from (28), in the harmonic gauge, is While if we hypothesize g ,ij and the second field equation of (28) becomes Then the general solution for g (2) ij from (28) is and the second line of (32) is mathematically satisfied. If we consider the trace of the second line of the set (32) we have a mathematical constraint for the gravitational potentials Φ, Ψ and we can affirm that only in GR the metric potentials are equals (or more generally their difference must be proportional to function |x| −1 ). 5 We choose a system of isotropic coordinates. Generally the set of coordinates (t, r, θ, φ) are called standard coordinates if the metric is expressed as ds 2 = g tt (t, r) dt 2 +g rr (t, r)dr 2 −r 2 dΩ while if one has ds 2 = g tt (t,

The pointlike solution
Let us consider a pointlike source with mass M . The energy-momentum tensor (17) is (we are not interesting to the internal structure) where ρ(x) is the mass density and u µ satisfies the condition g tt u t 2 = 1, u i = 0. Since (21) the expression (35) becomes The energy-momentum tensor satisfies the relations T where µ i . = |m i 2 |. The first choice in (37) corresponds to Yukawa-like behavior, while the second one to the oscillating case. However, the Green function for the negative squared mass is, in general, a linear combination of Cos and Sin with coefficients associated to the shift phase. For the sake of simplicity, we assumed here unitary coefficients. Both expressions are a generalization and/or an induced correction of the usual gravitational potential (∝ |x| −1 ), and when m i (27)) we recover the field equations of GR. Independently of algebraic sign of m i 2 , one can introduce two scale lengths µ i −1 . We note that in the case of f (R)-gravity, we obtain only one scale length (µ 1 −1 with f Y (0) = f Z (0) = 0) on the which the Ricci scalar evolves [55][56][57], but in f (X, Y, Z)-gravity we have an additional scale length µ 2 −1 on the which the Ricci tensor evolves.
If we choose m i 2 > 0 the curvature invariant X (2) (29) and the metric potentials Φ (30) and Ψ (33) are where r g = 2GM is the Schwarzschild radius and the relativistic invariant is The modified gravitational potential by f (R)-gravity is further modified by the presence of functions of R αβ R αβ and R αβγδ R αβγδ . The curvature invariant X (2) (the Ricci scalar) presents a massive propagation and when f (X, Y, Z) → f (R) we find the mass definition and propagation mode with m 2 disappear [55][56][57]. Exponential and oscillating behaviors of solutions (37) are compatible with respect to the previous obtained outcomes in the Newtonian limit of f (R)-gravity [50,51] . Indeed, in the case of a genuine f (R)-gravity also the mass definition (40) coincides again with the previous results. The choice of parameters (Yukawa-like case for the Green function) is made also in [72] and the constraint conditions on the derivatives of f obtained here (f Y (0)+4f Z (0) > 0 and 3f XX (0)+2f Y (0)+ 2f Z (0) < 0) are still compatible with respect to paper mentioned. Obviously the expressions Φ and Ψ in (38) satisfy the constraint condition (34).
The gravitational potential Φ, solution of Eqs. (59), has in general a Yukawa-like behavior depending on a characteristic length on which it evolves [55,76]. Then as it is evident the Gauss theorem is not valid 6 since the force law is not ∝ |x| −2 . The equivalence between a spherically symmetric distribution and point-like distribution is not valid and how the matter is distributed in the space is very important [55][56][57][58]. However the field equations (28) are linear and we can use the superposition principle. If we have a generic density mass ρ(x) the potential Φ becomes and in the case of a ball source with radius ξ the potential Φ outside the source is 6 It is worth noticing that also if the Gauss theorem does not hold, the Bianchi identities are always valid so the conservation laws are guaranteed.
is a geometric factor depending on the form of source and satisfying the condition lim x→0 F (x) = 1.
Besides the Birkhoff theorem results modified at Newtonian level: the solution can be only factorized by a space-depending function and an arbitrary time-depending function [56,76]. Furthermore the correction to the gravitational potential is depending on the only first two derivatives of f with rispect the X and the first derivatives with rispect the Y and Z in the point X = Y = Z = 0. This means that different analytical theories, from the third derivative perturbation terms on, admit the same Newtonian limit [55][56][57].
For a right physical interpretation of Φ in (38) we also impose the condition µ 1 − 4µ 2 < 0. This is a necessary condition for a correct interpretation of attractive gravitational potential, in addition to the condition µ 1 < µ 2 for a well with a negative minimum in |x| = 0. These conditions can be summarized in terms of scale lengths. In fact we find that the contribution induced by Ricci and Riemann square goes to zero faster than the contribution induced by Ricci scalar: In terms of Lagrangian it means that for any f (X, Y, Z)-gravity the correction to Hilbert-Einstein Lagrangian depends more on the contributions of Ricci scalar (f XX (0) = 0) than ones of Ricci and Riemann square ( Then, if we suppose f X (0) > 0 (for a right Newtonian limit), starting from the condition µ 1 < µ 2 we get a constraint on the derivatives of f with respect to curvature invariants In the case of f (R)-gravity (f Y (0) = f Z (0) = 0) we reobtain the same condition among the first and second derivatives of f [55][56][57][58].

f (X, Y, Z)-gravity and the quadratic Lagrangian
The outcome (38) can be obtained also by considering the so-called where a 1 , a 2 and a 3 are constants. In this case, [58] we find two characteristic lengths 2(3a 2 +a 3 ) 1/2 and the Newtonian limit of theory implies as solution Eqs. (38). We can state, then, the Newtonian limit of any f (X, Y, Z)-gravity can be reinterpreted by introducing the Quadratic Lagrangian and the coefficients have to satisfy the following relations A first considerations about (44) is regarding the characteristic lengths induced by f (X, Y, Z)-gravity. The second length µ 2 −1 is originated from the presence, in the Lagrangian, of Ricci and Riemann tensor square, but also a theory containing only Ricci tensor square could show the same outcome (it is successful replacing the coefficients a i of Quadratic Lagrangian or renaming the function f ). Obviously the same is valid also with the Riemann tensor square alone. Then a such modification of theory enables a massive propagation of Ricci Tensor and, as it is well known in the literature, a substitution of Ricci Scalar with any function of Ricci scalar enables a massive propagation of Ricci scalar. We can, then, affirm that an hypothesis of Lagrangian containing any function of only Ricci scalar and Ricci tensor square is not restrictive and only the experimental constraints can fix the arbitrary parameters. A second consideration is starting from the Gauss -Bonnet invariant defined by the relation G GB = X 2 − 4Y + Z [77]. In fact the induced field equations satisfy in four dimensions the following condition and by substituting them at Newtonian level ( tt ) in the equations (19) we find the field equations (ever at Newtonian Level) of Quadratic Lagrangian which can be recasted as The theory of gravity represented by the Lagrangian (46) is the more general theory considering all invariant curvatures, but we have to take into account a degeneracy. In fact we have different f (X, Y, Z)-gravity describing the same Newtonian Limit (see [55][56][57][58] for details). The solution of field equations or the observable quantities, as the galactic rotation curve, are parameterized by the derivatives of f , then we have different functions f (X, Y, Z) which admit the same physics in the weak field limit.

Stellar hydrostatic equilibrium in f (R)-gravity
Extended Theories of Gravity could have important applications at stellar level. Here we outline how stellar structure equations result modified if the underlying theory of gravity is not GR.

The Newtonian approach to the hydrostatic equilibrium
The condition of hydrostatic equilibrium for stellar structures in Newtonian dynamics is achieved by considering the equation [78] Together with the above equation, the Poisson equation gives the gravitational potential Φ as solution for a given matter density ρ. Since we are taking into account only static and stationary situations, here we consider only time-independent solutions 7 . In 7 The radius r is assumed as the spatial coordinate. It ranges from r = 0 at the stellar center to r = ξ at the surface of the star.
general, the temperature τ appears in Eqs. (47) and (48) the density satisfies an equation of state of the form ρ = ρ(p, τ ). In any case, we assume that there exists a polytropic relation between p and ρ of the form where K and γ are constant. The polytropic constant K is fixed and can be obtained as a combination of fundamental constants. However there are several realistic cases where K is not fixed and another equation for its evolution is needed. The constant γ is the polytropic exponent . Inserting the polytropic equation of state into Eq. (47), we obtain For γ = 1, the above equation can be integrated giving where we have chosen the integration constant to give Φ = 0 at surface (ρ = 0). The constant n is called the polytropic index and is defined as n = 1 γ−1 . Inserting the relation (51) into the Poisson equation, we obtain a differential equation for the gravitational potential Let us define now the dimensionless variables where the subscript c refers to the center of the star and the relation between ρ and Φ is given by Eq. (51). At the center (r = 0), we have z = 0, Φ = Φ c , ρ = ρ c and therefore w = 1. Then Eq. (52) can be written This is the standard Lané-Embden equation describing the hydrostatic equilibrium of stellar structures in the Newtonian theory [78].

Hydrostatic equilibrium in f (R)-gravity
Let us consider now how hydrostatic equilibrium for stellar structures can change extending the theory to f (R)-gravity. In this case the action (18) becomes The field equations are obtained by imposing f (X, Y, Z) → f (R) in the equations (19) and (20). We get where f = f X . Since we are interested in analyzing the modification of Lané-Embden equation, it is sufficient to solve only the field equation for the gravitational potential Φ. In order to achieve the Newtonian limit of the theory, the metric tensor g µν (21) can be approximated as follows and the field equations (26) become where for the sake of simplicity we set f (0) = 1. By the definition of mass (40) and R We note that for f (0) → 0 (or m → ∞) we have the standard Poisson equation: Φ = 4πGρ. This means that as soon as the second derivative of f is different from zero, deviations from the Newtonian limit of GR emerge. The equations (59) can be considered the modified Poisson equation for f (R)-gravity. They do not depend on gauge condition choice [58]. From the Bianchi identity, satisfied by the field equations (19), we have If the dependence on the temperature τ is negligible, i.e. ρ = ρ(p), this relation can be introduced into Eqs. (59), which become a system of three equations for p, Φ and R (2) and can be solved without the other structure equations. Let us suppose that matter satisfies still a polytropic equation p = K ρ γ . If we introduce the relation (51) into field equations (59) we obtain an integro-differential equation for the gravitational potential Φ, that is The integro-differential nature of Eq. (61) is the proof of the non-viability of Gauss theorem for f (R)-gravity. Adopting again the dimensionless variables is a characteristic length linked to stellar radius ξ, Eq. (61) becomes which is the modified Lané-Emden equation deduced from f (R)-gravity. Clearly the particular f (R)-model is specified by the parameters m and ξ 0 . If m → ∞ (i.e. f (R) → R), Eq. (64) becomes Eq. (54). We are only interested in solutions of Eq. (64) that are finite at the center, that is for z = 0. Since the center must be an equilibrium point, the gravitational acceleration |g| = dΦ/dr ∝ dw/dz must vanish for w (0) = 0. Let us assume we have solutions w(z) of Eq.(64) that fulfill the boundary conditions w(0) = 1 and w(ξ/ξ 0 ) = 0; then according to the choice (62), the radial distribution of density is given by and the pressure by For γ = 1 (or n = ∞) the integro-differential Eq. (64) is not correct. This means that the theory does not contain the case of isothermal sphere of ideal gas. In this case, the polytropic relation is p = K ρ. Putting this relation into Eq.(60) we have where the constant of integration is chosen in such a way that the gravitational potential is zero at the center. If we introduce Eq.(67) into Eqs. (59), we have Assuming the dimensionless variables z = |x| which is the modified "isothermal" Lané-Emden equation derived f (R)-gravity.

Solutions for the standard and modified Lané-Emden Equations
The task is now to solve the modified Lané-Emden equation and compare its solutions to those of standard Newtonian theory. Only for three values of n, the solutions of Eq.(54) have analytical expressions [78] We label these solution with GR since they agree with the Newtonian limit of GR. The surface of the polytrope of index n is defined by the value z = z (n) , where ρ = 0 and thus w = 0. For n = 0 and n = 1 the surface is reached for a finite value of z (n) . The case n = 5 yields a model of infinite radius. It can be shown that for n < 5 the radius of polytropic models is finite; for n > 5 they have infinite radius. From Eqs.(70) one finds z A general property of the solutions is that z (n) grows monotonically with the polytropic index n. In Fig. 1 we show the behavior of solutions w (n) GR for n = 0, 1, 5. Apart from the three cases where analytic solutions are known, the classical Lané-Emden Eq. (54) has to be be solved numerically, considering with the expression w (n) for the neighborhood of the center. Inserting Eq.(71) into Eq. (54) and by comparing coefficients one finds, at lowest orders, a classification of solutions by the index n, that is The case γ = 5/3 and n = 3/2 is the non-relativistic limit while the case γ = 4/3 and n = 3 is the relativistic limit of a completely degenerate gas. Also for modified Lané-Emden Eq. (64), we have an exact solution for n = 0. In fact, it is straightforward to find out where the boundary conditions w(0) = 1 and w (0) = 0 are satisfied. A comment on the GR limit (73) is necessary. In fact when we perform the limit m → ∞, we do not recover exactly w , since in the first equation of (59), when we perform f (R) → R, we have to eliminate the trace equation condition. In general, this means that the Newtonian limit and the Eddington parameterization of different relativistic theories of gravity cannot coincide with those of GR (see [16] for further details on this point).
The point Since the stellar radius ξ is given by definition By solving numerically the constraint 8 Eq. (75), we find the modified expression of the radius ξ. If m → ∞ we have the standard expression ξ = 3Φc 2πG valid for the Newtonian limit of GR. Besides, it is worth noticing that in the f (R)-gravity case, for n = 0, the radius is smaller than in GR. On the other hand, the gravitational potential −Φ gives rise to a deeper potential well than the corresponding Newtonian potential derived from GR [58].
In the case n = 1, Eq. (64) can be recast as follows wherew = z w. If we consider the solution of (76) as a small perturbation to the one of GR, we havẽ 8 In principle, there is a solution for any value of m.
The coefficient e −mξ < 1 is the parameter with respect to which we perturb Eq. (76). Besides these position ensure us that when m → ∞ the solution converge to something likew (1) GR (z). Substituting Eq. (77) in Eq.(76), we have and the solution is easily found Also in this case, for m → ∞, we do not recover exactly w GR (z). The reason is the same of previous n = 0 case [16]. Analytical solutions for other values of n are not available.
To conclude this section, we report the gravitational potential profile generated by a spherically symmetric source of uniform mass with radius ξ. By solving field Eqs. (59) inside the star and considering the boundary conditions w(0) = 1 and w (0) = 0, we get In the limit m → ∞, we recover the GR case w GR (z) = 1 − ξ 0 2 z 2 3ξ 2 . In Fig. 1 we show the behaviors of w (0) f (R) (z) and w (1) f (R) (z) with respect to the corresponding GR cases. Furthermore, we plot the potential generated by a uniform spherically symmetric mass distribution in GR and f (R)-gravity and the case w An important remark is in order at this point. In a general discussion of the problem, also time-dependent perturbations should be taken into account. In particular, quadratic gravity possess gosth-like instability related to the sign of coefficient in front of the term R 2 . In such cases, unstable or unphysical solutions come out and have to be discriminated with respect to physical cases. A discussion in this sense is in [38,39].

Rotation curves of galaxies
Another important astrophysical application of the above Newtonian limit of Fourth Order Gravity is related to the rotation curves of spiral galaxies. The approach can be considered also an alternative test for DM behavior.
In general, the motion of body embedded in the gravitational field is given by geodesic equation (5) and in the Newtonian limit, formally, we obtain the classical structure of the motion equation (7) The blue dashed-dotted line is the potential derived from GR (w GR (z)) and the red dashed-dotted line the potential derived from f (R)-gravity (w f (R) (z)) for a uniform spherically symmetric mass distribution. The assumed values are mξ = 1 and mξ 0 = .4. From a rapid inspection of these plots, the differences between GR and f (R) gravitational potentials are clear and the tendency is that at larger radius z they become more evident.
but the gravitational potential is given by (38) or generally by (41). The study of motion is very simple if we consider a particular symmetry of mass distribution ρ, otherwise the analytical solutions are not available. Our aim is to evaluate the corrections to the classical motion in the easiest situation: the circular motion. In this case we do not consider the radial and vertical motion. The condition of stationary motion on the circular orbit is where v c is the velocity. The distribution of mass can be modeled simply by introducing two sets of coordinates: the spherical coordinates (r, θ, φ) and the cylindrical coordinates (R, θ, z). An useful mathematical tool is the Gauss flux theorem for gravity: The gravitational flux through any closed surface is proportional to the enclosed mass. The law is expressed in terms of the gravitational field. The gravitational field g is defined so that the gravitational force experienced by a particle with mass m is F grav = m g. Since the Newtonian mechanics satisfies this theorem and, by thinking to a spherical system of mass distribution, we get, from (82), the equation where M (r) is the only mass enclosed in the sphere with radius r. The Green function of the f (X, Y, Z)-gravity ( = |x − x | −1 ), instead, does not satisfy the theorem [58]. In this case we must consider directly the gravitational potential (41) generated by using the superposition principle.
Apart the mathematical difficulties incoming from the research of gravitational potential for a given mass distribution, the non-validity of Gauss theorem implies, for example, that a sphere can not be reduced to a point (see solution (42)). In fact the gravitational potential generated by a ball (also with constant density) is depending also on the Fourier transform of ball [58]. However our aim is to consider not the simple case of motion of body in the vacuum, but the more interesting case of motion in the matter. So we must leave any possibility of idealization and consider directly the calculation of the potential (41).
Two remarks on the (41) are needed. The two corrections have different algebraic sign, and in particular the Yukawa correction with µ 1 implies a stronger gravitational force, while the second one (purely induced by Ricci and Riemann square) contributes with a repulsive force. By remembering that the motivations of extending the outcome of GR to new theories is supported by missing matter justifying the flat rotation curves of galaxies, the first correction is a nice candidate. A crucial point is given by the spatial range of correction. In fact the Yukawa corrections imply a massive propagation; then, more massive is the particle, shorter is the spatial range.
At last in Newtonian Mechanics the Gauss theorem gives us a spherically symmetric gravitational potential even if the spherically symmetric source is rotating. In GR as well as in Fourth Order Gravity, however, the rotating spherically symmetric source generates an axially symmetric space-time (the well known Kerr metric) and only if the source is at rest one has the space-time with the same symmetry (the well known Schwarzschild metric). Then the galaxy being a rotating system will generate an axially symmetric space-time while we are using the solution (41). This aspect is not contradictory because the solutions are calculated in the Newtonian limit (i.e. v 2 1) and under this assumption the Kerr metric collapses into Schwarzschild metric. In fact we have where Σ 2 = r 2 + η 2 cos 2 θ, H 2 = r 2 − r g r + η 2 , η = L/M and L is the angular momentum along the z-axis.

Mass models for galaxies
From the point of view of morphology, a galaxy can be modeled by considering at least two components: the bulge and the disk. Obviously the galaxy is a more complicated structure and there are others components, but for our aim this idealization is satisfactory. The bulge, generally, can be represented easily with cylindrical coordinates (but in a more crude idealization it is like a ball), while the disk has a radius bigger than the thickness. However we find that the rotation curve does not present the Keplerian behavior outside the matter, but the curve remain constant for any distance. Then we must formulate the existence of exotic matter that can justify the experimental observation. A simple discussion about the distribution of DM can be formulated by imposing the constant value of velocity in (83) for large distances. In fact we find A matter distribution as (85) has a problem when we want to calculate the total mass. In fact if we have ρ ∼ r −2 , the mass diverges. A such exotic behavior seems no-physical, but this outcome is only the consequence of constant rotation curve. In fact by increasing the distance also the mass must increase with the power law for any distance (83). However, since the Gauss theorem holds in GR, the matter outside the sphere of integration does not contribute to the gravitational flux and then we do not have difference with respect to the ordinary matter.
The same arguments are not valid in f (X, Y, Z)-gravity: the non-viability of Gauss theorem implies that the range of integration of DM could cover all range and also the matter outside is considered. A cut off is needed now. Here then for completeness we consider that the galaxy is composed by three components: the bulge, the disk and an alone of DM.
It should be noted that the spatial behaviors of DM (generally spherically symmetric) are made only a posteriori: the cornerstones of rotation curves study are the GR and the distribution of ordinary matter. Only after this assumption the distribution of DM is such as to justify the gap between the theoretical prediction and the experimental observation.
Before jumping to analysis of rotation curves we want to resume the principal spatial distributions of mass in the three galactic components. In literature there are many forms of density, but it is possible to resume them as follows.
More realistic models are the ones with mass density depending also on the z-coordinate for bulge and disk. Particularly one can consider the following choice In this choice, suggested by Dehnen & Binney [80], the bulge is described as a truncated power -law model (first line of (86) z ). β, γ, q and ξ t are the parameters. While for the disk, one adopted a double exponential where the total mass is M d = 2πξ d 2 Σ e ξ 0 ξ d with Σ = (48 ± 8) M /pc 2 and ξ 0 = 8.5 Kpc [81]. Finally in the case of DM the density profile is the NFW model [7,82] where g(x) = ln(1 + x) − x x+1 , M vir DM and ξ vir DM are the virial mass and virial radius and ξ s is a characteristic length. For the Milky Way ξ vir DM /ξ s = 10 ÷ 15. Leaving the axis-symmetry one can consider a more simple model: the spherical symmetry model. With this approach one has [83][84][85] where Γ(x) is the Gamma function, 0 ≤ γ < 3 is a free parameter and 0 ≤ α < 1 is the ratio of DM inside the sphere with radius ξ DM with respect to the total DM. The radius ξ DM and the mass M DM play conceptually the same role respectively of ξ vir DM and M vir DM . However, as before claimed, the hot point is the choice of DM model. Given all these models it is almost normal that there are many different estimates of DM. Therefore, the parameters of DM model may not be unique [87].

Rotation curves in f (X, Y, Z)-gravity
We are interested to evaluate the circular velocity (82) adopting the mass models (88). So the potential (41), by setting f X (0) = 1, becomes where Ξ is the distance on which we observe the rotation curve, K is the elliptic function and the modulus of distance is given by The constant a is a scale factor defined by the substitution R , r → a r , a R so all quantities are dimensionless. At last, by remembering r = √ R 2 + z 2 the circular speed (82) in the galactic plan is given by In the Figs. 2, 3, we report the spatial behaviors of rotation curve induced by the bulge and disk component. The behavior for any component is compared in the framework of GR, FOG, GR + DM and FOG + DM. The values of free parameters of model are in the first line in Table 1 and refereing to Milky Way. The values of scale lengths µ 1 , µ 2 are set at 10 −2 a −1 , 10 2 a −1 . In both components we note for R >> ξ b , ξ d the Keplerian behavior, while it is missing only when we consider also the DM component. The shape of the rotation curve is similar to ones obtained by varying the total mass and scale radius. For a given scale radius, the peak velocity varies proportionally to a square root of the mass. For a fixed total mass, the peak-velocity position moves inversely proportionally to the scale radius, or along a Keplerian line.
As it is known in literature f (X, Y, Z)-gravity, and in particular f (R)-gravity, mimics a partial contribution of DM. In fact the corrective term ∝ e −µ 1 |x| /|x| contributes to enhance the attraction and thus the rotation curve must increase to balance the force. In the case of the other term, we have a correction ∝ −e −µ 2 |x| /|x| that contributes, being repulsive, to decrease the velocity. However in both cases these terms are asymptotically null and f (X, Y, Z)-gravity and GR must lead to the same result.
Only with the addition of DM it is possible to raise the curve and have almost constant values. In the Fig. 4 we report the component of rotation curve induced by only auto-gravitating DM.  In Fig. 5 we show the global behavior (experimentally expected) of rotation curve compared with respect to the bulge, disk and DM component for the f (X, Y, Z)-gravity. While in Fig. 6 there is the global rotation curve in the framework of GR, FOG, GR + DM and FOG + DM. At last in Fig. 7 we replicate the outcome of Fig. 6 but we inserted the value µ 2 = 5 a −1 . In this case the rotation curve induced by f (X, Y, Z)-gravity allows lower values as previously we claimed.    Table 1. Parameters of models (88). The unity of mass is 10 10 M and a = 1 Kpc. From the experimental point of view we used an updated rotation curve of Milky Way by integrating the existing data from the literature, and plot them in the same scale [88]. The data used are available in a digitized from the URL http://www.ioa.s.u-tokyo.ac.jp/ sofue/mw/rc2009/. The unified rotation curve shows clearly the three dominant components: bulge, disk, and flat rotation due to the DM [89][90][91][92][93][94][95][96]. These data, finally, have been updated further by [97]. The whole set of data are plotted in Fig. 8 and on them the theoretical rotation curve induced by f (X, Y, Z)-gravity with DM has been superimposed. The values of best fit are shown in Table 1 with µ 1 = 10 −2 Kpc −1 and µ 2 = 10 2 Kpc −1 .
The same mass model (88) has been considered also for the galaxy NGC 3198. This galaxy has been chosen since the bulge is missing. Then we set M b = 0 in the (89). In Fig. 9 we show the experimental data [98] and the superposition of theoretical behavior. Also in this case we find a nice outcome for a new set of parameters shown in Table 1, while the values of µ 1 and µ 2 are the same of Milky Way.
The initial aim, i.e. to extend the GR to a new class of theories, as we claimed in the introduction, is to justify the rotation curve without the DM component. From the previous outcomes, we see that even if the f (X, Y, Z)-gravity, or better a f (R)-gravity, admits a stronger attractive force, it is unable to realize our aim. Also in this framework we need Dark Matter. Obviously we need a smaller amount of DM on the middle distances, but for large distances we have the same problems of GR.

Rotation curves in R n -gravity and f (X, Y, Z)-gravity
The problem of DM seems to have been solved in literature, in the framework of f (R)-gravity, by considering the Lagrangian L = R n with n Q [59,62]. In these papers the gravitational potential for a point-like source can be  Table 1. The values of "masses" are µ 1 = 10 −2 Kpc −1 and µ 2 = 10 2 Kpc −1 .  Table 1. The values of "masses" are µ 1 = 10 −2 Kpc −1 and µ 2 = 10 2 Kpc −1 .
where r c is a characteristic length and β is a dimensionless parameter. To recover the condition lim r→∞ Φ R n (r) = 0 one must have 0 ≤ β < 1. In the case β = 0 the GR is found. We comment about the physical behavior of potential (92) and we want to add some reflections considering the result of the rotation curve shown above. Before to analyze the mathematical properties of metric linked to potential (92), we want to show the different values of correction to the Newtonian potential. In Fig. 10 we report the radial behavior of the corrections to 1/r for the potentials (38) and (92) (to minimize the difference we considered only f (R)-gravity). From the plot we note a discrepancy between the two corrections. The correction by R − R 2 6 µ 1 2 -gravity acts over distances much smaller, while the correction induced by R n -gravity provides a potential nearly constant over large intervals and slowly goes to zero (∼ r β−1 ). For this aspect the potential (92) does not need the DM component. Then with a procedure of fine tuning of r c and β it was possible to justify the experimental rotation curve for a wide class of galaxies [62] when n = 3.5. This choice was possible because there must be a relationship β = β(n) so that the potential (92) was compatible with respect to the field equations. These are the positive aspects of the potential (92) used in [59,62]. We conclude this section by reviewing the fundamental weaknesses of R n -gravity.
• The potential (92) presents an analogous behavior of potential (38). In fact for r < r c one has 1−(r/rc) β 2 r > 0 then the correction is "repulsive" like one induced by Ricci tensor square, while for r > r c one has an attractive correction. Now by remembering the reason of extension of GR, now we have unlike a repulsive contribution for r < r c . If in f (X, Y, Z)-gravity we can delete the Ricci tensor square contribution and we have only the f (R)-gravity, in R n -gravity we must collapse only in GR.
• The potential (92) belongs to general class of solutions for R n -gravity classified by a perturbative method [58], but the solutions are n-independent. Obviously the general solutions (it would be hard challenge to find them) are n-dependent, but at first order with respect to the perturbative parameter and in the vacuum 9 the field equations are identically vanishing. So we say that the presence of matter has been not considered and the choice of arbitrary constant has been evaluated only by matching R n -gravity with GR in the limit β → 0. In fact by solving the field equations correctly in presence of matter (also with the point-like source) we would obtain solutions depending on the perturbative parameter and the technique is misplaced.
For these two aspects, but especially for the second point, R n -gravity does not admit the Newtonian limit if n = 1. The potential (92) does not follow a correct framework when extending the GR to the new theories we want to generalize the Newtonian potential. Generally all theories without Ricci scalar in the Lagrangian suffer from the same problem. For example also R 2 -gravity is in the same situation: is not possible to extend the solution in the matter [58]. Although we have solutions as 1/r with additional asymptotically flat terms, it is not automatic the assertion that these solutions are the Newtonian limit of theory.
Let us analyze now the mathematical properties of the metric trying to justify the difference of spatial behaviors in Fig. 10. To simplify the calculation we choose the set of the standard coordinates. The metric (21), from the expressions (38), becomes 10 where dΩ = dθ 2 + sin 2 θdφ 2 is the solid angle, while the element of distance linked to potential (92) can be written as follows where Φ SC R n (r) is the the potential (92) and Ψ SC R n (r) is the other potential missing in the paper [62]. The metrics (93) and (38) represent the same space-time at first order of r g /r. However in their analysis the knowledge of last potential is useless because its contribution in the geodesic motion is at fourth order. By following the paradigm of weak field limit at small velocity [55][56][57] for the R n -gravity we find 9 This parameter is generally c −2 , but the analysis is the same if we consider the dimensionless ratio r g /r. 10 The set of standard coordinates is defined by the condition to obtain the standard definition of the circumference with radius r. From the metric tensor (21) we must impose the condition 1 − 2Ψ(r) r 2 =r 2 for the new radial coordinate.
where K(β) and K X are constants depending, respectively, on the value of β and on the integral operation, while the Ricci scalar R could be an arbitrary function. In fact it needs some comment about the index n in the Ricci scalar. If n is a integer number, then the Ricci scalar can assume any value and can be also a generic space depending function. More attention is needed if n is a rational number. The field equations (56) take into account up to third derivatives with respect to the Ricci scalar, then we must ensure that the function f (R) and its derivatives are always well defined [58]. In this case for n < 3 the solution Ricci flat (R = 0) or space depending and asymptotically vanishing are excluded. Only solutions with constant values are allowed, but the algebraic sign is crucial. A such behavior is expected any time we have the condition lim R→0 f (R) = constant [58]. Now in all these considerations we do not recover the condition lim r→∞ Ψ SC R n (r) = 0: then we have a theory which does not provide the Minkowskian limit. It is using the first perturbative contribution of a metric component (providing the flatness at infinity), while other contributions in the remaining metric components (negligible in the Newtonian limit) do not cover the same asymptotic limit. Then only for n > 3 we can have the flatness at infinity.
Then we could say that for n > 3 the Minkowskian limit is recovered but a such perturbative approach can be performed only in the vacuum. The objection previously shown comes back. Up to third order (c −6 or r g 3 ) the geometrical side of field equation is identically null, but the matter side could not be null at first order (c −2 or r g ). The class of R n -gravity are examples of theories where the weak field limit procedure does not generate automatically the Minkowskian limit. In fact only if we consider theories satisfying the condition lim X→0 f (R) = 0 [58], their weak field limit is compatible with the request of asymptotically flatness. Moreover f (R)-gravity mimicking an additional source due to its scalar curvature [55,56,60] we would have a constant matter that pervades all space giving us a justification of more intense gravitational potential. In addition if lim R→0 f (R) = costant we do not have the Minkowskian limit, but we can interpret the apparent mass, only from the experimental point of view, as DM. These aspects, then, can be a mathematical motivation for different shape of point-like gravitational potential, but also source of further attention.

The Gravitational Lensing in f (X, Y, Z)-gravity
Gravitational lensing can be considered as one of the main probes for DM at galactic and extragalactic level. However, modified gravity could affect the lensing quantities vs DM. Let us develop this issue in the framework of Fourth Order Gravity.

The Point-Like Source
Let us start by defining the Lagrangian of a photon moving in the gravitational field with metric (93). It is where Ξ(r) . = 1 +Ξ(r) . = 1 + 1 3 e −µ 1 r − 4 3 e −µ 2 r , Λ(r) e −µ 2 r and the dot represents the derivatives with respect to the affine parameter λ. Since the variable θ does not have dynamics (θ = 0) we can choose for simplicity θ = π/2. By applying the Euler-Lagrangian equation to Lagrangian (96) for the cyclic variables t, φ we find two motion constants and respect to λ we find the "energy" of Lagrangian 11 By inserting the equations (97) into (98) we find a differential equation forṙ r + is the solution for leaving photon, whileṙ − is one for incoming photon. Let r 0 be a minimal distance from the lens center (Fig. 13). We must impose the conditionṙ ± (r 0 ) = 0 from the which we find Now the deflection angle α (Fig. 13) is defined by following relation where λ 0 is the value of λ corresponding to the minimal value (r 0 ) of radial coordinate r. By putting the expressions of J,φ andṙ + into (101) we get the deflection angle 11 The (96) is a quadratic form, so it corresponds to its Hamiltonian. which in the case r g /r 1 12 becomes where From the definition ofΞ andΛ we note that in the case f (X, Y, Z) → R we obtain F µ 1 , µ 2 (r 0 ) → 0. In a such way we extended and contemporarily recovered the outcome of GR. The analytical dependence of function F µ 1 , µ 2 (r 0 ) from the parameters µ 1 and µ 2 is given by evaluating the integral (104). A such as integral is not easily evaluable from the analytical point of view. However this aspect is not fundamental, since we can numerically appreciate the deviation from the outcome of GR. In fact in Fig. 11 we show the plot of deflection angle (103) by f (X, Y, Z)-gravity for a given set of values for µ 1 and µ 2 . The spatial behavior of α is ever the same if we do not modify µ 2 . This outcome is really a surprise: by the numerical evaluation of the function F µ 1 , µ 2 (r 0 ) one notes that the dependence of µ 1 is only formal. If we solve analytically the integral we must find a µ 1 independent function. However, this statement should not be justified only by numerical evaluation but it needs an analytical proof. For these reasons in the next section we reformulate the theory of Gravitational Lensing (GL) generated by a generic matter distribution and demonstrate that for f (R)-gravity one has the same outcome of GR. Figure 11. Comparison between the deflection angle of GR (solid line) and one of f (X, Y, Z)-gravity (dashed line) (103) for a fixed value µ 2 = 2 and any µ 1 .

Extended Sources
In this section we want to recast the framework of GL for a generic matter source distribution ρ(x) so the photon can undergo many deviations. In this case we leave the hypothesis that the flight of photon belongs always to the same plane, but we consider only the deflection angle as the angle between the directions of incoming and leaving photon. Finally we find the generalization of GL in f (X, Y, Z)-gravity including the previous outcome of deflecting point-like source (and resolving the integral (104)).
The starting point is the relativistic distance (39) where Φ is given by the superposition principle (41) and an analogous relation is valid also for Ψ. By introducing the four velocity u µ =ẋ µ = (u 0 , u) the flight of photon is regulated by the condition then u µ is given by In the Newtonian limit we find that the geodesic motion equation becomeṡ and by supposing |u| 2 = 1 we can recast the equation in a more known aspecṫ where ∇ ⊥ = ∇ − u |u| · ∇ u |u| is the two dimensions nabla operator orthogonal to direction of vector u. In GR we would had onlyu = −2 ∇ ⊥ Φ since we have Ψ = Φ. In fact the field equations (28) admit the pointlike solutions (38) and the condition (34) is satisfied [57]. We can affirm, then, that only in GR the metric potentials are equals (or more generally their difference must be proportional to function |x| −1 ). The constraint (34) has been found also many times in the context of cosmological perturbation theory [99][100][101][102][103].
The deflection angle (101) is now defined by equation where λ i and λ f are the initial and final value of affine parameter [104]. For a generic matter distribution we can not a priori claim that the deflection angle belongs to lens plane (as pointlike source), but we can only link the deflection angle to the difference between the initial and final velocity u. So we only analyze the directions of photon before and after the interaction with the gravitational mass. Then the (109) is placed by assuming α = ∆u = u i − u f . From the geodesic equation (108) the deflection angle becomes The formula (110) represents the generalization of deflection angle in the framework of GR. By considering the photon incoming along the z-axes we can set u i = (0, 0, 1). Moreover we decompose the general vector x R 3 in two components: ξ R 2 and z R. The differential operator now can be decomposed as follows ∇ = ∇ ⊥ +ẑ ∂ z = ∇ ξ +ẑ ∂ z , while the modulus of distance is = ∆( ξ, ξ , z, z ). Since the potentials Φ , Ψ 1, around the lens, the solution of (108) with the initial condition u i = (0, 0, 1) can be expressed as follows and we can substitute the integration with respect to the affine parameter λ with z. In fact we note and the deflection angle (110) becomes From the expression of potentials for a generic matter distribution (see potential (41) where f X (0) = 1 for example) we find the relations and the deflection angle (110) becomes In the case of hypothesis of thin lens belonging to plane (x, y) we can consider a weak dependence of modulus ∆( ξ, ξ , z, z ) into variable z so there is only a trivial error if we set z = 0. With this hypothesis the integral into z is incorporated by definition of two dimensional mass density Σ( ξ ) = dz ρ( ξ , z ). Since we are interested only to the GL performed by one lens we can extend the integration range of z between (−∞, ∞). Now the deflection angle is the following The last integral in (115) is vanishing because the integrating function is odd with respect to variable z. The expression (116) is the generalization of outcome (103) and mainly we found a correction term depending only on the µ 2 parameter. In the case of point-like source Σ( ξ ) = M δ (2) ( ξ ) we find and in the case of f (X, Y, Z) → f (R) (i.e. µ 2 → ∞ and F µ 2 ( ξ, ξ ) → 0) we recover the outcome of GR α = 2 r g ξ/| ξ| 2 . From the theory of GL in GR we know that the deflection angle 2 r g /r 0 is formally equal to 2 r g /| ξ| if we suppose r 0 = | ξ|. Besides both r 0 , | ξ| are not practically measurable, while it is possible to measure the so-called impact parameter b (see Fig. 13). But only in the first approximation these three quantities are equal.
In fact when the photon is far from the gravitational source we can parameterize the trajectory as follows and from the definition of angular momentum (97) in the case of t b we have By using the condition (100)ṙ ± (r 0 ) = 0 we find the relation among b and r 0 justifying then the position r 0 = | ξ| in the limit r g /r 1 (but also r g /r 0 1). In Fig. 12 we report the plot of deflection angle (118). The behaviors shown in figure are parameterized only by µ 2 and we note an equal behavior shown in Fig. 11. With the expression (118) Figure 12. Comparison between the deflection angle of GR (solid line) and of f (X, Y, Z)-gravity (dashed line) (118) for 0.2 < µ 2 < 2.
we have the analytical proof of statement at the end of previous section. In fact in the equation (118) we have not any information about the correction induced in the action (18) by a generic function of Ricci scalar (f XX = 0). This result is very important if we consider only the class of theories f (R)-gravity. In this case, since µ 2 → ∞, we found the same outcome of GR. From the behavior in Fig. 12 we note that the correction to outcome of GR is deeply different for r 0 → 0, while for r 0 → ∞ the behavior (118) approaches the outcome of GR, but the deviations are smaller. This difference is given by the repulsive correction to the gravitational potential (see metric (39)) induced by f (Y, Z). Only by leaving the thin lens hypothesis (the lens does not belong to plane z = 0) we can have the deflection angle depending by µ 1 (115). In fact in this case the third integral in (115) is not zero. Then in the case of thin lens we have a complete degeneracy of results: f (R)-gravity equivalent to the GR in the thin lens case (see also [68]). In order to find some differences, we must include the contributions generated by the squared Ricci tensor. However also in this case we do not get the right behavior: the deflection angle is smaller than the one in GR. Furthermore f (X, Y, Z)-gravity does not mimic the DM component by assuming the thin lens hypothesis.
A final remark is order at this point. The fact that f (R) gravity and GR give the same lensing effects is due to the fact that the corrections induced by f (R) on the g tt and g ij components of the metric (that is the potentials Φ and Ψ) are opposite and cancel out as demonstrated in details in Ref. [67,68]. This means that the only GR contribution remains in the deflection angle.

The lens equation
In order to demonstrate the effect of a deflecting mass, let us show the simplest GL configuration in Fig. 13. A point-like mass is located at a distance D OL from the observer O. The source is at distance D OS from the observer, and its true angular separation from the lens L is β, the separation which would be observed in the absence of lensing (r g = 0). The photon which passes the lens at distance r 0 ∼ b is deflected with an angle α.
Since the deflection angle (103) is equal to (118), for sake of simplicity we will use the "vectorial" expression. Then the expression (118), by considering the relation (121), becomes The condition that this photon reach the observer is obtained from the geometry of Fig. 13. In fact we find Here D LS is the distance of the source from the lens. In the simple case with a Euclidean background metric here, D LS = D OS − D OL ; however, since the GL occurs in the Universe on large scale, one must use a cosmological model [104]. Denoting the angular separation between the deflecting mass and the deflected photon as θ = b/D OL the lens equation for f (X, Y, Z, )-gravity is the following where θ E = 2 rg D LS D OL D OS is the Einstein angle and (125) Figure 13. The gravitational lensing geometric for a point-like source lens L at distance D OL from observer O. A source S at distance D OS from O has angular position β from the lens. A light ray (dashed line) from S which passes the lens at minimal distance r 0 is deflected by α; the observer sees an image of the source at angular position θ = b/D OL where b is the impact factor. D LS is the distance lens -source.
Since we have 0 < θ 2 F(θ) < 1 (Fig. 14) we can find a perturbative solution of (124) by starting from one in GR, θ GR . In fact by assuming θ = θ GR ± + θ * and neglecting θ * 2 F(θ * ) in (124) we find and in the case of β = 0 we find the modification to the Einstein ring In Fig. 15 we show the angular position of images with respect to the Einstein ring. Both the deflection angle and the position of images assume a smaller value than ones of GR. Then the corrections to the GR quantities are found only for the introduction in the action (18) of curvature invariants Y (or Z), while there are no modifications induced by adding a generic function of Ricci scalar X. The algebraic signs of terms concerning the parameter µ 2 are ever different with respect to the terms of GR in (38) and they can be interpreted as a "repulsive force" giving us a minor curvature of photon. The correction terms concerning the parameter µ 1 have opposite algebraic sign in the metric component g tt and g ij (38) and we lose their information in the deflection angle (113).
In both approaches we find the same outcomes µ 1 -independent because the matter source (in our case it is a point-like mass) is symmetric with respect to z-axes and we neglect the second integral in (115).  Obviously for a generic matter distribution the deflection angle is defined by (115) and the choice of second derivative of function of Ricci scalar is not arbitrary anymore.

Conformal transformations in the weak field limit
As a final issue, let us discuss the conformal transformations in the weak field limit approximation. It is important to stress that we developed all the above considerations in the Jordan frame in order to control the behavior of the standard matter that, in this frame, remains minimally coupled. In the Einstein frame, the matter is non-minimally coupled and this fact could give rise to some difficulties in order to control the dynamics. Clearly, this means that all the considered modifications insist on the geometric part. Controlling the relation between physical quantities passing from a frame to another could give hints to understand if conformal transformations are only a mathematical tool or have some physical meaning (see [21] for a detailed discussion).

Scalar tensor gravity in the Jordan frame
Let us start by considering the action of the scalar-tensor gravity in 4 dimensions. It is Let us note that the action is apparently more general than (128). In fact by substituiting F (φ) → φ, we obtain only a new definition of functions ω(φ) and V (φ) so the two formulations are essentially equivalent. The term L m is the minimally coupled ordinary matter contribution considered as a perfect fluid; ω(φ) is a function of the scalar field and V (φ) is its potential which specifies the dynamics. Actually if ω(φ) = ±1, 0 the nature and the dynamics of the scalar field is fixed. It can be a canonical scalar field, a phantom field or a field without dynamics (see e.g. [105,106] for details). In the metric approach, the field equations are obtained by varying the action (128) with respect to g µν and φ. The field equations are and the trace equation is Here we introduced V φ = dV dφ and ω φ (φ) = dω(φ) dφ . If we assume that the Lagrangian density L m of matter depends only on the metric components g µν and not on its derivatives, we obtain T µν = 1/2 L m g µν − δL m /δg µν . Let us consider a source with mass M . The energy-momentum tensor is the same of equation (35), (36). It is useful to get the expression of L m . In fact from the definition of energy momentum tensor we have and further from the mathematical properties of metric tensor we have then we find The variation of density is given by order to deal with standard self-gravitating systems, any theory of gravity has to be developed in its Newtonian or post-Newtonian limit depending on the order of approximation in terms of squared velocity v 2 [58].
In this context, also the scalar field φ is approximated as the Ricci scalar. In particular we get φ = φ (0) + φ (2) + . . . while the functions V (φ) and ω(φ) can be substituted by their corresponding Taylor series.
From the lowest order of field equations (130) we have and also in the scalar tensor gravity a missing cosmological component in the action (128) implies that the space-time is asymptotically Minkowskian; moreover the ground value of scalar field φ must be a stationary point of potential. In the Newtonian limit, we have These equations are not simply the merging of field equations of GR and a further massive scalar field, but come out to the fact that the scalar tensor gravity generates a coupled system of equations with respect to Ricci scalar R and scalar field φ. The gravitational potentials Φ, Ψ and the Ricci scalar R (2) are given by and supposing that 2 ω(φ (0) ) φ (0) − 3 = 0 we find for the scalar field φ (2) the Yukawa-like field equation where we introduced the mass definition It is important to stress that the potential Ψ can be found also as see for example [107]. By using the Fourier transformation, the solution of Eq. (139) has the following form The expressions (138) and (142) represent the most general solution of any scalar-tensor gravity in the Newtonian limit. Since the superposition principle is yet valid (as previously analyzed), it is sufficient to consider again the solutions generated by the pointlike source. Then the solutions are [57,58,107] φ (2) In the case V (φ) = 0, the scalar field is massless and ω(φ) = −ω 0 /φ, we obtain the well-known Brans-Dicke solutions [41] with Eddington's parameter γ = 1+ω 0 2+ω 0 [70] where the gravitational constat is defined as

Scalar tensor gravity in the Einstein frame
Let us now introduce the conformal transformation where Ω = Ω(x) is a nowhere vanishing, regular function, called a Weyl or conformal transformation. This transformation has been introduced to show that scalar-tensor theories are, in general, conformally equivalent to the Einstein theory plus minimally coupled scalar fields. However if standard matter is present, the conformal transformation generates the non-minimal coupling between the matter component and the scalar field. By applying the transformation (145), the action in (128) can be reformulated as follows in whichR is the Ricci scalar relative to the metricg µν and Ξ is a generic constant. The two actions (128) and (146) are mathematically equivalent. In fact the conformal transformation is given by imposing the condition The relations between the quantities in the two frames arẽ The field equations for the new fieldsg µν andφ are whereT µν and˜ are the re-definition of the quantities (35) The integration of field equations (149) is only formal because we do not know the analytical expression of the coupling function between the matter and the scalar fieldφ (see the third line of (148)). We can make some assumptions on the parameter Ξ and the functionω(φ) in the minimally coupled Lagrangian (146) and on the function ω(φ) in the nonminimally coupled Lagrangian (128). If we choosẽ ω(φ) = −1/2, Ξ = 1 and ω(φ) = −ω 0 /φ, the transformation between the scalar fields φ andφ is given by the first line in (148), that is where obviously ω 0 > −3/2 andφ 0 is an integration constant 14 . The potential W and the matter LagrangianL m are In both frames, the scalar fields are expressed as perturbative contributions on the cosmological background (φ (0) ,φ (0) ) with respect to the dimensionless quantity v 2 . Then also for the scalar fieldφ, we can consider the developφ =φ (0) +φ (2) + . . . . Such a develop can be applied to the transformation rule (151) and we obtaiñ Since we are interested in the Newtonian limit of field equations (149), we can assume, for the conformally transformed metricg µν , an expression as (39) but with some differences. In fact from the conformal transformation (145) and from the last line of (148), we havẽ then the conformally transformed metric becomes and the relation between the gravitational potentials in the two frames is Then the field equations (149) become 15 15 With the assumptions of the metric (155) the Ricci tensorR µν in the Newtonian limit has the form Φ φ (0) (a similar behaviour forR (2) ij ), where the Ricci scalar is scaled by the factor φ (0) 2 . The same scaling occurs for the Laplacian: where also in this case we have W (φ (0) ) = 0 and Wφ(φ (0) ) = 0. However these conditions are an obvious consequence of the conformal transformation of conditions V (φ (0) ) = 0 and V φ (φ (0) ) = 0. In fact we can figure out that V (φ) ∝ (φ − φ (0) ) 2 and then W (φ) ∝ eφ which, by using relations (153), satisfies the above conditions. Finally, we note that Wφφ(φ (0) ) = 2ω 0 +3 and by the definition of mass m φ 2 , given in Eq. (140), we obtain Wφφ(φ (0) ) = m φ 2 /φ (0) .
Finally, the energy-momentum tensorT µν is given by the following expressioñ whereg µνũ µũν = 1 thenũ 0 = φ (0) + 2Φ. In the Newtonian limit, we findT (2) 00 = ρ/φ (0) and It remains only to calculate the source term δL m /δφ of the scalar fieldφ (2) . From the third line of (148) and, by using the transformation rules (151), we find the coupling between the scalar field and the ordinary matter Then the system of Eqs. (157) becomes and their solutions in the case of pointlike source arẽ The difference in Eqs.(156) between the gravitational potentials is satisfied by using the expression of scalar field in the Jordan frame (first line of (143)) where, obviously, we set ω(φ) = −ω 0 /φ. In fact we find and an analogous relations is found also for the couples Ψ,Ψ. Furthermore we can check also the transformation rules (151) and (153) for the solutions (143) and (161) of the scalar fields φ,φ.
The redefinition of the gravitational constant G (as performed in the Jordan frame G → G * in the case of Brans-Dicke theory [41]) is not available when we are interested to compare the outcomes in both frame. In fact the couple of potentials Φ,Φ differs not only from the dynamical contribution of the scalar field (φ (2) ) but also from the definition of the gravitational constant. Furthermore, in the Einstein frame, the scalar fieldφ does not contribute (the coupling constant between R andφ is vanishing), then we find the same outcomes of GR with ordinary matter. However by supposing the Jordan frame as starting point and coming back via conformal transformation, we find that the gravitational constant is not invariant and depends on the background value of the scalar field in the Einstein frame, that is

The analogy between the scalar tensor gravity and f (R)-gravity
Recently, several authors claimed that higher-order theories of gravity and among them, f (R) gravity, are characterized by an ill defined behavior in the Newtonian regime. In particular, it is discussed that Newtonian corrections of the gravitational potential violate experimental constraints since these quantities can be recovered by a direct analogy with Brans-Dicke gravity simply supposing the Brans-Dicke characteristic parameter ω 0 vanishing (see [50][51][52] for a discussion). Actually, the calculations of the Newtonian limit of f (R)-gravity, directly performed in a rigorous manner, have showed that this is not the case [56][57][58]108] and it is possible to discuss also the analogy with Brans-Dicke gravity. The issue is easily overcome once the correct analogy between f (R)-gravity and the corresponding scalar-tensor framework is taken into account. It is important mentioning that several important results already achieved in the Newtonian regime, see e.g. [109,110], are confirmed by these rigorous approach and only the wrong analogies are ruled out.
In literature, it is shown that f (R) gravity models can be rewritten in term of a scalar-field Lagrangian non-minimally coupled with gravity but without kinetic term implying that the Brans-Dicke parameter is ω(φ) = 0. This fact is considered the reason for the ill-definition of the weak field limit that should be ω → ∞ inside the Solar System.
Let us deal with the f (R) gravity formalism in order to set correctly the problem. The field equations (56) can be recast in the framework of scalar-tensor gravity as son as we select a particular expression for the free parameters of the theory. The result is the so-called O'Hanlon theory [111] which can be written as The field equations are obtained by starting from ones of (130) By supposing that the Jacobian of the transformation φ = f R is non-vanishing, the two representations can be mapped one into the other considering the following equivalence From the definition of the mass (140) we have φ V φ (φ) − 2V (φ) = 3 m φ 2 φ (2) , then we have also f R R − 2f = 3 m φ 2 φ (2) and by performing the Newtonian limit on the function f [58], we get f R (0)R (2) = −3 m φ 2 φ (2) . The spatial evolution of Ricci scalar is obtained by solving the field equations (56) without using the conformal transformation [57,58]. The solution for the potentials Φ, Ψ are obtained simply by setting ω(φ) = 0 in Eqs. (143) and φ (0) = f R (0). In the case f (R) → R, from the second line of (165), V (φ) = 0 → m φ = 0 and the solutions (143) become the standard Schwarzschild solution in the Newtonian limit. Finally, we can consider a Taylor expansion 16 of the form f = f R (0) R (2) + f RR (0) 2 R (2) 2 so that the associated scalar field reads φ = f R (0) + f RR (0) R (2) . The relation between φ and R (2) is R (2) = φ−f R (0) f RR (0) while the self-interaction potential (second line of (165)) turns out the be V (φ) = − (φ−f R (0)) 2 2 f RR (0) satisfying the conditions V (f R (0)) = 0 and V φ (f R (0)) = 0. In relation to the definition of the scalar field, we can opportunely identify f R (0) with a constant value φ (0) = f R (0) which justifies the previous ansatz for matching solutions in the limit of GR. Furthermore, the mass of the scalar field can be expressed in term of the Lagrangian parameters as m φ 2 = 1 3 φ (0) V φφ (φ (0) ) = − f R (0) 3f RR (0) . Also in this case the value of mass is the same obtained by solving the problem without invoking the scalar tensor analogy [57,58]. However with this remark, it is clear the analogy between f (R)-gravity and a particular class of scalar tensor theories [111]. Finally, the above results have to be interpreted also with respect to results reported in [50][51][52]. In those works, the relation between f (R) gravity and scalar-tensor theories was used to compute the weak-field limit. Though f (R) gravity, in metric formalism, is related to the Brans-Dicke theory with ω 0 = 0, the equivalence is not complete because f (R) has a non-trivial potential, which is not present in the original Brans-Dicke theory. The role of the potential is to introduce a length scale (as we found here) and then the post-Newtonian parameter γ develops a spatial dependence. In Refs. [50][51][52] this spatial dependence is discussed despite the fact that ω 0 = 0. Note that formulas presented in [50][51][52] for the weak-field limit are in complete agreement with those presented here also if achieved with a different approach.

Discussion and Conclusions
The weak field limit of Extended Theories of Gravity has been discussed in view of some relevant astrophysical issues. In particular, we have considered the hydrostatic equilibrium of stars, the galactic rotation curves and the gravitational lensing. Finally we have analyzed the relations between the Jordan and Einstein frames in the same limit and focused our attention on the relations linking the gravitational potentials in both frames.
Specifically, we have considered theories containing generic functions of the Ricci scalar R, the squared Ricci tensor R αβ R αβ and the squared Riemann tensor R αβγδ R αβγδ . We obtain the analogy, in the Newtonian limit, with the so-called Quadratic Lagrangian containing the squared Ricci scalar and the squared Ricci tensor in addition to the linear term R. All contributions to the field equations related to the squared Riemann curvature invariant can be expressed by the other two curvature invariants (squared Ricci tensor and squared Ricci scalar) via the Gauss-Bonnet invariant. It is straightforward to show that the spherically symmetric solutions show a Yukawa-like dependence with two characteristic lengths. An important result is that, for generic Fourth Order Gravity models, two non-equivalent metric potentials come out. In the limit f → R, only the Newtonian potential is present as it has to be in GR. The presence of two gravitational potentials, together with the non-validity of the Birkhoff and Gauss theorems, are the main differences between Fourth Order Gravity and GR.
Coming to the astrophysical applications, we adopted a polytropic equation of state relating the mass density to the pressure and derived a modified Lané-Emden equation whose solutions can be compared to the analogous solutions coming from GR. As soon as one considers the limit f (R) → R, the standard hydrostatic equilibrium theory is fully recovered. Since the Gauss theorem is not valid in this context and the modified Lané-Emden equation is an integro-differential equation, the mass distribution plays a crucial role. The correlation between two points in the star is given by a Yukawa-like term of the corresponding Green function. These feature opens the possibility to address the structure of anomalous stars that, in standard stellar theory could not be consistently achieved [38].
As further astrophysical application, we have considered the rotation curves of spiral galaxies. Starting from the pointlike solutions and having formulated the expression of the rotation curves, we considered the principal galactic components (i.e. the bulge, the disk and the halo). The theoretical curves have been compared with the experimental data. The rotation curves have been evaluated by considering the bulge and the DM component spherically symmetric and a circular disk symmetric with respect to a plane where the radius is larger than the thickness. Also in the case of Fourth Order Gravity with standard matter, we find that the rotation curve has the Keplerian behavior and only if one adds some DM component, a reliable matching with experimental data is achieved. However the DM hypothesis gives rise to two serious problems: since the DM distribution is diverging when we consider the whole amount of mass, it is crucial the choice of the cut-off inside the integral. Another critical point is the choice of the mass model and the values of free parameters. These shortcomings can be consistently addressed by introducing a further scalar field as in [107]. Such scalar field, however, can be reinterpreted in terms of curvature invariants [107].
The gravitational lensing approach in Fourth Order Gravity has been pursued on two steps: firstly we considered a point-like source and the motion of photons. In the second step, we took into account the geodesic motion and reformulated the light deflection for a generic matter distribution. In the case of an axially symmetric matter distribution, we obtained the standard relation between the deflection angle and the orthogonal gradient of metric potentials. Otherwise, the angle is depending also on the travel direction of the photon. In particular, if there is a z-symmetry, the deflection angle does not depend explicitly on the parameters of Fourth Order Gravity. Starting from the definition of the masses, one have to note that the contribution of a generic function of the Ricci scalar is only in the missing parameter. Then f (R)-gravity admits the same geodesic trajectory of GR. If one wants to take into account corrections to GR, one needs to add a generic function of the squared Ricci tensor into the Hilbert-Einstein action. In this case, we find the deflection angle smaller than the one of GR. The mathematical motivation is a consequence of the algebraic sign of the parameter in front of the squared Ricci tensor. In fact it is different with respect to the GR term and we can interpret it as a "repulsive force" giving a lower curvature for the photon trajectory.
A similar result is found for the galactic rotation curve, where the contribution of the squared Ricci tensor gives a lower rotation velocity profile than the one derived from the weak field limit of GR. However, if we estimate the weight of the corrections induced by f (R)-gravity, we have a perfect agreement with the GR from the point of view of gravitational lensing. Only by adding f (R αβ R αβ ) in the action, we induce modifications in both gravitational lensing and galactic rotation profiles.
Finally, we have tackle the debate of selecting the physical frame by conformal transformations. Specifically, we have shown how the Newtonian limit behaves in the Jordan and in the Einstein frame. The general result is that Newtonian potentials, masses and other physical quantities can be compared in both frames once the perturbative analysis is correctly performed. The main result is that if such an analysis is carefully developed in a frame, the perturbative process can be controlled step by step leading to coherent results in both frames. In other words, also if the gauge invariance is broken, there is the possibility to control conformal quantities and fix the observables.