On effective models of regular black holes inspired by higher-derivative and nonlocal gravity

In this work we study static spherically symmetric solutions of effective field equations related to local and nonlocal higher-derivative 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 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 effective 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.


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 effective 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][5][6][7].
The 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.
The 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.
The underlying higher-derivative models we focus on belong to a large class of theories quadratic in curvatures, described by the action where G is the Newton constant and F 1,2 ( ) are form factors, analytic functions of the d'Alembert operator subjected to the constraint1 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: (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][9][10][11].
(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][13][14][15][16].
(iii) If F 1,2 (z) are polynomials of degree N 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 à la Lee-Wick [17][18][19][20][21][22].This subclass of polynomial-derivative gravity is also known as Lee-Wick gravity [20,21].
(iv) For nonlocal form factors of the type where 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][5][6][7].
Quadratic structures of the type R µναβ F 3 ( )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 affect 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 affected by the presence of these 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).
With 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 .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. Te 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].
Even 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].
To 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. 4The conclusion of this study can be summarized as follows: (I) All the solutions found are regular at r = 0.
(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.
Moreover, it was speculated whether some appealing properties of the regular solutions of fourth-derivative gravity could be shared by their higher-derivative counterparts, namely: (III) Nonlinearities are suppressed sufficiently close to r = 0; the solutions behave like the weak-field ones.
(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.
Last 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 where 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 The 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.Simply 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 where dΩ 2 = dθ 2 + sin 2 θdφ 2 and primes denote differentiation 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]) where In an equivalent form, ϕ(r) can be obtained by solving the equation where is 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 effective delta source (1.9) (see, e.g., [57,84]).Hence, in flat spacetime, the effect 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).The main idea, now, is to obtain a curved-space generalization of (1.8) to static spherically symmetric spacetimes in the form where T µ ν is an effective energy-momentum tensor reproducing the effects of the higher derivatives with being ρ eff 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 off, f (z) ≡ 1 and Eq.(1.10) results in the Schwarzschild solution.
The 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 ( ) together with the Einstein tensor in the equations of motion, namely, where 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 ( ) in the flat-spacetime form, from the inversion of the operator f ( ) one can recast (1.12) into (1.10)[7,85,[90][91][92][93][94][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 effective energymomentum tensor, i.e., ∇ µ T µ ν = 0 (see, e.g., [7]).This leads to the introduction of effective radial (p r ) and tangential (p θ ) pressures in the effective energy-momentum tensor, As 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].
In the present work, we study static spherically symmetric solutions to the system of equations (1.10), with the general line element where A(r) and B(r) are two functions to be determined.Differently from previous works [90,92], here we are interested in solutions with a nontrivial shift function, B(r) 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 g tt −g rr . 6 We 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 It 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, differs 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).
Finally, let us mention that the formulation of the problem in terms of the effective 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][103][104][105][106][107][108][109][110][111][112]).Therefore, the detailed study of the solutions of Einstein's equation (1.10) with the effective source (1.9) is relevant not only for higher-derivative gravities but also for alternative models of quantum gravity.
The 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 effective 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 effective field equations on general grounds, i.e., without particularizing to specific form factors.We also discuss the effect 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 effective 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 effective 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 ℏ = 1.

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 6 See [88,89,96] for further evidences that solutions in higher-derivative gravity models might require a nontrivial shift function. 7Several definitions of "spacetime regularity" exist in the literature, e.g., referring to the regularity of metric components or Christoffel symbols in a given coordinate chart, to the regularity of a set of curvature invariants, or to the geodesic completeness of the spacetime (see, e.g., [97][98][99][100] and references therein).As a matter of convention, in this work, we characterize a regular black hole by the absence of singularities in its Kretschmann scalar R µναβ R µναβ .For static and spherically symmetric geometries, this also implies in the regularity of all the curvature invariants constructed contracting an arbitrary number of Riemann and metric tensors [101].An extended definition of regularity, which also considers the absence of singularities in curvature invariants containing covariant derivatives of curvatures, will be discussed in Sec. 4. and Therefore, using the effective energy momentum tensor (1.13), the field equations (1.10) are equivalent to ) and The 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).The formal solution for A(r) can be directly obtained from Eq. (2.4).It only depends on the effective delta source, and it is given by where the mass function represents the effective mass inside a sphere of radius r.
Given the solution (2.7) for A(r), the function B(r) can be obtained through (2.5).The formal solution reads (2.9) The 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 effective pressures tend to zero sufficiently 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.Finally, since the solution for the function B(r) depends on the choice of the pressure components of the effective energymomentum tensor T µ ν , it is mandatory to comment on these other components.One of the simplest possible choices is to fix the effective radial pressure p r via an equation of state involving this component and the effective source ρ eff .In this case, both A(r) and B(r) are directly determined in terms of the source ρ eff .Then, the component p θ is given by the last equation of motion (2.6).This is completely equivalent to fixing p θ using the effective conservation equation ∇ µ T µ r = 0, namely, from which one can easily determine p θ once ρ eff and p r (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 p r = p r (ρ eff ), their motivations, and the properties of the corresponding function B(r).
As a conclusion of this section, the solution of the effective field equations (1.10) can be expressed as (2.11)

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 and only keep quantities up to the first order in a(r), b(r) and their derivatives.Also, in the nonrelativistic limit, one sets p r = p θ = 0. Thus, the linearized Eqs.(2.4) and (2.5) are whose solutions are the functions entering the linearized version of the metric (1.14), i.e., Since 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 with m(r) defined in (2.8), while b(r) is easily obtained from Finally, since ρ eff , a(r), and b(r) are proportional to M, the nontrivial component (2.10) of the continuity equation, is already O(M 2 ), implying that it is verified within the approximation.The Newtonian potential ϕ(r) is defined through the relation thus from (3.4) we obtain Taking the derivative of the previous equation, we find where we used (3.2) and (3.3).Therefore, Substituting these expressions in (3.4), we obtain the Newtonian metric in the form (1.5), presented in the Introduction.

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 sufficiently small r is, indeed, well approximated by the Newtonian-limit solutions.

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 ) ) ) If the four functions K i (r) are regular, all the curvature invariants without covariant derivatives are bounded as well.
In an equivalent way, one can use the Kretschmann scalar, to 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 sufficient that all functions K i (r) are finite.Thus, if K is bounded, all the K i (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).
To analyze the behavior for small r, let us assume that A(r) and B(r) are analytic8 around r = 0 and apply the power series expansion, into the expressions for the functions K i (r): ) Therefore, the regularity of the curvature invariants requires or, in other words, In 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 effective sources of models that contain six or more metric derivatives.
To 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]: Definition 1.Given a continuous bounded function ξ : [0, +∞) → R, we say that ξ(r) is p-regular for an integer p 1 if ξ(r) is at least 2p-times continuously differentiable and the first p odd-order derivatives of ξ(r) vanish as r → 0, namely, If 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.
Indeed, if a function ξ : [0, +∞) → R is p-regular, this means that also the map Ξ : R 3 → R given by Ξ( r) = ξ(| r|) is of class C 2p .On the other hand, if lim r→0 ξ (2p+1) (r) 0, then (switching to Cartesian coordinates) we have depending on whether x = 0 is approached from the left or the right.This means that Ξ is not (2p + 1)-times differentiable, 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. General 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 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 2n+1 for some n ∈ N, then there might exist curvature-derivative scalars with 2n 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.
A simple useful example is provided by the scalar R, which for the metric (1.14) reads 4) . (4.10) Once 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 Hence, the scalar R diverges like 1/r unless 5a 3 + 3b 3 = 0; in particular, it is regular if 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 effective 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.

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.
First, let us suppose that the spacetime is asymptotically flat such that we can write with 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. The next step is to apply this procedure and calculate the linearized version of the curvature tensor, (R µν αβ ) lin .Using (4.1), its components read ) ) ) Notice 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.
The nonlinear part of R µν αβ is then proportional to the differences which turn out to be proportional to the derivatives of B(r): ) ) In the same spirit as before, using the power series representation (4.3) for A(r) and B(r) we find, around r = 0, Therefore, if the curvature invariants are regular, i.e., if the relations in (4.5) are true, we get showing 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 K i , it follows that, for small r, such scalars behave approximately like the linearized ones.Furthermore, 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 effective pressures tend to zero sufficiently 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 p r vanishes sufficiently fast as r → 0.
It 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 In what concerns the curvature-derivative invariants, the situation is slightly different.Using once more R as an example and proceeding in the same way as before, we split (4.10) in terms of the linear and nonlinear parts, and using the series representation (4.3) for regular metrics (i.e., with Therefore, if R is finite, that is it, if 5a 3 + 3b 3 = 0, we see that its linear and nonlinear parts approach a constant when r → 0. The behavior of ( R) non−lin marks a difference between curvature-derivative and curvature invariants since the nonlinear part of the latter goes to zero when r → 0. Note, however, that if ( R) lin is finite at r = 0, the entire R is also finite (because its nonlinear part cannot generate a singularity).
It is also interesting to notice that, because of the covariant derivatives, ( R) non−lin contains terms that do not depend on the function B(r) [see Eq. (4.20)].This is different from what happens in the case of scalars without covariant derivatives.Therefore, even if B(r) = 0, R can receive a contribution from the nonlinear part of the metric.

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 effective 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 effective 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.
To start, let us recall some relevant results on the smeared delta source ρ eff (r); afterward, we analyze the consequences of these properties to the mass function (2.8) and, finally, to A(r).The properties of the effective 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.
(1a) If f (−k 2 ) → 1 when k 2 → 0 (in other words, if general relativity is recovered in the IR limit), then lim r→∞ ρ eff (r) = 0. (5.1) In particular, this is related to the solution being asymptotically flat.
Under the assumption (1a) one can also prove that: (2a) If there exists a k 0 such that f (−k 2 ) grows as k 2 for k > k 0 , then ρ eff (r) is integrable but it is not bounded.In fact, it diverges at the origin, This 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].
(3a) On the other hand, if f (−k 2 ) grows at least as fast as k 4 for sufficiently large k, the effective source is integrable and finite everywhere, i.e., ρ eff (r) is at least 0-regular.In this situation, ρ eff (r) has an absolute maximum at r = 0, max {ρ eff (r)} = lim r→0 ρ eff (r).
(5.3)This 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).
(4a) If the function f (−k 2 ) asymptotically grows at least as fast as k 4+2N for an integer N 1, then the effective source ρ eff (r) is (at least) N-regular (see Definition 1 in Sec.4.1).In other words, it is 2N times differentiable and their first N odd-order derivatives vanish at the origin, (5.4) (5a) As a corollary of the previous statement, if f (−k 2 ) asymptotically grows faster than any polynomial, then ρ eff (r) is ∞regular, i.e., it is an even function of r.This happens, e.g., in the nonlocal models with exponential form factors, From the above results it follows that, for any higher-derivative model form factor, the mass function satisfies [57] lim r→∞ m(r) = M. (5.5) Therefore, for a distant observer, the solution (2.7) for A(r) reduces to Schwarzschild10 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 Notice that this proves the regularity condition of Eq. (4.6).Being a continuously differentiable function, and given (5.7), it follows that A(r) is bounded for any higher-derivative model effective source.Thus, there is a certain critical mass M c such that M < M c implies A(r) > 0 for all r.
The 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 (∇r) 2 = A(r), (5.8) that 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.
A direct way of verifying the existence of the critical mass M c for effective sources is to consider the unitary mass function, defined by which 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 the function A(r) is strictly positive, and the metric does not have horizons.The critical mass is, therefore, To summarize, if f (−k 2 ) grows at least as fast as k 2 for sufficiently 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.
In 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.
In addition, it is worth noticing that if m(r) 1 for all r, then m(r) M and the position of the event horizon r + is bounded by the Schwarzschild radius, r + r s = 2GM. (5.12) In general, if ρ eff (r) is positive, the inequality (5.12) is always true.This statement can be proved by observing that if ρ eff (r) is positive, then the mass function m(r) is monotone.In fact, whence m(r) is a monotonically increasing function.Then, from Eq. (5.5) it follows m(r) M. On the other hand, if there is a region where ρ eff (r) < 0, then m(r) must have oscillations, and it may assume values larger than M (see [7,91,92,114] for explicit examples).
The value of the critical mass M c is related to the massive parameters of the model (see Sec. 7 for explicit examples).If M ≫ M c and those parameters are such that the smearing of the source is negligible at the scale 2GM, 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.
It 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 effective 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 effective source (1.9).An intriguing question is whether such features of the effective 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][87][88][89].Nevertheless, we remark that even for general functions A(r), asymptotically flat, everywhere regular geometries must have an even number of horizons [86].
To 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.
Applying Taylor's theorem and taking into account the results (3a) and (4a) on the behavior of the effective delta source near r = 0, it follows that: If the function f (−k 2 ) asymptotically grows at least as fast as k 4+2N for an integer N 0, then A(r) is at least (N + 1)-regular, i.e., it is 2N + 2 times differentiable and d 2n+1 dr 2n+1 A(r) r=0 = 0 for n = 0, . . ., N. (5.14) Notice 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].
Finally, as a corollary, if f (−k 2 ) asymptotically grows faster than any polynomial, A(r) is an even function-which also implies that the effective mass m(r) is an odd function.This can be proved, alternatively, from the result (5a) on the effective delta source: where 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.

