In this work we study static spherically symmetric solutions of effective field equations related to local and nonlocal higherderivative gravity models based on their associated effective delta sources. This procedure has been applied to generate modifications of the Schwarzschild geometry in several contexts e.g. modified gravity string theory noncommutative geometry generalized uncertainty principle scenarios but a general analysis of the possible equations of state and their influence on the solutions was still lacking. Here we aim to fill this gap in the literature and investigate whether these metrics might be able to reproduce features of the solutions of higherderivative gravity models. In particular we present an equation of state such that the solution matches the Newtonianlimit one in both regimes of large and small r. A significant part of the work is dedicated to studying the curvature regularity of the solutions and the comparison with the linearized solutions. Explicit metrics are presented for effective sources originating from local and nonlocal models. The results obtained here might be regarded as possible links between the previous research on linearized higherderivative gravity and the solutions of the nonlinear complete field equations which remain unknown at the moment.
['Regular black hole', 'Higher-derivative gravity', 'Nonlocal gravity', 'Effective source formalism', 'Models of quantum gravity', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Theory']
On effective models of regular black holes inspired by higherderivative and nonlocal gravity
{'On e ff ective models of regular black holes inspired by higher-derivative and nonlocal gravity': "Tib'erio de Paula Netto , 1, ∗ Breno L. Giacchini , 2, 1, † Nicol'o Burzill'a , 3, 4, ‡ and Leonardo Modesto 5, 4, §  \n1 Departamento de F'ısica, ICE, Universidade Federal de Juiz de Fora, Juiz de Fora, 36036-900, MG, Brazil  \n2 Institute of Theoretical Physics, Faculty of Mathematics and Physics,  \nCharles University, V Holeˇsoviˇck'ach 2, 180 00 Prague 8, Czech Republic  \n3 INFN Sezione di Roma Tor Vergata, 00133 Roma, Italy  \n4 Department of Physics, Southern University of Science and Technology, Shenzhen 518055, China  \n5 Dipartimento di Fisica, Universit'a di Cagliari, Cittadella Universitaria, 09042 Monserrato, Italy  \nIn this work we study static spherically symmetric solutions of e ff ective field equations related to local and nonlocal higher-derivative gravity models, based on their associated e ff ective delta sources. This procedure has been applied to generate modifications of the Schwarzschild geometry in several contexts (e.g., modified gravity, string theory, noncommutative geometry, generalized uncertainty principle scenarios), but a general analysis of the possible equations of state and their influence on the solutions was still lacking. Here, we aim to fill this gap in the literature and investigate whether these metrics might be able to reproduce features of the solutions of higher-derivative gravity models. In particular, we present an equation of state such that the solution matches the Newtonian-limit one in both regimes of large and small r . A significant part of the work is dedicated to studying the curvature regularity of the solutions and the comparison with the linearized solutions. Explicit metrics are presented for e ff ective sources originating from local and nonlocal models. The results obtained here might be regarded as possible links between the previous research on linearized higher-derivative gravity and the solutions of the nonlinear complete field equations, which remain unknown at the moment.", '1. INTRODUCTION': "Higher-derivative gravitational theories attract considerable attention due to their numerous applications in quantum gravity. Indeed, the simple integration of vacuum bubble diagrams of matter fields in curved spacetime produces curvature-squared terms in the e ff ective action [1], while the inclusion of the same type of structures in the classical gravitational action of general relativity renders the quantum theory renormalizable [2]. Adding more derivatives to the action enhances the ultraviolet (UV) behavior of the quantum theory even more. In fact, terms quadratic in curvatures and polynomial in derivatives can make the theory superrenormalizable [3]. Theories with an infinite number of derivatives may also be (super)renormalizable depending on the choice of the action nonlocal form factor [4-7].  \nThe literature on models with more than four metric derivatives in the gravitational action has increased substantially in the last decade. This is motivated by the possibility of solving the contradiction between renormalizability and unitarity present in the original formulation of fourth-derivative gravity. Even though the quantum properties of these models are well understood nowadays, the same cannot be said about their classical solutions.  \nThe aim of the present work is to systematically study a procedure for generating modifications of the Schwarzschild geometry inspired by higher-derivative gravity models and discuss the main features of the resulting metrics. Although these metrics are not solutions to the complete field equations of those models, they might reproduce characteristic properties of the fundamental solutions.  \nThe underlying higher-derivative models we focus on belong to a large class of theories quadratic in curvatures, described by the action  \nS = 1 16 π G ∫ d 4 x √ -g [ R + R µν F 1( /square ) R µν + RF 2( /square ) R ] , (1.1)  \nwhere G is the Newton constant and F 1 , 2( /square ) are form factors, analytic functions of the d'Alembert operator subjected to the constraint 1 F 1 /nequal -3 F 2. In what follows we list some popular choices of form factors and their respective quantum properties regarding renormalizability 2 and the presence of ghost-like degrees of freedom:  \n- (i) F 1( z ) = F 2( z ) ≡ 0 corresponds to the trivial case where we recover general relativity, which is a unitary but nonrenormalizable theory [8-11].\n- (ii) If F 1 , 2( z ) are nonzero constants we meet the fourth-derivative theory, also known as Stelle gravity [2]. This model is renormalizable, but it is nonunitary due to the presence of a spin-2 massive ghost (if quantized using the standard techniques of quantum field theory). Recent proposals for solving the problem of the ghost in the context of fourth-derivative gravity can be found, e.g., in [12-16].\n- (iii) If F 1 , 2( z ) are polynomials of degree N /greaterorequalslant 1 the theory is superrenormalizable [3]. The unitarity problem in these models can be tackled by choosing form factors such that all ghost-like degrees of freedom have complex masses and performing the quantization ' a la Lee-Wick [17-22]. This subclass of polynomial-derivative gravity is also known as Lee-Wick gravity [20, 21].\n- (iv) For nonlocal form factors of the type  \nF 2( /square ) = e H 2( /square ) -1 /square , F 1( /square ) = -e H 1( /square ) -1 /square -3 F 2( /square ) , (1.2)  \nwhere H 1 , 2( z ) are entire functions, the propagator of the theory only contains the massless pole corresponding to the graviton. Thus, these models are ghost-free at tree level and, depending on the choice of the functions H 1 , 2( z ), they can be (super)renormalizable [4-7].  \nQuadratic structures of the type R µναβ F 3( /square ) R µναβ and O ( R 3 ··· ) terms could also be included in the action (1.1) without substantially changing the general results concerning (super)renormalizability. Since the latter type of terms only a ff ect the interaction vertices (and not the propagator), provided that they do not contain more derivatives than the curvature-quadratic part, they cannot spoil the renormalizability. On the contrary, by a judicious choice of these structures, the theories described in (iii) and (iv) may even become finite at quantum level [3, 23]. As for the term quadratic in the Riemann tensor, it modifies both propagator and vertices in a very similar way to the structures in (1.1), and the analysis of renormalizability can be performed in an analogous manner (see, e.g., [3]). Nonetheless, the space of classical solutions can be significantly a ff ected by the presence of  \nthese two kinds of terms in the action. For the sake of simplicity, we do not consider such terms in the present work; instead, we focus only on a particular class of the model (1.1).  \nThe systematic search and study of spherically symmetric solutions in theories quadratic in curvature dates back to the 1970s [24-27] (see also [28, 29] and references therein for earlier considerations). Some authors have argued that the Schwarzschild solution is the only static, spherically symmetric, and asymptotically flat solution of fourth-derivative gravity in four dimensions (see, e.g., [30, 31]). However, in [32] it was demonstrated the existence of other solutions that deviate from the Schwarzschild one. After that, the properties of static spherically symmetric solutions in Stelle gravity became an active research topic [33-52]. The main results can be summarized very briefly as follows: In the complete fourth-derivative gravity, the solution that couples to a delta source contains a naked singularity; other solutions, corresponding to regular metrics, could describe wormholes. Event horizons may be present only in the particular case of Einstein-Weyl gravity.  \nWith respect to the models with more than four derivatives in the action or with nonlocal form factors, most of the solutions found in the literature are in the linear approximation [53-76]. 3 The study of the theory in the Newtonian limit is certainly interesting, for it gives the large-distance behavior of the metric-which should coincide with the asymptotic limit of the nonapproximated solution. The small-distance behavior of the linearized solutions has also been the subject of intensive studies, and it was shown that the modified Newtonian potential is bounded and the curvature scalars are regular in most of the higherderivative theories. To be more precise, all models of the type (1.1) with nonzero F 1( z ) and F 2( z ) have a finite Newtonian potential [53, 54]; while if F 1 , 2( z ) are nontrivial polynomials (i.e., the action has six or more derivatives) the model admits a regular Newtonian-limit solution, without singularities in the curvature invariants [55]. The generalization of these statements to nonlocal theories can be found in [56]. Basically, what is important for the absence of singularities is the UV behavior of the form factors [56, 57]. If H 1 , 2( z ) in (1.2) is such that F 1 , 2( z ) behave like a polynomial of degree N at high energies (like in the cases of Kuz'min and Tomboulis form factors [5, 6]), the nonlocal theory reproduces the regularity properties of the corresponding (local) polynomial model. If F 1 , 2( z ) grow faster than any polynomial (for large z ), not only the curvature invariants are regular at r = 0, but so are all the invariants that are polynomial in derivatives of curvature tensors (collectively known as curvature-derivative invariants [56, 83, 84]). For the detailed consideration of curvature-derivative invariants, see [56, 83].  \nEven though the linear version of the theory does not seem to be the appropriate setting to investigate whether higher derivatives could resolve the singularity, it has been speculated that the solution of the full nonlinear theory might be regular as well. It was even ventilated that close to r = 0 the gravitational field may be weak enough so that one could trust the linearized solution in a small region around the origin [85]-or, under certain circumstances, from far away up to the origin [76].  \nTo the best of our knowledge, static spherically symmetric solutions of the full nonlinear equations of motion of curvaturequadratic gravity theories with more than four derivatives of the metric have been considered only in [86], based on the Frobenius method and numerical analysis of the field equations in models with up to ten derivatives. 4 The conclusion of this study can be summarized as follows:  \n- (I) All the solutions found are regular at r = 0.\n- (II) There is a critical mass M c for the formation of a black hole. If M < M c the solution does not have an event horizon. On the other hand, for M > M c the solution possesses more than one horizon.  \nMoreover, it was speculated whether some appealing properties of the regular solutions of fourth-derivative gravity could be shared by their higher-derivative counterparts, namely:  \n- (III) Nonlinearities are suppressed su ffi ciently close to r = 0; the solutions behave like the weak-field ones.\n- (IV) The gravitational field is strong near the position of the horizons; in these regions, there is a transition between the linear and nonlinear regimes.  \nLast but not least, owing to the complexity of the equations of motion of higher-derivative gravity, the works [7, 90-92] considered a truncated form of the field equations associated to specific models defined by the action  \nS = 1 16 π G ∫ d 4 x √ -g [ R + G µν F ( /square ) R µν ] , (1.3)  \nwhere G µν is the Einstein tensor. The action (1.3) corresponds to a particular class of theories of the type (1.1), namely, the ones whose form factors satisfy the relation  \nF 1( /square ) = -2 F 2( /square ) ≡ F ( /square ) . (1.4)  \nwhere  \nf ( z ) = 1 + zF ( z ) . (1.7)  \nIn an equivalent form, ϕ ( r ) can be obtained by solving the equation  \n∆ ϕ = 4 π G ρ e ff , (1.8)  \nwhere  \nρ e ff ( r ) = M 2 π 2 r ∫ ∞ 0 d k k sin( kr ) f ( -k 2 ) (1.9)  \nis the smeared delta source , which follows from the inversion of the operator f ( ∆ ) in (1.6). 5 Therefore, Eq. (1.6) can be cast as a standard Poisson equation (1.8) with the e ff ective delta source (1.9) (see, e.g., [57, 84]). Hence, in flat spacetime, the e ff ect of the higher derivatives may be regarded as inducing a smearing of the original delta source, through the nonconstant function f ( z ) in the integrand of (1.9).  \nThe main idea, now, is to obtain a curved-space generalization of (1.8) to static spherically symmetric spacetimes in the form  \nG µ ν = 8 π G ˜ T µ ν, (1.10)  \nwhere ˜ T µ ν is an e ff ective energy-momentum tensor reproducing the e ff ects of the higher derivatives with  \n˜ T t t = -ρ e ff , (1.11)  \nbeing ρ e ff the smeared delta source (1.9). In this way, like Eq. (1.8) yields the modified Newtonian potential, the metrics which follow from (1.10) can be interpreted as higher-derivative modifications of Schwarzschild, through the function f ( z ). In fact, if the higher derivatives are switched o ff , f ( z ) ≡ 1 and Eq. (1.10) results in the Schwarzschild solution.  \nThe formula (1.10) can also be seen as a truncation of the full equations of motion for the theory (1.3) sourced by a pointlike source. Indeed, owing to the specific form of the higher-derivative structure in (1.3), it is possible to factor the operator f ( /square ) together with the Einstein tensor in the equations of motion, namely,  \nf ( /square ) G µ ν + O ( R 2 ··· ) = 8 π GT µ ν, (1.12)  \nwhere T µ ν is the matter usual energy-momentum tensor. Therefore, by discarding the terms O ( R 2 ··· ) of higher order in curvatures in (1.12) and taking the d'Alembertian in f ( /square ) in the flat-spacetime form, from the inversion of the operator f ( /square ) one can recast (1.12) into (1.10) [7, 85, 90-95]. Although frequently used in the literature, this procedure of truncating the d'Alembertian (to obtain the invertible form that does not depend on the metric) is seldom discussed. In this respect, it was argued that these truncations of the original equation (1.12) might be compensated by imposing the conservation of the e ff ective energymomentum tensor, i.e., ∇ µ ˜ T µ ν = 0 (see, e.g., [7]). This leads to the introduction of e ff ective radial ( pr ) and tangential ( p θ ) pressures in the e ff ective energy-momentum tensor,  \n˜ T µ ν = diag( -ρ e ff , pr , p θ, p θ ) . (1.13)  \nThe truncation of the equations of motion is sometimes referred to as the propagator approximation , or smeared delta source approximation . Since the present work is devoted precisely to solutions of this type, the remaining part of this introductory section describes the procedure and the underlying assumptions.  \nSimply put, the procedure consists in a curved-spacetime generalization of the method for calculating the Newtonian potential of a pointlike source in higher-derivative models. It can be shown (see, e.g., [84] or Sec. 3 below) that in the Newtonian limit, the solution of the equations of motion associated to (1.3) only depends on one potential , ϕ ( r ). In fact, the linear limit metric written in Schwarzschild coordinates is given by  \nd s 2 = -[ 1 + 2 ϕ ( r ) ] d t 2 + [ 1 + 2 r ϕ ' ( r ) ] d r 2 + r 2 d Ω 2 , (1.5)  \nwhere d Ω 2 = d θ 2 + sin 2 θ d φ 2 and primes denote di ff erentiation with respect to the coordinate r . The potential for a pointlike source is thus the solution of the modified Poisson equation (see, e.g., [55, 84])  \nf ( ∆ ) ∆ ϕ = 4 π G ρ, ρ ( /vector r ) = M δ (3) ( /vector r ) , (1.6)  \nAs we show in the present work, although the spacetime configurations constructed in this framework are not likely the solutions of the complete field equations, they reproduce relevant aspects of the findings of [86].  \nIn the present work, we study static spherically symmetric solutions to the system of equations (1.10), with the general line element  \nd s 2 = -A ( r ) e B ( r ) d t 2 + d r 2 A ( r ) + r 2 d Ω 2 , (1.14)  \nwhere A ( r ) and B ( r ) are two functions to be determined. Di ff erently from previous works [90, 92], here we are interested in solutions with a nontrivial shift function, B ( r ) /nequal 0. The motivation comes from the observation that the identity ϕ ( r ) = -r ϕ ' ( r ) (which is valid in Einstein gravity) for the Newtonian potential does not hold in general for higher-derivative models, thus to reproduce the weak-field behavior (1.5) we need a metric with gtt /nequal -g rr . 6  \nWe mainly focus on form factors that grow at least as fast as f ( z ) ∼ z 2 for large z . This corresponds to the smeared delta sources associated to models whose actions have at least six derivatives of the metric tensor [see (1.3) and (1.7)]. In the main part of this paper, we show that under some reasonable assumptions, the modified Schwarzschild metric that follows from the simple equation (1.10) can reproduce the properties (I)-(IV) of the static spherically symmetric solutions of [86]. By imposing some physical requirements on the spatial part of ˜ T µ ν , we also show that for all form factors such that f ( z ) ∼ z 2 or faster for large z , the solutions of (1.10) have curvature invariants without singularities. This result is at the core of a procedure to generate a large amount of regular black hole metrics. 7  \nIt is worth pointing out that the metrics we obtain here, regarded as higher-derivative modifications of Schwarzschild, are similar in spirit to the nonlocal modifications of the Kerr and Schwarzschild metrics recently studied in [82]. In fact, they are not solutions of higher-derivative gravity models, but encode higher-derivative modifications at equations whose original forms would result in the Schwarzschild metric. Our procedure, however, di ff ers from the one of [82]; for instance, their modification of the Schwarzschild metric takes the form A ( r ) = 1 + 2 ϕ ( r ) and B ≡ 0 (in our notation).  \nFinally, let us mention that the formulation of the problem in terms of the e ff ective field equations (1.10) with smeared delta sources is standard in several other approaches to quantum gravity such as noncommutative geometry, generalized uncertainty principle models, and string theory (see, e.g., [69, 102-112]). Therefore, the detailed study of the solutions of Einstein's equation (1.10) with the e ff ective source (1.9) is relevant not only for higher-derivative gravities but also for alternative models of quantum gravity.  \nThe paper is organized as follows. In Sec. 2 we present the equations for the metric components A ( r ) and B ( r ), which follows from the e ff ective field equations (1.10), and obtain their formal solutions depending on a generic form factor f ( z ). Then, a brief review of the linearized field equations and their solutions is carried out in Sec. 3. In Sec. 4 we discuss the curvature regularity of the solutions of the e ff ective field equations on general grounds, i.e., without particularizing to specific form factors. We also discuss the e ff ect of the nonlinearity of the field equations on the regularity of the solutions. From Sec. 5, the considerations start to become model-dependent. In that section, we study the general properties of the solution for A ( r ), showing how the form factor f ( z ) can define the regularity of A ( r ) and the possible structure of horizons of the metric. The similar analysis for the function B ( r ) is carried out in Sec. 6, where we discuss several choices of equations of state for the e ff ective pressures and their influence on the solutions. In particular, in Sec. 6.4 we introduce an equation of state such that the solution of the e ff ective field equations interpolates between a regular core, for small r , and the Newtonian-limit solution at large r . Examples where the scenarios described in the previous sections can be explicitly visualized are presented in Sec. 7, which contains one case of nonlocal form factor and three cases of local higher-derivative gravity with increasing complexity (up to the case of the most general polynomial form factor, including an arbitrary number of complex and multiple roots). Finally, Sec. 8 contains a summary of the results and our conclusion. Throughout this paper, we use the same sign conventions of [84], and we adopt the unit system such that c = 1 and /planckover2pi = 1.", '2. EFFECTIVE FIELD EQUATIONS: THE GENERAL SOLUTION FOR A ( r ) AND B ( r )': "As mentioned in the Introduction, our discussion of the higher-derivative modifications of the Schwarzschild metric is based on the system of equations (1.10). On the left-hand side we have the Einstein tensor, whose nonzero components evaluated for the generic metric (1.14) are  \nG t t = A ' r + A r 2 -1 r 2 , (2.1)  \nand  \nG θ θ = G φ φ = A ' r + 1 2 AB ' r + 3 4 A ' B ' + 1 4 AB ' 2 + 1 2 A '' + 1 2 AB '' . (2.3)  \nTherefore, using the e ff ective energy momentum tensor (1.13), the field equations (1.10) are equivalent to  \nA ' r + A r 2 -1 r 2 = -8 π G ρ e ff , (2.4)  \nAB ' r = 8 π G ( ρ e ff + pr ) (2.5)  \nand  \nA ' r + 1 2 AB ' r + 3 4 A ' B ' + 1 4 AB ' 2 + 1 2 A '' + 1 2 AB '' = 8 π Gp θ. (2.6)  \nThe first and the last equations are just the tt and θθ components of (1.10), while (2.5) is obtained from the rr equation by subtracting (2.1) from (2.2).  \nThe formal solution for A ( r ) can be directly obtained from Eq. (2.4). It only depends on the e ff ective delta source, and it is given by  \nA ( r ) = 1 -2 Gm ( r ) r , (2.7)  \nwhere the mass function  \nm ( r ) = 4 π ∫ r 0 d x x 2 ρ e ff ( x ) (2.8)  \nrepresents the e ff ective mass inside a sphere of radius r .  \nGiven the solution (2.7) for A ( r ), the function B ( r ) can be obtained through (2.5). The formal solution reads  \nB ( r ) = B 0 + 8 π G ∫ r r 0 d x x [ ρ e ff ( x ) + pr ( x ) ] A ( x ) . (2.9)  \nThe constant B 0 = B ( r 0) is responsible for a nontrivial time delay between the position r and an observer at r 0 (for more details, see the discussion in [97]). Of course, one can always fix B 0 by imposing the asymptotic condition lim r →∞ B ( r ) = 0, so that the metric reproduces the Minkowski geometry at infinity-provided that the e ff ective pressures tend to zero su ffi ciently fast for large r . The constant B 0 can also be changed by a redefinition of the time coordinate t , and for these reasons hereafter we shall take its value according to the convenience without further elaborations.  \nFinally, since the solution for the function B ( r ) depends on the choice of the pressure components of the e ff ective energymomentum tensor ˜ T µ ν , it is mandatory to comment on these other components. One of the simplest possible choices is to fix the e ff ective radial pressure pr via an equation of state involving this component and the e ff ective source ρ e ff . In this case, both A ( r ) and B ( r ) are directly determined in terms of the source ρ e ff . Then, the component p θ is given by the last equation of motion (2.6). This is completely equivalent to fixing p θ using the e ff ective conservation equation ∇ µ ˜ T µ r = 0, namely,  \np ' r = -1 2 ( A ' A + B ' ) ( pr + ρ e ff ) -2 r ( pr -p θ ) , (2.10)  \nfrom which one can easily determine p θ once ρ e ff and pr (and, consequently, A and B ) are specified. Indeed, if this equation is satisfied, the tt and rr components of (1.10) imply in (2.6). In Sec. 6, we discuss possible choices of equations of state in the form pr = pr ( ρ e ff ), their motivations, and the properties of the corresponding function B ( r ).  \nAs a conclusion of this section, the solution of the e ff ective field equations (1.10) can be expressed as  \nd s 2 = -( 1 -2 Gm ( r ) r ) exp ( B 0 + 8 π G ∫ r r 0 d x x [ ρ e ff ( x ) + pr ( x ) ] A ( x ) ) d t 2 + ( 1 -2 Gm ( r ) r ) -1 d r 2 + r 2 d Ω 2 . (2.11)  \nG r r = G t t + AB ' r , (2.2)", '3. NEWTONIAN LIMIT': "Here we make a brief digression to show how to reproduce the Newtonian-limit metric (1.5) using the field equations (2.4)(2.6). In the linear approximation we write  \nA ( r ) -→ 1 + a ( r ) , B ( r ) -→ b ( r ) (3.1)  \n↦  \n↦  \nand only keep quantities up to the first order in a ( r ), b ( r ) and their derivatives. Also, in the nonrelativistic limit, one sets pr = p θ = 0. Thus, the linearized Eqs. (2.4) and (2.5) are  \na ( r ) r 2 + a ' ( r ) r = -8 π G ρ e ff ( r ) , (3.2)  \nb ' ( r ) r = 8 π G ρ e ff ( r ) , (3.3)  \nwhose solutions are the functions entering the linearized version of the metric (1.14), i.e.,  \nd s 2 = -[1 + a ( r ) + b ( r )]d t 2 + [1 -a ( r )]d r 2 + r 2 d Ω 2 . (3.4)  \nSince Eq. (2.4) is already linear in the metric function A ( r ), there is actually no approximation involved in (3.2), which is only a rewriting of the former in terms of the function a ( r ) = A ( r ) -1. Its solution is  \na ( r ) = -2 Gm ( r ) r , (3.5)  \nwith m ( r ) defined in (2.8), while b ( r ) is easily obtained from (3.3),  \nb ( r ) = b 0 + 8 π G ∫ r r 0 d x x ρ e ff ( x ) . (3.6)  \nFinally, since ρ e ff , a ( r ), and b ( r ) are proportional to M , the nontrivial component (2.10) of the continuity equation,  \n∇ µ ˜ T µ r = 1 2 ( A ' A + B ' ) ρ e ff = 1 2 ( a ' + b ' ) ρ e ff , (3.7)  \nis already O ( M 2 ), implying that it is verified within the approximation.  \nThe Newtonian potential ϕ ( r ) is defined through the relation  \nthus from (3.4) we obtain  \ngtt = -[ 1 + 2 ϕ ( r ) ] , (3.8)  \nϕ ( r ) = a ( r ) + b ( r ) 2 . (3.9)  \nTaking the derivative of the previous equation, we find  \nϕ ' ( r ) = a ' ( r ) + b ' ( r ) 2 = -a ( r ) 2 r , (3.10)  \nwhere we used (3.2) and (3.3). Therefore,  \na ( r ) = -2 ϕ ' ( r ) r and b ( r ) = 2 [ ϕ ( r ) + ϕ ' ( r ) r ] . (3.11)  \nSubstituting these expressions in (3.4), we obtain the Newtonian metric in the form (1.5), presented in the Introduction.", '4. CURVATURE REGULARITY AND SUPPRESSION OF NONLINEARITIES NEAR r = 0': 'Before discussing the general properties of the solutions for A ( r ) and B ( r ) obtained in Sec. 2, it is instructive to review what conditions such functions should satisfy to regularize the curvature and curvature-derivative invariants. We also elaborate more on the appealing idea that nonlinearities are suppressed close to r = 0. We show that under some general assumptions on A ( r ) and B ( r ), the behavior of the curvature invariants for su ffi ciently small r is, indeed, well approximated by the Newtonian-limit solutions.', '4.1. Regularity of curvature invariants': "The regularity analysis of the curvature invariants constructed only by contracting an arbitrary number of Riemann and metric tensors can be greatly simplified by noticing that any of such scalars can be expressed as a combination of the components R µν αβ (i.e., as a sum of products of these components only, without the need of using the metric tensor) [101]. For a static spherically symmetric metric, R µν αβ has only four independent components, which in the case of a metric in the form (1.14) are given by  \nK 1 ≡ R tr tr = -A '' 2 -1 2 AB '' -3 4 A ' B ' -1 4 AB ' 2 , (4.1a)  \nK 2 ≡ R t θ t θ = R t φ t φ = -A ' + AB ' 2 r , (4.1b)  \nK 3 ≡ R r θ r θ = R r φ r φ = -A ' 2 r , (4.1c)  \nK 4 ≡ R θφ θφ = -A -1 r 2 . (4.1d)  \nIf the four functions Ki ( r ) are regular, all the curvature invariants without covariant derivatives are bounded as well.  \nIn an equivalent way, one can use the Kretschmann scalar,  \nK ≡ R µναβ R µναβ = 4 K 2 1 + 8 K 2 2 + 8 K 2 3 + 4 K 2 4 , (4.2)  \nto investigate the regularity of a whole set of curvature invariants. Indeed, since it is the sum of squares, for its regularity it is necessary and su ffi cient that all functions Ki ( r ) are finite. Thus, if K is bounded, all the Ki ( r ) are bounded as well and, as a consequence, so is any scalar formed only by curvature contractions [101] (e.g., the scalar curvature R , the square of the Ricci tensor R µν R µν , the square of the Weyl tensor C µναβ C µναβ or any polynomial scalar constructed with them). However, from the practical point of view, in order to verify the regularity of a given metric, it is simpler to work with the short expressions in (4.1) rather than the one in (4.2).  \nTo analyze the behavior for small r , let us assume that A ( r ) and B ( r ) are analytic 8 around r = 0 and apply the power series expansion,  \nA ( r ) = ∞ ∑ n = 0 anr n , B ( r ) = ∞ ∑ n = 0 bnr n , (4.3)  \ninto the expressions for the functions Ki ( r ):  \nK 1 = -a 2 -a 0 b 2 -b 1(3 a 1 + a 0 b 1) 4 + O ( r ) , (4.4a)  \nK 2 = -a 1 + a 0 b 1 2 r -a 2 -a 1 b 1 2 -a 0 b 2 + O ( r ) , (4.4b)  \nK 3 = -a 1 2 r -a 2 + O ( r ) , (4.4c)  \nK 4 = -1 -a 0 r 2 -a 1 r -a 2 + O ( r ) . (4.4d)  \nTherefore, the regularity of the curvature invariants requires  \na 0 = 1 and a 1 = b 1 = 0 , (4.5)  \nor, in other words,  \nA (0) = 1 , (4.6)  \nA ' (0) = B ' (0) = 0 . (4.7)  \nIn the next sections we show that these conditions are fulfilled for the general solution (2.11) in any theory whose form factor is such that f ( -k 2 ) ∼ k 4 or faster for large k -for example, for e ff ective sources of models that contain six or more metric derivatives.  \nTo close this section, we briefly comment on the regularity of curvature-derivative invariants. As highlighted in Refs. [56, 83], scalars that are polynomial in curvature tensors and their covariant derivatives can display a singular behavior at r = 0 if the regular metric of the form (1.14) contains odd powers of r in the Taylor expansion (4.3) of its components. Thus, some of the considerations can be stated in a simple form by using the definition of order of regularity of a function, as introduced in [56]:  \nDefinition 1. Given a continuous bounded function ξ : [0 , + ∞ ) → R , we say that ξ ( r ) is p-regular for an integer p /greaterorequalslant 1 if ξ ( r ) is at least 2 p -times continuously di ff erentiable and the first p odd-order derivatives of ξ ( r ) vanish as r → 0, namely,  \n0 /lessorequalslant n /lessorequalslant p -1 = ⇒ lim r → 0 ξ (2 n + 1) ( r ) = 0 . (4.8)  \nIf these conditions are satisfied, we shall also say that p is the order of regularity of ξ . This extends the definition of a regular, continuous function (which can also be called 0-regular) in a way that is convenient for studying the regularity of the spacetime.  \nIndeed, if a function ξ : [0 , + ∞ ) → R is p -regular, this means that also the map Ξ : R 3 → R given by Ξ ( /vector r ) = ξ ( | /vector r | ) is of class C 2 p . On the other hand, if lim r → 0 ξ (2 p + 1) ( r ) /nequal 0, then (switching to Cartesian coordinates) we have  \nlim x → 0 ∂ 2 p + 1 x Ξ ( x , 0 , 0) = ± ξ (2 p + 1) (0) , (4.9)  \ndepending on whether x = 0 is approached from the left or the right. This means that Ξ is not (2 p + 1)-times di ff erentiable, and if we take one of the metric functions A ( r ) or B ( r ) for ξ , this implies that the metric is not smooth at r = 0.  \nGeneral evidences of a relation between the order of the first odd power of r in the series and the minimal number of covariant derivatives in a divergent curvature-derivative scalar were presented in the works [56, 83]. (See [56] for a complete treatment of the problem at linear level, and [83] for the detailed consideration of invariants of the type /square n R at nonlinear level and further examples.) Taking into account these results, if the first odd power of r in the series expansion of the functions A ( r ) and B ( r ) is at order r 2 n + 1 for some n ∈ N , then there might exist curvature-derivative scalars with 2 n covariant derivatives of curvature tensors that diverge at r = 0. All local curvature-derivative scalars are expected to display a nonsingular behavior only if A ( r ) and B ( r ) are ∞ -regular, i.e., if the metric components are even functions .  \nA simple useful example is provided by the scalar /square R , which for the metric (1.14) reads  \n/square R = 1 4 r 4 ( 16 A -16 A 2 -16 rA ' + 32 rAA ' + 8 r 2 A ' 2 -8 rAB ' + 8 rA 2 B ' + 12 r 2 AA ' B ' -8 r 3 A ' 2 B ' + 4 r 2 A 2 B ' 2 -8 r 3 AA ' B ' 2 -2 r 4 A ' 2 B ' 2 -r 4 AA ' B ' 3 -8 r 2 AA '' -16 r 3 A ' A '' -28 r 3 AB ' A '' -6 r 4 A ' B ' A '' -5 r 4 AB ' 2 A '' -44 r 3 AA ' B '' -10 r 4 A ' 2 B '' -12 r 3 A 2 B ' B '' -17 r 4 AA ' B ' B '' -2 r 4 A 2 B ' 2 B '' -16 r 4 AA '' B '' -4 r 4 A 2 B '' 2 -24 r 3 AA (3) -4 r 4 A ' A (3) -8 r 4 AB ' A (3) -16 r 3 A 2 B (3) -18 r 4 AA ' B (3) -6 r 4 A 2 B ' B (3) -4 r 4 AA (4) -4 r 4 A 2 B (4) ) . (4.10)  \nOnce more, to investigate the behavior of the invariant near r = 0 we use the power series representation (4.3), already assuming a 0 = 1 and a 1 = b 1 = 0. The result is  \n/square R = -8(5 a 3 + 3 b 3) r + O ( r 0 ) . (4.11)  \nHence, the scalar /square R diverges like 1 / r unless 5 a 3 + 3 b 3 = 0; in particular, it is regular if a 3 = b 3 = 0-in other words, if A ( r ) and B ( r ) are 2-regular, as noted in [56, 83]. In the following sections, we show that there exist equations of state for the e ff ective pressures such that this condition is fulfilled in any theory whose form factor is such that f ( -k 2 ) ∼ k 6 or faster for large k [see, e.g., Eqs. (5.14) and (6.3) below], i.e., for smeared delta sources of models with eight or more derivatives. We remark, however, that the regularity of curvature-derivative scalars, which depend on higher-order derivatives of the metric, can be very sensitive to the choice of equation of state; an example of this is presented in Appendix A.", '4.2. Suppression of nonlinearities': "As stated in the Introduction, some authors have argued that close to r = 0, the gravitational field may be very weak, so the linearized solution reproduces the fate of the singularities. Now, we elaborate more on this idea, showing that this can be the case under some assumptions on the form of the metric components.  \nFirst, let us suppose that the spacetime is asymptotically flat such that we can write  \nA ( r ) = 1 + a ( r ) , (4.12)  \nwith lim r →∞ a ( r ) = 0. We also assume that a ( r ), B ( r ), and their derivatives are proportional to a parameter such that the term 'linearization' is understood as the evaluation of quantities at leading order in these functions, since they are already linear in the parameter (i.e., there is no linearization inside these functions). In the explicit examples we consider in Sec. 7, this parameter can be taken, e.g., as the mass 9 M .  \nThe next step is to apply this procedure and calculate the linearized version of the curvature tensor, ( R µν αβ )lin. Using (4.1), its components read  \nK lin 1 = -a '' + B '' 2 , (4.13a)  \nK lin 2 = -a ' + B ' 2 r , (4.13b)  \nK lin 3 = K 3 , (4.13c)  \nK lin 4 = K 4 . (4.13d)  \nNotice that since K 3 and K 4 are already linear in A ( r ) and do not depend on B ( r ), their linearized expressions coincide with the original ones.  \nThe nonlinear part of R µν αβ is then proportional to the di ff erences  \nK non -lin i ≡ Ki -K lin i , (4.14)  \nwhich turn out to be proportional to the derivatives of B ( r ):  \nK non -lin 1 = 1 4 [ 3 a ' B ' + (1 + a ) B ' 2 + 2 aB '' ] , (4.15a)  \nK non -lin 2 = aB ' 2 r , (4.15b)  \nK non -lin 3 = K non -lin 4 = 0 . (4.15c)  \nIn the same spirit as before, using the power series representation (4.3) for A ( r ) and B ( r ) we find, around r = 0,  \nK non -lin 1 = 1 4 b 1(3 a 1 + a 0 b 1) + ( a 0 -1) b 2 + O ( r ) , (4.16a)  \nK non -lin 2 = ( a 0 -1) b 1 2 r + a 1 b 1 2 + ( a 0 -1) b 2 + O ( r ) . (4.16b)  \nTherefore, if the curvature invariants are regular, i.e., if the relations in (4.5) are true, we get  \nK non -lin i = O ( r ) , (4.17)  \nshowing that for small r the nonlinearities in R µν αβ are suppressed. Since all the invariants polynomial in curvatures and metric tensors can be build using combinations of the functions Ki , it follows that, for small r , such scalars behave approximately like the linearized ones.  \nFurthermore, if the functions a ( r ) and B ( r ) approach the Newtonian solution for small r , not only are the nonlinearities suppressed, but the value of the curvature invariants evaluated at r = 0 is equal to the Newtonian one. This happens, indeed, with the general solution (2.11) if the e ff ective pressures tend to zero su ffi ciently fast for small r . As shown in Sec. 3, the solution a ( r ) of the Newtonian equation of motion (3.2) is the same one entering A ( r ) = 1 + a ( r ) that solves the nonlinear equations of motion (1.10); while for B ( r ), given that A (0) = 1, it is always possible to find a small enough r such that 1 / A ( r ) ≈ 1 and, as a result, the solution (2.9) approaches (3.6) if pr vanishes su ffi ciently fast as r → 0.  \nIt is also worthwhile to note that if B ( r ) = 0, then K non -lin i = 0 always, even for singular metrics . This explains why in general relativity the Kretschmann invariant evaluated using either the Schwarzschild solution or the Newtonian-limit metric gives the same result  \nK Sch. = 48 G 2 M 2 r 6 . (4.18)  \na ( r ) = -2 Mr 2 r 3 + M 2 L , B ( r ) = 0 .  \nIn what concerns the curvature-derivative invariants, the situation is slightly di ff erent. Using once more /square R as an example and proceeding in the same way as before, we split (4.10) in terms of the linear and nonlinear parts,  \n( /square R )lin = -4 a r 4 + 4 a ' r 3 -2 a '' r 2 -6 a (3) + 4 B (3) r -a (4) -B (4) , (4.19)  \n( /square R )non -lin ≡ /square R -( /square R )lin = a ( -4 a r 4 + 8 a ' r 3 -2 a '' r 2 -6 a (3) r -a (4) ) + a ' ( 2 a ' r 2 -4 a '' r -a (3) ) -B ' [ 2 a ' 2 r + 1 2 a ' 2 B ' + 3 2 a ' a '' -(1 + a ) ( 2 a r 3 + (1 + a ) B ' + 3 a ' r 2 -2 a ' B ' -7 a '' r -1 4 a ' B ' 2 -5 4 B ' a '' -2 a (3) )] -B '' [ 5 2 a ' 2 + (1 + a ) ( 11 a ' + 3(1 + a ) B ' r + 4 a '' + 17 4 a ' B ' + 1 2 (1 + a )( B ' 2 + 2 B '' ) )] -B (3) [ 4 a (2 + a ) r + 9 2 aa ' + 3 2 (1 + a ) 2 B ' ] -B (4) a (2 + a ) , (4.20)  \nand using the series representation (4.3) for regular metrics (i.e., with a 0 = 1, a 1 = b 1 = 0) we get  \n( /square R )lin = -8(5 a 3 + 3 b 3) r -60(3 a 4 + 2 b 4) + O ( r ) , (4.21)  \n( /square R )non -lin = -12 b 2 (6 a 2 + b 2) + O ( r ) . (4.22)  \nTherefore, if /square R is finite, that is it, if 5 a 3 + 3 b 3 = 0, we see that its linear and nonlinear parts approach a constant when r → 0. The behavior of ( /square R )non -lin marks a di ff erence between curvature-derivative and curvature invariants since the nonlinear part of the latter goes to zero when r → 0. Note, however, that if ( /square R )lin is finite at r = 0, the entire /square R is also finite (because its nonlinear part cannot generate a singularity).  \nIt is also interesting to notice that, because of the covariant derivatives, ( /square R )non -lin contains terms that do not depend on the function B ( r ) [see Eq. (4.20)]. This is di ff erent from what happens in the case of scalars without covariant derivatives. Therefore, even if B ( r ) = 0, /square R can receive a contribution from the nonlinear part of the metric.", '5. GENERAL PROPERTIES OF A ( r )': "In this section, we discuss some properties of (2.7) at a general level, without particularizing to any gravitational model e ff ective source. As discussed in the Introduction, by just considering the asymptotic behavior of the action form factor F ( z ) [or, equivalently, f ( z ) in (1.7)], we can explain many important physical properties of the solutions that hold true for large classes of higher-derivative e ff ective sources. For instance, we are particularly interested in the absence (or the presence) of an event horizon, of singularities in the curvature invariants, as well as in the behavior of the solution for large r . All the results presented in general form in this section can be seen realized in concrete examples in Sec. 7, where we choose some representative form factors and obtain the explicit expressions.  \nTo start, let us recall some relevant results on the smeared delta source ρ e ff ( r ); afterward, we analyze the consequences of these properties to the mass function (2.8) and, finally, to A ( r ). The properties of the e ff ective sources in higher-derivative gravity theories have been intensively studied in the works [56, 57] (see also [84] for an introduction and review). For this reason, we do not repeat their derivation here; we summarize the results as follows.  \n(1a) If f ( -k 2 ) → 1 when k 2 → 0 (in other words, if general relativity is recovered in the IR limit), then  \nlim r →∞ ρ e ff ( r ) = 0 . (5.1)  \nIn particular, this is related to the solution being asymptotically flat.  \nUnder the assumption (1a) one can also prove that:  \n- (2a) If there exists a k 0 such that f ( -k 2 ) grows as k 2 for k > k 0, then ρ e ff ( r ) is integrable but it is not bounded. In fact, it diverges at the origin,  \nρ e ff ( r ) ∼ r → 0 1 r . (5.2)  \nThis case corresponds to the smeared delta sources for Stelle gravity and nonlocal theories that behave as the fourthderivative theory in the UV. This type of source was also considered in the study of black holes in models with generalized uncertainty principle [112].  \n(3a) On the other hand, if f ( -k 2 ) grows at least as fast as k 4 for su ffi ciently large k , the e ff ective source is integrable and finite everywhere, i.e., ρ e ff ( r ) is at least 0-regular. In this situation, ρ e ff ( r ) has an absolute maximum at r = 0,  \nmax { ρ e ff ( r ) } = lim r → 0 ρ e ff ( r ) . (5.3)  \nThis is the case of local models with six or more metric derivatives and most nonlocal models (namely, those whose propagator is suppressed in the UV at least as fast as in the sixth-derivative gravity).  \n(4a) If the function f ( -k 2 ) asymptotically grows at least as fast as k 4 + 2 N for an integer N /greaterorequalslant 1, then the e ff ective source ρ e ff ( r ) is (at least) N -regular (see Definition 1 in Sec. 4.1). In other words, it is 2 N times di ff erentiable and their first N odd-order derivatives vanish at the origin,  \nd 2 n + 1 d r 2 n + 1 ρ e ff ( r ) ∣ ∣ ∣ ∣ r = 0 = 0 for n = 0 , · · · , N -1 . (5.4)  \n(5a) As a corollary of the previous statement, if f ( -k 2 ) asymptotically grows faster than any polynomial, then ρ e ff ( r ) is ∞ -regular, i.e., it is an even function of r . This happens, e.g., in the nonlocal models with exponential form factors, f ( -k 2 ) = e ( k 2 /µ 2 ) /lscript for /lscript ∈ N .  \nFrom the above results it follows that, for any higher-derivative model form factor, the mass function satisfies [57]  \nlim r →∞ m ( r ) = M . (5.5)  \nTherefore, for a distant observer, the solution (2.7) for A ( r ) reduces to Schwarzschild 10 with mass M . On the other hand, for small r , the limits (5.2) and (5.3) together with the definition (2.8) imply that  \nm ( r ) ∼ r → 0 { r 2 , if f ( -k 2 ) ∼ k 2 for large k , r 3 , if f ( -k 2 ) ∼ k 4 (or faster) for large k . (5.6)  \nAs consequence, combining (5.5) and (5.6) with (2.7) we get  \nlim r → 0 A ( r ) = lim r →∞ A ( r ) = 1 . (5.7)  \nNotice that this proves the regularity condition of Eq. (4.6). Being a continuously di ff erentiable function, and given (5.7), it follows that A ( r ) is bounded for any higher-derivative model e ff ective source. Thus, there is a certain critical mass M c such that M < M c implies A ( r ) > 0 for all r .  \nThe existence of M c means that there is a mass gap for the solution to be a black hole. Indeed, the presence of an event horizon is related to the invariant  \n( ∇ r ) 2 = A ( r ) , (5.8)  \nthat only depends on the function A ( r ). The event horizon is given by the largest root of A ( r ) = 0, and the metric does not describe a black hole if A ( r ) is strictly positive.  \nA direct way of verifying the existence of the critical mass M c for e ff ective sources is to consider the unitary mass function, defined by  \n˜ m ( r ) = m ( r ) M , (5.9)  \nwhich does not depend on M [see Eq. (1.9)] and it is such that lim r →∞ ˜ m ( r ) = 1. Notice that the function ˜ m ( r ) / r must assume a maximum value at some point r ∗ (since it tends to zero in the limits of large and small r , see Eq. (5.7), and there is a region where it is positive). Thus, if M is such that  \n2 GM ˜ m ( r ∗ ) r ∗ < 1 , (5.10)  \nthe function A ( r ) is strictly positive, and the metric does not have horizons. The critical mass is, therefore,  \nM c = r ∗ 2 G ˜ m ( r ∗ ) . (5.11)  \nTo summarize, if f ( -k 2 ) grows at least as fast as k 2 for su ffi ciently large k , that is it, for smeared delta sources of any higherderivative gravity theory , there exists a mass gap for a solution describing a black hole.  \nIn what concerns the solutions that do describe black holes (i.e., if M > M c), we notice that, because of (5.7), the function A ( r ) changes sign an even number of times. This means that the solution also contains at least one inner horizon besides the event horizon. The exact number of horizons depends on the function f ( -k 2 ) and on the value of M , but it is always an even number. In the limiting case M = M c, a pair of horizons merge into one extremal horizon.  \nIn addition, it is worth noticing that if ˜ m ( r ) /lessorequalslant 1 for all r , then m ( r ) /lessorequalslant M and the position of the event horizon r + is bounded by the Schwarzschild radius,  \nr + /lessorequalslant r s = 2 GM . (5.12)  \nIn general, if ρ e ff ( r ) is positive, the inequality (5.12) is always true. This statement can be proved by observing that if ρ e ff ( r ) is positive, then the mass function m ( r ) is monotone. In fact,  \nρ e ff ( r ) /greaterorequalslant 0 = ⇒ m ' ( r ) = 4 π r 2 ρ e ff ( r ) /greaterorequalslant 0 , (5.13)  \nwhence m ( r ) is a monotonically increasing function. Then, from Eq. (5.5) it follows m ( r ) /lessorequalslant M . On the other hand, if there is a region where ρ e ff ( r ) < 0, then m ( r ) must have oscillations, and it may assume values larger than M (see [7, 91, 92, 114] for explicit examples).  \nThe value of the critical mass M c is related to the massive parameters of the model (see Sec. 7 for explicit examples). If M /greatermuch M c and those parameters are such that the smearing of the source is negligible at the scale 2 GM , then we expect the solution to be very close to the Schwarzschild solution at the outer horizon. In particular, in the context of Eq. (5.12), r + approaches r s from below.  \nIt is important to keep in mind that the linearity of the function A ( r ) in the mass parameter M plays a central role when deriving the above conclusions regarding the number of horizons, oscillations of the e ff ective mass function, and the existence of a black hole mass gap. This underlying assumption can be traced back to Eq. (2.4) and the form of the e ff ective source (1.9). An intriguing question is whether such features of the e ff ective solutions are also present in the solutions of the complete field equations of the original higher-derivative / nonlocal gravity. The results on exact solutions of sixth- and higher-derivative gravities so far available in the literature are not enough to clarify the situation, for they mainly concern either the weak-field regime outside a horizon or local aspects of the solutions [86-89]. Nevertheless, we remark that even for general functions A ( r ), asymptotically flat, everywhere regular geometries must have an even number of horizons [86].  \nTo close this section on general features of the function A ( r ), we present some of the properties of its derivatives, which are related to the occurrence of singularities in the curvature invariants at r = 0, as discussed in Sec. 4.  \nApplying Taylor's theorem and taking into account the results (3a) and (4a) on the behavior of the e ff ective delta source near r = 0, it follows that: If the function f ( -k 2 ) asymptotically grows at least as fast as k 4 + 2 N for an integer N /greaterorequalslant 0, then A ( r ) is at least ( N + 1)-regular, i.e., it is 2 N + 2 times di ff erentiable and  \nd 2 n + 1 d r 2 n + 1 A ( r ) ∣ ∣ ∣ ∣ r = 0 = 0 for n = 0 , . . . , N . (5.14)  \nNotice that this result also covers the cases of the nonlocal form factors by Kuz'min and Tomboulis [5, 6], which behave like polynomials in the UV. In particular, if f ( -k 2 ) asymptotically grows at least as fast as k 4 (like in any polynomial model), we get A ' (0) = 0, which is a necessary condition for the absence of curvature singularities at r = 0 [see the discussion involving Eq. 4.7].  \nFinally, as a corollary, if f ( -k 2 ) asymptotically grows faster than any polynomial, A ( r ) is an even function-which also implies that the e ff ective mass m ( r ) is an odd function. This can be proved, alternatively, from the result (5a) on the e ff ective delta source:  \nρ e ff ( -r ) = ρ e ff ( r ) = ⇒ m ( -r ) = -4 π ∫ r 0 d x x 2 ρ e ff ( -x ) = -m ( r ) , (5.15)  \n↦  \nwhere we made the change of variable x →-x in the integration; thus, A ( -r ) = A ( r ). As recently discussed in [56, 83], if A ( r ) and B ( r ) are even analytic functions, not only the curvature invariants are finite at r = 0, but also all the curvature-derivative invariants. The analogous result for B ( r ) is going to be derived in the next section.  \n↦", '6. GENERAL PROPERTIES OF B ( r )': 'In this section, we analyze several choices for fixing the e ff ective pressures and underlying motivations. For each, we present the general properties of the function B ( r ) associated with it.', '6.1. Case of null pressures': "Considering  \nthe solution (2.9) for B ( r ) becomes  \nB ( r ) = B 0 + 8 π G ∫ r r 0 d x x ρ e ff ( x ) A ( x ) . (6.2)  \nThe properties of the function B ( r ) defined in (6.2) can be directly deduced from the results concerning the e ff ective source (mentioned in Sec. 5 and proved in [56, 57]) and the results obtained for A ( r ) in the previous section. Namely, if A ( r ) > 0 ∀ r then B ( r ) is bounded for any higher-derivative model. Moreover, if f ( -k 2 ) ∼ k 2 for large k , then B ( r ) is di ff erentiable but it is not 1-regular, as B ' (0) /nequal 0; while if f ( -k 2 ) grows at least as fast as k 4 + 2 N for a certain N /greaterorequalslant 0 and su ffi ciently large k , B ( r ) is at least ( N + 1)-regular, hence,  \nd 2 n + 1 d r 2 n + 1 B ( r ) ∣ ∣ ∣ ∣ r = 0 = 0 for n = 0 , · · · , N . (6.3)  \nThe solution (2.7) together with (6.2) corresponds to the case where Eq. (1.10) serves as a 'small-curvature approximation' of the full nonlinear equations of motion for the theory (1.3) sourced by a pointlike mass,  \nIn particular, for any model with six or more derivatives, i.e., if f ( -k 2 ) ∼ k 4 or faster for large k , B ' (0) = 0 . Together with (5.14), this last result is responsible for canceling the curvature singularity at r = 0 for models with f ( -k 2 ) growing at least as fast as k 4 [see the discussion involving Eq. 4.7]. Moreover, if f ( -k 2 ) grows faster than any polynomial, there exists a neighborhood of r = 0 where B ( r ) is an even function. This can be proved by setting r 0 = 0 in (6.2), making the change of integration variable x →-x , and using that for this type of form factor ρ e ff ( r ) and A ( r ) are also even functions.  \nT µν ( /vector r ) = M δ t µ δ t ν δ (3) ( /vector r ) . (6.4)  \nIndeed, if the original source is a pointlike particle, the inversion of the operator f ( z ) in (1.12), ignoring the terms O ( R 2 ··· ) of higher order in curvatures, automatically gives the e ff ective energy-momentum tensor (1.13) with pr = p θ = 0. However, a careful reader will notice that in this situation the e ff ective energy-momentum tensor ˜ T µ ν is not conserved [or, equivalently, the equation G θ θ = 0 is not satisfied, see (2.6)]. In other words, this choice of e ff ective pressures makes the system (2.4)-(2.6) over-determined. Approximate solutions can occur, nevertheless: The conservation equation (2.10) might hold only under the approximation , namely, in regions where the spacetime curvature is small-that is, in the whole spacetime (if M < Mc ) or in the limits of short and large distances (if M > Mc ), as we discuss in what follows.  \nFirstly, for M < Mc , the situation is similar to what happens in the Newtonian approximation (see the discussion in Sec. 3) since the conservation equation (2.10) for pr = p θ = 0 reduces to  \n-1 2 ( A ' A + B ' ) ρ e ff = O ( M 2 ) . (6.5)  \nSimilarly, using (2.7) and (6.2) in (2.3), it is possible to prove (after some manipulations) that  \nG θ θ = O ( M 2 ) , (6.6)  \nwhich means that the field equations are satisfied in a good approximation in the whole spacetime if the mass M is su ffi ciently small. How small the mass needs to be depends on the critical mass Mc , because for M < Mc the functions A ( r ) and B ( r ) are bounded, and so are the curvatures. Indeed, in this case, the spacetime curvature is limited by a certain value and might not grow too much. This will become more evident in the specific examples of Sec. 7.  \nOn the other hand, for M > Mc the solution (6.2) is ill-defined at the position of the horizons, given the occurrence of the term 1 / A ( r ) in the integrand. The dominance of the singular structure 1 / A ( r ) in (6.2) characterizes the high-curvature regions of the spacetime, where the 'small-curvature approximation' breaks down. Hence, the expression (2.9) can be used as long as the  \npr = p θ = 0 , (6.1)  \ninterval [ r 0 , r ] does not cross any horizon. Assuming that A ( r ) and B ( r ) are analytic functions, 11 far from the horizons and near r = 0 we can use their power series representation (4.3), which substituted into (2.3) give  \nG θ θ = 1 -a 0 r 2 -a 0 b 1 + 2 a 1 2 r + 1 4 ( a 0 b 2 1 + a 1 b 1 ) + O ( r ) . (6.7)  \nBut, from Eqs. (5.7), (5.14) and (6.3) we have a 0 = 1 and a 1 = b 1 = 0, such that  \nG θ θ = O ( r ) . (6.8)  \nThis means that the e ff ective field equations are satisfied in a good approximation for a su ffi ciently small r around the point r = 0 for any value of the mass M . Note that this is true only for the e ff ective sources of models with more than four derivatives of the metric in the action. In fourth-derivative gravity, we have a 1 , b 1 /nequal 0; therefore, the e ff ective equations (1.10) will not reproduce the properties described in [32, 33] of Stelle gravity for small r . Finally, the angular equations are also approximately satisfied for large r because lim r →∞ ρ e ff = 0 implies lim r →∞ G θ θ = 0 through the functions A ( r ) and B ( r ).  \n6.2.  \nCase  \npr  \n=  \nαρ  \ne  \nff  \nIf one sets pr = αρ e ff (6.9)  \nfor a certain constant α , but leaves the tangential pressure to be determined by the conservation equation (2.10), the result for B ( r ) is  \nB ( r ) = B 0 + (1 + α )8 π G ∫ r r 0 d x x ρ e ff ( x ) A ( x ) , (6.10)  \nwhile p θ is given by  \np θ = 1 + α 4 ( A ' A + B ' ) r ρ e ff + α 2 ( 2 ρ e ff + r ρ ' e ff ) . (6.11)  \nTherefore, the functions A ( r ) and B ( r ) can be obtained from the equations (2.4) and (2.5) alone, and the quantity p θ is constructed afterward. (It is straightforward to verify that, with p θ defined above, the remaining equation (2.6) is also satisfied.) Note that the general properties of (6.10) are the same as described in the previous Sec. 6.1 because this solution only di ff ers with respect to (6.2) by a multiplicative factor in front of the integral.  \nA drawback of this choice of equation of state is the occurrence of the term 1 / A ( r ) in the expression of the tangential pressure (6.11), for p θ ( r ) is likely to diverge if the metric contains horizons. For example, with the simplest choice 12 α = 0, using the field equations one gets  \nA ' ( r ) A ( r ) + B ' ( r ) = 2 Gm ( r ) r 2 A ( r ) . (6.12)  \nThus, if the metric has a horizon at r = rH , it follows 2 Gm ( rH ) = rH > 0, showing that the term in (6.12) certainly diverges at the horizon. The only possibility for the tangential pressure [see (6.11)]  \np θ ( r ) = Gm ( r ) ρ e ff ( r ) 2 rA ( r ) (6.13)  \nto be regular is if it happens that ρ e ff ( r ) vanishes at rH . In general, the singularity of the tangential pressure at the horizons is related to the singularity of B ' ( r ) at r = rH , see Eq. (6.10). These singularities are physical and cannot be removed by changing coordinates, for they also manifest in the curvature invariants.  \nOne possible way to avoid the obstacles caused by the divergences at the horizons is to choose α = -1. Indeed, for this particular value of the parameter, the e ff ective tangential pressure [see (6.11)],  \np θ = -ρ e ff -r ρ ' e ff 2 , (6.14)  \nis regular for any higher-derivative gravity model such that f ( -k 2 ) grows at least as fast as k 4 for su ffi ciently large k .  \nHowever, the right-hand side of the equation (2.5) vanishes with this choice, so that B = 0 [see also Eq. (6.10)]. Therefore, the metric (1.14) assumes the Schwarzschild-like form  \nd s 2 = -A ( r )d t 2 + d r 2 A ( r ) + r 2 d Ω 2 . (6.15)  \nAs mentioned in the Introduction, the disadvantage of (6.15) is that, since B = 0, these types of solutions cannot reproduce the asymptotic Newtonian-limit behavior of the higher-derivative gravities.  \nAnother starting point for obtaining this solution can be imposing (6.15) as Ansatz for the metric (which makes the system (2.4)-(2.6) consistent) and then finding  \npr = -ρ e ff (6.16)  \nas the solution of (2.5), and (6.14) as solution of (2.6). This approach was adopted in the frameworks of noncommutative gravity, generalized uncertainty principle, and other UV completions of gravity, e.g., in Refs. [90, 109-111], as well as in some local and nonlocal higher-derivative models [7, 91, 92].", '6.3. General case of regular pressures': "As it was evident in the preceding discussion, the main di ffi culty in defining B ( r ) is the presence of A ( r ) in Eq. (2.5) if the metric has horizons. We also saw that the choice of radial pressure pr = -ρ e ff was capable of curing pressure divergences at the horizons, with the side e ff ect that B ≡ 0. Moreover, it was argued that the curvature singularity at the horizons was related to a divergence in the e ff ective pressures. So, another possibility of regularizing these horizon physical singularities, but with a nontrivial B ( r ), is with an equation of state in the form  \npr ( r ) = A ( r ) σ ( r ) -ρ e ff ( r ) , (6.17)  \nwhere σ ( x ) is an Nb -regular continuous function (for an integer Nb /greaterorequalslant 0) such that the function ξ ( x ) = x σ ( x ) is integrable on [0 , + ∞ ).  \nIndeed, this choice eliminates the term 1 / A ( r ) in Eq. (2.9), so that the solution for B ( r ) becomes  \nB ( r ) = 8 π G ∫ r ∞ d x x σ ( x ) . (6.18)  \nBy definition, this integral is well defined for all values of r , and for convenience, we set r 0 = ∞ . With the assumptions about σ ( x ), it follows that B ( r ) is ( Nb + 1)-regular.  \nTo verify that the equation of state (6.17) makes the e ff ective pressures regular, first we notice that if ρ e ff ( r ) is Na -regular and σ ( x ) is Nb -regular, then A ( r ) is ( Na + 1)-regular (see Sec. 5) and pr ( r ) is at least N -regular, with N ≡ min { Na , Nb } . In what concerns the tangential pressure defined by the conservation equation (2.10),  \np θ = -ρ e ff + A σ + r 4 ( 3 A ' σ + AB ' σ + 2 A σ ' -2 ρ ' e ff ) , (6.19)  \nit is straightforward to verify that the order of regularity of each of the terms in the right-hand side of (6.19) is either Na or min { Na + 1 , Nb } ; therefore, p θ ( r ) is at least N -regular, just like pr ( r ).  \nOf course, in principle, there is a huge freedom to choose regular functions σ ( r ) capable of yielding regular black hole solutions. For example, Ref. [108] in the context of noncommutative geometry used σ ( r ) ∝ r ρ e ff ( r ). However, since for e ff ective sources coming from higher-derivative gravity this function is only 0-regular, B ( r ) cannot be more than 1-regular-regardless of the order of regularity of ρ e ff . Further discussion about this choice of σ ( r ) is carried out in Appendix A. On the other hand, the multiplication for higher powers of r can enhance the regularity of B ( r ); for instance, if ρ e ff ( r ) is Na -regular, the choice σ ( r ) ∝ r 2 n ρ e ff ( r ) (for a positive integer n ) makes B ( r ) to be ( Na + n + 1)-regular. We shall consider another form for the function σ ( r ) in what follows.  \n6.4. Case pr = [ A ( r ) -1] ρ e ff  \nAs a particular case of the equation of state presented in the previous subsection, let us consider a function σ ( r ) such that the resultant pressure components have the following physical properties:  \n- (1b) The e ff ective pressures vanish asymptotically for large r .\n- (2b) The e ff ective pressures are finite at the horizons.\n- (3b) The e ff ective pressures vanish at r = 0.  \nWhile (1b) and (2b) are natural in the framework of Sec. 6.3, (3b) actually imposes restrictions on the form of σ ( r ). 13 The physical motivation for these requirements is to guarantee that the solution is asymptotically flat, everywhere regular and that it matches the smallr behavior of the solution obtained in Sec. 6.1 (which approximates the solution whose original source is a pointlike mass).  \nIn the most interesting case (for this work) of a regular ρ e ff , a simple equation of state satisfying the conditions (1b)-(3b) is  \npr ( r ) = [ A ( r ) -1] ρ e ff ( r ) , (6.20)  \nwhich follows from the choice σ ( r ) = ρ e ff ( r ) in the context of Sec. 6.3. Accordingly, the solution (6.18) for B ( r ) reads  \nB ( r ) = 8 π G ∫ r ∞ d x x ρ e ff ( x ) , (6.21)  \nwhile the tangential pressure (6.19) is  \np θ = [ 3 4 A ' r + A ( 1 + 2 G π r 2 ρ e ff ) -1 ] ρ e ff + A -1 2 r ρ ' e ff . (6.22)  \nIt is straightforward to verify that (6.20) and (6.22) indeed vanish as r → 0 if ρ e ff is regular. 14 Also, from the general results of Sec. 6.3, we know that if ρ e ff is N -regular, pr and p θ are bounded and (at least) N -regular. A finer result is obtained in the case ρ e ff ( r ) is analytic 15 at r = 0-in that case pr and p θ are actually ( N + 1)-regular and, as a consequence, near the origin they are O ( r 2 ) even for N = 0. This statement can be proved by direct calculation, using the Taylor series expansion of Eqs. (6.20) and (6.22) [with (2.7) and (2.8)] and noticing the cancellation of the term of order r 2 N + 1 .  \nFrom the features of ρ e ff mentioned in Sec. 5, it follows that the function B ( r ) in (6.21) is bounded for any higher-derivative model. Moreover, from the considerations of Sec. 6.3, if ρ e ff is N -regular (or, equivalently, if f ( -k 2 ) grows at least as fast as k 4 + 2 N for a certain N /greaterorequalslant 0 and su ffi ciently large k ), B ( r ) is ( N + 1)-regular. In particular, the identity (6.3) also holds, showing that the behavior of B ( r ) in (6.21) near r = 0 is very similar to the one of (6.2), as expected.  \nOther choices of equations of state can also reproduce the same aforementioned physical e ff ect on the solutions. Indeed, this procedure has considerable ambiguity; for example, any equation of state in the form  \npr ( r ) = [ A /lscript ( r ) -1 ] ρ e ff ( r ) , /lscript = 1 , 2 , . . . , (6.23)  \nfollowing from σ ( r ) = A /lscript -1 ( r ) ρ e ff ( r ), would also lead to the above properties (1b)-(3b). In this work, for the sake of simplicity, we adopt (6.20).  \n∗ ∗ ∗  \nAs a conclusion of this section, here we presented a recipe, based on Eq. (1.10), of how to construct black hole solutions that are regular everywhere. First, choose a form factor that behaves as f ( -k 2 ) ∼ k 4 (or faster) for large k ; then, choose e ff ective pressures satisfying the conditions (1b) and (2b)-for example, through the function σ ( r ) defined in Sec. 6.3.  \nIn addition, if pr and p θ vanish su ffi ciently fast for large r and at the origin [condition (3b)], those solutions reproduce the asymptotic behavior of the corresponding linearized higher-derivative model in the limits r → ∞ and r → 0. At the same time, near the horizons the metric is regularized by the regular pressures. One might even argue that these nontrivial pressure components can mimic the e ff ect of the omitted higher-order terms when passing from Eq. (1.12) to (1.10). In our opinion, however, a proper address of this issue ultimately involves the study of the complete equations of motion, which exceeds the scope of this work.  \nFinally, for the equation of state (6.20) the metric is given by  \nd s 2 = -( 1 -2 Gm ( r ) r ) exp ( 8 π G ∫ r ∞ d x x ρ e ff ( x ) ) d t 2 + ( 1 -2 Gm ( r ) r ) -1 d r 2 + r 2 d Ω 2 , (6.24)  \nwith m ( r ) defined in (2.8), as usual.", '7. SELECTED EXAMPLES': 'In the previous sections we presented general results, which depended only on the asymptotic behavior of the form factor F ( z ) of the model (1.3). This section aims to provide concrete examples where the results above can be explicitly verified and other particular features of the solutions discussed. In Sec. 7.1 we present an example of nonlocal ghost-free gravity, while in Secs. 7.2 and 7.3 we discuss the cases of local higher-derivative gravity with simple poles. The general case of local higher-derivative gravity models, including degenerate poles, is considered in Sec. 7.4. A detailed analysis of the sixth-derivative Lee-Wick gravity is carried out in a separate publication [114].', '7.1. Nonlocal ghost-free gravity': 'One of the simplest nonlocal higher-derivative models is defined by the form factor [7, 69-72]  \nF ( z ) = e -z /µ 2 -1 z , (7.1)  \nfor which  \nf ( z ) = e -z /µ 2 , (7.2)  \nwhere µ > 0 has dimension of mass and defines the nonlocality scale. In this case, performing the integration (1.9), the smeared e ff ective delta source has a Gaussian profile [69],  \nρ e ff ( r ) = M µ 3 8 π 3 / 2 e -1 4 µ 2 r 2 , (7.3)  \nand the mass function is given by  \nm ( r ) = M [ erf ( µ r 2 ) -µ r √ π e -1 4 µ 2 r 2 ] . (7.4)  \nIn Fig. 1 we show the behavior of the e ff ective source and mass function, where the properties described in Sec. 5 can be explicitly verified.  \nGiven (7.4), the function A ( r ) in (2.7) reads  \nA ( r ) = 1 -2 GM r erf ( µ r 2 ) + 2 GM µ √ π e -1 4 µ 2 r 2 . (7.5)  \nIn order to evaluate the black hole mass gap, it is convenient to rewrite Eq. (7.5) in terms of the dimensionless function Z ( x ), namely,  \nA ( r ) = 1 -2 GM µ Z ( µ r ) , (7.6)  \nFIG. 1: Plot of the dimensionless forms of the e ff ective delta source (7.3) (left panel) and the mass function (7.4) (right panel) in terms of the variable x = µ r . It is immediate to verify that the graphs are in agreement with Eqs. (5.1) and (5.3) for the e ff ective source, and (5.5) and (5.6) for the mass function. Also, because ρ e ff ( r ) takes only positive values, the mass function is monotonic and m ( r ) /lessorequalslant M ; see the discussion around (5.13).  \n<!-- image -->  \nwhere  \nZ ( x ) = 1 x erf ( x 2 ) -e -x 2 4 √ π . (7.7)  \nThe plot of (7.7) is shown in Fig. 2, where one can verify that Z ( x ) is bounded, as expected from the general discussion in Sec. 5. Indeed, lim x → 0 Z ( x ) = lim x →∞ Z ( x ) = 0, and it has an absolute maximum Z max = 0 . 263 at the point x ∗ = 3 . 02. Hence, the critical mass (5.11) in this case is given by  \nM c = 1 2 G µ Z max = 1 . 9 M 2 P µ , (7.8)  \nwhere we used G = 1 / M 2 P .  \nAccordingly, if M < M c, A ( r ) is always positive and the metric has no horizon. On the other hand, for M > M c, there exists a region where A ( r ) takes negative values. Since the function Z ( x ) in (7.7) only has one local maximum, the equation A ( r ) = 0 can have at most two roots, r -and r + , with r + /greaterorequalslant r -. The values r ± are the positions of the event horizon and an inner horizon. In Fig. 3 we show the plot of A ( r ) for the two possible scenarios. Moreover, because m ( r ) /lessorequalslant M , the position of the event horizon is bounded by the Schwarzschild radius, i.e., r + /lessorequalslant r s = 2 GM .  \nFIG. 2: Plot of Z ( x ) in Eq. (7.7).  \n<!-- image -->  \nFIG. 3: Plot of A ( r ) for µ = M P = 1 with M = 2 . 5 > M c (solid line, two horizons), M = 1 . 9 = M c (dashed line, one extremal horizon) and M = 1 . 4 < M c (dotted line, no horizon).  \n<!-- image -->', '7.1.1. Solution for B ( r ) with pr = p θ = 0': "Firstly, let us analyze the framework of Sec. 6.1, i.e., with the e ff ective pressures pr = p θ = 0, corresponding to the case in which (1.10) is the small-curvature approximation for the full equations of motion of the theory (1.3), sourced by a pointlike mass. The solution for B ( r ) is obtained by substituting (7.3) and (7.5) in the integral (6.2), which can be evaluated numerically (if the integration interval does not cross a horizon). However, in order to study the curvature tensor components (4.1) it su ffi ces to consider the quantity  \nB ' ( r ) = 8 π G r ρ e ff ( r ) A ( r ) , (7.9)  \nbecause the curvature components (4.1) depend only on the derivatives of B ( r ).  \nBy using (7.3) and (7.5) into (7.9) we can evaluate each component in (4.1). For simplicity, in what follows we focus only on the scalar curvature,  \nR = 2 K 1 + 4 K 2 + 4 K 3 + 2 K 4 , (7.10)  \nbut similar results also hold for the other invariants, such as the Kretschmann scalar (4.2). The explicit expression for (7.10) is  \nR ( r ) = S ( r ) A ( r ) , (7.11)  \nFIG. 4: Comparison between R (solid line) and R lin (dashed line), respectively given by Eqs. (7.11) and (7.15). The panel (a) shows the case in which the mass is considerably smaller than the critical mass: R is positive and achieves its maximum at r = 0; the two solutions have the same qualitative behavior. For larger values of the mass (but still smaller than M c), R may achieve arbitrarily large negative values [panel (b)]. The maximum of | R | can occur at r /nequal 0, which shows a significant di ff erence between the two solutions; this marks the approximation breakdown. The panel (c) illustrates the typical case when M > M c, for which R diverges at the horizons.  \n<!-- image -->  \n≈  \nwhere  \nS ( r ) = GM µ 3 √ π [ 1 -3 GM r erf ( µ r 2 ) ] e -1 4 µ 2 r 2 + 3 G 2 M 2 µ 4 π e -1 2 µ 2 r 2 . (7.12)  \nIt is easy to verify that lim r →∞ R ( r ) = 0 and  \nlim r → 0 R ( r ) = GM µ 3 √ π . (7.13)  \nAt intermediate scales, although S ( r ) is a continuous, analytic even function of r , owning to the term 1 / A ( r ) in (7.11), the behavior of R ( r ) depends on the value of M compared to the critical mass M c.  \nIn Fig. 4 we show the typical behavior of (7.11) for M < M c and M > M c. We also compare both graphs with the Newtonianlimit solution (1.5), with the potential associated to the form factor (7.2) [7, 69-72],  \nϕ ( r ) = -GM r erf ( µ r 2 ) , (7.14)  \nfor which the linearized Ricci scalar reads  \nR lin = G µ 3 M √ π e -1 4 µ 2 r 2 . (7.15)  \nAs one can see from the graphs, in both cases, the linearized and nonlinearized solutions match pretty well for very small and large r , i.e., for r /lessmuch 1 /µ and r /greatermuch 1 /µ . Also, if M is su ffi ciently small, the two solutions do not deviate too much; we elaborate more on this in what follows.  \nAs discussed in Sec. 6, if M < M c we have A ( r ) > 0 everywhere and the scalar curvature (7.11) is bounded. To be more precise, if M < 0 . 66 M c the function S ( r ) is positive and R achieves the maximum value at r = 0, see Eq. (7.13) and Fig. 4a. For larger values of M within the range 0 . 66 M c < M < M c, the negative term in (7.12) becomes more relevant. It dominates for 0 . 95 M c < M < M c (Fig. 4b)-in this range, | R ( r ) | can assume arbitrarily large values as M approaches the critical mass. In fact, for M > M c the equation A ( r ) = 0 has two roots, and the scalar curvature diverges at the horizon positions r ± (see Fig. 4c). Therefore, combining these results with the numeric factors in Eq. (7.13) one can show that, even though R is bounded for any value of M smaller than M c, the small-curvature condition stated in the form | R | /µ 2 < 1 actually requires a slightly more stringent constraint, namely, M < 0 . 93 M c.  \nAnalogous results can be formulated for the other curvature-squared invariants, R 2 µν and R 2 µναβ . This explicitly shows that there exists a value for the mass M < M c such that the small-curvature condition R 2 /µ 4 < 1 holds ( R is a generic curvature tensor). The O ( R 2 ) corrections in (1.12) might be irrelevant in these circumstances. We emphasize, however, that the approximation holds for su ffi ciently small r even when M > M c because the nonlinearities in curvatures are suppressed in this limit, as discussed in general grounds in Sec. 4.  \nAs explained before, since the form factor (7.2) grows faster than any polynomial, the metric components A ( r ) and B ( r ) are even functions. This makes all curvature invariants to be singularity free at the origin, even those constructed with derivatives of the curvature tensors. Below, we present some values of these invariants at r = 0 [see also (7.13)]:  \nR µναβ R µναβ (0) = 5 G 2 M 2 µ 6 3 π , R µν R µν (0) = G 2 M 2 µ 6 π , (7.16)  \n/square R (0) = -GM µ 5 2 π ( 3 √ π + 2 G µ M ) , /square 2 R (0) = GM µ 7 12 π 3 / 2 ( 45 π + 138 √ π G µ M -52 G 2 µ 2 M 2 ) . (7.17)  \nIt is remarkable that with the very simple equation (1.10), one can obtain the same qualitative results of [86]. Namely, the curvatures are regular at r = 0 and weak around the origin, there exists a mass gap for the black hole formation, and the curvatures abruptly grow close to the horizon positions (if there are horizons).  \n7.1.2. Solution for B ( r ) with pr = [ A ( r ) -1] ρ e ff  \nLet us now obtain the solution associated with the equation of state (6.20), of Sec. 6.4, which regularizes the divergences of the pressures at the horizons. The quadrature (6.21) for B ( r ) can be explicitly evaluated in this case,  \nB ( r ) = -2 G µ M √ π e -1 4 µ 2 r 2 . (7.18)  \nSo, the metric is given by  \nd s 2 = -[ 1 -2 GM r erf ( µ r 2 ) + 2 GM µ √ π e -1 4 µ 2 r 2 ] exp ( -2 G µ M √ π e -1 4 µ 2 r 2 ) d t 2 + [ 1 -2 GM r erf ( µ r 2 ) + 2 GM µ √ π e -1 4 µ 2 r 2 ] -1 d r 2 + r 2 d Ω 2 , (7.19)  \nwhich is consistent with the Newtonian approximation described in Sec. 3. Indeed, expanding the temporal part of (7.19) in powers of M , at leading order we get gtt = -[1 + 2 ϕ ( r )], with ϕ ( r ) given by (7.14).  \nIn consonance with the discussion of Sec. 4, in this case the curvature invariants are regular everywhere. Indeed, in Fig. 5 we show the typical behavior of the scalar curvature R for M > M c; comparing with Fig. 4c, we notice that this choice of e ff ective pressures regularizes the singularities at the horizons. Moreover, the values of the curvature invariants at r = 0 are the same as in Eqs. (7.13) and (7.16), because both equations of state lead to solutions that have the same behavior for su ffi ciently small r , as explained in Sec. 6.4. (However, curvature-derivative invariants calculated at the two solutions might di ff er, as they involve higher-order metric derivatives.) We do not plot the case where M < M c since it is very similar to the one shown in Fig. 4a.  \nFinally, we point out that the metric (7.19) di ff ers from the one obtained in Ref. [108] inspired by noncommutative geometry. Although both come from the same Gaussian e ff ective source (7.3), the underlying equations of state are di ff erent, which leads to di ff erent shift functions B ( r ). In Appendix A, we carry out a more detailed comparison of the two solutions.  \nFIG. 5: Behavior of the Ricci scalar for the metric (7.19) when M > Mc . The scalar curvature is finite everywhere, even at the horizons r ± .  \n<!-- image -->", '7.2. Sixth-derivative gravity with real poles': "As a second example, let us consider the case where the form factor f ( z ) is a polynomial of degree two with simple real roots m 2 1 and m 2 2 , namely,  \nf ( z ) = ( m 2 1 -z )( m 2 2 -z ) m 2 1 m 2 2 , (7.20)  \nFIG. 6: Graph of Z ( x ) in Eq. (7.25). The function is positive for every x > 0 and it has a maximum Z max = 0 . 298 at x ∗ = 1 . 79.  \n<!-- image -->  \nwhich corresponds to a sixth-derivative gravitational action. The e ff ective delta source (1.9) for this model was calculated in [57] and it reads  \nρ e ff ( r ) = Mm 2 1 m 2 2 4 π ( m 2 2 -m 2 1 ) e -m 1 r -e -m 2 r r . (7.21)  \nOrdering the masses such that m 1 < m 2, it is easy to see that ρ e ff ( r ) > 0 for all r . According to the discussion involving Eq. (5.13), it follows that the mass function,  \nm ( r ) = M       1 + m 2 1 (1 + m 2 r ) e -m 2 r -m 2 2 (1 + m 1 r ) e -m 1 r m 2 2 -m 2 1       , (7.22)  \nis monotone and also positive; this can be explicitly verified by checking that m ' ( r ) > 0 and m (0) = 0.  \nIt is worth mentioning that the limit m 2 → m 1 of the expressions (7.21) and (7.22) is smooth, and it results in the formulas for the case of sixth-derivative gravity with degenerate masses, namely,  \nρ e ff ( r ) ∣ ∣ ∣ m 2 = m 1 = Mm 3 1 8 π e -m 1 r , m ( r ) ∣ ∣ ∣ m 2 = m 1 = M [ 1 -( 1 + m 1 r + 1 2 m 2 1 r 2 ) e -m 1 r ] . (7.23)  \nThe expression for A ( r ) is obtained from combining (2.7) and (7.22), and it can be cast in the form 16  \nSince the UV behavior of the form factor does not depend on the multiplicity of the roots, the e ff ective sources in (7.23) and (7.21) have the same order of regularity, as it can be easily checked by comparing their Taylor series. For further discussion regarding the limit m 2 → m 1 in linearized sixth-derivative gravity, see [67].  \nA ( r ) = 1 -2 GM m 1 m 2 m 2 2 -m 2 1 [ m 2 Z ( m 1 r ) -m 1 Z ( m 2 r )] , (7.24)  \nwhere we defined the dimensionless function  \nZ ( x ) = 1 x [1 -(1 + x ) e -x ] . (7.25)  \nSince Z ( x ) is bounded (see Fig. 6), A ( r ) is also bounded. Moreover, the function Z ( x ) is positive for any x > 0 and takes the maximal value Z max = 0 . 298 for x ∗ = 1 . 79. Therefore,  \nmax r { m 2 Z ( m 1 r ) -m 1 Z ( m 2 r ) } /lessorequalslant max r { m 2 | Z ( m 1 r ) | + m 1 | Z ( m 2 r ) |} /lessorequalslant ( m 1 + m 2) Z max , (7.26)  \nA ( r ) ∣ ∣ ∣ m 2 = m 1 = 1 -2 GM r [ 1 -( 1 + m 1 r + 1 2 m 2 1 r 2 ) e -m 1 r ] .  \nFIG. 7: Numerical solution of Eq. (7.30), which shows that it has only one solution for each q /greaterorequalslant 1. In the limit q → 1 the equation is solved for y = 3 . 38, while for large q the solution for y tends to 1.79.  \n<!-- image -->  \nwhich leads to the conclusion that A ( r ) is positive if M is small enough such that  \n2 GM m 1 m 2 m 2 -m 1 Z max < 1 . (7.27)  \nThis yields a lower-bound estimate to the critical mass M c of Eq. (5.11),  \nM c /greaterorequalslant 1 . 68 M 2 P ˜ µ , ˜ µ ≡ m 1 m 2 m 2 -m 1 , (7.28)  \nwhere again we used G = 1 / M 2 P . We remark that M c might actually be considerably bigger than the estimate (7.28), which is only a lower bound. Nevertheless, if M < 1 . 68 M 2 P / ˜ µ the solution describes a horizonless compact object regardless of the choice of the equation of state in the framework of Sec. 6.3.  \nOn the other hand, for M > M c, the solution has exactly two horizons. To prove this statement, it su ffi ces to show that the function A ( r ) in Eq. (7.24) has only one local minimum. In particular, this is true if the equation A ' ( r ) = 0 has only one solution 17 for a given pair of masses m 1 and m 2. If we define the dimensionless variables  \nq ≡ m 2 m 1 and y ≡ m 1 r , (7.29)  \nthe equation A ' ( r ) = 0 can be cast as  \n1 -q 2 + q 2 ( 1 + y + y 2 ) e -y -[ 1 + qy (1 + qy ) ] e -qy 1 -q 2 = 0 . (7.30)  \nThe numerical solution of this equation in the qy -plane is shown in Fig. 7, where one can see that it has only one solution for each value of q . In other words, given a pair ( m 1 , q )-or, equivalently, ( m 1 , m 2)-the function A ( r ) has only one minimum.  \nFigure 7 can also be regarded as the graph of a function y min( q ), which is related to the position r ∗ of the minimum 18 of A ( r ), namely,  \nr ∗ = y min( q ) m 1 . (7.31)  \nTherefore, in Fig. 7 we can identify two special distinct regimes, namely, q → 1 and q → ∞ . The former corresponds to the degenerate case m 1 = m 2, for which y min assumes its maximal value y min(1) = 3 . 38. As q increases, the function y min( q )  \nmonotonically decreases, and it tends to the constant value x ∗ = 1 . 79 [the same position of the extremum of Z ( x ), see Eq. (7.25)]. This latter regime corresponds to the case in which m 2 /greatermuch m 1, i.e., the massive mode m 2 is too heavy and has a very short range; consequently, for larger distances, the model behaves similarly to the fourth-derivative gravity with only the scale m 1. Both situations are analogous to what was described in [67] about the weak-field limit of sixth-derivative gravity. Since the outer (event) horizon r + is bounded from below by r ∗ , and by the Schwarzschild radius from above [see Eq. (5.12)], we conclude that  \n1 . 79 m 1 /lessorequalslant r + /lessorequalslant 2 GM . (7.32)  \nLet us stress that the case of sixth-derivative gravity in which the masses m 1 /nequal m 2 in (7.20) are complex conjugate is completely di ff erent from the one studied here, for the e ff ective source is not positive and the mass function has oscillations. For instance, this allows the metric to have more than two horizons, as we showed in detail in [114].  \nIn what concerns the nontrivial possibilities for B ( r ), let us only discuss the solution with regular pressures in the framework of Sec. 6.4, i.e., with the equation of state (6.20) (the case of null pressures is very similar to Sec. 7.1.1). In this case, Eq. (6.21) yields 19  \nB ( r ) = -2 GMm 1 m 2 m 2 2 -m 2 1 ( m 2 e -m 1 r -m 1 e -m 2 r ) . (7.33)  \nCombining (7.24) and (7.33) we obtain the explicit expression for the metric,  \nd s 2 = -     1 -2 GM m 1 m 2 m 2 2 -m 2 1 [ m 2 Z ( m 1 r ) -m 1 Z ( m 2 r )]      exp       -2 GMm 1 m 2 m 2 2 -m 2 1 ( m 2 e -m 1 r -m 1 e -m 2 r )       d t 2 +      1 -2 GM m 1 m 2 m 2 2 -m 2 1 [ m 2 Z ( m 1 r ) -m 1 Z ( m 2 r )]      -1 d r 2 + r 2 d Ω 2 . (7.34)  \nIt is easy to check that to the leading order in M , the above metric reduces to the Newtonian-limit form (1.5), where  \nϕ ( r ) = -GM r       1 -m 2 2 e -m 1 r -m 2 1 e -m 2 r m 2 2 -m 2 1       (7.35)  \nis precisely the modified Newtonian potential [53, 67], in agreement with the general discussion of Sec. 3.", '7.3. Polynomial gravity with simple real poles': "As a more general but less explicit example, let us consider the case where the form factor f ( z ) is a generic polynomial of degree N /greaterorequalslant 2 with real distinct roots, i.e.,  \nf ( z ) = 1 + zF ( z ) = N ∏ i = 1 m 2 i -z m 2 i , (7.36)  \nwith mi /nequal mj if i /nequal j . (The more general case in which the roots can be complex and / or degenerate is discussed in Sec. 7.4 below.) The expression for the e ff ective source can be obtained by expanding the integrand of (1.9) in partial fractions, using  \n1 f ( z ) = N ∑ i = 1 Ci m 2 i m 2 i -z , (7.37)  \nwhere the coe ffi cients Ci are given by  \nCi = N ∏ j = 1 j /nequal i m 2 j m 2 j -m 2 i (7.38)  \nB ( r ) ∣ ∣ ∣ m 2 = m 1 = -GMm 1(1 + m 1 r ) e -m 1 r .  \nand satisfy the identities [56, 115, 116]  \nN ∑ i = 1 Ci m 2 k i =        1 , for k = 0, 0 , for k = 1 , 2 , . . . , N -1 . (7.39)  \nSubstituting (7.37) into (1.9) it follows [57]  \nρ e ff ( r ) = M 4 π r N ∑ i = 1 Ci m 2 i e -mir . (7.40)  \nHence, the e ff ective mass function (2.8) reads  \nm ( r ) = M        1 -N ∑ i = 1 Ci (1 + mir ) e -mir        . (7.41)  \nWe did not encounter any e ff ective source (7.40) that assumes negative values; therefore, we conjecture that ρ e ff ( r ) above is strictly positive, which makes m ( r ) to be a monotonic function and the horizons to be bounded by the Schwarzschild radius, according to the discussion involving Eq. (5.13). Let us remark once more that the situation with complex quantities m 2 i is entirely di ff erent, and in that case ρ e ff ( r ) can achieve negative values, with important implications for the black hole solutions [114].  \nThe combination of (2.7) and (7.41) yields  \nA ( r ) = 1 -2 GM r        1 -N ∑ i = 1 Ci (1 + mir ) e -mir        . (7.42)  \nFrom the identity (7.39) with k = 0, it is easy to see that the property (5.7) is satisfied, in particular, lim r → 0 A ( r ) = 1. The equation (7.42) can be rewritten in terms of the same dimensionless function Z ( x ) in (7.25),  \nA ( r ) = 1 -2 GM N ∑ i = 1 Ci mi Z ( mir ) . (7.43)  \nSince Z ( x ) is bounded by 0 /lessorequalslant Z ( x ) /lessorequalslant Z max = 0 . 298, A ( r ) is also bounded. Therefore, the same reasoning we adopted to estimate the mass gap in Sec. 7.2 leads to the conclusion that A ( r ) is positive provided that M is small enough such that  \n2 GM ˜ µ Z max < 1 , (7.44)  \nwhere  \n˜ µ ≡ N ∑ i = 1 | Ci | mi . (7.45)  \nThis yields the lower-bound estimate to the critical mass M c in Eq. (5.11),  \nM c /greaterorequalslant 1 . 68 M 2 P ˜ µ . (7.46)  \nNotice that Eqs. (7.46) and (7.28) are essentially the same, the only di ff erence being the mass scale ˜ µ . This is expected, for the previous section's model is the particular case N = 2 of the model considered here.  \nFinally, the equation of state (6.20) in the framework of Sec. 6.4 results in the solution (6.21) for B ( r ), namely,  \nB ( r ) = -2 GM N ∑ i = 1 Ci mi e -mir . (7.47)  \nThe complete metric, thus, reads  \nd s 2 = -       1 -2 GM r [ 1 -N ∑ i = 1 Ci (1 + mir ) e -mir ]        exp        -2 GM N ∑ i = 1 Ci mi e -mir        d t 2 +        1 -2 GM r [ 1 -N ∑ i = 1 Ci (1 + mir ) e -mir ]        -1 d r 2 + r 2 d Ω 2 . (7.48)  \nIn agreement with the general discussion of Sec. 3, in the leading order in M this metric reduces to the Newtonian form (1.5) with  \nϕ ( r ) = -GM r        1 + N ∑ i = 1 Ci e -mir        (7.49)  \nbeing exactly the modified Newtonian potential evaluated in [53].", '7.4. General polynomial gravity': 'Finally, let us consider the most general case of polynomial gravity, where we allow the mass parameters mi to be complex and / or degenerate. The form factor reads  \nf ( z ) = n ∏ i = 1       m 2 i -z m 2 i       α i , (7.50)  \nwhere we assume that the polynomial f ( z ) has degree N /greaterorequalslant 2 and it has n /greaterorequalslant 1 distinct roots m 2 1 , m 2 2 , . . . , m 2 n , each of which has a multiplicity α 1 , α 2 , . . . , α n , such that  \nN = n ∑ i = 1 α i . (7.51)  \nSimilarly to the previous example, we can use the partial fraction decomposition to write  \n1 f ( z ) = n ∑ i = 1 α i ∑ j = 1 C i , j       m 2 i m 2 i -z       j , (7.52)  \nwhere the coe ffi cients can be determined by the residue method (see, e.g., [117]) and are given by  \nC i , j = ( -1) α i -j m 2 j i ( α i -j )! d α i -j d z α i -j ( m 2 i -z ) α i f ( z ) ∣ ∣ ∣ ∣ ∣ z = m 2 i . (7.53)  \n∣ These coe ffi cients satisfy identities that generalize the ones in (7.39), in particular,  \nn ∑ i = 1 α i ∑ j = 1 C i , j = 1 and n ∑ i = 1 C i , 1 m 2 i = 0 . (7.54)  \nCalculating the e ff ective delta source ρ e ff from (1.9) using (7.52), we find [57]  \nρ e ff ( r ) = M 4 π 3 / 2 n ∑ i = 1 α i ∑ j = 1 C i , j m 2 j i ( j -1)! ( r 2 mi ) j -3 2 K j -3 2 ( mir ) , (7.55)  \nwhere K ν is the modified Bessel function of the second kind [117]. As shown in [57], the latter identity in (7.54) is responsible for the 0-regularity of the above e ff ective source.  \nIt is also useful to write explicitly the expression of the modified Newtonian potential ϕ ( r ), obtained from the modified Poisson equation (1.6) and calculated in [55],  \nϕ ( r ) = -GM r         1 -2 √ π n ∑ i = 1 α i ∑ j = 1 A i , j ( j -1)! ( mir 2 ) j -1 2 K j -3 2 ( mir )         , (7.56)  \nwhere the coe ffi cients A i , j are defined such that 20  \n1 z f ( z ) = 1 z + n ∑ i = 1 α i ∑ j = 1 A i , j m 2 i       m 2 i m 2 i -z       j , (7.57)  \nnamely,  \nA i , j = ( -1) α i -j m 2( j -1) i ( α i -j )! d α i -j d z α i -j ( m 2 i -z ) α i z f ( z ) ∣ ∣ ∣ ∣ ∣ ∣ z = m 2 i . (7.58)  \nHence, comparing (7.52) and (7.57) one can show that  \nA i , j = α i ∑ k = j C i , k . (7.59)  \nAt this point, we can use relations (3.5), (3.11) and (7.56) to obtain compact expressions for the mass function and the function A ( r ),  \nm ( r ) = M         1 -4 √ π n ∑ i = 1 α i ∑ j = 1 A i , j ( j -1)! ( mir 2 ) j + 1 2 K j -5 2 ( mir )         , (7.60)  \nA ( r ) = 1 -2 GM r + 8 GM √ π r n ∑ i = 1 α i ∑ j = 1 A i , j ( j -1)! ( mir 2 ) j + 1 2 K j -5 2 ( mir ) . (7.61)  \nOne can easily check that the formulas (7.56), (7.60), and (7.61) reduce to the simple pole ones of Sec. 7.3 if we set α i = 1 ∀ i . For instance, in this case, (7.53) coincides with (7.38).  \nFinally let us consider the solution described in Sec. 6.4, with an e ff ective pressure defined by the equation of state (6.20). From Eq. (6.21), we obtain the nontrivial shift function  \nB ( r ) = -4 GM √ π n ∑ i = 1 α i ∑ j = 1 C i , j mi ( j -1)! ( mir 2 ) j -1 2 K j -1 2 ( mir ) . (7.62)  \nPutting together the expressions for A ( r ) and B ( r ) we can write the metric (6.15) in explicit form,  \nd s 2 = -[ 1 -2 GM r + 8 GM √ π r n ∑ i = 1 α i ∑ j = 1 A i , j ( j -1)! ( mir 2 ) j + 1 2 K j -5 2 ( mir ) ] exp [ -4 GM √ π n ∑ i = 1 α i ∑ j = 1 C i , jmi ( j -1)! ( mir 2 ) j -1 2 K j -1 2 ( mir ) ] d t 2 + [ 1 -2 GM r + 8 GM √ π r n ∑ i = 1 α i ∑ j = 1 A i , j ( j -1)! ( mir 2 ) j + 1 2 K j -5 2 ( mir ) ] -1 d r 2 + r 2 d Ω 2 . (7.63)', '8. SUMMARYAND CONCLUSIONS': 'In the present work, we studied static spherically symmetric solutions of the e ff ective field equations (1.10) sourced by an e ff ective energy-momentum tensor containing the smeared delta source. In particular, we discussed the role played by the e ff ective pressure components, several possible ways of defining them, and their consequences to the solution. Throughout the analysis, we kept the e ff ective delta source as general as possible, relating our results to certain generic properties of such sources; this guarantees the generality of the considerations and their applicability to a broad family of models defined by form factors f ( -k 2 ) that grow at least as fast as k 4 . A summary of the results can be formulated as follows.  \nIf pr = p θ = 0, the system of field equations is over-determined; nevertheless, it can admit an approximate solution in the regimes of small and large values of r . The relevance of this solution is its correspondence to the case in which the original source is a Dirac delta function. Indeed, recalling that the equivalent solution in general relativity would be the Schwarzschild one [118], it follows that (under the approximations considered) the solution with pr = p θ = 0 can be regarded as an analogous of Schwarzschild for nontrivial form factors.  \nOn the other hand, the e ff ective field equations can be solved in the whole spacetime if one considers nontrivial e ff ective pressures pr and p θ . However, there is considerable arbitrariness in the choice of pressure components if one considers equations of state without a guiding physical motivation. In Sec. 6, we discussed some possibilities of equations of state and the solutions associated. In particular, we argued that the choice pr = -ρ e ff makes the shift function B ( r ) be identically zero and, as a consequence, for large r , the resultant metric does not recover the modified Newtonian potential. We also showed that the Newtonian-limit behavior can be matched with the equation of state pr = [ A ( r ) -1] ρ e ff . The explicit solution of the e ff ective  \nTABLE I: Solutions obtained in this work for several form factors f ( z ); the shift function B ( r ) presented follows from the equations of state of Sec. 6.4, namely, pr = [ A ( r ) -1] ρ e ff . The function Z ( x ) is defined in Eq. (7.25), while the coe ffi cients Ci , C i , j , and A i , j are defined, respectively, in Eqs. (7.38), (7.53), and (7.59). The last column refers to the section where the solution is discussed.  \nfield equations with this equation of state was obtained for e ff ective sources originating from local and nonlocal higher-derivative gravity models, as summarized in Table I.  \nWe remark that all the solutions found, regardless of the choice of the e ff ective pressures, present a black hole mass gap. In other words, the metric does not have an apparent horizon if the mass is smaller than a critical value.  \nA significant part of the work was dedicated to studying the regularity of the solutions and how they compare with the Newtonian limit metric. We showed that if the form factor f ( -k 2 ) of the model grows at least as fast as k 4 , there are choices of equations of state that result in the curvature invariants being singularity-free. The same result is generalized to guarantee the regularity of sets of curvature invariants containing covariant derivatives; however, in this case, the behavior of the form factor should be improved accordingly. The whole set of invariants that are polynomial in curvature tensors and their covariant derivatives can be regular if f ( -k 2 ) grows faster than any polynomial (i.e., in some nonlocal gravity models).  \nRegarding the comparison with the nonrelativistic limit, we showed that the solutions related to the equation of state pr = [ A ( r ) -1] ρ e ff match the Newtonian-limit solutions not only in the regime of large r but also for small r (if the e ff ective source is regular). This means that the curvature invariants show the same leading behavior in these regions; however, higher-order derivative invariants may di ff er at r = 0.  \nWe highlight that even though our discussion was centered on smeared sources from higher-derivative gravity models, the same e ff ective delta sources can originate from other quantum gravity approaches, e.g., noncommutative geometry, generalized uncertainty principle models, and string theory [69, 90, 102-112]. Hence, the general results derived here can also be applied in other frameworks.  \nIt has yet to be determined if the singularities of the classical solutions of general relativity can be resolved at a more fundamental quantum level. Nevertheless, our results may suggest that they can be regularized by higher-derivative structures in the gravitational action. Indeed, the higher-derivative modifications of the Schwarzschild metric presented here are nonsingular and display the same qualitative behavior of the solutions reported in [86]. Although a rigorous analysis of the complete field equations in the case of models with six or more derivatives is still pending, this work can motivate and be useful for further research on this important topic.', 'ACKNOWLEDGMENTS': "This work was supported by the Basic Research Program of the Science, Technology, and Innovation Commission of Shenzhen Municipality (grant no. JCYJ20180302174206969). B.L.G. is supported by Primus grant PRIMUS / 23 / SCI / 005 from Charles University.  \nFor the sake of completeness, here we discuss the case of the equation of state that comes from the choice σ ( r ) ∝ r ρ e ff ( r ) in the framework of Sec. 6.3. This example is interesting for at least three reasons: i ) The pressure components are regular, but they do not vanish as r → 0, ii ) the equation of state introduces an odd power of r in the e ff ective pressures, and iii ) the same equation of state was used in Ref. [108] in the context of noncommutative geometry with the Gaussian e ff ective source (7.3), motivating a comparison of that solution with our metric (7.19).  \nHence, let us choose the equation of state  \npr = [ µ r 2 A ( r ) -1 ] ρ e ff , (A1)  \nwhere µ is a parameter with dimension of mass, which could be linked to the massive parameters of the higher-derivative model. In this case, the conservation equation (2.10) defines  \np θ = µ r 8 { 3 r ρ e ff A ' + [ (6 + rB ' ) ρ e ff + 2 r ρ ' e ff ] A } -1 2 ( 2 ρ e ff + r ρ ' e ff ) . (A2)  \nBoth pressures satisfy the criteria (1b) and (2b) of Sec. 6.4-but not (3b).  \nAn advantage of using (A1) is the simplicity of the expression for B ( r ), since it makes Eq. (2.5) become the same equation defining the e ff ective mass m ( r ). Therefore, together with the condition of matching Minkowski at infinity, it gives  \nB ( r ) = G µ [ m ( r ) -M ] . (A3)  \nNotice, however, that since the e ff ective pressures do not vanish at the origin, the behavior of this solution near r = 0 might be very di ff erent from the ones presented in the main part of the paper (see Sec. 6). As a simple example, recall that the mass function m ( r ) is an odd function for nonlocal models defined by form factors that grow faster than any polynomial [see Eq. (5.15) and the related discussion]; hence, aside from a constant, the series representation of (A3) has only odd powers of r , whereas, for the same models , the solutions (6.2), (6.10), and (6.21) have only even powers. Let us remind that the absence of odd powers in the series expansion of the metric is crucial to the regularity of scalars containing derivatives of curvature tensors [56, 83], as discussed also in Sec. 4.  \nThe presence of odd powers of r in the series of (A3) is generally expected for higher-derivative models. In fact, from the smallr behavior of the mass function in Eq. (5.6) we obtain the expansion of (A3) around r = 0,  \nB ( r ) = G µ M [ -1 + 4 πρ e ff (0) 3 M r 3 ] + O ( r 4 ) , (A4)  \nif f ( -k 2 ) ∼ k 4 (or faster) for large k . Thus, for the higher-derivative gravity models considered in this paper, one can always guarantee the presence of the term r 3 in the series expansion of B ( r ) if the equation of state (A1) is used. Moreover, the term linear in r is always absent, which ensures the regularity of curvature scalars without covariant derivatives (see Sec. 4).  \nTaking into account the relation between the asymptotic behavior of f ( -k 2 ) and the lowest odd power of r in the series of A ( r ) [see Eq. (5.14)], the Eq. (4.10) guarantees that the scalar /square R is singular at r = 0 for any model with f ( -k 2 ) ∼ k 6 (or faster) for large k (e.g., eighth- and higher-derivative gravity), if the equation of state (A1) is used. In fact, the explicit calculation gives  \n/square R = -32 π G µρ e ff (0) r + O ( r 0 ) . (A5)  \nIn the terms of [83], this happens because, while the first odd power in the series expansion of grr is at least O ( r 5 ) [see Eq. (5.14)], for gtt it is for sure O ( r 3 ), see Eq. (A4).  \nOnly in the case of f ( -k 2 ) ∼ k 4 (like in sixth-derivative gravity) /square R might be regular [83]. This occurs because A ( r ) contains a term a 3 r 3 in the series expansion [ a 3 = -2 π G ρ ' e ff (0) /nequal 0], which opens a possibility for /square R be regular if  \n-40 a 3 -32 π G µρ e ff (0) = 0 , (A6)  \nas it can be deduced from the comparison of (A5) and (4.10). One can also verify that under this exceptional condition, the other scalars with two covariant derivatives (such as ( ∇ λ R µναβ ) 2 and R µν /square R µν ) are also regular.  \nLast but not least, in what concerns the comparison with the Newtonian limit, we notice that the higher powers of r entering the equations of state considered in this Appendix make the e ff ective pressures vanish more slowly than the ones of Sec. 6.4, as r increases. Thus, the solution in the regime of large r might likely converge to the Newtonian-limit metric in a more slow way.  \n<!-- image -->  \nFIG. 8: Kretschmann scalar K and the ratio K / K Sch. with K evaluated at the metrics (A7) (solid line), (7.19) (dashed line), and at the Newtonian-limit metric (3.4) (dotted line) with potential (7.14).  \n<!-- image -->  \nAs a concrete example, let us consider the case of the Gaussian e ff ective source (7.3). The main choice of equation of state used in this paper yields the metric (7.19), while the choice (A1) results in  \nd s 2 = -[ 1 -2 GM r erf ( µ r 2 ) + 2 GM µ √ π e -1 4 µ 2 r 2 ] exp { -G µ M [ 1 -erf ( µ r 2 ) + µ r √ π e -1 4 µ 2 r 2 ]} d t 2 + [ 1 -2 GM r erf ( µ r 2 ) + 2 GM µ √ π e -1 4 µ 2 r 2 ] -1 d r 2 + r 2 d Ω 2 . (A7)  \nAs anticipated, the metric (A7) coincides with the one obtained in [108] through the replacement µ = 1 / √ θ and a straightforward manipulation involving the lower incomplete gamma function γ ( a , x ) using the formula  \nγ ( 3 2 , µ 2 r 2 4 ) = 1 2 [ √ π erf ( µ r 2 ) -µ r e -1 4 µ 2 r 2 ] . (A8)  \nIn Fig. 8a, we plot the Kretschmann scalar for the metrics (7.19), (A7), and for the Newtonian solution with the potential (7.14). It shows that the solution (A7) di ff ers significantly from the others as r → 0. In Fig. 8b, we compare the Kretschmann scalar evaluated at all the three solutions mentioned above with the one for the Schwarzschild solution, Eq. (4.18). The three solutions satisfy K / K Sch. → 1 as r →∞ , i.e., in the IR they tend to the Schwarzschild solution, but the presence of the higher powers of r in the equation of state (A1) makes the solution (A7) converge slightly more slowly.  \n- [1] R. Utiyama and B. S. DeWitt, Renormalization of a classical gravitational field interacting with quantized matter fields , J. Math. Phys. 3 , 608 (1962).\n- [2] K. S. Stelle, Renormalization of Higher Derivative Quantum Gravity , Phys. Rev. D 16 , 953 (1977).\n- [3] M. Asorey, J. L. L'opez and I. L. Shapiro, Some remarks on high derivative quantum gravity , Int. J. Mod. Phys. A 12 , 5711 (1997), arXiv:hep-th / 9610006.\n- [4] N. V. Krasnikov, Nonlocal Gauge Theories , Theor. Math. Phys. 73 , 1184 (1987) [Teor. Mat. Fiz. 73 , 235 (1987)].\n- [5] Yu. V. Kuz'min, Finite nonlocal gravity , Sov. J. Nucl. Phys. 50 , 1011 (1989) [Yad. Fiz. 50 , 1630 (1989)].\n- [6] E. T. Tomboulis, Superrenormalizable gauge and gravitational theories , arXiv:hep-th / 9702146.\n- [7] L. Modesto, Super-renormalizable Quantum Gravity , Phys. Rev. D 86 , 044005 (2012), arXiv:1107.2403.\n- [8] G. 't Hooft and M. Veltman, One loop divergencies in the theory of gravitation , Ann. Inst. H. Poincar'e Phys. Theor. A 20 , 69 (1974).\n- [9] S. Deser and P. van Nieuwenhuisen, One-loop divergences of quantized Einstein-Maxwell fields , Phys. Rev. D 10 , 401 (1974).\n- [10] M. H. Goro ff and A. Sagnotti, Quantum gravity at two loops , Phys. Lett. B160 , 81 (1985).\n- [11] M. H. Goro ff and A. Sagnotti, The ultraviolet behavior of Einstein gravity, Nucl. Phys. B266 , 709 (1986).\n- [12] C. M. Bender and P. D. Mannheim, No-ghost theorem for the fourth-order derivative Pais-Uhlenbeck oscillator model , Phys. Rev. Lett. 100 , 110402 (2008), arXiv:0706.0207.\n- [13] C. M. Bender and P. D. Mannheim, Exactly solvable PT-symmetric Hamiltonian having no Hermitian counterpart , Phys. Rev. D 78 , 025022 (2008), arXiv:0804.4190.\n- [14] A. Salvio and A. Strumia, Quantum mechanics of 4-derivative theories , Eur. Phys. J. C 76 , 227 (2016), arXiv:1512.01237.\n- [15] D. Anselmi, On the quantum field theory of the gravitational interactions , J. High Energy Phys. 06 , 086 (2017), arXiv:1704.07728.  \n- [16] J. F. Donoghue and G. Menezes, Unitarity, stability and loops of unstable ghosts , Phys. Rev. D 100 , 105006 (2019), arXiv:1908.02416.\n- [17] T. D. Lee and G. C. Wick, Negative metric and the unitarity of the S Matrix , Nucl. Phys. B 9 , 209 (1969).\n- [18] T. D. Lee and G. C. Wick, Finite Theory of Quantum Electrodynamics , Phys. Rev. D 2 , 1033 (1970).\n- [19] R. E. Cutkosky, P. V. Landsho ff , D. I. Olive and J. C. Polkinghorne, A non-analytic S matrix , Nucl. Phys. B 12 , 281 (1969).\n- [20] L. Modesto and I. L. Shapiro, Superrenormalizable quantum gravity with complex ghosts , Phys. Lett. B 755 , 279 (2016), arXiv:1512.07600.\n- [21] L. Modesto, Super-renormalizable or finite Lee-Wick quantum gravity , Nucl. Phys. B 909 , 584 (2016), arXiv:1602.02421.\n- [22] D. Anselmi and M. Piva, Perturbative unitarity of Lee-Wick quantum field theory , Phys. Rev. D 96 , 045009 (2017), arXiv:1703.05563.\n- [23] L. Modesto and L. Rachwał, Super-renormalizable and finite gravitational theories , Nucl. Phys. B 889 228 (2014), arXiv:1407.8036.\n- [24] V. A. Wynne and G. H. Derrick, A Theory of Gravitation Incorporating the Quadratic Action Principle of Relativity , Nuovo Cimento 15B , 181 (1973).\n- [25] K. S. Stelle, Classical Gravity with Higher Derivatives , Gen. Rel. Grav. 9 , 353 (1978).\n- [26] V. P. Frolov and G. A. Vilkovisky, Spherically Symmetric Collapse in Quantum Gravity , Phys. Lett. B 106 , 307 (1981).\n- [27] V. P. Frolov and G. A. Vilkovisky, Quantum gravity removes classical singularities and shortens the life of black holes , in: Proceedings of the Second Marcel Grossmann Meeting on the Recent Developments of General Relativity , edited by R. Ru ffi ni (North-Holland Publishing Company, 1982), pp. 455.\n- [28] E. Pechlaner and R. Sexl, On Quadratic Lagrangians in General Relativity , Commun. Math. Phys. 2 , 165 (1966).\n- [29] P. Havas, On Theories of Gravitation with Higher-Order Field Equations , Gen. Rel. Grav. 8 , 631 (1977).\n- [30] V. P. Frolov and I. L. Shapiro, Black Holes in Higher Dimensional Gravity Theory with Quadratic in Curvature Corrections , Phys. Rev. D 80 , 044034 (2009), arXiv:0907.1411.\n- [31] W. Nelson, Static Solutions for 4th order gravity , Phys. Rev. D 82 , 104026 (2010), arXiv:1010.3986.\n- [32] H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Black Holes in Higher-Derivative Gravity , Phys. Rev. Lett. 114 , 171601 (2015), arXiv:1502.01028.\n- [33] H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Spherically Symmetric Solutions in Higher-Derivative Gravity , Phys. Rev. D 92 , 124019 (2015), arXiv:1508.00010.  \n[34] Y. F. Cai, G. Cheng, J. Liu, M. Wang and H. Zhang,  \nFeatures and stability analysis of non-Schwarzschild black hole in quadratic gravity  \n,  \nJ. High Energy Phys.  \n01  \n108 (2016), arXiv:1508.04776.  \n- [35] X. H. Feng and H. Lu, Higher-Derivative Gravity with Non-minimally Coupled Maxwell Field , Eur. Phys. J. C 76 178 (2016), arXiv:1512.09153.\n- [36] B. Holdom and J. Ren, QCD analogy for quantum gravity , Phys. Rev. D 93 , 124030 (2016), arXiv:1512.05305.\n- [37] K. Lin, W. L. Qian, A. B. Pavan and E. Abdalla, (Anti-) de Sitter Electrically Charged Black Hole Solutions in Higher-Derivative Gravity , Europhys. Lett. 114 60006 (2016), arXiv:1607.04473.\n- [38] B. Holdom and J. Ren, Not quite a black hole , Phys. Rev. D 95 , 084034 (2017), arXiv:1612.04889.\n- [39] H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Lichnerowicz Modes and Black Hole Families in Ricci Quadratic Gravity , Phys. Rev. D 96 046006 (2017), arXiv:1704.05493.\n- [40] K. Kokkotas, R. A. Konoplya and A. Zhidenko, Non-Schwarzschild black-hole metric in four dimensional higher derivative gravity: analytical approximation , Phys. Rev. D 96 064007 (2017), arXiv:1705.09875.\n- [41] K. Goldstein and J. J. Mashiyane, Ine ff ective Higher Derivative Black Hole Hair , Phys. Rev. D 97 024015 (2018), arXiv:1703.02803.\n- [42] J. Podolsk'y, R. ˇ Svarc, V. Pravda and A. Pravdov'a, Explicit black hole solutions in higher-derivative gravity , Phys. Rev. D 98 021502 (2018), arXiv:1806.08209.\n- [43] R. ˇ Svarc, J. Podolsk'y, V. Pravda and A. Pravdov'a, Exact black holes in quadratic gravity with any cosmological constant , Phys. Rev. Lett. 121 231104 (2018), arXiv:1806.09516.\n- [44] A. Bonanno and S. Silveravalle, Characterizing black hole metrics in quadratic gravity , Phys. Rev. D 99 , 101501 (2019), arXiv:1903.08759.\n- [45] R. A. Konoplya and A. F. Zinhailo, Hawking radiation of non-Schwarzschild black holes in higher derivative gravity: a crucial role of grey-body factors , Phys. Rev. D 99 104060 (2019), arXiv:1904.05341.\n- [46] B. Holdom, A ghost and a naked singularity; facing our demons , arXiv:1905.08849.\n- [47] A. Salvio and H. Veermae, Horizonless ultracompact objects and dark matter in quadratic gravity , J. Cosmol. Astropart. Phys. 02 , 018 (2020), arXiv:1912.13333.\n- [48] U. Aydemir and J. Ren, Dark sector production and baryogenesis from not quite black holes , Chin. Phys. C 45 , 075103 (2021), arXiv:2011.13154.\n- [49] A. Bonanno and S. Silveravalle, The gravitational field of a star in quadratic gravity , J. Cosmol. Astropart. Phys. 08 , 050 (2021), arXiv:2106.00558.\n- [50] B. Holdom, 2-2-holes simplified , Phys. Lett. B 830 , 137142 (2022), arXiv:2202.08442.\n- [51] J. Daas, K. Kuijpers, F. Saueressig, M. F. Wondrak and H. Falcke, Probing quadratic gravity with the Event Horizon Telescope , Astron. Astrophys. 673 , A53 (2023), arXiv:2204.08480.\n- [52] S. Silveravalle and A. Zuccotti, Phase diagram of Einstein-Weyl gravity , Phys. Rev. D 107 , 6 (2023), arXiv:2210.13877.\n- [53] L. Modesto, T. de Paula Netto and I. L. Shapiro, On Newtonian singularities in higher derivative gravity models , J. High Energy Phys. 1504 , 098 (2015), arXiv:1412.0740.\n- [54] B. L. Giacchini, On the cancellation of Newtonian singularities in higher-derivative gravity , Phys. Lett. B 766 , 306 (2017), arXiv:1609.05432.\n- [55] B. L. Giacchini and T. de Paula Netto, Weak-field limit and regular solutions in polynomial higher-derivative gravities , Eur. Phys. J. C 79 , 217 (2019), arXiv:1806.05664.  \n- [56] N. Burzill'a, B. L. Giacchini, T. de Paula Netto and L. Modesto, Higher-order regularity in local and nonlocal quantum gravity , Eur. Phys. J. C 81 , 462 (2021), arXiv:2012.11829.\n- [57] B. L. Giacchini and T. de Paula Netto, E ff ective delta sources and regularity in higher-derivative and ghost-free gravity , J. Cosmol. Astropart. Phys. 1907 , 013 (2019), arXiv:1809.05907.\n- [58] N. Burzill'a, B. L. Giacchini, T. de Paula Netto and L. Modesto, Newtonian potential in higher-derivative quantum gravity , Phys. Rev. D 103 , 064080 (2021), arXiv:2012.06254.\n- [59] V. P. Frolov, A. Zelnikov and T. de Paula Netto, Spherical collapse of small masses in the ghost-free gravity , J. High Energy Phys. 1506 , 107 (2015), arXiv:1504.00412.\n- [60] V. P. Frolov, Mass-gap for black hole formation in higher derivative and ghost free gravity , Phys. Rev. Lett. 115 , 051102 (2015), arXiv:1505.00492.\n- [61] J. Boos, V. P. Frolov and A. Zelnikov, Gravitational field of static p-branes in linearized ghost-free gravity , Phys. Rev. D 97 , 084021 (2018), arXiv:1802.09573.\n- [62] J. Boos, J. Pinedo Soto and V. P. Frolov, Ultrarelativistic spinning objects in nonlocal ghost-free gravity , Phys. Rev. D 101 , 124065 (2020), arXiv:2004.07420.\n- [63] L. Buoninfante and B. L. Giacchini, Light bending by a slowly rotating source in quadratic theories of gravity , Phys. Rev. D 102 , 024020 (2020), arXiv:2005.05355.\n- [64] L. Buoninfante, A. S. Cornell, G. Harmsen, A. S. Koshelev, G. Lambiase, J. Marto and A. Mazumdar, Towards nonsingular rotating compact object in ghost-free infinite derivative gravity , Phys. Rev. D 98 , 084041 (2018), arXiv:1807.08896.\n- [65] I. Kol'aˇr and A. Mazumdar, NUT charge in linearized infinite derivative gravity , Phys. Rev. D 101 , 124005 (2020), arXiv:2004.07613.\n- [66] C. Heredia, I. Kol'aˇr, J. Llosa, F. J. M. Torralba and A. Mazumdar, Infinite-derivative linearized gravity in convolutional form , Class. Quant. Grav. 39 , 085001 (2022), arXiv:2112.05397.\n- [67] A. Accioly, B. L. Giacchini and I. L. Shapiro, Low-energy e ff ects in a higher-derivative gravity model with real and complex massive poles , Phys. Rev. D 96 , 104004 (2017), arXiv:1610.05260.\n- [68] V. P. Frolov and A. Zelnikov, Head-on collision of ultrarelativistic particles in ghost-free theories of gravity , Phys. Rev. D 93 , 064048 (2016), arXiv:1509.03336.\n- [69] A. A. Tseytlin, On singularities of spherically symmetric backgrounds in string theory , Phys. Lett. B 363 , 223 (1995), arXiv:hep-th / 9509050.\n- [70] W. Siegel, Stringy gravity at short distances , arXiv:hep-th / 0309093.\n- [71] T. Biswas, A. Mazumdar and W. Siegel, Bouncing universes in string-inspired gravity , J. Cosmol. Astropart. Phys. 03 , 009 (2006), arXiv:hep-th / 0508194.\n- [72] T. Biswas, E. Gerwick, T. Koivisto and A. Mazumdar, Towards singularity and ghost free theories of gravity , Phys. Rev. Lett. 108 , 031101 (2012), arXiv:1110.5249.\n- [73] J. Edholm, A. S. Koshelev and A. Mazumdar, Behavior of the Newtonian potential for ghost-free gravity and singularity-free gravity , Phys. Rev. D 94 , 104033 (2016), arXiv:1604.01989.\n- [74] E. Kilicarslan, Weak Field Limit of Infinite Derivative Gravity , Phys. Rev. D 98 , 064048 (2018), arXiv:1808.00266.\n- [75] L. Buoninfante, A. S. Koshelev, G. Lambiase and A. Mazumdar, Classical properties of non-local, ghost- and singularity-free gravity , J. Cosmol. Astropart. Phys. 09 , 034 (2018), arXiv:1802.00399.\n- [76] L. Buoninfante, A. S. Koshelev, G. Lambiase, J. Marto and A. Mazumdar, Conformally-flat, non-singular static metric in infinite derivative gravity , J. Cosmol. Astropart. Phys. 06 , 014 (2018), arXiv:1804.08195.\n- [77] E. Kilicarslan, pp-waves as exact solutions to ghost-free infinite derivative gravity , Phys. Rev. D 99 , 124048 (2019), arXiv:1903.04283.\n- [78] S. Dengiz, E. Kilicarslan, I. Kol'aˇr and A. Mazumdar, Impulsive waves in ghost free infinite derivative gravity in anti-de Sitter spacetime , Phys. Rev. D 102 , 044016 (2020), arXiv:2006.07650.\n- [79] I. Kol'aˇr, T. M'alek and A. Mazumdar, Exact solutions of nonlocal gravity in a class of almost universal spacetimes , Phys. Rev. D 103 , 124067 (2021), arXiv:2103.08555.\n- [80] I. Kol'aˇr, T. M'alek, S. Dengiz and E. Kilicarslan, Exact gyratons in higher and infinite derivative gravity , Phys. Rev. D 105 , 044018 (2022), arXiv:2107.11884.\n- [81] I. Kol'aˇr and T. M'alek, Infinite derivative gravity resolves nonscalar curvature singularities , Phys. Rev. D 108 , L081502 (2023), arXiv:2307.07285.\n- [82] V. P. Frolov and J. Pinedo Soto, Nonlocal modification of the Kerr metric , Symmetry 15 , 1771 (2023), arXiv:2308.00114.\n- [83] B. L. Giacchini, T. de Paula Netto and L. Modesto, Action principle selection of regular black holes , Phys. Rev. D 104 , 084072 (2021), arXiv:2105.00300.\n- [84] B. L. Giacchini and T. de Paula Netto, Regular Black Holes from Higher-Derivative E ff ective Delta Sources , in: Regular Black Holes: Towards a New Paradigm of Gravitational Collapse , edited by C. Bambi (Springer, Singapore, 2023), arXiv:2307.12357.\n- [85] Y. D. Li, L. Modesto and L. Rachwał, Exact solutions and spacetime singularities in nonlocal gravity , J. High Energy Phys. 12 , 173 (2015), arXiv:1506.08619.\n- [86] B. Holdom, On the fate of singularities and horizons in higher derivative gravity , Phys. Rev. D 66 , 084010 (2002), arXiv:hep-th / 0206219.\n- [87] J. M. Pawlowski and J. Trankle, E ff ective action and black hole solutions in asymptotically safe quantum gravity , arXiv:2309.17043.\n- [88] J. Daas, C. Laporte, F. Saueressig and T. van Dijk, Rethinking the E ff ective Field Theory formulation of Gravity , arXiv:2405.12685.\n- [89] B. L. Giacchini and I. Kol'aˇr, Toward regular black holes in sixth-derivative gravity , arXiv:2406.00997.\n- [90] L. Modesto, J. W. Mo ff at and P. Nicolini, Black holes in an ultraviolet complete quantum gravity , Phys. Lett. B 695 , 397 (2011), arXiv:1010.0680.\n- [91] Y. Zhang, Y. Zhu, L. Modesto and C. Bambi, Can static regular black holes form from gravitational collapse? , Eur. Phys. J. C 75 , 96 (2015), arXiv:1404.4770.  \n- [92] C. Bambi, L. Modesto and Y. Wang, Lee-Wick black holes , Phys. Lett. B 764 , 306 (2017), arXiv:1611.03650.\n- [93] C. Bambi, D. Malafarina and L. Modesto, Terminating black holes in asymptotically free quantum gravity , Eur. Phys. J. C 74 , 2767 (2014), arXiv:1306.1668.\n- [94] C. Bambi, D. Malafarina and L. Modesto, Black supernovae and black holes in non-local gravity , J. High Energy Phys. 04 , 147 (2016), arXiv:1603.09592.\n- [95] L. Modesto, Super-renormalizable Multidimensional Quantum Gravity , Astron. Rev. 8 , 4 (2013), arXiv:1202.3151.\n- [96] P. Bueno and P. A. Cano, On black holes in higher-derivative gravities , Class. Quant. Grav. 34 , 175008 (2017), arXiv:1703.04625.\n- [97] V. P. Frolov, Notes on nonsingular models of black holes , Phys. Rev. D 94 , 104056 (2016), arXiv:1609.01758.\n- [98] C. Bejarano, G. J. Olmo and D. Rubiera-Garcia, What is a singular black hole beyond General Relativity? , Phys. Rev. D 95 , 064043 (2017), arXiv:1702.01292.\n- [99] R. Carballo-Rubio, F. Di Filippo, S. Liberati and M. Visser, Geodesically complete black holes , Phys. Rev. D 101 , 084047 (2020) arXiv:1911.11200.\n- [100] T. Berry, A. Simpson and M. Visser, General class of 'quantum deformed' regular black holes , Universe 7 , 165 (2021), arXiv:2102.02471.\n- [101] K. A. Bronnikov and S. G. Rubin, Black Holes, Cosmology, and Extra Dimensions (World Scientific, Singapore, 2012).\n- [102] M. Cadoni, M. Oi and A. P. Sanna, E ff ective models of nonsingular quantum black holes , Phys. Rev. D 106 , 024030 (2022), arXiv:2204.09444.\n- [103] M. Cadoni and A. P. Sanna, Nonsingular black holes from conformal symmetries , Class. Quant. Grav. 40 , 145012 (2023) arXiv:2302.06401.\n- [104] A. Akil, M. Cadoni, L. Modesto, M. Oi and A. P. Sanna, Semiclassical spacetimes at super-Planckian scales from delocalized sources , Phys. Rev. D 108 , 044051 (2023), arXiv:2211.01657.\n- [105] P. Nicolini, A. Smailagic and E. Spallucci, Noncommutative geometry inspired Schwarzschild black hole , Phys. Lett. B 632 , 547 (2006), gr-qc / 0510112.  \n[106] J. R. Mureika and E. Spallucci,  \nVector unparticle enhanced black holes: exact solutions and thermodynamics  \n, Phys. Lett. B  \n693  \n, 129  \n(2010), arXiv:1006.4556.  \n- [107] P. Gaete, J. A. Helayel-Neto and E. Spallucci, Un-graviton corrections to the Schwarzschild black hole , Phys. Lett. B 693 , 155 (2010), arXiv:1005.0234.\n- [108] P. Nicolini and E. Spallucci, Noncommutative geometry inspired wormholes and dirty black holes , Class. Quant. Grav. 27 , 015010 (2010), arXiv:0902.4654.\n- [109] P. Nicolini, Nonlocal and generalized uncertainty principle black holes , arXiv:1202.2102.\n- [110] M. Isi, J. Mureika and P. Nicolini, Self-completeness and the generalized uncertainty principle , J. High Energy Phys. 11 , 139 (2013), arXiv:1310.8153.  \n[111] P. Nicolini, E. Spallucci and M. F. Wondrak,  \nQuantum Corrected Black Holes from String T-Duality  \n, Phys. Lett. B  \n797  \n, 134888 (2019),  \narXiv:1902.11242.  \n- [112] M. Knipfer, S. Koppel, J. Mureika and P. Nicolini, Generalized Uncertainty Principle and Black Holes in Higher Dimensional SelfComplete Gravity , J. Cosmol. Astropart. Phys. 08 , 008 (2019), arXiv:1905.03233.\n- [113] S. A. Hayward, Formation and evaporation of regular black holes , Phys. Rev. Lett. 96 , 031103 (2006), arXiv:gr-qc / 0506126.\n- [114] N. Burzill'a, B. L. Giacchini, T. de Paula Netto and L. Modesto, Regular multi-horizon Lee-Wick black holes , J. Cosmol. Astropart. Phys. 11 , 067 (2023), arXiv:2308.12810.\n- [115] I. Quandt and H. J. Schmidt, The Newtonian limit of fourth and higher order gravity , Astron. Nachr. 312 , 97 (1991), arXiv:gr-qc / 0109005.\n- [116] A. Pais and G. E. Uhlenbeck, On Field theories with nonlocalized action , Phys. Rev. 79 , 145 (1950).\n- [117] I.S. Gradshteyn and I.M. Ryzhik, Tables of integrals, series and products , 7th ed., edited by A. Je ff rey and D. Zwillinger (Academic Press, San Diego, 2007).\n- [118] H. Balasin and H. Nachbagauer, The energy-momentum tensor of a black hole, or what curves the Schwarzschild geometry? , Class. Quant. Grav. 10 , 2271 (1993), arXiv:gr-qc / 9305009."}