Black holes, cosmological solutions, future singularities, and their thermodynamical properties in modified gravity theories

Along this review, we focus on the study of several properties of modified gravity theories, in particular on black-hole solutions and its comparison with those solutions in General Relativity, and on Friedmann-Lemaitre-Robertson-Walker metrics. The thermodynamical properties of fourth order gravity theories are also a subject of this investigation with special attention on local and global stability of paradigmatic f(R) models. In addition, we revise some attempts to extend the Cardy-Verlinde formula, including modified gravity, where a relation between entropy bounds is obtained. Moreover, a deep study on cosmological singularities, which appear as a real possibility for some kind of modified gravity theories, is performed, and the validity of the entropy bounds is studied.


I. INTRODUCTION
General Relativity (GR) has been the most successful gravitational theory of the last century, fully accepted as a theory that describes the macroscopic geometrical properties of space-time. For an isotropic and homogeneous geometry, GR leads to Friedmann equations which describe in an appropriate way the cosmological evolution with radiation and then matter dominated epochs. Nevertheless, the development of observational cosmology in the last decades with experiments of increasing precision like supernovae observations [1][2][3] has revealed that the Universe is in a stage of accelerated expansion. GR provided with usual matter sources is not able to explain this phenomenon. Moreover, GR does not account either for the cosmological era known as inflation [4], believed to have taken place before the radiation stage and that could alleviate some problems of standard cosmology like the horizon and the flatness problem [5]. In addition, GR with usual baryonic matter cannot explain the observed matter density determined by fitting the standard ΛCDM model to the WMAP7 data [6], the latest measurements from the BAO (Baryon Acoustic Oscillations) in the distribution of galaxies [7] and the Hubble constant (H 0 ) measurement [8]. Thus, GR requires the introduction of an extra component called dark matter that accounts for about 20% of the energy content of our Universe.
A more puzzling problem is associated to the present accelerated expansion of the Universe. There are also a large amount of different explanations. One of them, assuming the validity of GR, postulates the existence of an extra cosmic fluid, the dark energy (DE), whose state equation p = ω DE ρ (where p and ρ are the pressure and the energy density of the fluid) demand ω DE < −1/3 in order to provide an accelerated cosmic expansion [9][10][11]. The cosmological constant is the simplest model of DE, corresponding to an equation of state ω DE = −1. However, if we assume that the cosmological constant represents the quantum vacuum energy, its value seems to be many orders of magnitude bigger than the observed one [12].
Thus, alternative explanations to the cosmological constant have been largely studied in the last years.
These theories essentially modify GR by considering actions different from the Einstein-Hilbert one [13][14][15][16][17][18][19][20][21][22][23][24]. Examples are Lovelock theories, free of ghosts and whose field equations contain second derivatives of the metric at most; string theory inspired models, which include a Gauss-Bonnet term in the Lagrangian (see [25] and references therein); scalar-tensor theories like Brans-Dicke one, in which gravitational interaction is mediated by both a scalar field and GR tensor field; or the so-called f (R) theories (see [26][27][28][29][30][31][32][33] for a recent and extensive review). In this investigation we shall restrict ourselves to f (R) theories in the metric formalism (where the connection depends on the metric, so the present fields in the gravitational sector of the action come only from the metric tensor) in the Jordan frame. In this frame, the gravitational Lagrangian is given by R + f (R), where f (R) is an arbitrary function universe related with the entropy proportional to the area of the size of the universe, a kind of holographic principle. Finally, the formula may be rewritten as a dynamical entropy bound from which a number of entropy bounds, proposed earlier, follow. However, such a relation with the 2d CFT has been attempted to be extended for general perfect fluids (see [99]), and more general scenarios including modified gravity [100], but the equivalence could not be extended. Furthermore, it was found that the bound on the entropy may be violated in some particular scenarios, mainly when the universe moves close to a future singularity. Moreover, in spite of that in general the bound is violated sufficiently close to the singularity, even when some quantum effects are included, the violation can be also produced much before the singularity occurs [100].
Hence, all of these suggest a deep connection between gravitational physics and thermodynamics, which is still an open issue nowadays, and could help to improve our knowledge on the way to reconstruct consistent gravitational theories that may solve some of the most important problems on theoretical physics. Along this work, some of these results are reviewed, and we may suggest new scenarios where gravitational physics may be tested The present review is organized as follows: In Section II we revised the rudiments of the f (R) modified gravity theories that will be used throughout the paper. Then, Section III is devoted to study static and spherically symmetric configurations in these theories, with special focus on electromagnetic solutions and perturbative approach. The following Section IV deals with the Kerr-Newman black holes in f (R) theories. Section V encompasses a general revision of thermodynamics of f (R) black holes for both anti-de Sitter and Kerr-Newman scenarios. The thermodynamical stability is then studied for two paradigmatic examples of f (R) models in Section VI. Section VII deals with the reconstruction of F (R) actions for some particular cosmological solutions and how this kind of theories are capable to reproduce the whole cosmological evolution. In Section VIII, we review some of the thermodynamical properties of FLRW equations, and its derivation from the first law of thermodynamics, which can be easily extended to F (R) gravity. Section IX is devoted to some extensions of the Cardy-Verlinde formula to more complex configurations of FLRW universes, including modified gravity. In Section X, future singularities are studied in the context of entropy bounds, and the possible violation of the proposed bounds on the energy and entropy of the universe. Finally, we conclude the paper by giving our conclusions in Section XI.

II. F (R) THEORIES OF GRAVITY
In this section, F (R) gravity is introduced, which basically consists in an extension of the Einstein-Hilbert action by considering more complex gravitational actions consisting on functions of the Ricci scalar, and which can be rewritten in terms of a Brans-Dicke-like theory as shown below. The action for F (R)-gravity is given by: Here L m denotes the Lagrangian of some kind of matter, and κ 2 = 8πG is the gravitational coupling constant. In the present work we employ the natural units system in which = c = 1. Note also that our definition for the Riemann tensor is: Field equations are obtained by varying the action Equation (1) with respect to g µν , where T δg µν is the energy-momentum tensor and prime denotes here derivatives with respect to R. Nevertheless, in general it is very difficult to get exact solutions directly from field Equation (3), although some particular solutions can be obtained, as shown in the next sections. However, the action Equation (1) is equivalent to a kind of scalar-tensor theory with a null kinetic term, which can be used to reconstruct the appropriate action for some particular solutions, specially for cosmological solutions, as shown in Section VII. Then, the action Equation (1) is rewritten as follows [101]: By the variation of the action with respect to the metric tensor g µν , the field equation is obtained: g µν (P (φ)R + Q(φ)) + P (φ)R µν + g µν P (φ) − ∇ µ ∇ ν P (φ) = κ 2 T (m) µν (5) In addition, the action Equation (4) gives an extra equation for the scalar field φ, obtained directly from the action by varying it with respect to φ: P (φ)R + Q (φ) = 0 (6) where the primes denote derivatives with respect to φ. This equation can be solved such that the scalar field is a function of the Ricci scalar R, φ = φ(R), and then, by replacing this result in the action Equation (4), the action Equation (1) is recovered: With the aim of proposing a realistic alternative to GR, a possible modification consists of adding a function of the scalar curvature, f (R) to the Einstein-Hilbert (EH) Lagrangian. Therefore the gravitational sector of the action Equation (1) becomes: By performing variations with respect to the metric, the modified Einstein Equation (3) turns out to be: (1 + f R )R µν − 1 2 (R + f (R))g µν + D µν f R = 8πG T µν (9) where f R ≡ df (R)/dR and D µν ≡ g µν − ∇ µ ∇ ν with ≡ ∇ α ∇ α and ∇ is the usual covariant derivative. These equations may be writtenà la Einstein by isolating on the l.h.s. the Einstein tensor and the f (R) contribution on the r.h.s. as follows [102]: We can also find the expression for the scalar curvature by contracting Equation (9) with g µν , which gives: Note that, unlike GR where R and T are related algebraically, for a general f (R) those two quantities are dynamically related. Thus, in principle, non-null curvature solutions are possible in vacuum for certain f (R) models.
In addition, one may study the possible corrections at local scales of the Newtonian law and the violations of the equivalence principle by using again an auxiliary scalar field A and applying a conformal transformation. Then, we may rewrite the action Equation (8) as: Then, by varying this action with respect to A, it gets R = A, and the action Equation (8) is recovered. By the conformal transformationg µν = e σ g µν , where σ = − ln(1 + f (A)), the action Equation (12) is rewritten in the so-called Einstein frame: where the tilde represents that the scalar curvature is evaluated with respect to the transformed metricg µν , and the potential V (σ) is given by: Here, one has to solve the equation σ = − ln(1 + f (A)) in order to get A = A(σ). This conformal transformation affects the matter sector in the gravitational action by a term proportional to e 2σ , which can imply large corrections to the Newtonian law unless the potential Equation (14) is constructed in order to avoid them at local scales by the so-called chameleon mechanism [103,104], which implies that the mass of the scalar field σ should be enough large at local scales in order to avoid large corrections to the Newtonian law [105]: Hence, if the mass is much larger at local scales than the curvature, the possible corrections to the Newtonian law will be negligible. In this sense, several models that accomplish this rule have been proposed (see [106][107][108][109][110]), which are also capable to reproduce the cosmological history.

III. STATIC AND SPHERICALLY SYMMETRIC BLACK HOLES IN f (R) GRAVITIES
Let us study in this section the f (R) solutions under staticity and spherical symmetry assumptions. The most general D ≥ 4 dimensional metric external to static and spherically symmetric-therefore non-rotating-configurations can be written as (see [111]): (16) or alternatively: where dΩ 2 D−2 is the metric on the S D−2 sphere and identification λ(r) = e −2Φ(r) A(r) and µ(r) = A(r) can be straightforwardly established if required. Since the metric is static, the scalar curvature R in D dimensions depends only on r and it is given, for the metric parametrization Equation (16), by: where the prime denotes derivative with respect to r coordinate.

