Can static regular black holes form from gravitational collapse?

Starting from the Oppenheimer-Snyder model, we know how in classical general relativity the gravitational collapse of matter form a black hole with a central spacetime singularity. It is widely believed that the singularity must be removed by quantum gravity effects. Some static quantum-inspired singularity-free black hole solutions have been proposed in the literature, but when one considers simple examples of gravitational collapse the classical singularity is replaced by a bounce, after which the collapsing matter expands for ever. We may expect 3 possible explanations: $i)$ the static regular black hole solutions are not physical, in the sense that they cannot be realized in Nature, $ii)$ the final product of the collapse is not unique, but it depends on the initial conditions, or $iii)$ boundary effects play an important role and our simple models miss important physics. In the latter case, after proper adjustment, the bouncing solution would approach the static one. We argue that the"correct answer"may be related to the appearance of a ghost state in de Sitter spacetimes with super Planckian mass. Our black holes have indeed a de Sitter core and the ghost would make these configurations unstable. Therefore we believe that these black hole static solutions represent the transient phase of a gravitational collapse, but never survive as asymptotic states.


I. INTRODUCTION
In classical general relativity, under the main assumptions of the validity of the strong energy condition and of the existence of global hyperbolicity, the collapse of matter inevitable produces a singularity of the spacetime [1]. At a singularity, predictability is lost and standard physics breaks down. According to the weak cosmic censorship conjecture, spacetime singularities formed from collapse must be hidden behind an event horizon and the final product of the collapse must be a black hole [2]. In 4-dimensional general relativity, the only uncharged black hole solution is the Kerr metric [3], which reduces to the Schwarzschild solution in the spherically symmetric case. The Oppenheimer-Snyder model is the simplest fully analytic example of gravitational collapse, describing the contraction of a homogeneous spherically symmetric cloud of dust [4]. It clearly shows how the collapse produces a spacetime singularity and the final product is a Schwarzschild black hole.
In analogy with the appearance of divergent quantities in other classical theories, it is widely believed that spacetime singularities are a symptom of the limitations of classical general relativity, to be removed by quantum gravity effects. While we do not have yet any robust and reliable theory of quantum gravity, the resolution of spacetime singularities has been investigated in many quantum gravity inspired models. Very different approaches have studied corrections to the Schwarzschild/Kerr solution, finding black hole metrics in which the curvature invariants are always finite [5] 1 .
In the same spirit, one can study the modifications to the Oppenheimer-Snyder solution and to other models of collapse. In this case, the singularity is replaced by a bounce, after which the cloud starts expanding [6]. It is therefore disappointing that the quantum-gravity corrected model of collapse does not reproduce to the quantum-gravity corrected black hole solution.
In this paper, we want to investigate this apparent contradictory result. First, we determine both the quantum-gravity corrected static black hole metric and the quantum-gravity corrected homogeneous collapse solution within the same theoretical framework, since the ones reported in the literature come from different models. We find that the problem indeed exists. Second, we try to figure out the possible reason. One possibility is that the static regular black hole spacetimes are ad hoc solutions, but they cannot be created in a collapse and therefore they are physically irrelevant. The collapse always produces an object that bounces. Another possible explanation is that the final product of the collapse depends on the initial conditions. The collapse of a homogeneous cloud creates an object that bounces, while with other initial conditions (not known at present) the final product is a static regular black hole. Lastly, it is possible that the simple homogeneous collapse oversimplifies the model, ingoing and outgoing energy fluxes between the interior and the exterior solutions are important, and, after proper readjustment that seems to be difficult to have under control within an analytic approach, the collapsing model approaches the static regular black hole solution. Our quantum-gravity inspired theories are unitary, super-renormalizable or finite at quantum level, and there are no extra degrees of freedom at perturbative level around flat spacetime. This should rule out the possibility that the explanation of our puzzle is due to the fact that these models may not be consistent descriptions of quantum gravity. However, these theories display a ghost state in de Sitter spacetime when the cosmological constant exceeds the square of the Planck mass. This fact may be responsible for our finding and answers the question in the title of this paper. Our black holes have indeed a de Sitter core with an effective cosmological constant larger than the square of the Planck mass when the black hole mass exceeds the Planck mass. The presence of a ghost makes the solutions unstable and therefore they cannot be the final product of the gravitational collapse.
The content of the paper is as follows. In the next section, we briefly review the classical homogeneous and spherically symmetric collapse model. In Section III, we derive the spherically symmetric black hole solutions in a super-renormalizable and asymptotically-free theory of gravity with the family of form factors proposed by Krasnikov [7] and Tomboulis [8]. In Section IV, we study the spherically symmetric homogeneous collapse in the same models. Summary and conclusions are in Section V.