GENERAL PROPERTIES OF B(r)
In this section, we analyze several choices for fixing the effective pressures and underlying motivations.For each, we present the general properties of the function B(r) associated with it.

Case of null pressures
Considering p r = p θ = 0, ( the solution (2.9) for B(r) becomes 2) The properties of the function B(r) defined in (6.2) can be directly deduced from the results concerning the effective 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 differentiable but it is not 1-regular, as B ′ (0) 0; while if f (−k 2 ) grows at least as fast as k 4+2N for a certain N 0 and sufficiently large k, B(r) is at least (N + 1)-regular, hence, In 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 ρ eff (r) and A(r) are also even functions.The 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, Indeed, 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 effective energy-momentum tensor (1.13) with p r = p θ = 0.However, a careful reader will notice that in this situation the effective 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 effective 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 < M c ) or in the limits of short and large distances (if M > M c ), as we discuss in what follows.
Firstly, for M < M c , the situation is similar to what happens in the Newtonian approximation (see the discussion in Sec. 3) since the conservation equation (2.10) for p r = p θ = 0 reduces to − 1 2 Similarly, using (2.7) and (6.2) in (2.3), it is possible to prove (after some manipulations) that which means that the field equations are satisfied in a good approximation in the whole spacetime if the mass M is sufficiently small.How small the mass needs to be depends on the critical mass M c , because for M < M c 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.
On the other hand, for M > M c 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 interval [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 But, from Eqs. (5.7), (5.14) and ( 6.3) we have a 0 = 1 and a 1 = b 1 = 0, such that This means that the effective field equations are satisfied in a good approximation for a sufficiently small r around the point r = 0 for any value of the mass M. Note that this is true only for the effective sources of models with more than four derivatives of the metric in the action.In fourth-derivative gravity, we have a 1 , b 1 0; therefore, the effective 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→∞ ρ eff = 0 implies lim r→∞ G θ θ = 0 through the functions A(r) and B(r).

Case p r = α ρ eff
If one sets p r = α ρ eff (6.9) for a certain constant α, but leaves the tangential pressure to be determined by the conservation equation (2.10), the result for B(r) is while p θ is given by Therefore, 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 differs with respect to (6.2) by a multiplicative factor in front of the integral.A 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 choice12 α = 0, using the field equations one gets Thus, if the metric has a horizon at r = r H , it follows 2Gm(r H ) = r H > 0, showing that the term in (6.12) certainly diverges at the horizon.The only possibility for the tangential pressure [see (6.11)] to be regular is if it happens that ρ eff (r) vanishes at r H .In general, the singularity of the tangential pressure at the horizons is related to the singularity of B ′ (r) at r = r H , see Eq. (6.10).These singularities are physical and cannot be removed by changing coordinates, for they also manifest in the curvature invariants.
One 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 effective tangential pressure [see (6.11)], is regular for any higher-derivative gravity model such that f (−k 2 ) grows at least as fast as k 4 for sufficiently large k.However, 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 As 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.

General case of regular pressures
As it was evident in the preceding discussion, the main difficulty 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 p r = −ρ eff was capable of curing pressure divergences at the horizons, with the side effect that B ≡ 0.Moreover, it was argued that the curvature singularity at the horizons was related to a divergence in the effective 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 p r (r) = A(r)σ(r) − ρ eff (r), (6.17)where σ(x) is an N b -regular continuous function (for an integer N b 0) such that the function ξ(x) = xσ(x) is integrable on [0, +∞).Indeed, this choice eliminates the term 1/A(r) in Eq. (2.9), so that the solution for B(r) becomes By 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 (N b + 1)-regular.
To verify that the equation of state (6.17) makes the effective pressures regular, first we notice that if ρ eff (r) is N a -regular and σ(x) is N b -regular, then A(r) is (N a + 1)-regular (see Sec. 5) and p r (r) is at least N-regular, with N ≡ min{N a , N b }.In what concerns the tangential pressure defined by the conservation equation (2.10), it is straightforward to verify that the order of regularity of each of the terms in the right-hand side of (6.19) is either N a or min{N a + 1, N b }; therefore, p θ (r) is at least N-regular, just like p r (r).Of 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ρ eff (r).However, since for effective 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 ρ eff .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 ρ eff (r) is N a -regular, the choice σ(r) ∝ r 2n ρ eff (r) (for a positive integer n) makes B(r) to be (N a + n + 1)-regular.We shall consider another form for the function σ(r) in what follows.

Case p r = [A(r) − 1] ρ eff
As 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: (1b) The effective pressures vanish asymptotically for large r.
(2b) The effective pressures are finite at the horizons.
(3b) The effective pressures vanish at r = 0.While (1b) and (2b) are natural in the framework of Sec.6.3, (3b) actually imposes restrictions on the form of σ(r). 13The physical motivation for these requirements is to guarantee that the solution is asymptotically flat, everywhere regular and that it matches the small-r behavior of the solution obtained in Sec.6.1 (which approximates the solution whose original source is a pointlike mass).
In the most interesting case (for this work) of a regular ρ eff , a simple equation of state satisfying the conditions (1b)-(3b) is (6.20) which follows from the choice σ(r) = ρ eff (r) in the context of Sec.6.3.Accordingly, the solution (6.18) for B(r) reads while the tangential pressure (6.19) is It is straightforward to verify that (6.20) and (6.22) indeed vanish as r → 0 if ρ eff is regular. 14Also, from the general results of Sec.6.3, we know that if ρ eff is N-regular, p r and p θ are bounded and (at least) N-regular.A finer result is obtained in the case ρ eff (r) is analytic 15 at r = 0-in that case p r 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 2N+1 .
From the features of ρ eff 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 ρ eff is N-regular (or, equivalently, if f (−k 2 ) grows at least as fast as k 4+2N for a certain N 0 and sufficiently 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.
Other choices of equations of state can also reproduce the same aforementioned physical effect on the solutions.Indeed, this procedure has considerable ambiguity; for example, any equation of state in the form p r (r) = A ℓ (r) − 1 ρ eff (r), ℓ = 1, 2, . . ., (6.23) following from σ(r) = A ℓ−1 (r)ρ eff (r), would also lead to the above properties (1b)-(3b).In this work, for the sake of simplicity, we adopt (6.20).* * * As 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 effective pressures satisfying the conditions (1b) and (2b)-for example, through the function σ(r) defined in Sec.6.3.
In addition, if p r and p θ vanish sufficiently 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 effect 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.
Finally, for the equation of state (6.20) the metric is given by with m(r) defined in (2.8), as usual.

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].

Nonlocal ghost-free gravity
One of the simplest nonlocal higher-derivative models is defined by the form factor [7,[69][70][71][72]] for which where µ > 0 has dimension of mass and defines the nonlocality scale.In this case, performing the integration (1.9), the smeared effective delta source has a Gaussian profile [69], and the mass function is given by In Fig. 1 we show the behavior of the effective source and mass function, where the properties described in Sec. 5 can be explicitly verified.Given (7.4), the function A(r) in (2.7) reads In 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, x where The 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 where we used G = 1/M 2 P .Accordingly, 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 + 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) M, the position of the event horizon is bounded by the Schwarzschild radius, i.e., r + r s = 2GM.Firstly, let us analyze the framework of Sec.6.1, i.e., with the effective pressures p r = 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 suffices to consider the quantity because the curvature components (4.1) depend only on the derivatives of B(r).By 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, but similar results also hold for the other invariants, such as the Kretschmann scalar (4.2).The explicit expression for (7.10) is FIG. 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 0, which shows a significant difference 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. where 12) It is easy to verify that lim r→∞ R(r) = 0 and At 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 .In 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][70][71][72], for which the linearized Ricci scalar reads As 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 ≪ 1/µ and r ≫ 1/µ.Also, if M is sufficiently small, the two solutions do not deviate too much; we elaborate more on this in what follows.
As 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.66M 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.66M c < M < M c , the negative term in (7.12) becomes more relevant.It dominates for 0.95M 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.93M c .
Analogous 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 sufficiently 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.
As 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)]: ) It 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).Let 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, So, the metric is given by which 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 g tt = −[1 + 2ϕ(r)], with ϕ(r) given by (7.14).
In 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 effective 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 sufficiently small r, as explained in Sec.6.4.(However, curvature-derivative invariants calculated at the two solutions might differ, 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.
Finally, we point out that the metric (7.19) differs from the one obtained in Ref. [108] inspired by noncommutative geometry.Although both come from the same Gaussian effective source (7.3), the underlying equations of state are different, which leads to different shift functions B(r).In Appendix A, we carry out a more detailed comparison of the two solutions.) when M > M c .The scalar curvature is finite everywhere, even at the horizons r ± .

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, which corresponds to a sixth-derivative gravitational action.The effective delta source (1.9) for this model was calculated in [57] and it reads Ordering the masses such that m 1 < m 2 , it is easy to see that ρ eff (r) > 0 for all r.According to the discussion involving Eq. ( 5.13), it follows that the mass function, is monotone and also positive; this can be explicitly verified by checking that m ′ (r) > 0 and m(0) = 0.It 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, Since the UV behavior of the form factor does not depend on the multiplicity of the roots, the effective 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].The expression for A(r) is obtained from combining (2.7) and (7.22), and it can be cast in the form 16 where we defined the dimensionless function Since 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, max ) 16 In particular, in the degenerate case m 2 = m 1 we obtain  (7.30), which shows that it has only one solution for each q 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.
which leads to the conclusion that This yields a lower-bound estimate to the critical mass M c of Eq. (5.11), where 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.68M 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.
On the other hand, for M > M c , the solution has exactly two horizons.To prove this statement, it suffices 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 17for a given pair of masses m 1 and m 2 .If we define the dimensionless variables q ≡ m 2 m 1 and y ≡ m 1 r, (7.29) the equation A ′ (r) = 0 can be cast as 1 − q 2 + q 2 1 + y + y 2 e −y − 1 + qy(1 + qy) e −qy 1 − q 2 = 0. (7.30) The 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.Figure 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, Therefore, 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) monotonically 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 ≫ 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 1.79 m 1 r + 2GM.(7.32)Let us stress that the case of sixth-derivative gravity in which the masses m 1 m 2 in (7.20) are complex conjugate is completely different from the one studied here, for the effective 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].
In 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