A. Spherically Symmetric and Static Constant Curvature Solutions: Generalities
At this stage it is interesting to ask about which are the most general static and spherically symmetric metrics with constant scalar curvature R 0 .
This curvature can be found solving the Equation (18) R = R 0 . Thus, it is straightforward to check that, for a constant Φ(r) = Φ 0 , the general solution is [112,113]: with a 1 and a 2 being arbitrary integration constants. In fact, for the particular case D = 4, R 0 = 0 and Φ 0 = 0, the metric can be written exclusively in terms of the function: By establishing the identifications a 1 = −2G N M and a 2 = Q 2 , this solution corresponds to a Reissner-Nordström solution, i.e., a charged massive BH solution with mass M and charge Q. Further comments about this result will be made in Section III B.

B. Spherically Symmetric and Static Constant Curvature Solutions
Let us study here which are the vacuum constant curvature solutions for f (R) theories provided that the metric under consideration is of the form Equation (16)-or alternatively Equation (17). By inserting the metric Equation (16) into the general f (R) gravitational action S g in Equation (8), and making variations with respect to the metric functions, A(r) and Φ(r), the equations of motion become: and 22) or alternatively in by using the metric parametrization given by Equation (17) the equations of motion now are expressed as: where f , f and f denote derivatives of f (R) with respect to the curvature R. In the following calculations in this section, we use the set of Equations (21) and (22). In order to simplify the difficulty of the above equations, let us first consider here the case of constant scalar curvature R = R 0 solutions. Then the equations of motion Equations (21) and (22) reduce to: and As commented in Section II, the constant curvature solutions of f (R) theories in vacuum satisfy: provided that 2(1 + f (R 0 )) = D. Thus from Equation (25) one concludes (excluding pathologic cases) that Φ (r) = 0 and then Equation (26) becomes: This is a f (R)-independent, linear second order inhomogeneous differential equation. This statement must be understood in the sense that Equation (28) does not involve explicitly any f (R) parameter. Nonetheless, it is obvious from Equation (27) that the value of R 0 does depend upon the f (R) model under study. Therefore the metric solution A(r) will also inherit information about the model under consideration. Then, Equation (28) can be easily integrated to give the general solution: which depends on two arbitrary constants C 1 and C 2 . However, two important remarks need to be done at this stage: • Firstly, by comparison with Equation (19), one can see that the term with the power r 2−D is absent. This fact will be studied in Section III C.
• Secondly, this solution has no constant curvature in the general case since, as we found above, the constant curvature requirement demands C 1 = 1. This issue just requires a constant fixing (or equivalently a time reparametrization) and does not affect the solution.
Then, for negative R 0 , the found solution in Equation (29) with C 1 = 1 is basically the D-dimensional generalization obtained by Witten [71] of the BH in AdS space-time solution considered by Hawking and Page [70]. With the natural choice Φ 0 = 0 the solution can be written as: where where µ D−2 is the area of the D − 2 sphere, l 2 ≡ −D(D − 1)/R 0 is the asymptotic AdS space scale squared and M is the mass parameter usually found in the literature. Thus, the main conclusion of this analysis [112,113] is that the only static and spherically symmetric vacuum solutions with constant (negative) curvature of any f (R) gravity is just the Hawking-Page BH in AdS space. However this kind of solution is not the most general static and spherically symmetric metric with constant curvature, as can be seen by comparison with the solutions found in Equation (19). Therefore there exist constant curvature BH solutions that cannot be obtained as vacuum solutions of any f (R) model. In particular, the term with r 2−D power cannot be mimicked by any f (R) model. As we shall show immediately, the most general solution as given by Equation (19) can be described as a charged BH solution in a D = 4 f (R)-Maxwell theory.

C. f (R) Solutions Combined with Electromagnetism
Let us briefly mention here for completeness with the previous results the case of charged BH in f (R) theories that we shall study in detail in Section IV. We shall limit ourselves to the D = 4 case, since otherwise (D = 4) the r 2−D term appearing in Equation (19) cannot be originated by the Maxwell action to be presented below. Let us consider here a f (R) gravitational Lagrangian supplemented with the usual Maxwell action for standard Electromagnetism. Therefore, the considered action represents a generalization of the Einstein-Maxwell action: where as usual Considering an electromagnetic potential of the form A µ = (V (r), 0) and the static spherically symmetric metric Equation (17), one finds that the solution with constant curvature R 0 reads [112,113]: Notice that unlike the EH case, the contribution of the BH charge to the metric tensor is corrected by a (1 + f (R 0 )) −1 factor. This fact is of particular importance in order to constrain f (R) models as astrophysically viable: the sign of the factor 1 + f (R 0 ) needs to be positive in order to guarantee that the solution Equation (33) keeps the same sign as in the Reissner-Nordström solution. Otherwise, for a charge Q the third term in the solution Equation (33) will be negative and therefore this fact becomes experimentally observable (affecting for instance to the geodesic movements around a charged BH). In fact, the requirement 1 + f (R 0 ) > 0 has been usually accepted [114] as one of the gravitational viability conditions for f (R) models and is also related, as will be shown here, with thermodynamical viability.

D. Perturbations Around Schwarzschild-(anti)-de Sitter Solutions
In the previous section we have revised the static spherically symmetric solutions with constant curvature. In EH theory, the constant curvature ansatz Equation (27) provides the most general static and spherical symmetric solution straightforwardly. Without entering into details, this is the well-known Birkhoff-Jebsen theorem result [115,116]. However, it is not transparent for this to also be the case in f (R) theories, i.e., that the most general static and spherically symmetric solution for f (R) theories possesses constant curvature. Furthermore, it seems that a generic solution would behave differently in the Jordan and Einstein frames, both related by a conformal transformation, at least in the perturbative analysis (see [117,118]).
In this section we revise the basic rudiments of the perturbative analysis of the problem originally studied in [112,113], only in the Jordan frame. In that analysis, it was assumed that the modified gravitational Lagrangian is a small perturbation around the EH Lagrangian. Therefore, the arbitrary f (R) function is of the form: where α 1 is a dimensionless parameter and g(R) is assumed to be analytic in α. This last assumption is based on the fact that the scalar curvature R is also α-dependent as will be seen in expressions Equation (35). Therefore g(R) may also be expanded in powers of α. By assumption and for the sake of simplicity, g(R) is assumed analytic in the expansion in powers of α. Now, assuming that the λ(r) and µ(r) functions appearing in the metric Equation (17) are also analytical in α, they can be written as follows: where {λ 0 (r), µ 0 (r)} are the unperturbed solutions for the EH action with cosmological constant, i.e., S(A)dS solution, given by: which are the standard BH solutions in a D-dimensional Schwarzschild-(A)dS space-time. The factor C 2 can be chosen by performing a coordinate t reparametrization so that both functions in Equation (36) can be identified. For the moment, the background solutions are kept as given in Equation (36) since the possibility of getting λ(r) = µ(r) in the perturbative expansion will be discussed later on. The perturbative procedure is as follows: by inserting Equations (34) and (35) into the general Einstein modified Equations (23) and (24), equations for each order in α parameter powers are obtained. Please refer to [112,113] for details.
Notice that from the obtained results up to second order in α, the corresponding metric has constant scalar curvature for any value of the integration constants. As a matter of fact, this metric is nothing but the standard Schwarzschild-(A)dS geometry, and can be easily recast in the usual form by making a trivial reparametrization in the time coordinate.
Therefore, this analysis proved that at least up to second order, the only static, spherically symmetric solutions which are analytical in α are the standard Schwarzschild-AdS space-times irrespectively of the f (R) model under consideration. This result is valid provided that both the solutions and the gravitational Lagrangian are analytic in the α parameter and that g(R) corresponds to a small deviation from the usual EH gravitational Lagrangian. To conclude Sections III and III D, we can summarize by saying that in the context of f (R) gravities the only spherically symmetric and static solutions of negative constant curvature are the standard BH in AdS space. The same result applies in the general case (without imposing constant curvature) in perturbation theory up to second order. However, the possibility of having static and spherically symmetric solutions with non-constant curvature cannot be excluded in the case of f (R) functions, which are not analytical in α parameter.

IV. KERR-NEWMAN BLACK HOLES IN f (R) THEORIES
In this section, we shall focus our attention in obtaining constant curvature R 0 vacuum solutions for fields generated by massive charged objects in the frame of f (R) gravity theories. Hence, the appropriate action to derive the field equations will be now Equation (32) and D = 4 will be assumed throughout this section. The obtained equations (in G N = c = = k B = 1 units) are: In the studied case D = 4, the trace of the previous equation leads again to expression Equation (27) due to the conformal character of the Faraday tensor, i.e., F µ µ = 0. For these f (R) field equations the axisymmetric, stationary and constant curvature R 0 solution that describes a BH with mass, electric charge and angular momentum is analogous to the one found by Carter and published for the first time in 1973 [119]. In Boyer-Lindquist coordinates, the metric takes the form [120,121]: with ∆ r := r 2 + a 2 1 − R 0 12 r 2 − 2M r + Q 2 (1 + f (R 0 )) ρ 2 := r 2 + a 2 cos 2 θ ; ∆ θ := 1 + R 0 12 a 2 cos 2 θ ; Ξ := 1 + R 0 12 a 2 (39) where M , a and Q denote the mass, spin and electric charge parameters respectively. Notice that, unlike in the GR case, the contribution of the charge of the BH to the metric is corrected by a (1 + f (R 0 )) −1/2 factor. This feature was already obtained for Reissner-Nordström BH presented in expression Equation (33) and originally in [112,113]. On the other hand, the required potential vector and electromagnetic field tensor in Equation (37) solutions for metric Equation (38) are respectively: To lighten notation, from now on let us use the notationQ 2 ≡ Q 2 / (1 + f (R 0 )) as a normalized electric charge parameter of the BH.