II. BLACK HOLES AND GRAVITATIONAL COLLAPSE IN CLASSICAL GENERAL RELATIVITY
The most general spherically symmetric metric describing a collapsing cloud of matter in comoving coordinates is given by where dΩ 2 represents the line element on the unit twosphere and ν, R, and G are functions of t and r. The energy-momentum tensor is given by T µ ν = diag(−ρ(r, t), p r (r, t), p θ (r, t), p θ (r, t)) , (2) and the Einstein equations reduces to where the ′ denotes a derivative with respect to r, and the˙denotes a derivative with respect to t. The function F (r, t) is called Misner-Sharp mass and is defined by The whole system has a gauge degree of freedom that can be fixed by setting the scale at a certain time. One usually sets the area radius R(r, t) to be equal to the comoving radius r at the initial time t i = 0, i.e. R(r, 0) = r. We can then introduce a scale factor a R(r, t) = ra , that will go from 1, at the initial time, to 0, at the time of the formation of the singularity. The condition to describe collapse isȧ < 0. For a homogeneous perfect fluid, p r = p θ = p(t). The simplest case is the gravitational collapse of a cloud of dust, p = 0, which is the well-known Oppenheimer-Snyder model [4]. From Eq. (3), we see that F is proportional to the amount of matter enclosed within the shell labeled by r at the time t. For dust, from Eq. (4) it follows that F is independent of t, so there are no inflows and outflows through any spherically symmetric shell of radial coordinate r. If r b denotes the comoving radial coordinate of the boundary of the cloud, F (r b ) = 2M Sch , where M Sch is the Schwarzschild mass of the vacuum exterior. Let us note that in the general case of a perfect fluid that is not true, and for a non-vanishing pressure the homogeneous spherically symmetric interior must be matched with a non-vacuum Vaidya exterior spacetime. Eq. (5) reduces to ν ′ = 0 and one can always choose the time gauge in such a way that ν = 0. From Eq. (6), we find that G is independent of t and we can write G = 1 + f (r). In the homogeneous marginally bound case (representing particles that fall from infinity with zero initial velocity), f = 0 and therefore G = 1.
In the homogeneous spherically symmetric gravitational collapse of a cloud of dust, one finds that the energy density is given by where ρ 0 is the energy density at the initial time t i = 0, and the scale factor is The model has a strong curvature singularity for a = 0, which occurs at the time t s = 2/ √ 3ρ 0 , as can be seen from the divergence of the Kretschmann scalar The boundary of the cloud collapses along the curve R(r b , t) = r b a(t) and the whole cloud becomes trapped inside the event horizon at the time t tr < t s for which For p = 0, the exterior solution is not Schwarzschild, but the Vaidya spacetime, because there is a nonvanishing flux through the boundary r b and F (r b , t) does depend on time. In the radiation case, one finds where ρ 0 is the energy density at the initial time t i = 0 and the singularity occurs at the time t s = 3/ρ 0 /2. Like in the dust model, the final outcome is a Schwarzschild black hole. For the generic case of a perfect fluid with equation of state p = ωρ, the scale factor is (for ω = −1) where t s is the time of the formation of the singularity The time of the formation of the event horizon is still given by R(r b , t tr ) = F (r b , t tr ) and occurs before the time t s .

