ar X iv :2 00 7. 04 06 5v 2 [ gr -q c] 2 1 Se p 20 20 Junction conditions in Palatini f(R) gravity Gonzalo J. Olmo1, 2, ∗ and Diego Rubiera-Garcia3, † 1Departamento de F́ısica Teórica and IFIC, Centro Mixto Universidad de Valencia - CSIC. Universidad de Valencia, Burjassot-46100, Valencia, Spain 2Departamento de F́ısica, Universidade Federal da Paráıba, 58051-900 João Pessoa, Paráıba, Brazil 3Departamento de F́ısica Teórica and IPARCOS, Universidad Complutense de Madrid, E-28040 Madrid, Spain (Dated: September 22, 2020) We work out the junction conditions for f(R) gravity formulated in metric-affine (Palatini) spaces using a tensor distributional approach. These conditions are needed for building consistent models of gravitating bodies with an interior and exterior regions matched at some hypersurface. Some of these conditions depart from the standard Darmois-Israel ones of General Relativity and from their metric f(R) counterparts. In particular, we find that the trace of the stress-energy momentum tensor in the bulk must be continuous across the matching hypersurface, though its normal derivative need not to. We illustrate the relevance of these conditions by considering the properties of stellar surfaces in polytropic models, showing that the range of equations of state with potentially pathological effects is shifted beyond the domain of physical interest. This confirms, in particular, that neutron stars and white dwarfs can be safely modelled within the Palatini f(R) framework. I. INTRODUCTION Many models of gravitating bodies of both theoreti- cal and observational interest involve the consideration of two different patches of space-time, glued at some hy- persurface separating the interior from the exterior (see e.g. [1] and references therein). Such models can suc- cessfully yield both configurations with an interior filled with some matter distribution separated from an exte- rior vacuum solution, as in stellar bodies [2, 3], or two matter-filled regions such as in branes/domain walls [4–7] or thin-shells [8–10]. For these constructions to be well defined from a mathematical point of view, a number of conditions must hold at the matching hypersurface, which within the context of General Relativity (GR) are known as the Darmois-Israel matching (or junction) con- ditions [11, 12]. These conditions involve the continuity of the first fundamental form across the matching hyper- surface, and a number of conditions relating the stress- energy tensor with the discontinuity in the extrinsic cur- vature on that hypersurface. Given the physical relevance of such models and condi- tions for the modeling of different bodies, a timely ques- tion is to what extent they can be applied to other the- ories of gravity beyond GR (see e.g. [13–16] for some re- views). The consideration of these theories are motivated by a variety of reasons in view of the challenges posed by the exploration of the strong-field regime of the gravita- tional interaction in the multimessenger astronomy era [17], and by the difficulties of the early and late-time cosmological models [18]. Typically such theories intro- duce additional difficulties in the form of higher-order equations of motion and much more involved dynamics, ∗ gonzalo.olmo@uv.es † drubiera@ucm.es which naturally requires an upgrade of the junction con- ditions [19–24]. This introduces additional constraints, which troubles the construction of such matched mod- els. Perhaps the simplest family of extensions of GR are those based on a function of the curvature scalar, the so- called f(R) gravity, which has attracted a great deal of attention in the literature [25]. In this case, it has been proven that many known matched solutions in GR do not survive when considered in the f(R) scenario [26]. Thus, the analysis of the junction conditions beyond GR is of utmost importance in order to consistently explore the phenomenology of such models. The main aim of this paper is to work out the junc- tion conditions corresponding to the Palatini formula- tion of f(R) theories. In the Palatini approach, metric and affine connection are regarded as independent de- grees of freedom, and the field equations are obtained from independent variations of the action with respect to both of them [28]. As opposed to their metric coun- terpart, the resulting theory has second-order field equa- tions, with the new gravitational dynamics being due to new nonlinear couplings engendered by the matter fields, and no extra propagating degrees of freedom are intro- duced. Moreover, Palatini f(R) models may pass the post-Newtonian (solar system) tests [29], are compatible with current gravitational wave observations, and have a rich phenomenology in black hole [30], stellar [31], and cosmological [32] scenarios. The equations of motion and junction conditions corre- sponding to Palatini f(R) theories turn out to be clearly different from those corresponding to GR and to the met- ric formulation of f(R) theories. Differences arise even in the conservation equations on the matching hypersurface, being the same in GR and in metric f(R) [26] but differ- ent in Palatini f(R) gravity. Our results will also prove useful, in particular, to better understand an issue that has undermined progress in the study of stellar models in Palatini f(R) theories in the last years. We are refer- http://arxiv.org/abs/2007.04065v2 mailto:gonzalo.olmo@uv.es mailto:drubiera@ucm.es 2 ring to the potential emergence of curvature divergences near the surface of polytropic stellar models [33–36]. As we will see, a consistent distributional treatment of the field equations puts forward that such divergences are an artifact of the approach followed in [33–36] and that the equations on the matching surface are free from such pathologies for polytropic equations of state of physical interest. II. JUNCTION CONDITIONS FOR f(R) GRAVITY IN PALATINI APPROACH A. Action and field equations The action of f(R) gravity can be written as S = 1 2κ2 ∫ d4x √−gf(R) + ∫ d4x √−gLm(gµν , ψm) , (1) where κ2 is Newton’s constant in suitable units, g is the determinant of the space-time metric gµν , f(R) is some function of the curvature scalar R ≡ gµνRµν(Γ), with Γ ≡ Γλ µν the affine connection, and Lm denotes the mat- ter Lagrangian with ψm representing a collection of mat- ter fields. In the metric formalism the connection is taken to be given by the Levi-Civita one, that is, the one satisfying ∇Γ µ( √−ggαβ) = 0, and the variation of the action (1) with respect to the metric yields the field equations fRGµν(g)− RfR − f 2 gµν − fRR(∇µ∇νR− gµν∇α∇αR) − fRRR(∇µR∇νR− gµν∇α∇αR) = κ2Tµν ,(2) where fR ≡ df/dR, and Tµν = 2√ −g δSm δgµν is the stress- energy tensor. Tracing in this equation with gµν yields 3✷fR +RfR − 2f = κ2T , (3) where T is the trace of the stress-energy tensor. As is well known, the presence of the derivative operator ✷ yields a dynamical character to the (scalar) field φ ≡ fR, which can be read off as an additional propagating degree of freedom. On the other hand, in the Palatini formulation, the field equations follow by independent variations of the action (1) with respect to the metric and the connection. After solving the connection equation, the resulting equa- tions for the metric can be conveniently cast as [27, 28] fRGµν(g) + RfR − f 2 gµν − [∇µ∇νfR − gµν✷fR] (4) + 3 2fR [∇µfR∇νfR − 1 2 gµν(∂fR) 2] = κ2Tµν , with a similar notation as in the metric case. Tracing (4) with gµν yields RfR − 2f = κ2T . (5) This equation implies that R ≡ R(T ) which, instead of a differential equation, is simply an algebraic relation link- ing curvature and matter. This fact introduces a radical difference in the properties and dynamics of the theory with respect to the metric formulation and, therefore, a distributional analysis of the Palatini f(R) equations must take care of that fact. Let us point out that the above equations (2) and (4) for the metric and the Palatini versions of f(R) theories can be seen as particular cases of Brans-Dicke scalar- tensor theories with a potential. In particular, the trace equations (3) and (5) can be seen as specific cases of the general equation (3 + 2ω)✷φ+ 2V (φ) − φVφ = κ2T , (6) with the metric case corresponding to ω = 0 while Pala- tini having ω = −3/2. The scalar field is identified as φ ≡ fR and its potential is given by V (φ) = φR(φ)−f(R(φ)). Note that the scalar field is non-dynamical in the Palatini case, while in any other case it does propagate. In the fol- lowing, we will consider the distributional description of Eq.(4), which is equivalent to the ω = −3/2 Brans-Dicke theory. B. Distributional stress-energy tensor We consider two smooth four-dimensional manifolds (M±, gµν), and let us denote by V ± two bounded regions living in M± with boundaries Σ±. These regions are matched (“glued”) at a time-like hypersurface Σ with the natural identification of their boundaries as Σ+ = Σ−, and with the space-time metric gµν assumed to be well defined through the entire manifold and, in particular, to be continuous (but not necessarily differentiable) at Σ, which is a basic requirement in the junction condi- tions formalism [37]. Since there may be discontinuities in several geometric quantities across Σ, the suitable tool to deal with this scenario is that of tensorial distribu- tions, namely, tensor fields with compact support on the manifold. One of the main goals of the junction condi- tions is to identify the singular parts of the curvature and matter tensorial distributions on Σ and find the allowed discontinuities of these quantities across it. For the sake of our analysis we borrow, to a large extent, the notation introduced in Ref.[38], and we refer the reader there for further details on the mathematical setup employed here. Given the particular way the matter fields feed the new dynamics in Palatini f(R) gravity, as follows from Eq.(5), we start our analysis by introducing a distribution for the stress-energy tensor as Tµν = T+ µνθ + T− µν(1 − θ) + τµνδ Σ . (7) Hereafter underbars indicate distributions; T± µν are the stress-energy tensors in V ±, respectively; θ is the scalar distribution defined by the (locally integrable) Heaviside function, the latter taking the value of 1 in V +, 0 in 3 V −, and any intermediate reference value in Σ; δΣ is a scalar Delta-type distribution with support on Σ acting upon any test function X as < δ,X >≡ ∫ ΣX ; and τµν accounts for the singular part of the stress-energy tensor on Σ. Similarly, the distributional form of the trace of the stress-energy tensor reads T = T+θ + T−(1 − θ) + τδΣ , (8) where τ ≡ τµµ is the trace of the singular part of the stress-energy tensor. To cast the field equations (4) in distributional form, we have to deal with several potential difficulties in the transition from regular tensorial functions to tensorial distributions. Let us isolate first the term ∇µfR∇νfR = f2 RRR 2 T∇µT∇νT , (9) where RT ≡ dR/dT = κ2/(RfRR − fR). From (8), the distributional form of the covariant derivative of the trace of the stress-energy tensor reads ∇µT = ∇µT +θ +∇µT −(1− θ) + nµ[T ]δ Σ +∇µ(τδ Σ) . (10) Here nµ is the unit vector normal to Σ defined via ∇µθ = nµδ Σ, and brackets represent a discontinuity (“jumps”) across Σ in the quantity contained there. Thus [T ] ≡ T+ − T− is the jump in the trace of the stress- energy tensor. Given that the term ∇µT∇νT in (9) leads to products of the form θ · δ and δ · δ, which are not well defined in distributional theory, in order to have consis- tent equations in distributional sense we must require the conditions [T ] = 0 (11) τ = 0 . (12) The first of these conditions simply indicates the need for the continuity of the trace of the stress-energy ten- sor across Σ, which was already expected on grounds of continuity and standard differentiability of the tensorial equations. The second condition, τ = 0, is new and en- tirely due to the distributional nature of the quantities defined on the hypersurface. In the context of GR, it would imply a vanishing brane tension. In the modified gravity context, its implications are still to be under- stood. Note that these restrictions involve the traces of the stress-energy tensor and its singular part on Σ, but tell us nothing about the full structure of these objects. Bearing in mind (11), let us consider next the contri- bution to the field equations (4) given by the term (in distributional sense) ∇µ∇νT = ∇µ∇νT +θ+∇µ∇νT −(1−θ)+nµnνbδ Σ , (13) where we have defined the scalar quantity b ≡ nα[∇αT ] , (14) which is the discontinuity in the covariant derivative of the trace in the normal direction to the hypersurface Σ. This expression shall be needed later. C. Decomposition of the Einstein tensor on Σ The distributional form of the Einstein tensor takes the form [37] Gµν = G+ µνθ +G− µν(1 − θ) + Gµνδ Σ , (15) where Gµν represents the singular part of the Einstein tensor on Σ. Plugging this expression in the field equa- tions (4), and taking into consideration the constraint (11) on the stress-energy tensor, the singular part of the field equations (4) can be written as fRΣ Gµν + fRRRT |Σ bhµν = κ2τµν , (16) where the subindex Σ denotes quantities evaluated on Σ, while hµν is the projector on Σ, namely hµν ≡ gµν − nµnν , (17) usually referred to as the first fundamental form, and which must be continuous across Σ [37]. Given that Gµν can be expressed as [38] Gµν = −[Kµν ] + hµν [K ρ ρ ] , (18) where K± µν ≡ hρβh σ µ∇± ρ nσ , (19) is the second fundamental form on V ±, respectively, we can trace in Eq.(16) with the metric gµν and using (12) together with the fact that (18) implies that Gρ ρ = 2[Kρ ρ ], one arrives to [Kρ ρ ] = − 3fRRRT 2fR ∣ ∣ ∣ ∣ Σ b . (20) Plugging this relation back in (16) one finally gets − [Kµν ] + 1 3 hµν [K ρ ρ ] = κ2 τµν fRΣ , (21) which is clearly consistent with the tracelessness of τµν in (12). Note that these results largely depart from the corresponding expressions in GR, where the singu- lar part of the Einstein equations Gµν = κ2τµν does not introduce the constraint (12), since the latter is associ- ated to the contribution (9), which is vanishing when f(R) = R− 2Λ (the Einstein-Hilbert action with cosmo- logical constant term Λ). Thus, in such a case, instead of (21) one has −[Kµν ] + hµν [K ρ ρ ] = κ2τµν , whose trace reads 2[Kρ ρ ] = κ2τ and, therefore, the brane tension in GR is non-vanishing in general. As opposed to that, we have just seen that in Palatini f(R) gravity the brane tension vanishes but [Kρ ρ ] 6= 0. In the metric version of f(R), however, one typically has [Kρ ρ ] = 0 [26]. It is important to note that the condition (12) imposes a con- straint between the energy density and the two principal pressures on Σ, which effectively reduces the number of degrees of freedom of the matter on the hypersurface Σ. This is expected to have important implications for the resolution of the conservation equations. 4 D. Energy conservation equations As shown in [38], the Bianchi identities hold in the distributional sense, that is ∇µG µ ν = 0. This can be explicitly written as two sets of equations: (K+ ρσ +K− ρσ)Gρσ = 2nρnσ[Rρσ ]− [R] (22) DρGρν = −nρhσν [Rρσ ] , (23) where Dρ ≡ hρ α∇α denotes the covariant derivative on Σ. On the other hand, given that (5) leads to R = R(T ) and that [T ] = 0, it follows that [R] = 0. Our next goal is to extract additional information from the two equations (22) and (23) about the behavior of the stress-energy tensor on the hypersurface Σ. To this end, let us start by taking a covariant derivative upon Eq.(16) to find κ2Dρτ ρ ν = fRRRT |Σ ((DρTΣ)Gρ ν +Dνb) + fRΣ DρGρ ν , + (fRRRR 2 T + fRRRTT ) ∣ ∣ Σ b(DνTΣ) , (24) where we have used the relation DνfR = fRRRT |Σ (DνTΣ) with TΣ denoting the value of T on Σ. Let us deal with the different terms in this expression, starting with the one in [Rµν ], as follows from the replacement of (23) in (24). Projecting on the field equations (4) with nρhσν we find fRn ρhσνRρσ − nρhσν∇ρ∇σfR + 3 2fR (nρ∇ρfR)fRRRT (DνTσ) = κ2nρhσνTρσ .(25) The jump of this expression across Σ becomes fRn ρhσν [Rρσ ] = κ2nρhσν [Tρσ] + nρhσν [∇ρ∇σfR] − 3f2 RR 2fR R2 T b(DνTΣ) . (26) Bearing in mind that ∇ρ∇σfR = ∇ρ(fRRRT∇σT ) (27) = (fRRRR 2 T + fRRRTT )∇ρT∇σT + fRRRT∇ρ∇σT , then its jump across Σ can be written as nρhσν [∇ρ∇σfR] = (fRRRR 2 T + fRRRTT )b(DνTΣ) + fRRRTn ρhσν [∇ρ∇σT ] . (28) To deal with the last term in this expression we introduce, following [26], a decomposition of the form [∇ρ∇σT ] = Anρnσ + nρ((Dσb)− [Kρσ](D ρTΣ)) (29) + nσ((Dρb)− [Kρν ](D ρTΣ)) + b 2 (K+ ρσ +K− ρσ) , where A ≡ nµnν [∇µ∇νR]. This way, Eq.(28) reads nρhσν [∇ρ∇σT ] = (Dνb)− [Kαν ](D αTΣ) . (30) We are almost there. Combining Eqs.(26), (28) and (30) with (24), and using (23) and (20) we find, after a bit of algebra, that (24) simplifies down to Dρτρν = −nρhσν [Tρσ] , (31) which coincides with the expression found in GR and in metric f(R) gravity. Let us now extract information from Eq.(22). The most important point here is to compute the quantity on its left-hand side using the field equations (4). Pro- jecting with fRn ρnσ and comparing the result on both sides of Σ recalling that gµν must be continuous across Σ, one finds fRn ρnσ[Gρν ] + hρσ[∇ρ∇σfR]− 3 2fR hρσ[∇ρfR∇σfR] + 3 4fR [(∇µfR) 2] = κ2nρnσ[Tρσ] . (32) Let us isolate the various contributions to this equation. First we have the term hρσ[∇ρ∇σfR] = (fRRRR 2 T + fRRRTT )[(DρTΣ) 2] + fRRRTh ρσ[∇ρ∇σT ] = b 2 fRRRTh ρσ(K+ ρσ +K− ρσ) , (33) where in the last line we have explicitly implemented the fact that, since T must be continuous across Σ, then (DρTΣ) = 0. The second piece of (32) that we need to consider is hρσ[∇ρfR∇σfR] = f2 RRR 2 Th ρσ[∇ρT∇σT ] = f2 RRR 2 T [(DρTΣ) 2] = 0 , (34) the last equality due to the same considerations as before. The third relevant piece of (32) is [(∇µfR) 2] = (nµnν + hµν)∇µfR∇νfR = R2 T f 2 RR[b 2] , (35) where we have defined [b2] ≡ [(nα∇αT ) 2] = 2bnα∂αTΣ (see Appendix D.2. of Ref.[39] for details). With all these elements, Eq.(32) reads fRn ρnσ[Gρν ] = κ2nρnσ[Tρσ]− fRRRTh ρσ b 2 (K+ ρσ +K− ρσ) − 3R2 T f 2 RR fR [b2] . (36) Plugging this expression in the right-hand-side of Eq.(32) we finally find (K+ ρσ +K− ρσ)τ ρσ = 2nρnσ[Tρσ]− 3R2 T f 2 RR fR [b2] . (37) Likewise the analysis of the second fundamental form above, this result recovers the one of GR only if fRR = 0, that is, f(R) = R− 2Λ. To summarize the main results of this section, the junc- tion conditions in Palatini f(R) gravity are: the conti- nuity of the first fundamental form (17) across Σ; the continuity of the trace of the stress-energy tensor across Σ, namely Eq.(11); the vanishing of the trace of its sin- gular part on Σ, namely Eq.(12); and Eq.(21) for the dis- continuity of the trace of the second fundamental form, which we recall it departs from the GR and the metric f(R) cases [26]. In cases allowing for shells or branes, the energy content on Σ is determined by Eqs.(21), (31) and (37). 5 III. APPLICATION: POLYTROPIC STELLAR MODELS For a specific implementation of this framework, we will now apply the formulas derived in the previous sec- tion to address an issue raised in [33–35] that affects mod- els of stellar structure in Palatini f(R) gravity. Using the standard tensorial approach, it turns out that when one attempts to match a polytropic stellar model with equa- tion of state ρ(P ) = ( P K )1/γ + P γ − 1 , (38) to an external Schwarzschild solution, certain difficulties arise. Curvature divergences, in particular, develop at the matching surface, the latter characterized by vanish- ing pressure. Knowing that this happens for a range of polytropic indices of physical interest, namely, 3/2 < γ < 2, which includes the relevant case γ = 5/3 that describes a gas of non-relativistic degenerate fermions, one might be tempted to conclude that the viability of Palatini f(R) models could be in danger. Yet one should bear in mind that polytropic models are just crude statistical approxi- mations of some stars. On the other hand, given that the divergent terms are directly related to what happens at the matching hypersurface, it seems appropriate to use the more rigorous approach of tensorial distributions de- veloped above to reassess the problem. Under this light, it seems relevant to see what Eq.(16), which describes the dynamics at the matching hypersurface, has to say about the behavior of curvature there. A glance at Eq.(16) puts forward that the discontinuity in the derivative of the trace of the stress-energy tensor in the direction normal to Σ is the key element that controls the modified dynamics on Σ. In fact, given that the trace T must be continuous across the surface (as follows from (11)) and that R = R(T ) must also be (due to (5)) the coefficient fRΣ must be smooth and generically expected to be close to unity in regions of low energy density such as the surface of a star. This is the case, for instance, of the quadratic model f(R) = R+aR2, for which Eq.(5) yields R(T ) = −κ2T and we have fR = 1 as T → 0. Thus, whenever b ≡ nα[∇αT ] is not zero, the nonlinearities of the f(R) theory will contribute to the surface dynamics. Continuing with the quadratic model, it is easy to see that the term fRRRT = −2aκ2 is just a constant and, therefore, b is the element that may bring in new effects. Considering a polytrope with equation of state (38) it follows that b = nr +∂rT + − nr −∂rT − = −nr −∂rT − ≈ (3 − ρP )Pr. According to the analysis of [36], one finds that ρP ∼ P (1−γ)/γ and Pr ∼ P 1/γ , which leads to b ∼ P (2−γ)/γ . It turns out that this quantity only diverges if γ > 2, shifting the problematic range of γ beyond the domain of direct physical interest. For γ < 2, it is evident that b = 0 at the surface. This makes Eq.(16) formally identical to its GR counterpart, as one would expect in any model which only departs from GR at high curvature/densities. The conservation equations (31) and (37) also degenerate into their GR counterparts. We have just seen that a careful analysis of the junc- tion conditions using distribution theory has been able to get rid of a problem that arises when a standard ten- sorial approach is considered. In the usual approach, in order to have a well defined geometry, the trace of the matter stress-energy tensor must be continuous and dif- ferentiable because, otherwise, the second derivatives of the matter fields that appear in the field equations would lead to divergences. In the distributional approach, how- ever, we have seen that those requirements are slightly re- laxed at a matching hypersurface because one is forced to have a continuous T but a normal discontinuity in ∇µT is allowed and still provides sensible results. A glance at the analysis of [36] shows that the divergent terms that arise at the surface of polytropic models are entirely due to the evaluation of the term ∇µ∇νT at the match- ing surface, where P = 0. Specifically, ∇µ∇νT gener- ates two1 terms of the form ρPPrr and ρPPP 2 r which go as ∼ P (3−2γ)/γ and diverge for γ > 3/2 when P → 0. In the distributional approach, however, the terms with ∇µ∇νT provide a contribution on Σ proportional to b, which measures the discontinuity in ∇µT . This prop- erty is precisely what cures the singular terms on the matching surface, since the seemingly pathological con- tributions ρPPrr and ρPPP 2 r are now replaced by ρPPr. Thus, a smooth matching surface between an external Schwarzschild solution and an internal polytropic fluid in Palatini f(R) gravity is possible whenever γ < 2 when a proper distributional treatment of the matching surface is considered. IV. CONCLUSION In this work we have derived the junction conditions corresponding to Palatini f(R) theories using distribu- tion theory. We have found that the equations that re- late the extrinsic curvature to the matter sources on the matching surface differ from those found in GR and in the metric version of f(R) theories. This is a non-trivial result which has many repercussions for the building of matched solutions of any kind within these theories, and which must be carefully implemented in order not to get to misleading results. As an specific example, we have considered the rele- vant case of polytropic stellar models. In such a case, a naive approach to match the internal fluid solution to the external Schwarzschild solution, as followed in [33– 36], leads to the unacceptable emergence of curvature divergences at the matching surface for a range of poly- tropic indices of physical interest. Here we have shown that a more rigorous treatment using distribution theory 1 We note that the divergence due to ρPPrr was overlooked in [36]. 6 modifies this conclusion, shifting the range of potentially pathological equations of state beyond the domain of rel- evance for physical applications. The junction conditions derived here, therefore, should be applied whenever one attempts to build stellar models that match vacuum ex- ternal solutions in any Palatini f(R) theory allowing, in particular, to consistently model neutron stars and white dwarfs, among others. Further applications of these junctions conditions be- yond stellar bodies include, for instance, the rigorous construction of thin-shell wormholes and the consistency of braneworld models in Palatini f(R) gravity. More- over, the approach followed here using distribution the- ory could be enlarged to include other Palatini theories of gravity beyond the f(R) case, involving further con- tractions of the Ricci/Riemann tensors such as, for in- stance, in the class of the so-called Ricci-based gravities [40], which are ghost-free, consistent extensions of GR. In these more general theories additional contributions of the stress-energy tensor (and not just its trace, as in the present case) are expected to yield novelties in the shape of the junction conditions as compared to GR, to Palatini f(R) theories, and to their metric counter- parts. This would allow, for instance, to reassess the claim of Ref.[41] on the existence of surface singulari- ties in Eddington-inspired Born-Infeld gravity (see [42] for a description of these theories, and [43, 44] for spe- cific applications) in pretty much the same way as for the Palatini f(R) case studied here. Work along these lines is currently underway. ACKNOWLEDGMENTS GJO is funded by the Ramon y Cajal contract RYC- 2013-13019 (Spain). DRG is funded by the Atracción de Talento Investigador programme of the Comunidad de Madrid (Spain) No. 2018-T1/TIC-10431, and acknowl- edges further support from the Ministerio de Ciencia, In- novación y Universidades (Spain) project No. PID2019- 108485GB-I00/AEI/10.13039/501100011033, and the Fundação para a Ciência e a Tecnologia (FCT, Portugal) research projects Nos. PTDC/FIS-OUT/29048/2017 and PTDC/FIS-PAR/31938/2017. This work is sup- ported by the Spanish project FIS2017-84440-C2- 1-P (MINECO/FEDER, EU), the project H2020- MSCA-RISE-2017 Grant FunFiCO-777740, the project SEJI/2017/042 (Generalitat Valenciana), the Con- solider Program CPANPHY-1205388, the Severo Ochoa grant SEV-2014-0398 (Spain) and the Edital 006/2018 PRONEX (FAPESQ-PB/CNPQ, Brazil) Grant No. 0015/2019. This article is based upon work from COST Actions CA15117 and CA18108, supported by COST (European Cooperation in Science and Technology). [1] H. Stephani et al., “Exact solutions of Einstein’s field equations” (CUP, 2003). [2] N. K. Glendenning, “Compact stars”, 2nd Ed. (Springer, 2001). [3] G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, arXiv:1912.05202 [gr-qc]. [4] L. Randall and R. Sundrum, Phys. Rev. Lett. 83 (1999) 3370. [5] L. Randall and R. Sundrum, Phys. Rev. Lett. 83 (1999) 4690. [6] G. R. Dvali, G. Gabadadze and M. Porrati, Phys. Lett. B 485 (2000) 208. [7] A. Vilenkin, Phys. Rept. 121 (1985) 263. [8] N. M. Garcia, F. S. N. Lobo and M. Visser, Phys. Rev. D 86 (2012) 044026. [9] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela and P. Pani, Phys. Rev. D 94 (2016) 084031. [10] F. S. Lobo, A. Simpson and M. Visser, Phys. Rev. D 101 (2020) 124035. [11] W. Israel, Nuovo Cim. B 44 (1966) 1. [12] G. Darmois, Mémorial des Sciences Mathématiques, Fas- cicule 25 (GauthierVillars, Paris, 1927). [13] S. Capozziello and M. De Laurentis, Phys. Rept. 509 (2011) 167. [14] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, Phys. Rept. 513 (2012) 1. [15] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692 (2017) 1. [16] L. Heisenberg, Phys. Rept. 796 (2019) 1. [17] E. Berti et al., Class. Quant. Grav. 32 (2015) 243001. [18] P. Bull et al., Phys. Dark Univ. 12 (2016) 56. [19] S. C. Davis, Phys. Rev. D 67 (2003) 024030. [20] E. Gravanis and S. Willison, Phys. Lett. B 562 (2003) 118. [21] N. Deruelle, M. Sasaki and Y. Sendouda, Prog. Theor. Phys. 119 (2008) 237. [22] A. Padilla and V. Sivanesan, JHEP 1208 (2012) 122. [23] S. Vignolo, R. Cianci and S. Carloni, Class. Quant. Grav. 35 (2018) 095014. [24] A. de la Cruz-Dombriz, P. K. S. Dunsby and D. Saez- Gomez, JCAP 1412 (2014) 048. [25] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13 (2010) 3. [26] J. M. M. Senovilla, Phys. Rev. D 88 (2013) 064015. [27] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82 (2010) 451. [28] G. J. Olmo, Int. J. Mod. Phys. D 20 (2011) 413. [29] G. J. Olmo, Phys. Rev. D 72 (2005) 083505. [30] C. Bambi, A. Cardenas-Avendano, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 93 (2016) 064016. [31] G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, Phys. Rev. D 100 (2019) 044020. [32] V. M. Enckell, K. Enqvist, S. Rasanen and L. P. Wahlman, JCAP 1902 (2019) 022. [33] E. Barausse, T. P. Sotiriou and J. C. Miller, Class. Quant. Grav. 25 (2008) 062001. [34] E. Barausse, T. P. Sotiriou and J. C. Miller, Class. Quant. Grav. 25 (2008) 105008. http://arxiv.org/abs/1912.05202 7 [35] E. Barausse, T. P. Sotiriou and J. C. Miller, EAS Publ. Ser. 30 (2008) 189. [36] G. J. Olmo, Phys. Rev. D 78 (2008) 104026. [37] C. J. S. Clarke and T. Dray, Class. Quant. Grav. 4 (1987) 265. [38] M. Mars and J. M. M. Senovilla, Class. Quant. Grav. 10 (1993) 1865. [39] B. Reina, J. M. M. Senovilla and R. Vera, Class. Quant. Grav. 33 (2016) 105008. [40] V. I. Afonso, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 97 (2018) 021503. [41] P. Pani and T. P. Sotiriou, Phys. Rev. Lett. 109 (2012) 251102. [42] J. Beltran Jimenez, L. Heisenberg, G. J. Olmo and D. Rubiera-Garcia, Phys. Rept. 727 (2018) 1. [43] Y. Tavakoli, C. Escamilla-Rivera and J. C. Fabris, An- nalen Phys. 529 (2017) 1600415. [44] M. Banados and P. G. Ferreira, Phys. Rev. Lett. 105 (2010) 011101.