A. Event Horizons
Let us now briefly study the horizon structure of these BH: according to the event horizon definition, i.e., g rr = 0, these hypersurfaces are found as the roots of the equation ∆ r = 0, i.e., This is an algebraic fourth order equation that, following [120,121], can be rewritten as: where • r − is always a negative solution with no physical meaning, • r int and r ext are the interior and exterior horizon respectively, and • r cosm represents-provided that it arises-the cosmological event horizon for observers between r ext and r cosm .
The latter divides the accessible region from the hidden region for an exterior observer. The Equation (42) can be solved by using L. Ferrari's method for quartic equations. The existence of real solutions for this equation is given by a factor h, which we shall denote as horizon parameter: In Figure 1 we show the zeros of the fourth order polynomial ∆ r that lead to different astrophysical configurations, depending on the values and signs of R 0 and h.
Let us stress at this stage that, from a certain positive value of the curvature R crit 0 onward, the h factor goes to zero for two values of a. They can be understood as follows: • Upper spin bound, a = a max for which the BH turns extremal-the interior and exterior horizons have merged into a single horizon with a null surface gravity. This is the usual configuration for the BH to become extremal.
• Lower spin bound, a = a min , below which the BH turns into a marginal extremal BH. This value can be understood as the cosmological limit for which a BH preserves its exterior horizon without being "torn apart" due to the relative recession speed between two radially separated points induced by the cosmic expansion [120,121].
In Figure 2 we revisit the original plot shown in [120,121], where for certain values of the electric charge parameter Q the allowed ranges of the spin a parameter values are presented.

V. BLACK HOLES THERMODYNAMICS IN f (R) THEORIES
As mentioned in the Introduction, the study of BH thermodynamics has drawn a lot of attention when studying these configurations in alternative theories of gravity. We shall divide this section in two main parts: the study of thermodynamics in AdS case and then the study of KN case as the natural generalization of the former. This approach will enable to introduce first the rudiments of the Euclidean method so the reader can familiarize with it, and then extend the formulation for more general objects such as the KN BH configurations. Finally, we will provide some examples by considering f (R) models that illustrate the procedure for both configurations. In the AdS case, the dimension of the space-time will be arbitrary and then fixed to several values in order to illustrate the found results. On the contrary, for the KN study, dimension will be fixed from the very beginning to D = 4 for the sake of simplicity.

A. BH Thermodynamics for AdS Configuration
In order to study the different thermodynamic quantities for the f (R) BH in AdS space, the temperature is usually considered as the basic quantity from which the rest of the thermodynamic quantities are derived. In principle, there are two different approaches of introducing this quantity for the AdS solutions. Firstly, one of the possible definitions is the one provided by the Euclidean quantum gravity [122][123][124][125] that was revised in [112,113] yielding: Another possible definition of temperature was firstly proposed in [126] stating that temperature can be given in terms of the horizon gravity K as T K ≡ K 4π . Then, by reminding the metric definition Equation (16), it is straightforward to obtain that T K = T E . In conclusion, both definitions lead to the same result for this kind of solution. Notice also that in any case the temperature depends only on the behavior of the metric near the horizon and it is independent of the underlying gravitational action. This last statement can be understood by remarking that either different actions or different field equations giving rise to the same solutions of the form Equation (16) will also present the same temperature. This is not the case for other thermodynamic quantities as will be seen later. Taking into account the results in previous sections, the usual analyses for f (R) theories [112,113] have focused on the constant curvature AdS BH solutions with Φ = 0 as a natural choice, and A(r) as given by expression Equation (30). Then, both definitions of temperature lead to: Notice that for a given dimension D, the temperature is a function of r H only, i.e., it depends only on the BH size. In the limit r H going to zero the temperature diverges as T ∼ 1/r H and for r H going to infinite T grows linearly with r H . Consequently T has a minimum corresponding to a temperature: The existence of this minimum was established in [70] for D = 4 by Hawking and Page. More recently, Witten extended this result to higher dimensions [71]. This minimum will be relevant in order to set the regions with different thermodynamic behaviors and stability properties. Explicit solutions for A(r H ) = 0 in Equation (30) can be found for D = 4, 5 and were found in [112,113].
Another important temperature value is the one corresponding to a horizon radius r H = l. Let us denote it as T 1 and from Equation (45) it is given by: Notice that for D > 2 we have T 0 < T 1 .
In order to compute the remaining thermodynamic quantities, the standard procedure is to consider the Ddimensional Euclidean action defined as follows: When the previous expression is evaluated on some metric with a periodic Euclidean time of period β, it equals β times the free energy F associated to this metric. The computation by Hawking and Page [70], generalized to higher dimensions by Witten [71], was extended to f (R) theories in [72] and [112,113]. In these references the difference of this action, when it is evaluated on the BH and on the AdS metric, was computed leading to: where R 0 = −D(D − 1)/l 2 . From this expression it is straightforward to obtain the free energy as F = ∆S E /β. One can also see that provided that the condition −(R 0 + f (R 0 )) > 0 is fulfilled-which is the usual case in EH gravity-we have F > 0 for r H < l and F < 0 for r H > l.
Once that both temperature and free energy have been calculated, the total thermodynamical energy may now be obtained as: where M is the mass defined in Equation (31). This is one of the possible definitions for the BH energy for f (R) theories but other definitions are available in the literature, see for instance [127] for a more general discussion.
For the EH action with non-vanishing cosmological constant, we have f (R) = −(D − 2)Λ D and then it is immediate to find E = M . However this is not the case for general f (R) gravitational Lagrangians as seen in Equation (50). Notice that positive energy in AdS space-time requires Now the entropy S can be obtained from the well-known relation S ≡ βE − βF . Then one gets: Notice that once again positive entropy requires R 0 + f (R 0 ) < 0. For the EH action non-vanishing cosmological constant, we have R 0 + f (R 0 ) = −2(D − 1)/l 2 and then the famous Hawking-Bekenstein result [128] is recovered: Finally we can compute the heat capacity C which can be written as: Then it is easy to find: For the particular case of the EH action with cosmological constant we find: In the Schwarzschild limit, i.e., l going to infinity, this formula turns into: which is the well-known result for standard Schwarzschild BH, i.e., C < 0. Provided that the condition (R 0 + f (R 0 )) < 0 is accomplished, the f (R) general case lead, like in the EH case, to C > 0 for r H > r H0 (the large BH region) and C < 0 for r H < r H0 (the small BH region). For r H ∼ r H0 (T close to T 0 ) C is divergent. Notice that in EH gravity, C < 0 necessarily implies F > 0 since T 0 < T 1 .
In summary, f (R) theories accomplishing the condition: provide a scenario analogous to the one described in full detail by Hawking and Page in [70] for the EH case. This analysis was originally shown in [112,113] and then widely used in literature.
In principle, particular f (R) models violating the condition R 0 + f (R 0 ) < 0 can be naively considered. However, in such cases, the mass and the entropy would be negative and therefore the AdS BH solutions for these models would be unphysical. Therefore R 0 + f (R 0 ) < 0 can be regarded as a necessary condition for f (R) theories in order to support the usual AdS BH solutions. Using Equation (27), this condition is fully equivalent to 1 + f (R 0 ) > 0. This last condition has a clear physical interpretation in f (R) gravity theories (see [114] and references therein). Indeed, it can be interpreted as the condition for the effective Newton constant G ef f = G D /(1 + f (R 0 )) to be positive. It can also be interpreted from the quantum point of view as the condition which prevents the graviton from becoming a ghost. Therefore the requirement Equation (57) can be also regarded as a thermodynamical viability condition for f (R) gravity theories when AdS configurations are studied [112,113].
Once that the main thermodynamical quantities have been derived, let us now revise the local and global stabilities phenomenology: For T < T 0 , the only possible state of thermal equilibrium in an AdS space is pure radiation with negative free energy and hence there is no stable BH solutions. For T > T 0 there exist two possible BH solutions; the small (and light) BH and the large (heavy) BH. The small one has negative heat capacity and positive free energy as the standard Schwarzschild BH. Therefore it is unstable under Hawking radiation decay. For the large BH there appear two possibilities; if T 0 < T < T 1 then both the heat capacity and the free energy are positive and the BH will decay by tunneling into radiation, but if T > T 1 then the heat capacity is still positive but the free energy becomes negative. In this case the free energy of the heavy BH will be less than that of pure radiation. Then pure radiation will tend to tunnel or to collapse to the BH configuration in equilibrium with thermal radiation. This phenomenology will be applied to different f (R) models in Section VI.