B(r)
Combining (7.24) and (7.33) we obtain the explicit expression for the metric, It is easy to check that to the leading order in M, the above metric reduces to the Newtonian-limit form (1.5), where is precisely the modified Newtonian potential [53,67], in agreement with the general discussion of Sec. 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 2 with real distinct roots, i.e., with m i m j if i 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 effective source can be obtained by expanding the integrand of (1.9) in partial fractions, using where the coefficients C i are given by ) 19 In the degenerate scenario m 2 = m 1 this expression reduces to and satisfy the identities [56,115,116] Substituting (7.37) into (1.9) it follows [57] ρ Hence, the effective mass function (2.8) reads We did not encounter any effective source (7.40) that assumes negative values; therefore, we conjecture that ρ eff (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 different, and in that case ρ eff (r) can achieve negative values, with important implications for the black hole solutions [114].
The combination of (2.7) and (7.41) yields From 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), Since Z(x) is bounded by 0 Z(x) Z max = 0.298, A(r) is also bounded.Therefore, the same reasoning we adopted to estimate the mass gap in Sec.The complete metric, thus, reads (7.48) In agreement with the general discussion of Sec. 3, in the leading order in M this metric reduces to the Newtonian form (1.5) with being exactly the modified Newtonian potential evaluated in [53].