III. QUANTUM-GRAVITY INSPIRED BLACK HOLES
We start from the classical Lagrangian of the renormalized theory in [9], where ✷ Λ = p µ p µ /Λ 2 and we use the signature (−, +, +, +). The main properties of our theoretical framework are discussed in Appendix A, where we show that the theories studied in this paper are unitary, superrenormalizable or finite at the quantum level, and there are no extra degrees of freedom (ghosts or tachyons) in flat spacetime. We note that we are not considering "local higher derivative theories of gravity", but "weakly non-local theories of gravity". Here with non-locality we mean that we have an operator with an infinite number of derivatives, while in a local theory the number of these derivatives would be finite. Weakly because it is only the whole sum that makes the theory non-local. However, the non-locality is not enough to have a good theory. We need that the propagator is the standard one times an entire function without zeros, singularities or poles in the whole complex plane. In this case, the theory does not have ghosts by construction, because the residue of the propagator at the pole is the same as the one of general relativity. The regularization of the solutions is thus due to the choice of the form factor and to the absence of interactions at high energy or, in other words, to asymptotic freedom. More details can be found in the last paper in [6]. The equations of motion for the theory up to terms quadratic in the curvature are The right hand side can be considered as an effective energy-momentum term, defined by S µν = V (−✷ Λ )T µν .
Within this approximation, the left hand side is compatible with the Bianchi identity, so the effective energymomentum tensor is conserved Now we want to solve the field equations for a static and spherically symmetric source. From the static property, the four-velocity is u µ = (u 0 , 0); that is, only the timelike component is non-zero, and u 0 = (−g 00 ) −1/2 . For simplicity, we consider a point source. The T 0 0 component is given by . Spherical or Cartesian-like coordinates are adopted to make calculations easier. In the spherically symmetric case, the metric is assumed to have the Schwarzschild form, where F (r) is and m(r) is the mass enclosed within the radius r. The effective energy-momentum tensor is defined by The S 0 0 component can be rewritten as where r = | x|. This is the representation for effective energy density, and, once the form factor is specified, it can be numerically solved to get ρ e (r). The radial integral contains the term sinkr kr . If we expand the sine function, we have Independently of the choice of the form factor V (z), the leading term is a constant, which means that at r = 0 we will always have a positive effective energy density proportional to the mass M and fully determined by V (z). As long as the convergence velocity of V (z) is larger than k −3 , we are able to get a finite effective energy density. In classical general relativity, V (z) = 1 and the result is not finite. The covariant conservation and the additional condition g tt = −g −1 rr completely determines S µν . The Einstein equations give From Eq. (25), we find the mass enclosed in the radius r At r = 0, m = 0 and F = 1. For r → ∞, m is a constant and F (r) → 1. We have thus solutions with two or more horizons. Moreover, at r = 0 a constant energy density gives a de Sitter spacetime, independently of the choice of the form factor. Here the effective cosmological constant is of order κ 2 M Λ 3 and this, as argued at the end of this paper, may be the key-ingredient to address our question. We now consider two specific form factors, proposed respectively by Krasnikov [7] and by Tomboulis [8]: Here p γ+1 (z) is a polynomial of order γ + 1. The superrenormalizability of the theory requires γ ≥ 3 and in what follows we will only consider the minimal renormalizable theory with γ = 3. In the low energy limit, z ≡ −✷ Λ → 0, and to recover general relativity we need V (z) → 1. So p γ+1 (0) = 0. Moreover, we should expect deviations from general relativity when z = 0, so for any z > 0 we have p γ+1 (z) = 0. This argument is in accordance with the restriction for p γ+1 (z). We can therefore consider three cases: p γ+1 (z) = z 4 , and p γ+1 (z) = z 4 ± 6z 3 + 10z 2 . The latter two cases are taken as a generalization of the first one. Fig. (1) shows ρ e (r) obtained from the numerical integration of Eq. (23) with the Krasnikov form factor and n = 1, 5, 10. The left panel is for M = 10, the right panel for M = 1. Since ρ e (r) ∝ M , if we change M we change the scale, without altering the shape of ρ e (r). For n ≥ 2, ρ e is not monotonic and can assume negative values. Because of that, it is possible to have more than two horizons. These plots show also that the effective energy density ρ e approaches a finite value for r → 0. Fig. (2) shows the corresponding F (r) functions. For M = 10, we find at least two horizons, while for M = 1 there is no horizons. For a large n, the oscillations of ρ e are stronger and therefore it is possible to form more than two horizons. However, this multi-horizon situation only exists when M ∼ 10 and n is large. For instance, we found that there are just two horizons when M = 100 and n = 10.
Figs. 3 and 4 are for the Tomboulis form factor. Fig. 3 shows the effective energy density ρ e (r): like for the Krasnikov form factor, ρ e is finite for r = 0 and it has an oscillatory behavior near ρ e = 0, so that it can be negative at some radii. Fig. 4 shows the function F (r). For M = 10 we have two horizons, while there is no horizons when M = 1. The choice of p γ+1 (z) mainly affects the region r < 2M .

