Spherically Symmetric Black Hole Formation in Painlevé-Gullstrand Coordinates
Jonathan Ziprick and Gabor Kunstatter
Department of Physics and Astronomy
University of Manitoba
Winnipeg, Manitoba, Canada R3T 2N2
Department of Physics and Winnipeg Institute of Theoretical Physics
University of Winnipeg
Winnipeg, Manitoba, Canada R3B 2E9
We perform a numerical study of black hole formation from the spherically symmetric collapse of a massless scalar field. The calculations are done in Painlevé-Gullstrand (PG) coordinates that extend across apparent horizons and allow the numerical evolution to proceed until the onset of singularity formation. We generate spacetime maps of the collapse and illustrate the evolution of apparent horizons and trapping surfaces for various initial data. We also study the critical behaviour and find the expected Choptuik scaling with universal values for the critical exponent and echoing period consistent with the accepted values of and , respectively. The subcritical curvature scaling exhibits the expected oscillatory behaviour but the form of the periodic oscillations in the supercritical mass scaling relation, while universal with respect to initial PG data, is non-standard: it is non-sinusoidal with large amplitude cusps.
Black holes are among the most interesting solutions to Einstein’s field equations. They demonstrate in various ways the incompleteness of our classical understanding of gravity and provide hints toward developing a quantum theory. The relevant issues are the endpoint of Hawking radiation, the breakdown of general relativity at curvature singularities and the microscopic source of black hole entropy. Given the lack of experimental guidance, black holes provide a vital testing ground for new ideas in quantum gravity.
Despite a great deal of analytic and numerical work, relatively little is known about the dynamics of black hole formation except in special, highly symmetric situations, since the equations for collapsing matter fields tend to be very difficult to work with. Important pioneering work in this area was done by Choptuik  in the early 1990’s, who took advantage of modern computational methods to numerically create small mass black holes from the collapse of spherically symmetric massless scalar fields. In the small mass limit he observed unexpected critical behaviour that has since been confirmed in numerous papers for a variety of matter fields (see  for an in-depth review).
There are three fundamental properties of the critical collapse that are relevant here. First, there is a critical zero mass black hole solution that plays the role of an intermediate attractor for the dynamical evolution. In the case of massless scalar fields the critical solution exhibits discrete self similarity (DSS) described in terms of a generic set of field variables by
where is the critical solution, and are the Schwarzshild radius and time, is some integer and is a dimensionless number. To explain this more clearly, consider the profiles of the critical solution field variables frozen in time at . After a certain amount of time elapses, the same profiles will be found on a scale that is smaller by a factor of . After an additional time has passed, the same profile will again be found on a scale that is times smaller than at . This gives rise to the more descriptive name for DSS in gravitational collapse, scale echoing.
The existence of this intermediate attractor gives rise to scaling relationships for near critical solutions. For any one parameter family of data in the limit of small black holes one finds a mass scaling relation:
where is the black hole mass on formation, is a real parameter referred to as the critical exponent, is a periodic function, is any parameter in the initial data and is the threshold value of that parameter that separates the data into two categories: supercritical (leading to black hole formation) and subcritical (leading to matter dispersion). Two separate analytical studies (, ) have shown that DSS results in the period of being
Note that is often referred to as a “small oscillation” or “wiggle”, but its explicit functional form must be determined numerically .
Garfinkle and Duncan have shown that a similar scaling relation exists for the maximum curvature at the origin for subcritical solutions . For any one parameter family of data in the barely subcritical limit one finds:
where is the Ricci scalar at the origin and is a periodic component, again with period . This relation is similar to that of mass scaling except for the factor in front of being rather than . This derives from mass having units of length and curvature having units of length (with ).
The third aspect of critical collapse we note here is universality. DSS and critical scaling have been observed in the gravitational collapse of many different types of fields. The constants and are determined by the properties of the critical solution and are therefore universal: they vary for different kinds of matter but are independent of the form of initial data and the parameter in that initial data that is varied. Numerical and semi-analytical calculations for the massless scalar field have given and (, , , , , ). In all previous calculations, done in Schwarzschild or double null coordinates, the functions and were well approximated by a small amplitude sine function with period .
In this paper we study numerically black hole formation from the spherically symmetric collapse of a massless scalar field in Painlevé-Gullstrand coordinates (, )111These coordinates are sometimes referred to as “flat slice” coordinates because in vacuum the surfaces of constant time are flat.. PG coordinates were used by Wilczek and collaborators  to study Hawking radiation and more recently were proposed by Husain and Winkler  as a reasonable starting point to investigate the quantum dynamics of black holes. The key distinguishing feature of PG coordinates is that they are regular across apparent horizons. This allows us to evolve initial data until the onset of singularity formation when our chosen boundary conditions (zero mass density at the origin) are no longer valid.
We first of all take advantage of the PG coordinate system to run the code beyond initial horizon formation to study in detail the formation of apparent horizons and the evolution of the resulting trapping surfaces. We then study the critical behaviour of the system and find our data to be in good agreement with previous results for the critical exponent and echoing period. Interestingly, we find the maximum curvature scaling relation to show the usual small oscillations but discover a new form for the oscillations in the mass scaling relation. The form of these mass scaling oscillations is found to be universal in the sense that it is the same for all families of initial data (that we tested) specified in PG coordinates. They are however quite different from those observed in Schwarzschild or double null coordinates. The reason for this difference is as yet unclear, but it is almost certainly related to the fact that the location of the initial apparent horizon can depend on the choice of spacetime slicings.
While these results are interesting in themselves, our analysis is also important because it paves the way for studying the effects of quantum corrections on both the critical behaviour and the dynamics of trapping surfaces .
The paper is organized as follows: in Section 2 we lay out the foundation of our mathematics and develop the equations of motion used in our computer code. Section 3 describes our code giving relevant details of the numerical techniques and boundary conditions. Following this, we present the results of our simulations, discussing the scaling relationships and global properties of the spacetime (including trapping surfaces). We conclude with some remarks on these findings and consider future directions for this work.
2 Equations of Motion
The dynamical system that we wish to study is that of a collapsing spherically symmetric, massless scalar field in four spacetime dimensions. It turns out to be convenient to express these equations in a somewhat non-traditional parameterization, using the formalism of 2d dilaton gravity (see  for a general review of dilaton gravity and  for applications to black holes). In addition to providing a relatively simple form for the evolution equations in 4 dimensions, the formalism is generally applicable to a variety of spherically symmetric theories in four or more dimensions . We will describe the formalism in some detail because the gauge fixing is a crucial part of our analysis. We begin by stating the classical action and defining a metric. We then perform a canonical transformation and fix the gauge to yield dynamical equations for the scalar field in PG coordinates. Throughout this paper we use a prime to represent and an overdot for .
The action for a massless scalar field coupled to a dilaton in 2d is
where is Newton’s consant in 2d, is the metric determinant, is the Ricci scalar and is the dilaton coupling.
The relationship of this action to spherically symmetric gravity in spacetime dimensions is given as follows. The physical metric is:
where is the line element of the unit -sphere and
With the further substitutions
the dilaton action (2.1) is (up to boundary terms) precisely equal to that of a spherically symmetric massless scalar field minimally coupled to Einstein gravity in spacetime dimensions:
In the context of spherically symmetric gravity, performing the dimensional reduction in the action before varying yields the same dynamics as imposing spherical symmetric in the field equations themselves.
In order to make our coordinate choice explicit and rigorous, we now summarize the Hamiltonian analysis for the theory. We first define the metric in modified ADM form:
where , and are arbitrary functions of the spacetime coordinates.
The authors of  derive a partially gauge fixed Hamiltonian. For completeness, we review this procedure before applying a second gauge fixing condition that yields the completely gauge fixed metric in PG form.
The first gauge choice is
which yields the consistency condition on the lapse and shift functions:
The equations are put in a more transparent form using the canonical transformation:
where is the momentum conjugate to . is then conjugate to with their Poisson bracket being
The resulting Hamiltonian is
is called the mass function because it approaches a constant at spatial infinity where it is equal to the ADM mass of the solution. is the energy density of the scalar field.
We now completely fix the gauge by choosing
This condition, along with Eqs.(2.17) and (2.11), implies that . These gauge conditions produce the non-static generalization of PG coordinates, as can be seen by using (2.12) and the gauge conditions in the case of a vacuum to reduce the line element to PG form :
More importantly, the spatial slices will be seen to be regular across apparent horizons that form during the evolution.
We are now able to write the equations of motion for the scalar field in fully reduced form:
where and are the solutions to
The right hand side of (2.23) defines the instantaneous mass density of the configuration. It has the expected contribution from the energy density of the scalar and an additional coupling term that corresponds to the contribution from its self gravity.
The above equations need to be supplemented by boundary conditions for and . Without loss of generality we choose at . A change in this value corresponds to a trivial rescaling of the time coordinate. We also fix at , which guarantees that the metric is flat in the neighbourhood of the origin. In fact, a non-zero value of at the origin signals the formation of a singularity, so that this boundary condition would have to be relaxed/modified in order to integrate the equations past singularity formation.
3 Computational Methods
Since we work in four dimensions (n=2), we have , and . For convenience we choose , which from Eq.(2.4) corresponds to choosing the characteristic length scale to be . Note that the first gauge choice (2.11) implies that our spatial coordinate is the radius of invariant two spheres.
The first step in the iteration process is to specify initial and configurations. We work with two forms:
where , and are the parameters which may be varied to study critical collapse. We shall refer to these two forms of initial data as gaussian and tanh respectively. Note that since the mass density depends on , the tanh form has one mass peak while the gaussian data has two. For both cases we define an initial standing wave by choosing .
The most interesting behaviour occurs as matter approaches the origin on ever finer scales. In order to observe this small behavoiur while keeping a large enough lattice to contain all of the mass within the system (numerical errors occur when mass leaves the grid), we refine the -spacing near the origin so that it is several orders of magnitude smaller than the spacing near the outer endpoint. We use an adaptive time step refinement according to the minimum found across the spacial slice using the condition
where is the inverse of the local speed of an ingoing null geodesic. This provides stability by preventing information from moving over too many -points in a single time step. Having proportional to greatly increases computational times for finer -spacing, placing a practical lower limit on the size of black holes that we are able to create. However, this is not a significant concern as we are able to work on small enough scales to produce critical mass scaling.
Both time and space integrations are performed using Runge-Kutta methods accurate to fourth order in the local grid spacing. For derivatives, we experimented with finite differences and cubic splines finding the same qualitative behaviour using both, which gives some indication of reliability. In the end we chose to use finite differences exclusively since they require less computational time.
4.1 Black hole evolution
Since our coordinates are regular across apparent horizons, we are able to evolve well beyond initial black hole formation up to the point when a singularity begins to form and our boundary conditions at the origin are no longer valid. Spacetime diagrams of black hole formation and evolution are given in figures 1(a) and 1(b) for tanh and gaussian initial data respectively. Using this same data, we have also constructed animations of the mass density profile throughout this process, available for viewing at http://theoryx5.uwinnipeg.ca/users/jziprick/.
We note that apparent horizons form in pairs within constant PG time slices. This is to be expected from looking at the apparent horizon condition:
which is initially positive across the -axis. Minima begin to form as mass collapses inward. Black hole formation is indicated by the appearance of a single horizon where the function first dips to zero. This horizon splits into two as the minimum drops below zero; the outer horizon continues to grow as mass passes through it while the inner horizon moves toward the origin. Depending upon the initial mass distribution, subsequent horizon pairs may form as demonstrated in figure 1(b) which depicts the evolution of gaussian initial data. The horizon paths define a two dimensional trapping surface within which all null geodesics point towards the origin.
4.2 Critical behaviour
In our study of the critical behaviour, we varied the parameters , , and in the gaussian form of initial data, and in the tanh form. The critical values were found to an accuracy of using a binary search. We tried a variety of grid resolutions with -spacings near the origin ranging from to , which in turn affects the -spacings. The accuracy was roughly equal in all cases where black hole size was not significantly affected by the resolution (when ). However, finer grids allowed for data to be found for lower values of thus showing critical scaling over a wider range of a given parameter.
To demonstrate scale echoing, we present a plot of the scalar field at as a function of for the nearest solution to criticality that was numerically achieved using an -spacing near the origin of . is the time of black hole formation in the critical solution given by the value approached asymptotically in the small mass limit. Since we choose as our boundary condition at , is equivalent to the proper time at the origin. is seen to be periodic in with a period of found by averaging over the two oscillations between the minima at and , with an error determined by the uncertainty in . This result is in good agreement with the accepted value of .
Table 1 lists the critical exponent and period found for various lattice resolutions and families of initial data using both the mass and maximum curvature scaling relations. Our results in all cases are seen to agree with the accepted values of and .
|type of scaling||form of initial data||parameter varied||resolution near origin|
A plot of the maximum subcritical curvature scaling behaviour is given in figure 3(a). was determined using a least squares linear fit over the points where the solution is seen to follow the intermediate attractor. The error quoted indicates how this value changes depending on the range of points chosen. The period and associated error (given in the first row of table 1) were determined by fitting a sine function to the residual left after subtracting the linear fit from this data.
Plots of the supercritical mass scaling behaviour are presented in figures 3(b)–3(d). On each of these we have drawn a line that osculates the peaks of the curves. The values of were determined by using a least squares linear fit through the data points that are very near this line (within 5 % of )222Note that our results are insensitive to this choice at the given level of precision.. The error quoted is twice the standard deviation in the calculated slope. The periods are found by averaging the distance between the cusp-like valleys, with the error representing the uncertainty in locating the cusps. The data used in these calcultions consisted of the two periods corresponding to the lowest range of the parameter that produced black holes such that the size was not significantly affected by lattice resolution (for example, the periods between and in figure 3(c) were used to calculate the values given in the second row of table 1).
Previous studies of critical collapse (done in either Schwarzschild or null coordinates) have shown both scaling relationships to be a slight oscillation about a straight line that is well approximated by a sine function. While this is confirmed by our results for the subcritical maximum curvature scaling, it is clearly not the case for our supercritical mass scaling results where we find a periodic form with relatively large amplitude and a distinctly different functional form. Considering that our results agree with the known parameters of Choptuik scaling for varied initial data and mesh resolutions, and that our new-found periodic form is universal within the class of coordinate systems that we are using, it is unlikely that this periodic form is a numerical artifact.
According to Gundlach  the function that determines the small wiggle is universal. However, the location and time of formation of the initial apparent horizon can depend on the choice of spacetime slicing. Ours is the first calculation in PG coordinates. These are non-static coordinates for which data are specified on spacelike surfaces that are regular across future horizons. These unique features of PG coordinates are undoubtedly at the root of the new periodic form in the mass scaling relationship. It would be of great interest to continue the near critical evolution past initial horizon formation in order to study explicitly the scaling relation of the final horizon radius, which is independent of spacetime slicing. This is currently under investigation.
We have successfully demonstrated Choptuik scaling for a massless scalar field in PG coordinates and revealed an interesting large amplitude component to the oscillations in the mass scaling relation. Our gauge choice allows integration of the dynamical equations beyond horizon formation until the onset of singularity formation. This allowed us to study in detail the evolution of black holes showing explicitly that apparent horizons form in pairs with the number of pairs depending on the complexity of initial data.
Particularly exciting is the prospect of incorporating quantum corrections into the gravitational potential333GK is grateful to R.B. Mann and A. Buchel for raising the question about the possible effects of quantum gravity corrections on critical collapse.. Loop Quantum Gravity suggests that the underlying discreteness at the quantum level will give rise in the semiclassical limit to an effective short range repulsive component to the gravitational potential which in turn can lead to singularity avoidance [19, 20]. Husain  has done an investigation of the effect of such a correction on Choptuik scaling using double null coordinates which did not allow the code to be extended past horizon formation. Since we are able to investigate the dynamics across the entire spacetime we can investigate directly the effects of short range modifications to the gravitational potential near the onset of singularity formation. We have begun such an investigation  and preliminary results indicate dynamical singularity avoidance that allows us to calculate numerically the non-singular spacetime showing black hole formation and evaporation. These results will be presented in a future publication.
We thank Ramin Daghigh, David Garfinkle, Viqar Husain, Randy Kobes, Jorma Louko and Ari Peltola for helpful discussions and correspondence. The authors are grateful to Westgrid for providing the computer resources used for some of the simulations. GK thanks H. Maeda for useful discussions regarding trapping surfaces. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada.
-  M.W. Choptuik, “Universality and Scaling in Gravitational Collapse of a Massless Scalar Field”, Phys. Rev. Lett. 70(1) 9 (1993).
-  C. Gundlach, “Critical phenomena in gravitational collapse”, Phys. Rept. 376(6) 339 (2003).
-  C. Gundlach, “Understanding critical collapse of a scalar field”, Phys. Rev. D 55(2) 695 (1997).
-  S. Hod and T. Piran, “Fine structure of Choptuik’s mass-scaling relation”, Phys. Rev. D 55(2), 440 (1997).
-  D. Garfinkle, G. Duncan, “Scaling of curvature in sub-critical gravitational collapse”, Phys. Rev. D 58 064024 (1998).
-  D. Garfinkle, “Choptuik Scaling in Null Coordinates”, Phys. Rev. D 51, 5558 (1995).
-  R.S. Hamad e and J. M. Stewart, “The spherically symmetric collapse of a massless scalar field”, Class. Quantum Grav. 13, 497 (1996).
-  J. Bland, G. Kunstatter, M. Becker and B. Preston and V. Husain, “Dimension-Dependence of the Critical Exponent in Spherically Symmetric Gravitational Collapse”, Class. Quantum Grav. 22 5355 (2005).
-  P. Painlevé, “La mécanique classique et la théorie de la relativité”, C. R. Acad. Sci. (Paris) 173 677 (1921).
-  A. Gullstrand, “Allgemeine Lösung des statischen Eink rperproblems in der Einsteinschen Gravitationstheorie”, Arkiv. Mat. Astron. Fys. 16(8) 1 (1922).
-  P. Kraus and F. Wilcek, Nucl. Phys. B (1995); M. Parikh and F. Wiczek, Phys. Rev. Letts. (2000);
-  V. Husain and O. Winkler, “Flat Slice Hamiltonian Formalism for dynamical black holes”, Phys. Rev. D71 104001 (2005) [gr-qc/0503031]; “Quantum black holes from null expansion operators”, Class. Quant. Grav. 22 L135-L14 (2005) [gr-qc/0412039]; “Quantum resolution of black hole singularities”, Class. Quant. Grav. 22 L127-L134 (2005) [gr-qc/0410125].
-  G. Kunstatter and J. Ziprick, work in progress.
-  D. Grumiller, W. Kummer and D.V. Vassilevich, “Dilaton gravity in two dimensions”, Phys. Rept. 369 327 (2002).
-  G. Kunstatter, R. Petryk and S. Shelemy, “Hamiltonian Thermodynamics of Black Holes in Generic 2-D Dilaton Gravity”, Phys. Rev. D 57, 3537 (1998).
-  R.G. Daghigh, J. Gegenberg and G. Kunstatter, “A Partially Gauge Fixed Hamiltonian for Scalar Field Collapse”, Class. Quant. Grav. 24 2099 (2007).
-  G. Kunstatter and J. Louko, “Transgressing the horizons: time operator in two dimensional dilaton gravity”, Phys. Rev. D 75 024036 (2007).
-  D. Garfinkle, C. Cutler and G. Duncan, “Choptuik scaling in six dimensions”, Phys.Rev. D 60 104007 (1999).
-  A. Ashtekar and M. Bojowald, “Quantum Geometry and the Schwarzschild Singularity”, Class. Quant. Grav. 23, 391 (2006) [arXiv:gr-qc/0509075].
-  L. Modesto, “Loop Quantum Black Hole”, Class. Quant. Grav. 23, 5587 (2006) [arXiv:gr-qc/0509078]; “Black Hole Interior from Loop Quantum Gravity”, arXiv:gr-qc/0611043; C. G. Boehmer and K. Vandersloot, “Loop Quantum Dynamics of the Schwarzschild Interior”, Phys. Rev. D 76, 104030 (2007) [arXiv:0709.2129 [gr-qc]]; M. Campiglia, R. Gambini and J. Pullin, “Loop Quantization of Spherically Symmetric Midi-Superspaces : The Interior Problem”, AIP Conf. Proc. 977, 52 (2008) [arXiv:0712.0817 [gr-qc]]; R. Gambini and J. Pullin, “Black holes in Loop Quantum Gravity: the Complete Space-Time”, arXiv:0805.1187 [gr-qc].
-  V. Husain, “Critical behaviour in quantum gravitational collapse”, arXiv:0808.0949v1 (2008).