General polynomial gravity
Finally, let us consider the most general case of polynomial gravity, where we allow the mass parameters m i to be complex and/or degenerate.The form factor reads where we assume that the polynomial f (z) has degree N 2 and it has n 1 distinct roots m 2 1 , m 2 2 , . . ., m 2 n , each of which has a multiplicity α 1 , α 2 , . . ., α n , such that (7.51) Similarly to the previous example, we can use the partial fraction decomposition to write where the coefficients can be determined by the residue method (see, e.g., [117]) and are given by C i, j = (−1) α i − j m 2 j i (α i − j)! Calculating the effective delta source ρ eff from (1.9) using (7.52), we find [57] ρ where 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 effective source.It 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], where the coefficients A i, j are defined such that 20   1 z f (z) ) 20 Comparing to the notation used in [55], the correspondence between the coefficients A i, j (used there) and A i, j (used here) is through A i, j = ( j−1)!A i, j /m  field equations with this equation of state was obtained for effective sources originating from local and nonlocal higher-derivative gravity models, as summarized in Table I.
We remark that all the solutions found, regardless of the choice of the effective 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.
A 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).
Regarding the comparison with the nonrelativistic limit, we showed that the solutions related to the equation of state p r = [A(r) − 1] ρ eff match the Newtonian-limit solutions not only in the regime of large r but also for small r (if the effective source is regular).This means that the curvature invariants show the same leading behavior in these regions; however, higher-order derivative invariants may differ at r = 0.
It 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.(b) Ratio K/K Sch. .
As a concrete example, let us consider the case of the Gaussian effective source (7.3).The main choice of equation of state used in this paper yields the metric (7.19)