IV. NON-LOCAL GRAVITY INSPIRED COLLAPSE
Now we want to find the quantum-gravity corrected solution of the gravitational collapse for a homogeneous and spherically symmetric cloud. The scale factor a(t) is determined through the propagator approach. We first write the metric as a flat Minkowski background plus a fluctuation h µν , where η µν = diag(−1, 1, 1, 1). The conformal scale factor a(t) and the fluctuation h µν (t, x) are related by the following relations: After a gauge transformation, we can rewrite the fluctuation in the usual harmonic gauge The fluctuation now reads We can then switch to the standard gravitational "barred" fieldh ′ µν defined bȳ The classical solution for the homogeneous and spherically symmetric gravitational collapse is known. We can thus compute the Fourier transformh(E) defined in (36).
If we know the classical solution a(t), we can find the distributionh(E) that provides the correct solution for V (−p 2 /Λ 2 ) = 1. We can then use a different form factor to find the corresponding quantum-gravity corrected scale factor a(t). Here we want to find the solutions for the gravitational collapse of a cloud of dust and radiation with the form factors of Krasnikov and Tomboulis.  energy density (right panel) for the Krasnikov form factor with n = 2 and 6 in the case of a cloud of radiation, ω = 1/3. The dust case is shown in Fig. 6. Classically, the scale factor a(t) monotonically decreases and finally vanishes. The corresponding energy density ρ(t) therefore diverges at the time t = 0. In the quantum-gravity corrected picture, we find a bounce: a(t) reaches a minimum and then start increasing. Far from the bounce, the classical and the quantum solutions are similar, while the difference becomes important at high densities.
The calculations with the Tomboulis form factor turn out to be significantly more complicated. We thus use the following approximated form factor (44) As shown in Fig. (7), such an approximated form factor is very similar to the exact one, with a tiny deviation near z = ±1. While such an approximated form factor has a discontinuity at z = ±1, the Fourier transform is applicable as the limits at z = ±1 in both directions are well-defined. Here we have only studied the minimal renormalizable theory with γ = 3 and have chosen p γ+1 (z) = z 4 . The part |p γ+1 (z)| > 1 is trivial and we can obtain an analytic result. For the other part of the integral, we need to use a regularization prescription and separate the finite and the divergent parts. The detailed calculation is in Appendix B. Figs. 8 and 9 show, respectively, the radiation and dust case. The qualitative behavior of the scale factor a(t) and the energy density ρ(t) is the same as the one of the Krasnikov model.