B. BH Thermodynamics for KN Configuration
Following the reasoning at the beginning of this section, the research concerning the KN configuration will be also done for BH with a well-defined horizon structure for negative values of R 0 .
In order to study the different thermodynamical properties of Kerr-Newman BH in f (R) theories, the usual approach is to start studying the temperature of the exterior horizon r ext ≡ r ext R 0 , a,Q, M .
For that purpose, the Euclidean action method [122][123][124][125] is again considered. The procedure consists now of performing the change of variables t → −iτ , a → ia on the metric (38) in order to obtain the Euclidean section (see [120,121] for details). After some calculations, the inverse Hawking temperature now becomes: For null spin and charge, this expression becomes naturally Equation (45).
BH horizon temperature can also be obtained through Killing vectors, as is explained in [126] where temperature is defined as: T κ, KN ≡ κ/4 π, with κ the surface gravity defined now as χ µ ∇ µ χ ν = κ χ ν . It can be verified that κ is the same at any horizon point and consequently T κ, KN = T E, KN as obtained in [129]. Now that we know the expression for the temperature for KN configurations, the natural step in order to obtain the remaining thermodynamical quantities is to consider the corresponding extension of the Euclidean action with respect to the previously considered Equation (48). It now becomes: with Y the integration region. As is described in [70], and fully derived in [120,121] the result for ∆S E, KN becomes: where Φ e is the electric potential of the horizon as seen from infinity: and Q is the physical electric charge of the BH, obtained by integrating the flux of the electromagnetic field tensor at infinity, which happens to be: From expressions Equations (61) and (62), it can be seen that Φ e and Q definitions are f (R)-model independent.
The reason for this lies in the fact that these calculations involve the vector potential and the electromagnetic field tensor, which are independent of the considered f (R) model. Further analysis of the action revealed [120,121] that it became singular for h = 0, as could be expected from extremal BH, whose temperature T E, KN = 0 makes the β KN factor diverge. Since thermodynamical potentials are obtained by dividing the action by the β factor, they still remain well-defined at T E, KN = 0. By using the expression Equation (60) we can immediately obtain Helmholtz free energy F . For rotating configurations, the correct definition is now: where the term Ω H J comes from the required Legendre transformation to fix angular momentum, being J the angular momentum of the BH and Ω H the angular velocity of the rotating horizon. For the explicit expressions, we refer the reader to [120,121]. Thus, performing some calculations one gets: Once again, expression Equation (64) reduces to Equation (49) for a non-rotating uncharged BH. Provided that the condition 1 + f (R 0 ) > 0 is required to hold in order to obtain positive values of the mass, the inspection of the numerator of F renders that F > 0 for values of the horizon below r limit ext (with an associated mass M limit through Equation (41)), and F < 0 for larger values.
Using the appropriate thermodynamical relations [130], we can derive the entropy S of the KN BH, which reads: If the area A H of the exterior horizon is now computed, one gets: Thus, if expression Equation (66) is substituted into the entropy expression Equation (65), the entropy can be expressed as: Consequently 1 + f (R 0 ) > 0 is again a mandatory condition to obtain a positive entropy, as was shown in the AdS as well. Note that if for example we want to recover GR with a cosmological constant taking f (R 0 ) = −2Λ, this relation turns into: S = A H /4, i.e., the famous Bekenstein result [128].
Once that the temperature T and the entropy S of the KN BH are known, the natural step further can be taken and study the heat capacity at constant scalar curvature R 0 and at fixed spin a and chargeQ parameters. From the definition: we obtain the expression: with Den(Q, r ext , R 0 , a) ≡ Den defined as: At this stage, it is interesting to find out for which values of R 0 , a,Q and M the denominator of the thermal capacity becomes zero, i.e., the thermal capacity goes through an infinite discontinuity, i.e., the BH experiences a phase transition. Following the analysis presented in [120,121], two kinds of BH on this subject can be distinguished depending on the values of the a,Q and M parameters and scalar curvature R 0 : • fast BH, without phase transitions and always positive heat capacity C > 0.
• slow BH, presents two phase transitions for two determined values of r ext .
In Figure 3  As was remarked in [120,121], KN-AdS BH, unlike Schwarzschild-AdS BH case [70], are allowed for any value of the temperature T . Hence stability of each BH will be exclusively given by the corresponding values of heat capacity C and free energy F as functions of the a,Q, M and R 0 parameters, that ultimately define the BH. However, for a set of fixed values of a,Q and R 0 , the mass parameter must be bigger than a minimum M min (characterized by T = 0) to host BH configuration, otherwise radiation is the only possible equilibrium state up to such a minimum mass. For bigger masses, we shall distinguish between the fast and the slow BH.
• Fast BH, with bigger values of the spin and the electric charge than the slow ones, shows a heat capacity always positive and a positive free energy up to a M = M limit value, and negative onwards. Thus, this BH is unstable against tunneling decay into radiation for mass parameter values of M < M limit . For M > M limit , free energy becomes negative, therefore smaller than that of pure radiation, that will tend to collapse to the BH configuration in equilibrium with thermal radiation.
• Slow BH shows a more complex thermodynamics, being necessary to distinguish between four regions delimited by the mass parameter values: M min < M I < M II < M limit , as follows: -For M min < M < M I and for M II < M < M limit , both the heat capacity and the free energy are positive, which means that the BH is unstable and decays into radiation by tunneling. -If M I < M < M II , the heat capacity becomes negative but free energy remains positive, being therefore unstable and decays into pure thermal radiation or to larger values of mass. -Finally, for M > M limit the heat capacity is positive whereas the free energy is now negative, thus tending pure radiation to tunnel to the BH configuration in equilibrium with thermal radiation.
Let us conclude this section as we did the previous one, by emphasizing that, although not quantitatively, the thermodynamical behavior of KN-AdS BH in f (R) modified gravity theories is qualitatively analogous to that of GR [130]. This was an original important result presented for the first time in [120,121]. In that investigation, authors proved that some thermodynamical quantities like the physical mass M, the angular momentum J and the entropy S differ from those of GR by a (1 + f (R 0 )) factor. One might think that the free energy F and the heat capacity C also differ in that same multiplicative way, but because of the dependence with the normalized charge parameterQ, this same term (1 + f (R 0 )) leads to a more complex deviation from the GR behavior.

VI. THERMODYNAMICS IN AdS AND KERR-NEWMAN: PARTICULAR EXAMPLES A. Thermodynamics in AdS
In order to illustrate the calculations sketched in Section V A, let us consider some particular f (R) models in whose heat capacity C and the free energy F will be determined that were originally studied in [112,113]. As was pointed out above, those two quantities provide the local and global stability of AdS BH. For the particular models that are studied below, the constant scalar curvature R 0 could be calculated exactly by using Equation (27). For the sake of simplicity the D-dimensional Schwarzschild radius in Equation (31) was taken as R D−3 S ≡ 2. The models we have considered are the following.
This model has been widely studied because the αR 2 term with α > 0 can account for the accelerated expansion of the Universe [131]. This model can also explain the observed temperature anisotropies observed in the CMB, and can become a viable alternative to scalar field inflationary models, since reheating after inflation would have its origin on the production of particles during the oscillation phase of the Ricci scalar [132]. For β = 2, this model was recently pointed out [16,18] as a possible explanation to dark matter.
Substituting this model in Equation (27) for arbitrary dimension and by only considering non-vanishing curvature solutions, the found solution is: For D > 2, the condition (2β − D)α < 0 provides well-defined scalar curvatures R 0 . Two separated regions have thus to be studied: Region 1 {α < 0, β > D/2} and region 2 {α > 0, β < D/2}. For this model we also get: Notice that in region 1, 1 + f (R 0 ) > 0 for D > 2, since in this case β > 1 is straightforwardly accomplished. In region 2, we find that for D > 2, the requirement R 0 + f (R 0 ) < 0, i.e., 1 + f (R 0 ) > 0, fixes β < 1, since this is the most stringent constraint over the parameter β in this region. Therefore the physical space of parameters in region 2 is restricted to be {α > 0, β < 1}. In Figures 4 and 5, taken from [112,113], one can see the physical regions in the parameter space (α, β) corresponding to the different signs of (C, F ).   This model has been proposed in [106] as cosmologically viable and has attracted a lot of attention in order to constrain its validity in different astrophysical and cosmological scenarios. Throughout this study, let us consider n = 1 for the sake of simplicity. A vanishing curvature solution also appears in this model and the two nontrivial curvature solutions are given by [112,113]: where κ > 0 and κ > (8D − 16)/D 2 are required for real R 0 solutions. The corresponding 1 + f (R 0 ) values for Equation (72) are: Since 1 + f (R 0 ) > 0 is required, one has sign(R ± 0 ) = sign(αβ). In fact it turns out that 1 + f (R − 0 ) is not positive for any allowed value of κ and therefore this curvature solution R − 0 is excluded for our study (see [112,113] for details). On the other hand, 1 + f (R + 0 ) > 0 only requires κ > 0 for dimension D ≥ 4 and therefore f (R 0 ) < 0. Therefore only two accessible regions need to be studied: region 1 {α > 0, β < 0} and region 2, {α < 0, β > 0}.

