Weak-field limit and regular solutions in polynomial higher-derivative gravities

In the present work we show that, in the linear regime, gravity theories with more than four derivatives can have remarkable regularity properties if compared to their fourth-order counterparts. To this end, we derive the expressions for the metric potentials associated to a pointlike mass in a general higher-order gravity model in the Newtonian limit. It is shown that any polynomial model with at least six derivatives in both spin-2 and spin-0 sectors has regular curvature invariants. We also discuss the dynamical problem of the collapse of a small mass, considered as a spherical superposition of nonspinning gyratons. Similarly to the static case, for models with more than four derivatives the Kretschmann invariant is regular during the collapse of a thick null shell. We also verify the existence of the mass gap for the formation of mini black holes even if complex and/or degenerate poles are allowed, generalizing previous considerations on the subject and covering the case of Lee-Wick gravity. These interesting regularity properties of sixth- and higher-derivative models at the linear level reinforce the question of whether there can be nonsingular black holes in the full nonlinear model.


Introduction
The problem of quantizing gravity is a long-standing one and many conceptually different approaches have been used to tackle it during the last almost five decades. One of the elements considered in the path towards quantum gravity are the higher derivatives and the role they play in the ultraviolet (UV) regime, where classical and quantum singularities show up. Motivations for the introduction of curvature-squared terms in the action come already at semiclassical level, from the observation that the renormalization of quantum field theory on curved background requires such higher-derivative terms [1] (see also [2,3] for a review). Furthermore, even though general relativity (GR) is not perturbatively renormalizable, its fourth-order counterpart is [4]. Increasing the number of derivatives in the action can make the theory even more regular. For example, in the local theories with more than four derivatives it is possible to achieve superrenormalizability [5]. Indeed, the models with six derivatives have divergences only up to 3-loops, while in those with more than eight derivatives only 1-loop divergences remain.
The benefits that higher derivatives bring in what concerns renormalization, however, come together with a serious drawback regarding unitarity. Although it is possible to associate the new degrees of freedom of the theory to positive-norm states in the Hilbert space, some of them may carry negative energy [4,5]. These so-called ghost states introduce instabilities in the theory, with the possibility of a boundless vacuum decay via the emission of an arbitrary amount of energy in the form of gravitons. In such a scenario it makes sense to study classical and quantum aspects of models which can offer insight on how to deal with, e.g., the tension between renormalizability and unitarity, or the most appropriate form of treating (or avoiding) ghosts and related instabilities 1 .
In this regard, two models have been the subject of interesting investigations in recent years. The first one we mention is the Lee-Wick gravity [11,12]-see, e.g., [10,[13][14][15][16] for further developments and applications. This theory is defined by the Einstein-Hilbert action enlarged by curvature-squared terms which contain polynomial functions of the d'Alembert operator, such as R µν F 1 (✷)R µν and RF 2 (✷)R. A general action of this type can be called polynomial higher-derivative gravity and was introduced in [5]; the Lee-Wick gravity assumes, furthermore, that the polynomials F i are such that all the massive poles of the propagator which correspond to ghost modes are complex. Hence, the physical spectrum of the theory contains the usual massless graviton and, possibly, a healthy massive scalar particle (as the lightest scalar excitation is not a ghost [5]). The pairs of complex conjugate massive modes are understood as virtual particles only and should decay to healthy particles. It was claimed that the presence of those complex poles do not violate the unitarity of the S-matrix if the Lee-Wick quantization prescription is used [11,12]. Therefore, this could be a form of restoring the unitarity, weakening the tension between renormalizability and unitarity.
Another proposal for dealing with the problem of ghosts is to avoid them, at least at tree level, by replacing the polynomials F i of the action by nonpolynomial functions of the d'Alembertian, which makes the theory nonlocal [17][18][19][20] (see also the earlier works [21,22]). It is possible to choose these functions in such a manner that the theory propagator contains only the pole at k 2 = 0, which corresponds to the graviton 2 . Owed to the absence of ghosts, this theory is sometimes called ghost-free gravity. As pointed out in [43], however, quantum corrections may prompt the emergence of an infinite amount of complex ghost poles. Therefore, the study of the Lee-Wick gravity may be useful also to the better understanding and development of nonlocal UV extensions of GR.
The present work revisits two topics that have previously been investigated in the context of local and nonlocal higher-derivative gravity models, namely, the nonrelativistic limit [10,20,[44][45][46][47][48] and the collapse of small mass spherical shells [49,50]. Our focus is on general polynomial gravity, with a special attention given to the case of complex poles-Lee-Wick gravity-and also, for the sake of completeness, higher-order (degenerate) poles. In this sense, the results presented here both generalize and refine previous considerations on the aforementioned topics, as we describe in what follows.
The presence of higher derivatives in the gravitational action tends to ameliorate both classical and quantum divergences. The former can be viewed, e.g., on the Newtonian potential and on the effect of the gravitational collapse. The latter is related, as mentioned before, to the (super)renormalizability of the theory. Since 1977 it is known that the fourth-derivative gravity is renormalizable and has finite nonrelativistic potentials [4,51]. This relation was recently extended to superrenormalizable higher-order gravity theories with real poles, which were shown to have a finite potential too [45]. On the other hand, the introduction of higher derivatives only in the R 2 -sector of the theory results in a nonrenormalizable model with a divergent (modified) Newtonian potential [44]. These examples of simultaneous occurrence of classical and quantum singularities raised the question of whether there is a fundamental relation between them [43,45,52]. The negative to this conjecture was given in [46], where it was shown that the Newtonian singularity is canceled in all the polynomial gravity theories with at least one massive mode in each sector, which included Lee-Wick and also some nonrenormalizable models.
Nevertheless, the proof of the finiteness of the potential carried out in [46] was based on the calculation only of the terms which give divergent contributions to the potential, and on the demonstration of an algebraic relation between the poles of the propagator of the theory. In the Sec. 2 of the present work we derive the expression for the weak-field metric potentials to all orders in r-including the case of degenerate poles-and obtain an alternative verification of the cancellation of the Newtonian singularity. This simpler demonstration is based on partial fraction decomposition and on the use of the heat kernel method for deriving gravitational potentials introduced in [50].
Having the expression of the linearized metric for a pointlike source in a general local higher-derivative gravity, it is possible to go beyond the analysis of the finiteness of the potential and discuss the regularity of the curvature invariants. This is carried out also in Sec. 2, where we show that the metric is regular only if the model contains more than four derivatives in both the scalar and tensor sectors. This includes local superrenormalizable and Lee-Wick gravity models. Following the aforementioned parallel between quantum and classical singularities, one can say that GR is nonrenormalizable and has divergent Newtonian potential, fourth-order gravity is renormalizable and has a finite gravitational potential (but their curvature invariants in the weak-field approximation diverge), and the higher-order gravities which are superrenormalizable have a complete regular nonrelativistic limit, i.e., the metric potentials and the curvatures have no singularities.
In the Sec. 3, the static solution found in the preceding section is used to obtain the metric associated to a nonspinning gyraton, which is an approximation to an ultrarelativistic massive particle without angular momentum [53][54][55]. The procedure comprise applying a boost to the nonrelativistic metric and then taking the Penrose limit (see, e.g., [56]). This turns out be an intermediate step to the analysis of the collapsing null shells, which is discussed in Secs. 4 and 5. The field generated by a nonspinning gyraton was derived in the context of the nonlocal ghost-free gravity in [49], and for the polynomial gravity with simple poles in [50]. In the present work we show that this metric has the same small-distance behavior in all nontrivial polynomial gravity theories. To conclude this section, some particular explicit examples are presented for the cases of complex and degenerate poles.
The collapse of small mass shells is analyzed in Secs. 4 and 5 which, respectively, discusses the case of an infinitesimally thin shell and of a shell with a finite thickness. By small mass we mean that we work only with linearized equations for the gravitational field, in agreement to what was developed in the previous sections. The interest in this scenario is the possibility of formation of mini black holes, e.g., owed to the collision of ultrarelativistic particles [48]. The formalism we follow was introduced in detail in [49], where it was applied to the ghost-free gravity. It was later generalized in [50], where the case of polynomial models with simple poles was considered. Our extension of the latter work to general polynomial models verifies the conclusion that there exists a mass gap to the formation of mini black holes. The presence of a mass gap is typical in higherderivative gravity models, which is known since the 1980s [57], and means that a black hole can only be formed if its mass is larger than a certain value. This is in contrast to what happens in GR, where any mass can become a black hole, provided it is concentrated in a sufficiently small region.
Also, in Secs. 4 and 5 we discuss the emergence of singularities during the collapse of null shells within general polynomial gravities by analyzing the Kretschmann scalar R 2 µναβ . In particular, in Sec. 5 we show that the Kretschmann scalar for a collapsing thick null shell is regular for all models with more than four derivatives in the spin-2 sector. This characterizes the class of models for which R 2 µναβ can have the logarithmic singularities found in [50]. Actually, only the fourth-derivative model has such divergence. Further discussion concerning parallels between local and nonlocal higher-derivative gravity is carried out in Sec. 6, where we also draw our conclusions.

Newtonian limit
In the static weak-field approximation we consider metric fluctuations around Minkowski spacetime and work with the equations of motion at the linear level. The only relevant terms in the action which contribute for the linearized field equations are those of second order in the perturbation h µν . Consequently, at the Newtonian limit a general higher-derivative gravity model can be reduced to the action where κ = 8πG and F 1 and F 2 are functions of the d'Alembert operator. If F 1 and F 2 are nonzero polynomial functions, not necessarily of the same degree, we say it is a polynomial higher-derivative model 3 . Otherwise, we say the theory is nonlocal. Let us note that the term R µναβ F 3 (✷)R µναβ is irrelevant for our purposes since, by means of the Bianchi identities and integrations by parts, one can prove that (see, e.g., [5]) Hence, the effect of such Riemann-squared term can be reproduced, at the linear level, by a redefinition of the functions F 1 and F 2 .
Performing the expansion (1), the bilinear form of the action (2) reads where we introduced the condensate notations and The variational principle then yields the field equations, where T µν is the energy-momentum tensor sourcing the field. As far as we are interested in a pointlike source, we assume where ρ = m δ 3 (r) is the mass density. In this case the metric can be written in the isotropic form Here ϕ = ϕ(r) and ψ = ψ(r) are the Newtonian potentials and r = x 2 + y 2 + z 2 . The metric potentials can be obtained by solving which are, respectively, the 00-component and the trace of the equations of motion (7). Moreover, the substitution ✷ → ∆ was implemented, as the metric is static. Instead of solving the system above directly for ϕ and ψ, it is more convenient to work with their linear combination in the form of Once the equations are solved for χ and ω it is straightforward to obtain the original metric potentials via The reason for working with χ and ω is threefold: first, the field equations for these new potentials have a simple structure in terms of the functions a and c. In fact, Eqs. (10) and (11) are equivalent to a(∆)∆χ = κρ , Second, the functions a and b above correspond precisely to the terms which appear, respectively, in the sectors of spin-2 and spin-0 of the propagator associated to the theory (2) [58], where P (2) and P (0−s) are, respectively, the spin-2 and spin-0 projection operators (see, e.g., [2]; tensorial indices and the terms which are gauge-dependent were omitted for simplicity). Indeed, the roots of the equations a(−k 2 ) = 0 and b(−k 2 ) = 0 determine the massive poles of the propagator and, therefore, the (massive) spectrum of the model. In this spirit, Eq. (12) splits the metric potentials into the contributions owed to the spin-2 modes (through χ) and to the scalar modes (via ω). Based on this relation between the roots of the equations a(−k 2 ) = b(−k 2 ) = 0 and the poles of the propagator, throughout the present work we shall refer to these quantities as either "roots" or "poles". The third motivation for working with the special combination in the form χ = ϕ + ψ is that the potential χ turns out to be the relevant one for the collapse of the spherical null shell (see discussion in Sec. 3 and also in Ref. [50]). The situation resembles what occurs in the light bending [59]. Qualitatively, this happens because in the ultrarelativistic limit the interaction between particles and the gravitational field is similar to that of photons.

Heat kernel solution
Equations (14) and (15) have the very same structure, the only difference being the operator function. From now on we assume that F 1 and F 2 are polynomial functions, as our interest in this work is on higher-derivative polynomial gravity. Then, a, b and c are also polynomials and the equations for χ and ω are essentially the same, the only difference being the coefficients of the polynomials. Therefore, we explicitly work out the solution for (14) and, mutatis mutandis, write down the solution for (15). The solution for χ can be easily found by means of the heat kernel approach, based on the Laplace transformation, as carried out in [50].
Indeed, introducing the Green's function for the Eq. (14) we have the integral solution Let us now assume that the inverseĤ −1 (∆) of operator (19) can be written as the Laplace transform of some function f (s), that is, Then, the x-representation of the Green's functionĜ reduces to where is the heat kernel of the Laplacian. By choosing x = r and x ′ = 0 , formula (20) simplifies to Particularizing for the higher-derivative gravity model (2), according to the fundamental theorem of algebra we can write the polynomial function a(−ξ) in the factored form 5 is a root of the equation a(−ξ) = 0 and α i is its multiplicity. Notice that if N is the degree of a(∆)-i.e., if there are 2(N + 1) derivatives in the spin-2 sector-then N i=1 α i = N . Since our interest is in general polynomial models, we shall not make any initial restriction on the complex or real nature of the quantities m 2 i , nor on their multiplicity. The function f (s) for the general higher-derivative gravity can be promptly obtained by substituting (25) into (21) and inverting the Laplace transform using expansion in partial fractions [60]. The result is where the coefficients A i,j are obtained from the comparison with the partial fraction decomposition of H −1 (−ξ), namely, Also, for compactness of notation, it is useful to define the symbol A i,j for j > α i by The potential χ can thus be evaluated by substituting (26) into (24), which gives where we assume that Re m 2 i > 0 for the integrals to converge. Under the change of variables sm 2 i → s each of the above integrals becomes with the last integral being carried out along the line Γ = {w ∈ C : w = m 2 i t, t ∈ R + }. In the case of a real root m 2 i the integration remains along the positive real axis, while for complex roots the integration line undergoes a rotation in the complex plane, but its points still satisfy Re w > 0. However, the integrand h(s) on the r.h.s. of (29) is an analytical function with only a removable singularity at the origin, and which vanishes for |s| → ∞. Therefore, the integral of h(s) along the oriented contour We conclude that even in presence of complex roots m 2 i it is possible to perform the integration along the positive real axis. Then, where we chose the square root of m 2 i with positive real part and recognized in the integral a representation of the modified Bessel function of the second kind K ν [60]. Hence, the potential χ is given by In deriving this result it was assumed that Re m 2 i > 0 and Re m i > 0. The last assumption is physically justified by the requirement that the potential decays to zero at large distances, as well as to avoid tachyons on the model. The former assumption, however, is related to the heat kernel method used to solve (14) and the premise that the operatorĤ −1 has the form of (21). Actually, the solution (31) also holds for the cases in which the polynomial a(−ξ) has roots with | Im m i | > Re m i > 0, as the Bessel functions provide the analytical continuation of each term in (28) viewed as function of an arbitrary m 2 i with Re m i > 0. We point out that it is possible to obtain the potential (31), even though with a longer calculation, directly for the general case of | arg m i | < π/2 by means of the Fourier transform method and using the Basset's representation of the modified Bessel functions [61].
The case of GR (a ≡ 1) is a trivial example of the previous formulas, as f GR (s) = −1 and χ GR (r) = −2Gmr −1 . Another direct example is if a(−ξ) = 0 has only nondegenerate (ND) roots. Then α i = 1 for all i, and f (s) boils down to [50] while the potential is given by [50] Since the only assumption in finding the solution for χ was that it satisfied Eq. (14), one can write down the solution for ω which satisfies (15). Let N ′ be the degree of the polynomial b(−ξ) and let −m ′2 i (with i ∈ {1, 2, ..., N ′ }) be the roots of the equation b(−ξ) = 0, each of them with multiplicity α ′ i . Then, the formula for ω(r) can be obtained by simply making the substitution (χ, (27) and (31).
In view of (13), the modified Newtonian potential ϕ for a general higher-derivative polynomial gravity is given by while ψ reads As noted before, the quantities m i are the masses of the extra degrees of freedom with spin-2, while m ′ i are related to the scalar ones. Moreover, the potentials are real despite the possibility of complex poles in the propagator. The cancellation of the imaginary part takes place because K n (z) = K n (z) for n ∈ R, and Aī ,j = A i,j , where the subscript index i refers to the complex pole conjugate to m 2 i . The general potential (34) generalizes previous considerations found in the literature which took into account real massive poles only in the scalar sector [44], or simple real poles [45] and simple complex poles [46,50] in scalar and tensor sectors. As noticed in [10,44,46], it is possible to obtain the potential for the case of degenerate poles by considering limits of the potential with only simple poles. However, this procedure may be ambiguous when applied to poles with α i > 2. The formula (34) clarifies the situation, as it explicitly allows for arbitrary multiplicity.

Finiteness of the metric potentials
If both χ and ω are finite, so are the metric potentials ϕ and ψ. As noticed in [50], if the roots of a(−ξ) = 0 are all simple, then χ is finite. In what follows we use the general formula (31) to show that χ is finite for an arbitrary nontrivial polynomial a of the form (5). Using the similarity between the solution for χ and ω, it then follows that these properties are valid also for ω defined by a nontrivial b given by (16). As a conclusion, if a and b have degree of at least one, then the potentials ϕ and ψ are finite at r = 0. This can be viewed as an explicit verification of the result obtained in [46], where only the terms of order r −1 were evaluated and the presence of degenerate poles was dealt with by the procedure of taking limits.
To this end, let us rewrite (31) separating the terms for which j > 3/2: where the summation over j ≥ 2 is considered only if α i > 1. For j ≥ 2 and small r the functions K j− 3 2 (m i r) behave like r −j+3/2 . Hence, all the terms with j ≥ 2 are finite at r = 0. It remains to check if the terms with j = 1 manage to cancel the Newtonian singularity. Since the terms with j = 1 have the form Therefore, the potential (31) can be written as where χ 0 and χ 1 are constants. Using the identity i A i,1 = 1 (see Eq. (123) of the Appendix A), it follows that the Newtonian singularity at r = 0 is canceled by the higher-derivative correction terms, even in presence of complex and/or degenerate poles.
The same reasoning holds for the potential ω, and therefore the metric potentials ϕ and ψ are finite, verifying the result of [46]. The condition for the cancellation of the singularity of the potential is the presence of at least one massive mode in the spin-2 and in the spin-0 sectors. For example, if F 1 = 0 but F 2 = 0 then ω is finite but ϕ and ψ are not [44].

Regularity of the curvature invariants
As it is well known, the finiteness of the potential is not enough to guarantee the regularity of the solution, as the curvature can still be singular. For a general metric in the form (9), e.g., the Kretschmann invariant is given by which clearly diverges if ϕ ′ (0) and ψ ′ (0) are not zero.
In order to find more rigorously the conditions for having regular curvature invariants, let us assume that both metric potentials are finite and write In terms of ϕ n and ψ n the Kretschmann scalar reads Therefore, the invariant R 2 µναβ is regular if, and only if, ϕ 1 = ψ 1 = 0. Actually, this is the same condition for the regularity of the set of curvature invariants: 6 where C µναβ is the Weyl tensor and χ n = ϕ n + ψ n , ω n = ϕ n − 2ψ n , n = 0, 1, 2, . . . .
In this spirit, one may be tempted to ask whether the condition ϕ 1 = ψ 1 = 0 is common in higher-derivative gravity models. For example, it is known that the fourth-order model does not satisfy this condition when coupled to a δ-source [51] (see also [64][65][66] for a more recent discussion). However, for the sixth-order gravity with a pair of complex poles this condition holds [14], and general considerations supporting the conjecture that for theories with more than four derivatives the condition ϕ 1 = ψ 1 = 0 should hold was given in [64]. We here address a more direct answer to this question by explicitly showing the sufficient conditions for having a regular metric for linearized higher-derivative gravity coupled to a δ-source.
To this end, let us extend to order r our previous calculations which showed that the potentials χ and ω are finite in higher-derivative gravity. Using the general expression for the potential (36) and the series expansion of the modified Bessel functions [60] it is not difficult to verify that the terms which contribute to order r yield But the term in brackets inside the summation over j ≥ 3 is Thus, where we define In the Appendix A we show that if the polynomial a(−ξ) is of degree N > 1, then S 1 = S 2 (see Eq. (126))-recall that 2(N + 1) is the number of derivatives in the spin-2 sector of the action. It follows that for theories of order higher than four, the nonrelativistic potential χ is not only finite, but it is also regular, i.e., χ 1 = 0. On the other hand, for the case of N = 1 with a root at ξ = −m 2 1 one has the trivial result A 1 = 1 and S 2 = 0, which gives χ 1 = Gmm 2 1 . This reasoning can be immediately extended to the potential ω, for which We conclude that the condition for the regularity of the curvature invariants 7 is the presence of at least two massive modes (or one degenerate pole) in each of the spin-2 and the spin-0 sectors-which is equivalent to having a and b of degree higher than one 8 . In other words, all higher-derivative theories defined by nonconstant polynomials F 2 and F 1 = 3F 2 are regular at the Newtonian limit. In particular, this holds for the superrenormalizable local higher-derivative gravity models, including Lee-Wick models.
In this context, the only possibilities for having a nonregular solution for a point source in the Newtonian limit is to have F 1 (✷) = const. or F 1 (✷) = 3F 2 (✷). In the first case the spin-2 sector contains the massless pole corresponding to the graviton and, possibly, one massive (ghost) particle. In terms of the definition of π-regularity [50], which means that π 1 = 0 for a metric potential π(r), we can say such a solution is not χ-regular, but it could be ω-regular provided that F 2 (✷) ∼ ✷ p with p ≥ 1. For the second case, i.e., if F 1 (✷) = 3F 2 (✷), the solution is not ω-regular. Of course, for the solution to be regular it must be both χ-and ω-regular. In particular, Stelle's fourth order gravity is not regular at the Newtonian limit when coupled to a δ-source [51,[64][65][66], even though it can be ϕ-regular for particular choices of parameters, namely, if m ′ 1 = 2m 1 .

Small-r conformally flat solutions
In view of the Eq. (46), it follows that a χ-regular solution yields C 2 µναβ = 0 at r = 0. The components of the Weyl tensor read 7 It is also possible to verify that under these conditions all individual components of the curvature tensors remain finite. 8 We point out that the effect of the regularization of the curvature can be viewed also in the polynomial theories as a regularization of the source in the Poisson equation for the metric potentials [62].
As we can see, most of them do not contain terms with r −1 powers, which implies that if the potentials are finite in the origin, the same is true for the corresponding components. The exception is the only independent component, C trtr . In fact, Thus, we conclude that this component is finite in r = 0 only for χ-regular theories. In such a case, χ 1 = 0 and the components of the Weyl tensor tend to zero as r → 0, which means that the metric is approximately conformally flat near the origin.

Ultrarelativistic limit
Up to this point we restricted considerations to the Newtonian limit. In the following sections the weak-field potential χ will be used to discuss the emergence of a singularity in the collapse of null shells. As a first step towards the gravitational field of a collapsing shell, we shall obtain the field associated to an ultrarelativistic point-particle, which may be done by the following procedure. First, we perform a Lorentz transformation into Eq. (9), which yields the metric of a moving object with velocity β. Thereafter, we take the limit β → 1 while keeping the relativistic mass of the object fixed (Penrose limit), i.e., being M the mass of the ultrarelativistic particle and γ = (1 −β 2 ) −1/2 the Lorentz factor. The resultant metric corresponds to a nonspinning gyraton [56]. In order to apply this scheme to the solution found in the previous section, let us rewrite the metric (9) in the form where is the flat spacetime metric and is the perturbation. Now, consider a boost in the x-direction, Introducing the null coordinates v = t ′ + x ′ and u = t ′ − x ′ , Eqs. (61) and (62) read Therefore, after applying the boost to the metric (58) one gets and In the limit β → 1 the form of flat metric (65) remains the same, while the perturbation goes to This shows, as mentioned before, that the dominant contribution in the ultrarelativistic limit comes from the special combination of the metric potentials in the form of χ = ϕ+ψ.
Owed to this fact, in this section and in Secs. 4 and 5 we restrict considerations to the spin-2 sector of the theory. In this spirit, when we refer to, e.g., "theories with more than four derivatives" is must be understood that these derivatives are on the spin-2 sector. The function Φ can be evaluated through (67) by combining Eqs. (24) and (57) and recalling that Indeed, taking into account that r 2 = γ 2 u 2 + y 2 + z 2 after the boost, it follows which can be written as where we defined the function F : R → R via The integral (71) typically has an infrared divergence, owed to the massless nature of the graviton. To overcome this problem one can introduce an infrared cutoff Ω for large s. Any change in the cutoff parameter can be absorbed into a redefinition of the coordinates. In other words, this ambiguity just reflects the freedom in the gauge choice. Quantities with classical physical meaning, such as the curvature tensors, do not depend on Ω (for a more detailed exposition see, e.g., [49]). For example, in the case of GR f (s) = −1, so that Here E 1 (z) is the exponential integral function. As Ω is a huge arbitrary cutoff, we assume z ≪ Ω 2 and write where γ is the Euler-Mascheroni constant, and terms of order z/Ω 2 and higher were discarded. For the general higher-derivative model (2) the function f (s) is given by Eq. (26), which yields By applying the same arguments used in Sec. 2.1 it is possible to express the function F Ω in terms of modified Bessel functions of the second kind, Before we present some explicit calculations for the cases of complex and degenerate poles, let us show a general property of the function F Ω (z) within higher-derivative gravity, based on its definition in (75). On the one hand, Eq. (73) shows that in GR F GR Ω (z) ∼ ln z diverges as z → 0. On the other hand, in [50] it was shown that this divergence do not occur in the case of polynomial gravity with simple poles, because the leading terms of F Ω (z) for small z are linear in z or of the type z ln z. Now we prove that this feature is present also in the general polynomial theory. Indeed, for small arguments the modified Bessel functions of the second kind K n (z) (n ∈ N) can be expanded as where c i are constants and γ is the Euler-Mascheroni constant. Substituting these expressions in (75) and using (123) it follows (c ′ is a new constant) The constants S n are defined just like in (52), while S is given by (80) Note that in any higher-derivative gravity model the singular term ln z which stems in GR (see (73)) is canceled by a specific combination of the contribution owed to each massive mode through K 0 (m i √ z). This is a direct consequence of the cancellation of the Newtonian singularity discussed in Sec. 2.2 and in Refs. [45,46]. It is also useful to point out that, while the constant S ′ 1 is nontrivial for all higher-derivative polynomial models, the quantities S 2 and S ′ 2 only appear if there is at least one pole with multiplicity equal or larger than 2, and P 3 is relevant only for models with at least one pole for which α i ≥ 3-this justifies our choice for the subscript labels.

Particular cases and examples
To close this section let us consider some examples of the diversity of scenarios which occur in higher-derivative gravity. In particular, we present explicit calculations for the sixth-order gravity, which is the simplest model which admits complex poles or real degenerate poles. We shall return to these examples in the next section, when analyzing the gravitational field of collapsing null shells.

Fourth-order gravity
In the fourth-order gravity there is only one possible scenario: the equation a(−ξ) = 0 has one real simple root at ξ = −m 2 1 . Therefore, S 1 = m 2 1 and S = m 2 1 ln m 2 1 , so that As the other examples will show, and in consonance with the discussion in Sec. 2.3, this is the only case in which the small-z expansion of F (z) contains the term z ln z.

Models with more than four derivatives
For any model of order higher than four there is the identity S 1 = S 2 (see Eq. (126) of the Appendix A). Hence, Eq. (79) can be cast in a very simple form: This result is both a generalization and a simplification of the analogous expression derived in [50], as it accounts for the possibility of degenerate poles and also rules out the terms of the type z ln z-which turn out to be possible only in the fourth-derivative gravity.

Non-degenerate models
The case of nondegenerate roots was investigated in Ref. [50]. Here we show that our general considerations correctly reproduces this particular case. If all the roots of a(−ξ) = 0 are simple, then α i = 1 ∀ i and the general expression (75) for F (z) reduces to [50] Now, let us assume that N > 1 (the case of N = 1 was discussed in the Example 3.1.1). Inasmuch as all the roots are nondegenerate, it follows that S 2 = S ′ 2 = P 3 = 0, whence S = S ′ 1 . Therefore, for small z the function F (z) behaves like

Maximally degenerate models
We say the higher-derivative model of order N > 1 is maximally degenerate if the equation a(−ξ) = 0 has only one root at ξ = −m 2 1 , with multiplicity N . In such a case, the following relations are valid: Thus, for small z the function F (z) can be written as

6th-order gravity with simple poles
For a pair of simple poles m 2 1 and m 2 2 , Eqs. (52) and (80) yield S 2 = 0 and S ′ 1 = 0. If these poles are simple and real the function f (s) is given by which yields, for small z, In the case of two conjugate simple complex roots with m 1 = α + iβ and m 2 = α − iβ, it follows and 3.1.6 6th-order gravity with degenerate poles and As the particular case of the maximally degenerate model with N = 2, it holds S ′ 1 = S ′ 2 = m 2 1 ln m 2 1 , which gives S = 0. We note that Eq. (93) can be obtained from the analogous equations for simple poles by taking the limit m 2 → m 1 in (88)-or the limit β → 0 in (90). While this procedure of taking the limit is simple to carry out in the case of two roots (see, e.g., [10] for more examples), the situation might be not so clear if one is to consider a higher-order root. In such a case it is preferred to work with the general formula (75), or (82), as discussed in Sec. 2.

Thin null shell collapse
In this section we analyze the collapse of a null shell and the formation of mini black holes. Following Refs. [49,50], we first consider a shell with vanishing thickness (thin null shell). For this case the Kretschmann curvature invariant is still singular, but this singularity is consequence of the nonphysical approximation of a infinitesimally thin shell.
The field associated to a null shell with vanishing thickness (also called δ-shell) can be obtained, at the linearized level, by the superposition of an infinite amount of gyratons spherically distributed and which pass through one given point O [49], which we take as the origin of the coordinate system. This point is the vertex of the null cone representing the shell, so that for t < 0 the shell is collapsing towards the apex O and for t > 0 the shell proceeds its expansion after the collapse. It can be shown that, outside the shell, the averaged metric perturbation dh 2 resulting from this distribution of nonspinning gyratons is given by (see [49] for a detailed derivation of this result) where we use spherical coordinates, so that dΩ 2 = dθ 2 + sin θ 2 dφ 2 is the metric of the unit sphere and is the complete metric. Here F (z) is defined by Eq. (71), as given by the metric (67) associated to a single gyraton.

Apparent horizon
The formation of black holes is closely related to the invariant where f = ̺ 2 ≡ g θθ . Indeed, the points for which g = 0 correspond to an apparent horizon [67]. If it happens that g(t, r) is strictly positive then the collapsing shell generates no apparent horizon. For the general metric (95) the invariant g is given by [50] where If there is a positive constant C such that then g is positive anywhere provided that M < (2GC) −1 . Therefore, in order to show the existence of a mass gap to the formation of mini black holes in general higher-derivative gravity it suffices to verify that the function r −1 q(r 2 − t 2 ) is bounded. In Ref. [50] it was shown that for nondegenerate models there is the mass gap. In what follows we extend this result to the general polynomial model. For F (z) given by Eq. (75) we have As a finite sum of continuous functions defined for all z ∈ R + , q(z) is also continuous. Hence, if q(z) has any divergence it can only take place for large or small z. The former divergence does not occur, because the functions K j (z) exponentially decay as |z| → ∞, in such a way that q(z) → 1 as z → ∞. On the other hand, assuming N > 1, for small arguments one has whence q(z) → 0 as z → 0. Being the asymptotic limits finite, it follows that q(z) is bounded. Now let us analyze the function r −1 q(r 2 − t 2 ). The function r −1 is continuous, it vanishes for large r and only diverges as r → 0. In this regime, however, q(r 2 − t 2 ) dominates over r −1 , since |t| < r outside the shell implies in r 2 − t 2 < r 2 . Thus, Similar analysis can be applied for the case N = 1, with the same result [50]. We conclude that r −1 q(r 2 − t 2 ) is bounded for general higher-derivative polynomial gravity models, which implies in the existence of the mass gap for the formation of mini black holes. The size of the gap depends on the scale λ = max i {m −1 i } defined by the massive excitations of the model; such scale could be affected by a gravitational seesaw-like mechanism as discussed in [13] (see also [10,63] for experimental bounds on this scale).
To give an example of an explicit calculation, consider the sixth-order gravity with degenerate poles discussed in Sec. 3.1.6. The associated function q(z) is given by Following [50] we put β 2 ≡ 1 − t 2 r −2 and v ≡ m 1 βr, so that r −1 q(r 2 − t 2 ) = m 1 βV (v), with The function V (v) is positive and reaches its maximum of about 0.249 at v ≈ 2.324 (see Fig. 1). Thus, where we recalled that outside the shell the parameter β ranges in the interval (0, 1). Therefore, if M 2(Gm 1 ) −1 the collapse does not result in a black hole.

Kretschmann scalar
Even though there is a mass gap for the mini black hole formation in the general higher-derivative gravity, the Kretschmann curvature invariant R 2 µναβ is not regular at r = 0. Indeed, for a metric in the form (95) it is given by [49,50] and primes denote differentiation with respect to the argument z. For N > 1 and small arguments q(z) is given by (101), yielding where β 2 ≡ 1 − t 2 r −2 ranges between 0 and 1 outside the shell. As the collapse proceeds, the Kretschmann scalar diverges 9 for r → 0. This very same behavior occurs in the nonlocal ghost-free gravity [49], and it was previously verified to occur also in the particular case of polynomial models with simple poles [50]. Actually, in view of these results on similar models it is natural to expect the nonregularity of the Kretschmann invariant, as in these cases the function F (z) has the same linear dependence on z for small arguments. As pointed out in [49,50] this singularity of R 2 µναβ is associated to the nonphysical assumption of infinitesimally thin shell. The physical imploding shell must have finite thickness, which tends to regularize the curvature. In order to check this assumption, in the next section we analyze the case of the collapsing thick null shell.

Thick null shell collapse
In the linear regime one can build the metric associated to thick null shell by superposing a set of thin shells collapsing to the same spatial point O, which we take as origin of the coordinate system. Of course, there are infinite possibilities of distributing the total energy of the shell throughout its thickness. Since our goal here is to show that a nonsingular source regularizes the Kretschmann scalar in the polynomial gravity (and ameliorates the divergence for the fourth-order model), we choose the most simple density function by assuming that the density ρ(t) at r = 0 remains constant during all the collapse, being null before/after it. Such a definition of the energy flux passing at O suffices to determine the density profile of the shell, insomuch as each element of the fluid moves at the speed of the light and no self interaction is considered inside the shell. Therefore, for a shell with total mass M and thickness (or duration) τ , where we set t = 0 as the moment when half of the total mass crosses O. The corresponding metric perturbation can be obtained by averaging the metric (94) of the thin null shells with respect to the density ρ [49] The collapse of a thick null shell define specific spacetime domains (see, e.g., [49]). In the present work we restrict considerations to the domain I near t = r = 0, where (and when) the shell assumes its highest density-favoring the mini black hole formation and the emergence of singularities. This domain is characterized by the intersection of the in-coming and the out-coming fluxes of null fluid, and is formally defined by the locus of the spacetime points for which r + |t| < τ /2. Moreover, the metric is stationary inside I, for the energy density is constant. Taking into account that only the δ-layers which cross O at times t ′ ∈ (t − r, t + r) contribute to the field inside this domain, it is not difficult to verify that Eq. (111) yields [49] where we defined Particularizing this solution for gravity models with six or more derivatives in the action, we substitute the expression (82) for F (z) around z = 0. It follows The Kretschmann scalar associated to this solution is which is regular r = 0, as anticipated. It is worthwhile to mention that the nonsingularity of the source is not enough, by itself, to guarantee the regularity of the curvature. In fact, in general relativity F (z) ∼ ln z, which gives R 2 µναβ ∼ r −4 for the collapsing thick null shell. Also, the presence of the term z ln z in the small-z expansion of F (z) could yield logarithmic divergences in the Kretschmann scalar. Such singularity was considered in [50] as a possibility for general higher-derivative polynomial gravity (see Eq. (79)). Nonetheless, it only occurs for the fourth-order model, since for nontrivial polynomial theories there is the relation S 1 = S 2 which regularizes the nonrelativistic potential χ.
Explicitly, the Kretschmann scalar for a collapsing thick null shell in the fourthderivative gravity follows from (81) and reads [50] with c ≡ 2γ − 2 + ln m 2 1 . The origin of this singularity can be traced back to the nonrelativistic limit. Indeed, in [50] it was shown that, for polynomial theories with simple poles, the nonregularity of the potential χ implied in a singular Kretschmann scalar for the collapsing thick null shell.
We have seen that the divergences are softened when a δ-shell is substituted by a thick shell. It is therefore natural to expect the existence of a mass gap to the formation of mini black holes for a collapsing thick null shell too. For the sake of completeness, we calculate the invariant g(r) on the the domain I for the solution (114), which reads Since r < τ on I, it follows that Hence, for a given τ it is also possible to avoid the existence of an apparent horizon inside I provided that the mass M is sufficiently small.

Summary and discussion
Let us summarize the results obtained. We derived the solutions for the Newtonian potentials associated to a pointlike mass in a general polynomial higher-derivative gravity, i.e., allowing the presence of complex and degenerate poles (with arbitrary order) on the propagator. This includes the class of (super)renormalizable theories and, in particular, Lee-Wick gravity models. It was verified, in agreement to [46], that the metric potentials remains finite in r = 0 provided that there is at least one massive mode in each spin-2 and spin-0 sectors. This is not a sufficient condition, however, to ensure the regularity of the solution, because there can be singularities in the curvatures.
Indeed, since the 1970s it is known that Stelle's fourth-order gravity possesses curvature singularities in the linear regime [51]. On the other hand, there were evidences that such singularities would be regularized in the models which contain more than four derivatives in the action [14,64]. Using the expressions (34) and (35) derived in Sec. 2 for the nonrelativistic potentials ϕ and ψ, we have shown that in a generic polynomial higherderivative gravity with more than four derivatives in both scalar and spin-2 sectors the curvatures remain finite at the origin. This happens because, under these circumstances, the coefficient of the linear term in the series expansion of the metric potentials is always equal to zero.
In the ensuing part of the paper we considered the dynamical process of the spherically symmetric collapse of null shells in linearized higher-derivative polynomial gravities. Here we generalized the discussion carried out in [50] to include the possibility of degenerate poles. If one allows the shell to have a certain thickness, then the Kretschmann invariant could become finite during the collapse, which means that the metric is regular. This observation on the regularity of the metric of the thick shell is a refinement of the result derived in [50]. Indeed, the logarithmic divergences of the Kretschmann invariant which in principle may occur in polynomial theories are actually ruled out in most of the cases, due to a specific algebraic relation between the poles of the propagator. Only in the fourth-order gravity these logarithmic divergences are possible. Finally, we have shown that, like in the case of polynomial gravities with simple poles in the propagator [50], there exists a mass gap for the mini black hole formation also in the models with higher-order poles.
With the results obtained in the present work it is possible to observe some similarities between the nonlocal (ghost-free) gravity and the local (polynomial) model with more than four derivatives. First, in both theories there is the cancellation of the Newtonian singularity of the metric potentials associated to a δ-source [17,19,20,45,46]. Second, it is known that in the nonrelativistic limit the nonlocal ghost-free gravity has a regular solution for the field generated by a pointlike source [19,68]. Our results show that in a generic polynomial higher-derivative gravity with more than four derivatives in each sector the Newtonian limit is regular too. A third similarity is the regularity of the metric of the collapsing shell. In fact, if a thin shell is considered, nonpolynomial and nontrivial polynomial theories have a Kretschmann scalar which diverges quadratically for small r [49,50]. This is, however, the consequence of the nonphysical assumption of an infinitesimally thin shell. If the shell has some thickness, then in both theories the Kretschmann invariant becomes finite during the collapse. This happens, again, because the leading term in the expansion of Eq. (93) around z = 0 is the linear one, just like what occurs in the nonlocal ghost-free gravity (see [49,50]). Solely in the fourth-order gravity the divergences in the Kretschmann invariant are possible, a situation analogous to what happens in the Newtonian limit. Moreover, in both theories there is a mass gap for the mini black hole formation. Indeed, this feature is present in any higher-derivative model with an arbitrary number of derivatives in the spin-2 sector [49,50,57], since in these theories there is a new mass scale. These four parallels between polynomial and ghost-free gravity theories can be supportive of the view that the nonlocal models may be considered as the limit of a theory with an infinite amount of complex poles hidden at the infinity [43]. In this sense, it is useful to notice that many good regularity properties of the nonlocal gravity [19,68,69] can be achieved without the need of losing locality at the classical level, and may be common to models with at least six derivatives.
All the results which were mentioned above have been obtained in the linear approximation. The most interesting question is whether there can be nonsingular solutions in the full nonlinear regime of polynomial gravity theories. The first step in this direction was done within the fourth-order gravity in Ref. [51], where the asymptotic analysis of the static field equations near the origin was carried out via the Frobenius technique. It was shown the existence of three families of solutions: a set of nonsingular solutions, and two sets of singular ones-one of them containing the Schwarzschild solution. The presence of the Schwarzschild solution is expected, because by means of the Gauss-Bonnet relation where E = R 2 µναβ − 4R 2 µν + R 2 , it is possible to completely remove the Riemann-squared term of the action, and it is clear that any vacuum solution of the Einstein equations (R µν = 0) is also a solution of the fourth-order gravity [51,70]. Nonetheless, in this model the Schwarzschild solution is not coupled to a positive-definite matter source [51]. More recently, some new aspects of the nonlinear static spherically solutions in fourthorder gravity were considered in [64][65][66] by means of numerical methods. In particular, it was studied what happens when the asymptotic solutions in strong-field regime near r = 0 are linked with the weak-field solution at large r in the form of a combination of Newtonian and Yukawa potentials (such a potential is the particular case of our general result, Eqs. (34) and (35)). In summary, the result is that for a δ-like source the solution has no horizon and falls to a timelike singularity at r = 0. Actually, the presence of the singularity in this solution is expected in view of the fact that R 2 µναβ diverges yet at the linear regime. Moreover, the absence of horizon in the full fourth-order model is guaranteed by a general theorem [65,66,71]; and only the particular theories where the R 2 -term is excluded from the action could have horizons.
In what concerns the theories with derivatives higher than fourth, in Ref. [64] the asymptotic solutions near r = 0 were studied by the Frobenius series expansion method in models with up to 10 derivatives in the action. It was shown that there is no Schwarzschildlike solutions, or other ones with singularity. Only the nonsingular solutions remain in the static spherically symmetric case for sixth-and higher-order theories 10 . The nonexistence of the exact Schwarzschild solution is due to the absence of the Gauss-Bonnet relation for the higher-order terms. The analogue relation (3) is insufficient to eliminate the effect of the Riemann-squared terms in the nonlinear regime, since O(R 3 ) structures still remain. Also, the nature of these nonsingular solutions implies that the complete solutions with large r behavior given by Eqs. (34) and (35) must have no horizon or an even number of horizons. Another interesting result of Ref. [64] is the necessity of theories with six or more derivatives to the possible elimination of the de Sitter-like horizons.
The results of the present work, in light of [64], bring more motivations for further investigation of the spherically symmetric static solutions in the full nonlinear regime for the polynomial theories with more than four derivatives. It would also be interesting to know whether in these theories there is some kind of no-horizon theorem, and we expect to revisit this issue in the future. In case of a positive answer, the complicated numerical search of solutions might be simplified.

A Useful identities with the coefficients A i,j
Let a(z) be a polynomial function of degree N ≥ 1 which satisfies a(0) = 1. The quantities A i,j defined by (27) are related to the coefficients a i,j of the partial fraction expansion of In fact, a i,j = A i,j (j − 1)! and, in particular, A i,1 = a i,1 and A i,2 = a i,2 . Proceeding the regrouping of the r.h.s. into a single fraction one obtains 10 We point out, however, that the method based on the expansion in Frobenius series around r = 0 is not sufficient to rule out the existence of singularities, as there may be solutions with a violent singularity which does not admit such representation at the origin. Also, it is possible to have solutions with singularities at a finite radius.
Comparing the numerators of the fractions above order by order in ξ, one obtains for the highest order term (N = i α i ) The substitution of this result into (39) shows that the Newtonian singularity is canceled in general higher-derivative models. Now, let us assume that N ≥ 2. Comparing both sides of (121) for the term proportional to ξ N −1 one obtains Since, for a given i, and using (123), it follows that In terms of the definitions in the Eq. (52), the identity above reads S 2 = S 1 . We recall that this relation is valid only if N ≥ 2. The case N = 1 implies in S 1 = m 2 1 and S 2 = 0.