V. SUMMARY AND CONCLUSIONS
In the present paper, we have studied both the static black hole solution and the homogeneous spherically symmetric collapse of a cloud of matter in a superrenormalizable and asymptotically free theory of gravity. The spacetime singularity predicted in classical general relativity is removed in both the cases. In the literature there were so far some scattered results in different theoretical frameworks. Here we have studied this issue in more details within the Krasnikov and Tomboulis models.
Static and spherically symmetric singularity free black hole solutions have been obtained. At the origin, the effective energy density is always finite and positive, independently of the exact expression of the form factor V (z). In other words, these black holes have a de Sitter core in their interior, where the effective cosmological constant is of order κ 2 M Λ, κ 2 = 32πG N , M is the black hole mass, and Λ is the energy scale of the theory which is naturally to expect to be close to the Planck mass. The singularity of the spacetime is therefore avoided due to the repulsive behavior of the gravitational force. For a large family of form factors, the effective energy density can be negative in some regions, which eventually provides the possibility of having multi-horizon black holes. In the homogeneous and spherically symmetric collapse of a cloud of matter, the formation of the singularity is always replaced by a bounce. Far from the bounce, the collapse follows the classical solution, while it departs from it at high densities. Strictly speaking, asymptotic freedom is sufficient to remove the singularity, but the presence of a bounce requires also a repulsive character for gravitational field in the high energy regime.
In conclusion, we have provided some convincing examples that show how the final products of the quantumgravity corrected collapse solutions are not the quantumgravity corrected Schwarzschild black hole metrics. This is not the result that one would expect a priori. There may be 3 natural explanations.
1. Static regular black holes cannot be created in any physical process. In this case, even if they are solution of a theory, they are much less interesting than their classical counterparts that can be created in a collapse.
2. The final product of the gravitational collapse is not unique. The collapse of a homogeneous and spherically symmetric cloud of matter does not produce a static regular black hole, but the collapsing matter bounces and then expands. With different initial conditions, not known at the moment, static regular black holes may form.
3. The simple example of a homogeneous cloud of matter oversimplifies the picture and misses important physics. As discussed in Section II, in the classical dust case we have a homogeneous interior and a Schwarzschild exterior without ingoing or outgoing flux through any spherical shell of comoving radial coordinate r. However, that is not true in general, and the exterior spacetime is a generalized Vaidya solutions with ingoing or outgoing flux of energy. This means that the homogeneous solution is not stable and must evolve to an inhomogeneous model. While the bounce can still occur, after it the collapsing matter may not expand forever. The boundary effects are important and, after proper readjustment that can unlikely be described without a numerical strategy, the collapse approaches the static black hole solution.
The possibility 1 excludes the possibilities 2 and 3, but the latter may also coexist. Here we have focused on the asymptotically free gravity theory with the Lagrangian given in Eq. (17), and we have shown that the issue indeed exists. The theoretical model is not sick, and therefore we cannot attribute the problem to the fact that we  are considering a non-consistent quantum theory. However, it is easy to compute the propagator for this class of theories around de Sitter spacetime background and to show the presence of a ghost when the cosmological constant exceeds the square of the Planck mass (in preparation). This fact may explain the fact that our static black holes are not the final product of the gravitational collapse. These black holes have indeed a de Sitter core in which the effective cosmological constant is κ 2 M Λ. If the black hole mass M exceeds the Planck mass, there is a ghost and the black hole is unstable. Therefore the solutions here presented cannot be the final product, but only an intermediate phase of the gravitational collapse.

ACKNOWLEDGMENTS
This work was supported by the NSFC grant No. 11305038, the Shanghai Municipal Education Commission grant for Innovative Programs No. 14ZZ001, the Thousand Young Talents Program, and Fudan University.