B. KN-AdS Thermodynamics
In this section we will study the thermodynamics of the two f (R) models presented in the last section but when a Kerr-Newman configuration is considered.
First, as was shown above, the condition 1 + f (R 0 ) > 0 must be satisfied in order to obtained physical results, so first of all the viability of each model depending upon the values of the parameters that define them will be checked out. Thus, for each model, we will study the range of parameters that allows the existence of Kerr-Newman BH, i.e., we will graphically display a max and a min (if present) in terms of those parameters, and study the region of a confined between these two surfaces. This phenomenology was absent in the SAdS study since the configuration was spherically symmetric and therefore a = 0. Finally we will focus on the thermodynamical quantities that define BH stability depending again upon the model range of parameters. In order to sketch graphically the thermodynamical behavior, let us then consider the BH parameter values: M = 1, a = 0.4 andQ = 0.2. Let us remind at this stage that this analysis is restricted to R 0 < 0 condition, as explained in the beginning of Section V.
For the sake of simplicity, let us introduce the dimensionless variables notation [120,121] that we will use in the graphical visualization of the results to be presented: where the previous symbols keep their original meaning as defined before. The considered models are as follows.
Let us remind the solution Equation (70) found above for this model. The viability condition 1 + f (R 0 ) > 0 restricts the range of parameters that define this f (R) model.
In Figure 6, we show the range of the spin parameter a for which BH are allowed, depending on the parameters α and β and for charge parameters values Q = 0 and Q = 0.75. For the solution R 0 < 0 there exist two possible regions: • Region 3 {α < 0, β > 2}. For this region, the value of a max decreases suddenly from its normal value to 0 near the frontier of the region, i.e., α ≈ 0 and β ≈ 2; other values of α and β display a relatively low curvature and the existence of BH is assured.
• Region 4 {α > 0, β < 1}. This region only shows problems when 0 < β < 1, where a max → 0, but, as β becomes more negative, the surface of a max slowly acquires higher values, recovering its usual value for β = −2. We have graphically schematized in Figure 7 the possible thermodynamical configurations of a BH for Regions 3 and 4, i.e., R 0 < 0 following the results presented in [120,121]. As we can see in Region 3, only for values close to the region boundary (α ≈ 0 and β ≈ 2) is BH stable both locally and globally (C > 0 and F < 0). Also, as the absolute value of α increases, it is more probable that the BH shows global instability (C > 0 and F > 0) and finally global and local instability (C < 0 and F > 0). In Region 4, however, we find the opposite behavior: near the limit of the region values (α ≈ 0 and β ≈ 1) we have global and local instability, and, as α gets bigger, the BH becomes more stable both locally and globally. Let us consider for the sake of simplicity again n = 1, thus having a two-parameters model, since we can define γ = β/α and then obtain: Replacing the latter in Equation (27), one gets Equation (72), that for D = 4, and focusing on the negative constant curvature solution, yields: Keeping in mind that we have to satisfy 1 + f (R 0 ) > 0, κ happens to be restricted to values of κ > 1. On the other hand, computation of 1 + f (R 0 ) reveals that R 0 is only a valid solution for values of κ and γ in the Region 2: {κ > 1, γ < 0}, being as mentioned R 0 < 0. In Figure 8 we show the range of the spin parameter a for which BH are allowed, depending on the parameters κ and γ for certain values of the chargeQ. Region 2 displays an a max surface that descends slightly as the value of κ increases, and more significantly when γ ≈ 0, becoming 0 when this limit is reached. We graphically schematize in Figure 9 the different possible thermodynamical configurations for the region in which R 0 < 0, in this case the region 2. In this parameter space region, a combination of high values of κ and small absolute values of γ favor BH local and global stability, and, conversely, a combination of small values of κ and high absolute values of γ lead to BH instability, first globally and then also locally.

VII. COSMOLOGICAL SOLUTIONS IN MODIFIED GRAVITY
Modified gravity has become very popular in recent years because it is capable to reproduce the dark energy epoch, and even the inflationary phase, with no need to introduce an additional component in the universe. Hence, the reconstruction of appropriate theories that contain some particular cosmological solutions have involved an important effort (see [37, 38, 42-46, 101, 133-149]). The same is true for the reconstruction of so-called viable modified gravities, which can recover GR in some limit while they introduce additional effects at cosmological scales (see [106][107][108][109][110]). However, as in every dark energy model, the required fine tuning on the form of the action as on free parameters of the model have resulted in a large effort in searching additional and independent predictions in the frame of modified gravities, including other aspects different than cosmology as the ones studied in the above sections (see also the reviews [26][27][28][29][30][31][32][33] and references therein). In this section we are interested to make a short review on cosmological solutions in the frame of modified gravities, and in particular in F (R) gravity, while the thermodynamical aspects are studied in the following sections. As it is widely accepted, our universe seems to be well described by FLRW metrics: and specifically the universe seems to be spatially flat, k = 0. Hence, by introducing the metric in F (R) field Equation (3), modified FLRW equations are obtained: Recent astronomical data indicates that the well-known ΛCDM model where DE is represented by an effective cosmological constant fits quite well the observational data. Nevertheless, it is shown that F (R) gravity can even reproduce exactly the ΛCDM model in absence of a cosmological constant s shown below (see [37,150]). Let us consider ΛCDM cosmology: where we have introduced the number of e-foldings N = ln a a0 instead of the cosmological time t, and we have assumed a pressureless fluid. We can easily solve the continuity equation for ρ m : while the first FLRW equation in Equation (78) can be rewritten in terms of the number of e-foldings, which yields [148]: where G(N ) = H 2 . Hence, for ΛCDM cosmology Equation (79), the scalar curvature is given by R = 3G (N ) + 12G(N ) = 12H 2 0 + κ 2 ρ 0 a −3 0 e −3N , and the Equation (81) gives a hypergeometric equation [37]: whose solution is given by: where x = R Hence, the action Equation (83) is capable to reproduce ΛCDM model with no cosmological constant. Nevertheless, the action Equation (83) is composed by hypergeometric functions, which do not seem to provide a successful explanation for dark energy, since the theory turns out much complicated than the usual GR with a cosmological constant. However, other viable models of F (R) gravity capable to mimic a cosmological constant at the present time have been suggested in the literature (see [106][107][108][109][110]), which in general seem to cross the phantom barrier [151]. Then, let us now study a model where the dominant component is phantom-like. Such kind of system can be easily expressed in the standard General Relativity when a phantom fluid is considered, where the FLRW equation reads H 2 (t) = κ 2 3 ρ ph (the subscript ph denotes the phantom nature of the fluid). As the EoS for the fluid is given by p ph = w ph ρ ph with w ph < −1, by using the conservation equation, the solution yields a(t) = a 0 (t s − t) −H0 , where a 0 is a constant, and t s is the so-called Rip time. Then, the solution describes a universe that ends in a Big Rip singularity at the time t s . The same behavior can be achieved in F (R) theory with no need to introduce a phantom fluid. The expression for the Hubble parameter as a function of the number of e-foldings is given by H 2 (N ) = H 2 0 e 2N/H0 . Then, the Equation (81), with no matter contribution, takes the form: . This equation is an Euler equation whose solution yields: Thus, the phantom cosmology a(t) = a 0 (t s − t) −H0 can be also obtained in the frame of F (R) theory and no phantom fluid is needed. In addition, the equivalence between f (R) gravity and scalar-tensor theories Equation (4) can be used to easily obtain cosmological solutions in f (R) gravity. For the metric Equation (77), Friedmann equations read: We redefine the scalar field such that it is chosen to be the time coordinate φ = t. Hence, taking into account the Equations (87) and the continuity equation for the perfect fluid ρ m , the Hubble parameter may be calculated as a function of the scalar field φ, H = g(φ). By combining the Equation (87), the function Q(φ) is deleted, and it yields: By solving this equation for a given function P (φ), a cosmological solution H(t) is found, and the function Q(φ) is obtained by means of equations in Equation (87): If we neglect the contribution of matter, then the Equation (88) is a first order differential equation on g(φ), and it can be easily solved. The solution found is the following [147]: where C is an integration constant. Hence, we have shown that F (R) gravity represents a serious candidate to explain the phenomena of late-time acceleration in the expansion of the universe. In the following sections, we explore the thermodynamical properties of FLRW metrics in the frame of modified gravity as well as an attempt to generalize Cardy-Verlinde formula.

VIII. FIRST LAW OF THERMODYNAMICS AND FLRW EQUATIONS
In this section, we review some of the most important results obtained in the study of FLRW metrics and their thermodynamics, specifically the relation between the first law of thermodynamics and FLRW equations, which can be deduced by assuming the well known relation between the entropy and the horizon area, as well as the relation between the gravity surface and the temperature [82][83][84]. Such study has been extended not only to f (R) gravity, as is also shown below (see [85]), but also to many other different scenarios [86][87][88][89][90][91][92][93][94][95][96]. In fact, it was found in [80] and extended in [81] for f (R) gravities, a direct deduction of gravity field equations by assuming the relation of the horizon area and the entropy. Nevertheless, here we are only interested to show the deduction of FLRW equation from the thermodynamics defined in the corresponding horizons of a particular spacetime. Let us consider a n + 1 FLRW metric Equation (77), with arbitrary spatial curvature k = −1, 0, 1. This metric may be rewritten as: where the indexes a, b = 0, 1, being x 0 = t, x 1 = r, andr = a(t)r, while dΩ 2 n−1 expresses the line element of an n − 1 unit sphere. Then, the apparent horizon for this general kind of spacetimes is defined as h ab ∂ ar ∂ br = 0, whose radius yields:r As will be shown below, one can assign an entropy proportional to the area of the apparent horizon as well as an energy flux through the horizon. This is the reason for considering the apparent horizon instead of other horizons defined in FLRW spacetimes. We may define the work density and the energy flux vector as [82][83][84]: where the work w refers to the change experimented by the apparent horizon, while the vector Ψ a refers to the flux of energy through the horizon. Then, the variation of the internal energy can be written as: where A is the area and V is the volume of an n-dimensional space with radiusr. This equation represents the unified first law according to [152]. The first term in the r.h.s. in Equation (95) takes the form: ∇A +r n−2 ∇ Ẽ r n−2 (96) where κ is the surface gravity: Note that the energy is defined as: Hence at the apparent horizon, the second term in the r.h.s. in Equation (96) vanishes, as on the black hole event horizon, so that the expression Equation (96) can be identified with the thermodynamical relation dE = −T dS if the surface gravity is assumed to be the temperature κ/2π, and the entropy yields the usual expression, proportional to A/4G. Now if we assume a perfect fluid filling the universe, with an energy-momentum tensor given by T µν = (ρ + p)u µ u ν + pg µν , the energy flux vector yields: Then, the change of the internal energy through the apparent horizon Equation (96) becomes: Identifying the entropy and temperature associated to the apparent horizon, and using the first law of thermodynamics dE = −T dS, the second FLRW equation is obtained: By combining it with the continuity Equation (80), the first FLRW equation is also obtained. Hence, FLRW equations can be deduced by using only thermodynamical arguments and identifying the entropy of the apparent horizon and its temperature, which is an identity that supports the idea of the relation between gravitational field equations and thermodynamical ones, as suggested by Jacobson in [80].
In the same way, this analysis can be easily extended to f (R) gravity as well as to other modified theories of gravity. Specifically, for F (R) gravity, the entropy relation is defined (see previous sections) as follows: Hence, following the same steps performed above for GR, the FLRW equations for f (R) gravity Equation (78) are recovered [85]. Then, a strong relation between thermodynamics of the horizons and gravitational field equations is established. In the next section, another aspect of the relation of FLRW equations and thermodynamics is explored, in that case following the holographic principle through the entropy bounds and the Cardy-Verlinde formula.