(A7)
As 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 In 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) differs 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.
or faster) for large k .

MFIG. 1 :
FIG.1: Plot of the dimensionless forms of the effective 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 effective source, and (5.5) and (5.6) for the mass function.Also, because ρ eff (r) takes only positive values, the mass function is monotonic and m(r) M; see the discussion around (5.13).

FIG. 5 :
FIG.5: Behavior of the Ricci scalar for the metric(7.19)when M > M c .The scalar curvature is finite everywhere, even at the horizons r ± .

FIG. 6 :
FIG.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.

FIG. 7 :
FIG.7: Numerical solution of Eq.(7.30), which shows that it has only one solution for each q 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.

7 . 2 |C i |m i . ( 7 . 45 )
leads to the conclusion that A(r) is positive provided that M is small enough such that 2GM μZ max < 1This yields the lower-bound estimate to the critical mass M c in Eq.(5.11), .(7.46) and (7.28) are essentially the same, the only difference being the mass scale μ.This is expected, for the previous section's model is the particular case N = 2 of the model considered here.Finally, the equation of state(6.20)  in the framework of Sec.6.4 results in the solution (6.21) for B(r), namely,B(r) = −2GM N i=1C i m i e −m i r .(7.47)

53 )
These coefficients satisfy identities that generalize the ones in(7.39), in particular,

TABLE 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]ρ eff .The function Z(x) is defined in Eq. (7.25), while the coefficients C i , 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.
, while the choice (A1) results in