Appendix A: Theoretical framework
For the theories studied in this paper, the classical action is where G µν is the Einstein tensor and κ 2 = 32πG N . All the non-polynomiality is in the form factor V (−✷/Λ 2 ), which must be an entire function. Λ is the Lorentz invariant energy scale and it is not subject to infinite or finite (non analytic) renormalizations. The natural value of Λ is of order the Planck mass and in this case all the observational constraints are satisfied. Indeed, at classical level all the corrections to the Einstein-Hilbert action are suppressed by 1/Λ, and if the value of Λ is large the theory can reduce to general relativity at low energies. At quantum level, the introduction of non-local operators in the action could potentially lead to strong non-localities generated by the renormalization group flow towards the infrared, in disagreement with observations. This is not the case here, as a consequence of the Donoghue argument [10] and of the fact that our action only involves entire functions. The entire function V (−✷/Λ 2 ) must have no poles in the whole complex plane, in order to ensure unitarity, and must exhibit at least logarithmic behavior in the ultraviolet regime, to give super-renormalizablitity at the quantum level. The theory is uniquely specified once the form factor is fixed, because the latter does not receive any renormalization: the ultraviolet theory is dominated by the bare action (that is, counterterms are negligible). In this class of theories, we only have the graviton pole. Since V (z) is an entire function, there are no ghosts and no tachyons, independently of the number of time derivatives present in the action. This is the main reason to introduce a non-polynomial Lagrangian. At the phenomenological level the form factor could be experimentally constrained, for example measuring the corrections to the gravitational potential, or hypothetically measuring a cross section in a scattering process at high energy. Concerning the difficulties with particular form factors and non-local operators, we note that the class of operators introduced by Krasnikov and by Tomboulis are well defined in the Euclidean as well as in the Lorentzian case, because (k 2 E ) 2 = (k 2 ) 2 , where k E is the momentum in the Euclidean space [7]. More details on the ultraviolet and infrared properties of this class of theories can be found in the last paper in [6]. We note that here it is possible to anti-screen gravity in the UV without introducing extra degrees of freedom because the theory is characterized by an entire function that goes to zero in the UV. We present below two different ways to see this.
Tree-level unitarity -The general and clear way to address the unitarity problem in Lagrangian formalism can be summarized as follows (for more details, see [11]): 1) we calculate the propagator expanding the action to the second order in the graviton fluctuation, 2) we calculate the amplitude with general external energy tensor sources, and 3) we evaluate the residue at the poles. A general theory is well defined if "tachyons" and "ghosts" are absent, in which case the corresponding propagator has only first poles at k 2 − M 2 = 0 with real masses (no tachyons) and with positive residues (no ghosts). In our class of theories, we only have one pole in k 2 = 0 with positive residue.
When we introduce a general source, the linearized action including the gauge-fixing reads where P (2) µνρσ and P µνρσ are the graviton projectors and θ µν = η µν − k µ k ν /k 2 . The transition amplitude in momentum space is where g is an effective coupling constant. To make the analysis more explicit, we can expand the sources using the following set of independent vectors in the momentum space where ǫ i are unit vectors orthogonal to each other and to k. The symmetric stress-energy tensor reads The conditions k µ T µν = 0 and k µ k ν T µν = 0 place constrains on the coefficients a, b, d, e i , f i . Introducing the spinprojectors and the conservation of the stress-energy tensor k µ T µν = 0 in (A5), the amplitude results where T = T µ µ . Clearly, there is only the graviton pole in k 2 = 0 and the residue at k 2 = 0 is For D > 3, we find that Res (A) k 2 =0 > 0 (because H(0) = 0), which means that the theory is unitary. Instead, for D = 3 the graviton is not a dynamical degree of freedom and the amplitude is zero.
As an example of this quantum transition, we can consider the interaction of two static point particles. Here T µ ν = diag(ρ, 0, 0, 0) with ρ = M δ( x) and the amplitude (A8) simplifies to which is positive in D > 3 and zero for D = 3 since, again, there are no local degrees of freedom in D = 3.
Källén-Lehmann (KL) representation -Any action in the class of theories here presented is defined in terms of an entire function with no zeros, singularities, nether poles. Therefore the usual derivation of the KL decomposition goes straightforward. Let us construct step by step such representation for a weakly nonlocal prototype theory. It is easy to show that, Therefore a complete set of states satisfy, where we defined the state λ p by boosting of an arbitrary boost p a momentum zero state λ 0 and we summed over all p. Assuming that 0|φ(x)|0 = 0 and Lorentz invariance of φ(0), We thus get for x 0 > y 0 : where the contour C F is closed in the lower half p 0 -plane, while the integral along the real axes is evaluated displacing the poles for Re(p 0 ) > 0 in the lower half plane and for Re(p 0 ) < 0 in the upper half plane (in our gravitational theory m λ = 0). Analog relation is obtained for x 0 < y 0 , therefore: The matrix element is Lorentz invariant and only depends on the mass of the state m λ , therefore we can write (A15) making use of the definite positive spectral density function ρ(s) defined by, as follows 0|T(φ(x)φ(y))|0 = +∞ 0 ds ∆ F (x − y; s)ρ(s), where we defined the function in analogy with the Feynman propagator. The spectral representation sees all the poles, but in this case they are the same as in general relativity, because we have introduced a function with no poles or zeros. In the case of a local theory with a finite number of derivatives, this is impossible and the theory is sick. We expect a single particle state to contribute to ρ(s) as an isolated delta The combination of the two equations gives the result of I I = − 1 2 (I 1 + I 2 ) + D , where I 1 + I 2 = π θ=0 dθ e −t sin θ sin[(n − 1)θ − t cos θ] + e t sin θ sin[(n − 1)θ + t cos θ] , )Res e iEt E n + µ n .
For radiation, n = 2 and we have dθ e −t sin θ sin[θ − t cos θ] + e t sin θ sin[θ + t cos θ] + π µ , where the divergent part is D = π µ for µ → 0. The dust case is slightly different and it is more convenient to choose the exponential form for approximation of V (z), so the integral is where α > 0 is not an integer. With the transform for the integration variable E = ρe iθ , dE = iEdθ, and applying the contour integral we have The first two terms I 1 + I 2 can be simplified as I 1 + I 2 = 2 1 0 dρ 1 ρ α e − 1 2 ρ 16 e iρt e −iπα (i sin απ) .
The integration of the unit circle gives the finite part of the required integral, and the integration around the origin gives the divergent part of order α − 1, For dust case, we have α = 7/3 and the divergence is of order 4/3.