IX. GENERALIZATION OF CARDY-VERLINDE FORMULA
This section is devoted to consideration of the Cardy-Verlinde formula for general scenarios, and in particular for modified gravity.
The Cardy-Verlinde formula was established by Verlinde in [98] as a thermodynamical relation between FLRW equations with a conformal fluid (a radiation-like fluid) and the Cardy formula (see [97]), and it was extended for a general perfect fluid in [99].
In [100], attempts to reproduce the Cardy-Verlinde formula was suggested for different scenarios, including modified gravity. Other aspects of the Cardy-Verlinde formula such as the presence of a viscous fluid have been studied in [153][154][155]. In this section, we mainly review those results. We consider a (n + 1)-dimensional spacetime described by the FLRW metric Equation (77). By inserting the metric Equation (77) in the Einstein field equations, the FLRW equations for GR are derived: Here ρ i = E i /V and p i are the energy-density and pressure of the matter component i that fills the Universe. In this section, we consider only the k = 1 closed Universe. Moreover, we assume an equation of state (EoS) of the form p i = w i ρ i with w i constant for each fluid, and assume no interaction between the different components. Then, the continuity Equation (80) has the form nowρ i + nH (ρ i + p i ) = 0 and by solving it, we find that the i fluid depends on the scale factor as: Let us now review the case of [99], where just one fluid with EoS p = wρ and w = constant is considered. The total energy inside the comoving volume V , E = ρV , can be written as the sum of an extensive part E E and a subextensive part E C , called the Casimir energy (for other ways to account the Casimir energy in a FLRW universe, see [156,157]), and takes the form: Under a rescaling of the entropy (S → λS) and the volume (V → λV ), the extensive and subextensive parts of the total energy transform as: Hence, by assuming that the Universe satisfies the first law of thermodynamics, the term corresponding to the Casimir energy E C can be seen as a violation of the Euler identity according to the definition in [98]: Since the total energy behaves as E ∼ a −nw and by the definition Equation (106), the Casimir energy also goes as E C ∼ a −nw . The FLRW universe expands adiabatically, dS = 0, so the products E C a nw and E E a nw should be independent of the volume V , and be just a function of the entropy. Then, by the rescaling properties Equation (107), the extensive and subextensive part of the total energy can be written as functions of the entropy only [99]: Here α and β are undetermined constants. By combining these expressions with Equation (106), the entropy of the Universe is written as a function of the total energy E and the Casimir energy E C [99]: which for w = 1/n (radiation-like fluid) reduces to [98]: which has the same form as the Cardy formula given in [97]. The first FLRW Equation (104) can be rewritten as a relation between thermodynamics variables, and yields: It is easy to check that for the bound proposed in [98], E C ≤ E BH , the equation for the entropy Equation (111) coincides with the first FLRW Equation (112) when the bound is reached. We will see below that when there are several fluid components, the same kind of expression as in [98] cannot be found as well as its extension to modified gravity, nor is there the same correspondence with the FLRW equation when the bound is saturated.

A. Multicomponent Universe
If m fluids are considered with arbitrary EoS, p i = w i ρ i , the expression for the total entropy is simple to derive just by following the same method as above. The total entropy is given by the sum of the entropies for each fluid [100]: This expression cannot be reduced to one depending only on the total energy unless very special conditions on the nature of the fluids are assumed. Let us for simplicity assume that there are only two fluids with EoS given by p 1 = w 1 ρ 1 and p 2 = w 2 ρ 2 , w 1 and w 2 being constants. We can substitute the fluids by an effective fluid described by the EoS: and p eff = 1 2 (p 1 + p 2 ), ρ eff = 1 2 (ρ 1 + ρ 2 ). Then, by using the energy conservation Equation (80), we find ρ 1 ∼ (a/a 0 ) −n(1+w1) and ρ 2 ∼ (a/a 0 ) −n(1+w2) , where a 0 is assumed to be the value of the scale factor at the time t 0 . The effective EoS parameter w eff can be expressed as a function of the scale factor a(t): The total energy inside a volume V becomes: As the energy is proportional to two different powers of the scale factor a, it is not possible to write it as a function of the total entropy only. As a special case, if the EoS parameters are w 1 = w 2 = w eff , the formula for the entropy reduces to Equation (110), and coincides with the CV formula when w eff = 1/n.
As another case, one might consider that for some epoch of the cosmic history, w 1 w 2 . Taking also a >> a 0 , we could then approximate the total energy by the function E T ∝ a −nw2 . From Equation (106) the Casimir energy would also depend on the same power of a, E C ∝ a −nw2 . The expression Equation (110) is again recovered with w = w 2 .
Thus in general, when a multicomponent FLRW Universe is assumed, the formula for the total entropy does not resemble the Cardy formula, nor does it correspond to the FLRW equation when the Casimir bound is reached. It becomes possible to reconstruct the formula Equation (111), and establish the correspondence with the Cardy formula, only if we make specific choices for the EoS of the fluids.

B. Inhomogeneous EoS Fluid and Modified Gravity
Let us now explore the case of an n + 1-dimensional Universe filled with a fluid satisfying an inhomogeneous EoS. This kind of EoS, generalizing the perfect fluid model, has been considered in several papers as a way to describe effectively the DE (see [158,159]). We assume an EoS expressed as a function of the scale factor: This EoS fluid could be taken to correspond to modified gravity, or to bulk viscosity [158,159]. By introducing Equation (117) in the energy conservation Equation (80) we obtain: Here we have performed a variable change t = t(a) such that the prime over ρ denotes derivative with respect to the scale factor a. The general solution of this equation is: and C is an integration constant. As shown above, only for some special choices of the functions w(a) and g(a) can the formula Equation (111) be recovered. Let us assume, as an example, that w(a) = −1 and g(a) = −a m , with m being constant. Then, the energy density behaves as ρ ∝ a m . Hence, by following the same steps as described above, the extensive and subextensive energy go as a m+n , and by imposing conformal invariance and the rescaling properties Equation (107), we calculate the dependence on the entropy to be: The expression for the entropy is easily constructed by combining these two expressions and substituting the extensive part by the total energy. This gives us the same expression as in Equation (110) with w = −(n + m)/n. Note that for m = −(1 + n), the formula Equation (111) is recovered and also its correspondence with the CFT formula. However for a generic power m, the CV formula cannot be reconstructed, like the cases studied above. Only for some special choices does the correspondence work, leading to the identification between the FLRW equation and the Cardy formula.
In the same way, we may reconstruct Cardy-Verlinde formula for modified gravity, and specifically for f (R) gravity. Note that f (R) gravity can be seen as GR plus an inhomogeneous fluid. For a closed 3 + 1 FLRW Universe, the modified FLRW equations are expressed now as: where primes denote derivatives respect to R and dots with respect to t. These equations can be rewritten in order to be comparable with those of standard GR, similarly as in Equation (78), by shifting terms to the right side of Equation (121). For such a propose, the geometric terms can be presented as an effective energy-density ρ F (R) and a pressure p F (R) : Then, an EoS for the geometric terms can be defined as p F (R) = w F (R) ρ F (R) . We can define an effective energy-density ρ = ρ m /F (R) + ρ F (R) and pressure p = p m /F (R) + p F (R) . Hence, for some special cases the formula for the entropy developed above can be obtained in F (R)-gravity (for an early attempt deriving a CV formula in a specific version of F (R)-gravity, see [160]). For example, for an F (R) whose solution gives ρ ∝ a −3(1+w eff ) , the formula for the entropy Equation (110) is recovered although in general, as in the cases studied above, no such expression can be given. On the other hand, one could assume that the geometric terms do not contribute to the matter sector. Supposing a constant EoS matter fluid, the expression for the entropy is given by Equation (110), although the cosmic Cardy formula Equation (112) does not have the same form, and in virtue of the modified first FLRW Equation (121) the form of the Hubble entropy S H , the total energy E and the Bekenstein energy E BH will be very different. It is not easy to establish correspondence between two such approaches. Note that using the effective fluid representation the generalized CV formula may be constructed for any modified gravity.
Let us now consider de Sitter space solution in F (R)-gravity (for review of CV formula in dS or AdS spaces, see [161,162]). As was pointed in [69], almost every function F (R) admits a de Sitter solution. This can be easily seen from the first FLRW equation in Equation (121). A de Sitter solution is given by a constant Hubble parameter H(t) = H 0 ; then by inserting in Equation (121) we obtain the following algebraic equation: Here R 0 = 12H 2 0 and the contribution of matter is neglected. Then, for positive roots H 0 of this equation, the corresponding F (R) leads to the de Sitter solution which may describe inflation or DE. In this case the formula for the entropy Equation (110) can be reproduced for w = −1, and even the universal bound Equation (125) can hold by taking a critical size of the Universe. The formula that relates the cosmic bounds in Equation (128) is easily obtained also in F (R)-gravity for a de Sitter solution. In such a case one can identify: which corresponds to the first FLRW equation written as S 2 H + (S B − S BH ) 2 = S 2 B . Thus, one can conclude that dynamical entropy bounds are not violated for modified gravity with de Sitter solutions. Note that quantum gravity effects may be presented also as an effective fluid contribution. In case when de Sitter space turns out to be the solution, even with the account of quantum gravity the above results indicate that dynamical cosmological/entropy bounds are valid. In other words, the argument indicates the universality of dynamical bounds. It seems that their violation is caused only by future singularities if they are not cured by quantum gravity effects. Note that a large number of modified gravity theories do not contain future singularities; they are cured by higher derivatives terms.

X. ON THE COSMOLOGICAL BOUNDS NEAR FUTURE SINGULARITIES
In [98], Verlinde proposed a new universal bound on cosmology based on a restriction of the Casimir energy E C ; cf. his entropy formula Equation (111). This new bound postulated was: where E BH = n(n − 1) V 8πGa 2 . It was deduced by the fact that in the limit when the Universe passes between strongly and weakly self-gravitating regimes, the Bekenstein entropy S B = 2πa n E and the Bekenstein-Hawking entropy S BH = (n − 1) V 4Ga , which define each regime, are equal. This bound could be interpreted to mean that the Casimir energy never becomes able to reach sufficient energy, E BH , to form a black hole of the size of the Universe. It is easy to verify that the strong (Ha ≥ 1) and weak (Ha ≤ 1) self-gravity regimes have the following restrictions on the total energy: From here it is easy to calculate the bounds on the entropy of the Universe in the case when the Verlinde formula Equation (111) is valid; this is (as shown in the above sections) for an effective radiation dominated Universe w eff ∼ 1/n. The bounds for the entropy deduced in [98] for k = 1 are: where S B is the Bekenstein entropy defined above, and S H is the Hubble entropy given by Equation (112). Note that for the strong self-gravity regime, Ha ≥ 1, the energy range is E C ≤ E BH ≤ E. According to the formula Equation (111) the maximum entropy is reached when the bound is saturated, E C = E BH . Then S = S H , such that the FLRW equation coincides with the CV formula, thus indicating a connection with CFT. For the weak regime, Ha ≤ 1, the range of energies goes as E C ≤ E ≤ E BH and the maximum entropy is reached earlier, when E C = E, yielding the result S = S B . The entropy bounds can be extended to more general cases, corresponding to an arbitrary EoS parameter w. By taking the bound Equation (125) to be universally valid, one can easily deduce the new entropy bounds for each regime, from the expression of the entropy Equation (110). These new bounds, discussed in [99], differ from the ones given in Equation (127), but still establish a bound on the entropy as long as the bound on E C expressed in Equation (125) is taken to be valid. The entropy bounds can be related through the first FLRW equation, yielding the following quadratic expression (for k = 1): Here, we would like to review what happens to the bounds, particularly to the fundamental bound Equation (125), when the cosmic evolution is close to a future singularity (see [100]); then the effective fluid dominating the cosmic evolution could have an unusual EoS. As shown below, for some class of future singularities, such a bound could soften the singularities in order to avoid violation of the universal bound Equation (125). It could be interpreted to mean that quantum effects become important when the bound is reached. However, as the violation of the bound could happen long before the singularity even in the presence of quantum effects, it could be a signal of breaking of the universality of the bound Equation (125). Let us first of all give a list of the possible future cosmic singularities, which can be classified according to [163] as: • Type I ("Big Rip"): For t → t s , a → ∞ and ρ → ∞, |p| → ∞.
• Type IV: For t → t s , a → a s and ρ → ρ s , p → p s but higher derivatives of Hubble parameter diverge (see [164]).
Note that the above list was suggested for the case of a flat FLRW Universe (see [165][166][167][168][169]). As we consider in this section a closed Universe (k = 1), we should make an analysis to see if the list of singularities given above is also valid in this case. It is straightforward to see that all the singularities listed above can be reproduced for a particular choice of the effective EoS. To show how the cosmic bounds behave for each type of singularity, we could write an explicit solution of the FLRW equations, expressed as a function of time depending on free parameters that will be fixed for each kind of singularity. Then, the Hubble parameter may be written as follows: where m is a constant properly chosen for each type of singularity. Note that this is just a solution that ends in the singularities mentioned above, but there are other solutions which also reproduce such singularities. We will study how the cosmic bounds behave near each singularity listed above. As pointed out in [163,, around a singularity quantum effects could become important as the curvature of the Universe grows and diverges in some of the cases. In other words, approaching the finite-time future singularity the curvature grows and universe reminds the early universe where quantum gravity effects are dominant ones because of extreme conditions. Then one has to take into account the role of such quantum gravity effects which should define the behavior of the universe just before the singularity. Moreover, they may act so as to prevent the singularity occurrence. In a sense, one sees the return of quantum gravity era. However, the consistent quantum gravity theory does not exist so far. Then, in order to estimate the influence of quantum effects to universe near to singularity, one can use the effective action formulation. We will apply the effective action produced by conformal anomaly (equivalently, the effective fluid with pressure/energy-density corresponding to conformal anomaly ones) because of several reasons. It is known that at high energy region (large curvature) the conformal invariance is restored so one can neglect the masses. Moreover, one can use large N approximation to justify why large number of quantum fields may be considered as effective quantum gravity. Finally, in the account of quantum effects via conformal anomaly we keep explicitly the graviton (spin 2) contribution. The conformal anomaly T A has the following well-known form: Here we assume for simplicity a 3 + 1 dimensional spacetime. Then, F is the square of a 4D Weyl tensor and G is the Gauss-Bonnet invariant: The coefficients b and b in Equation (130) are described by the number of N scalars, N 1/2 spinors, N 1 vector fields, N 2 gravitons and N HD higher derivative conformal scalars. They can be written as: As b is arbitrary it can be shifted by a finite renormalization of the local counterterm. The conformal anomaly T A can be written as T A = −ρ A + 3p A , where ρ A and p A are the energy and pressure densities respectively. By using Equation (130) and the energy conservation equation ρ A + 3H(ρ A + p A ) = 0, one obtains the following expression for ρ A [163,191]: The quantum corrected FLRW equation is given by: We study now how the bounds behave around the singularity in the classical case when no quantum effects are added, and then include the conformal anomaly Equation (130) quantum effects in the FLRW equations. We will see that for some cases the violation of the cosmic bound can be avoided.

A. Big Rip Singularity
This type of singularity has been very well studied and has become very popular as it is a direct consequence in the majority of the cases when the effective EoS parameter is less than −1, the so-called phantom case .
Observations currently indicate that the phantom barrier could have already been or it will be crossed in the near future, so a lot of attention has been paid to this case. It can be characterized by the solution Equation (129) with m ≤ −1, and this yields the following dependence of the total energy density on the scale factor near the singularity, when a 1, for a closed Universe (k = 1): where we have chosen H 1 = 2/n|1 + w| with w < −1, m = −1 and H 0 = 0 for clarity. This solution drives the Universe to a Big Rip singularity for t → t s , where the scale factor diverges. If the singularity takes place, the bound Equation (125) has to be violated before this happens. This can be seen from Equation (135), as the Casimir energy behaves as E C ∝ a n|w| while the Bekenstein-Hawking energy goes as E BH ∝ a n−2 . Then, as w < −1, the Casimir energy grows faster than the BH energy. Therefore, close to the singularity where the scale factor becomes very big, the value of E C will be much larger than E B , thus violating the bound Equation (125). Following the postulate from [98] one could interpret the bound Equation (125) as the limit where General Relativity and Quantum Field Theory converge, such that when the bound is saturated quantum gravity effects should become important. QG corrections could help to avoid the violation of the bound and may be the Big Rip singularity occurrence. As this is just a postulate based on the CV formula, which is only valid for special cases as shown in the sections above, the bound on E C could not be valid for any kind of fluid.
Let us now include the conformal anomaly Equation (130) as a quantum effect that becomes important around the Big Rip. In such a case there is a phase transition and the Hubble evolution will be given by the solution of the FLRW Equation (134). Let us approximate to get some qualitative results, assuming 3 + 1 dimensions. Around t s the curvature is large, and |ρ A | >> (3/κ 2 )H 2 + k/a 2 . Then ρ ∼ −ρ A , and from Equation (133)we get: We assume that the energy density, which diverges in the classical case, behaves now as: where λ is some negative number. By using the energy conservation equationρ + 3H(1 + w)ρ = 0, the Hubble parameter goes as H ∼ 1/(t s − t). We can check if this assumption is correct in the presence of quantum effects by inserting both results in Equation (136). We get: Hence as b > 0 and b < 0, ρ becomes negative, which is an unphysical result. Thus ρ should not go to infinity in the presence of the quantum correction. This is the same result as obtained in [163] where numerical analysis showed that the singularity is moderated by the conformal anomaly, so that the violation of the bound that naturally occurs in the classical case can be avoided/postponed when quantum effects are included. Now we consider the case where f (R) has the form: when the curvature is small or large. Then if the matter has an EoS parameter w > −1, by solving (121) we find: Then there may appear a singularity at t = 0, which corresponds to the Big Bang singularity, or at t = t s , which corresponds to the Big Rip singularity. Since the Casimir energy behaves as E C ∼ a −nw but E BH ∼ a n−2 , only when t → 0, E C dominates in case that n > 2 and w ≥ 0 or in case that n ≥ 2 and w > 0. Even in the phantom phase where 2α n(1+w) < 0, the bound Equation (125) is not violated.

B. Sudden Singularity
This kind of singularity is also problematic with respect to the bounds. Nevertheless, as long as the energy density ρ does not diverge, the violation of the bound may be avoided for some special choices. The sudden singularity can be described by the solution Equation (129) with 0 < m < 1, and constants H 0,1 > 0. Then the scale factor goes as: which gives a(t) ∼ e h0t (de Sitter) close to t s . From the first FLRW equation the total energy density becomes: which tends to a constant ρ ∼ H 2 0 + e −2H0ts for t → t s . Then the Casimir energy grows as E C ∝ H 2 0 a n + a n−2 , while E BH ∝ a n−2 close to t s . The BH energy grows slower than the Casimir energy, and the bound is violated for a finite t. However, by a specific choice of the coefficients, the violation of the bound Equation (125) could be avoided. For H 0 = 0, and by some specific coefficients, the bound could be obeyed. In general, it is very possible that E C exceeds its bound. In the presence of quantum corrections, the singularity can be avoided but the bound can still be violated, depending on the free parameters for each model. We may assume that in the presence of the conformal anomaly for n = 3, the energy density grows as [163]: where ρ 0 and ρ 1 are constants, and λ is now a positive number. Then the divergences on the higher derivatives of the Hubble parameter can be avoided, as is shown in [163]. Nevertheless, E C still grows faster than E BH , such that the Universe has to be smaller than a critical size in order to hold the bound Equation (125) as is pointed in [99] for the case of a vacuum dominated universe.

C. Type III Singularity
This type of singularity is very similar to the Big Rip, in spite of the scale factor a(t) being finite at the singularity. The solution Equation (129) reproduces this singularity by taking −1 < m < 0. The scale factor goes as: where for simplicity we take H 0 = 0. Then, for t → t s , the scale factor a(t) → a s . To see how E C behaves near the singularity, let us write it in terms of the time instead of the scale factor: where m < 0. Hence, the Casimir energy diverges at the singularity, while E BH ∝ a n−2 s takes a finite value for the singularity time t s , so the bound is clearly violated long before the singularity. Then, in order to maintain the validity of the bound Equation (125), one might assume, as in the Big Rip case, that GR is not valid near or at the bound. Even if quantum effects are included, as was pointed in [163], for this type of singularity the energy density diverges more rapidly than in the classical case, so that the bound is also violated in the presence of quantum effects.

D. Type IV Singularity
For this singularity, the Hubble rate behaves as: Here H 1 (t) and H 2 (t) are regular function and do not vanish at t = t s . The constant α is not integer and larger than 1. Then the scale factor behaves as Near t = t s , the first term dominates and all quantities like ρ, p, and a etc. are finite and therefore the bound Equation (125) would not be violated near the singularity.

E. Big Bang Singularity
When the matter with w ≥ 0 coupled with gravity and dominates, the scale factor behaves as Then there appears a singularity at t = 0, which may be a Big Bang singularity. Although the Big Bang singularity is not a future singularity, we may consider the bound Equation (125) when t ∼ 0. Since n(1 + w) > 2, the energy density behaves as ρ ∼ a −n(1+w) and therefore the Casimir energy behaves as E C ∼ a −nw . On the other hand, we find E BH ∼ a n−2 . Then when n > 2 or when n ≥ 2 and w > 0, E C dominates when a → 0, that is, when t → 0, and the bound Equation (125) is violated. This tells us, as expected, that quantum effects become important in the early universe.
Thus, in the above we have explored what happens near the future cosmic singularities. We have seen that in general, and with some very special exceptions in the case of Type II and Type IV, the bound will be violated if one assumes the validity of GR close to the singularity. Even if quantum corrections are assumed, it seems that the bound will be violated, although in the Big Rip case the singularity may be avoided when quantum effects are incorporated. It is natural to suggest, in accordance with Verlinde, that the bound on the Casimir energy means a finite range for the validity of the classical theory. When this kind of theory becomes saturated, some other new quantum gravity effects have to be taken into account. Hence, the universality of the bound Equation (125) is not clear and may hold just for some specific cases, like the radiation dominated Universe. Furthermore, even in absence of a future singularity, a possibility allowed also for the case that w < −1 (see [192]), the bound may be violated, specially for the case of a "Little Rip" (see [193]), a stage of the universe when the expansion would be strong enough to break some binding systems such as Solar System. That scenario, absence of singular points, can be very possible for a lot of DE models, and also in modified gravity (see [149,194,195])

XI. CONCLUSIONS
In this work we have considered both static spherically symmetric solutions and massive charged and spinning solutions in f (R) theories of gravity in a vacuum situation within the metric formalism. The first case was analyzed in arbitrary dimensions, whereas for the second we focused on four dimensions.
We first discussed the constant curvature case (including charged black holes in four dimensions) in static and spherically symmetric configurations. Then we studied the general case without imposing, a priori, the condition of constant curvature. For the uncharged case, we tackled this calculation through a perturbative analysis around the EH case, which made it possible to study those solutions being regular in the perturbative parameter. We have reproduced the explicit expressions up to second order for the metric coefficients, which only give rise to constant curvature (Schwarzschild AdS) solutions as in the EH case [112,113].
With regards to the massive, charged and spinning BH for f (R) gravity, we have derived the metric tensor that describes such objects. We presented how this solution differs from that found by Carter [119] by a proper redefinition of the spin, electric charge and vacuum scalar curvature and therefore that this type of solutions are supported in this kind of modified gravity theories.
Further study of the Kerr-Newman metric allowed us to describe the different black-hole-type configurations derived from the presence of different horizons: usual, extremal and marginal extremal black holes, naked singularities and naked extremal singularities.
On the other hand, we have also calculated thermodynamical quantities for the AdS as well as for KN configurations. We paid special attention to the issue of the stability of these kinds of solutions. As a very remarkable result, we have presented that the condition for a f (R) theory of gravity to support these kinds of black holes is given by R 0 + f (R 0 ) < 0 where R 0 is the constant curvature of the AdS space-time. In fact, this condition also implies that the effective Newton constant is positive and consequently the graviton does not become a ghost. By employing the Euclidean action method, we revised how the mass, the energy or the entropy of these BH (AdS and KN) differ from those predicted in GR just by a multiplicative factor 1 + f (R 0 ). This factor has to be positive in order to guarantee a positive mass and entropy for these kinds of BH.
For the AdS configuration, we sketched the fact that the qualitative thermodynamic behavior of the AdS BH in f (R) gravity theories is the same as the one found by Hawking and Page. On the other hand, with regards to the KN configurations, the analysis of the BH heat capacity revealed how two different types of BH can be distinguished: fast and slow, showing the latter two phase transitions. With respect to the horizon structure in the KN configuration, we have shown how horizons can only exist for values of the spin lower than a maximum value a max , and that from a certain positive value of the curvature onward, only above a minimum value a min .
Then, for two paradigmatic f (R) models, we have investigated the stability of the different possible configurations that arise from the values of the free energy and the heat capacity, quantities that depend on the particular f (R) model under study. We have analyzed explicitly the rich thermodynamical phenomenology that characterizes these gravitational models with a complete set of different figures that were originally presented in [112,113] and [120,121]. For doing so, we have studied the parameter regions in which BH in such models are locally stable and globally preferred.
Furthermore, we have investigated the deep connection between FLRW equations and the first law of thermodynamics. The study of FLRW metrics is a very important issue as the universe seems to be well described by this kind of metrics, where eventually one should assume extra components to explain the cosmological evolution. However, as suggested above, modified gravity can even avoid the need to introduce dark energy, and even a kind of inflaton, to explain the whole cosmological evolution. Hence, FLRW metrics and their connection to thermodynamics turns out to be a fundamental issue. Then, as suggested in [82][83][84], FLRW equations can be derived by assuming the entropy relation with the area of the apparent horizon. Even more, such a relation can be extended to f (R) gravity when one assumes the modified relation for the entropy of the horizon.
In addition, we have also reviewed the possible extension of generalized CV formula for multicomponent fluids, generalized in the sense that an inhomogeneous EoS (including modified f (R) gravity) was assumed. For some special cases the formula is reduced to the standard CV formula expressing the correspondence with 2d CFT theory. The dynamical entropy bound for all above cases was found. The universality of dynamical entropy bound near all four types of the future singularity, as well as the initial Big Bang singularity, was investigated. Note that the study of future singularities is a crucial step since a lot of dark energy models contain some kind of singularity. In fact, so-called viable modified gravities usually cross the phantom barrier (w < −1) (see [196]), and although it seems that future singularities do not occur, the "Little Rip" event may take place (see [151]). It was proved that except for some special cases of Type II and Type IV singularity, the dynamical entropy bound is violated near the singularity. Taking into account quantum effects of conformally invariant matter does not improve the situation. Hence, the dynamical entropy bound may not be universal and its violation simply indicates that the situation will be changed with the introduction of quantum gravity effects. This deep connection between gravitation and thermodynamics in so many scenarios, as the ones reviewed here in extended gravitational theories and others (see [197]), may suggest that gravity could be just an emergent force as suggested in [198], which may imply the birth of new physics that may lead to the search of a more fundamental theory.
Experimental checks to test the validity of a particular f (R) model can be applied in multiple realms of astrophysics. Among others, future data from CMB tensor perturbations [199] as well as gravitational collapse predictions [200][201][202] and the growth of cosmological perturbations (and subsequent structure formation) [47][48][49][50] may constrain the available range of parameters for viable models [203,204]. In addition, the study of the weak lensing in f (R) gravity can also provide complementary constraints on the models and free parameters as pointed out in [205]. The general viability conditions for f (R) theories prior to model specifications are summarized in reference [114]. Furthermore, the study of star structure and evolution in modified gravity theories may lead additional and independent tests for this class of fourth order theories of gravity (see [206]). Finally, the aforementioned constraints may be enlarged not only by studying astrophysical BH stability but also if quantum gravity scale is accessible at particle accelerators. For example, it has been shown that LHC could produce about one microBH per second [207,208]. Stability and thermodynamical properties of the produced BH might shed some light about the underlying theory of gravity.
DSG acknowledges support from a postdoctoral contract from the University of the Basque Country (UPV/EHU). DSG is also supported by the research project FIS2010-15492, and also by the Basque Government through the special research action KATEA and UPV/EHU under program UFI 11/55.