New results on the dynamics of critical collapse

The critical phenomena in gravitational collapse discovered by Choptuik demonstrate the rich dynamics of the Einstein equations [1]. Consider the gravitational collapse of generic families of a massless scalar field, whose initial data are parameterized by p. The parameter p measures the strength of the gravitational interaction. Strong interactions (high p) result in black hole formation, whereas weak interactions (low p) disperse the matter field to infinity, and flat spacetime remains. By fine-tuning p to the threshold of black hole formation, $ p=p_* $, critical collapse occurs.

In supercritical collapse, a tiny black hole forms, the mass of which has the scaling relation $ m_{BH}\propto|p-p_*|^{\gamma} $, where $ \gamma\simeq 0.37 $. The critical collapse solution exhibits a universality feature, namely, the spacetime produced by different families of critical initial data approaches the same solution after a finite time. The solution also displays discrete self-similarity: it is invariant after rescaling the spacetime by a certain factor. Since this discovery, similar results have been obtained in many other models (see Ref. [2] for a review). The critical behavior in the 3D scalar collapse was studied in Ref. [3]. The critical phenomena in the collapse of electromagnetic waves were simulated in Refs. [4, 5]. The interplay between multiple near-critical fields in spherical scalar collapse was investigated in Ref. [6]. The critical phenomena in the bald/scalarized black hole phase transition occurring in Einstein-Maxwell-scalar theory were reported in Ref. [7]. In simulations, near the threshold, all the intermediate solutions are attracted to the critical solution and remain in this state for a long time. The system resembles the behavior of quasinormal modes. In the late stage, the intermediate solutions decay into Reissner-Nordström black holes for subcritical cases or scalarized charged black holes for supercritical cases, resembling the quasinormal modes.

Analytic interpretations are important for understanding the dynamics of gravitational collapse. In Refs. [8−10], critical collapse was treated as an eigenvalue problem. By imposing discrete self-similarity, the global structure of the critical collapse spacetime was constructed with the pseudo-Fourier method. The rescaling factor Δ became an eigenvalue and was solved with high precision. The scaling law of the black hole mass in supercritical collapse was recovered analytically via the perturbation approach in Ref. [9]. Critical collapse was analyzed using a renormalization group method in Refs. [11, 12]. In Ref. [13], using an explicit approximate solution, a true solution was shown to exist. In Ref. [14], using a typical log-periodic formula in discrete scale invariance systems, the authors obtained an approximate analytic solution for the spacetime near the center. Approximate analytic expressions for the metric functions and matter field near the central singularity in black hole formation were obtained in Refs. [15, 16]. In Ref. [17], the equations for the matter field in critical collapse were analyzed with certain terms in the equations dropped. Furthermore, approximate expressions for certain combinations of the metric functions and derivatives of the scalar field were obtained.

In this study, considering the significance of analytic results, we obtain approximate analytic expressions for the metric functions and matter field in the large-radius region using numerical data. We also investigate the dynamics in the central region. We find that owing to the boundary conditions at the center, the equation of motion for the scalar field in the central region is reduced to the flat-spacetime form.

This paper is organized as follows. In Sec. II, we describe the methodology for simulating critical collapse. In Secs. III and IV, we study the dynamics in the large-radius and central regions, respectively. Finally, the results are summarized in Sec. V.

We simulate the critical collapse of a spherically symmetric massless scalar field ϕ in polar coordinates,

The energy-momentum tensor for the scalar field is $ T_{\mu\nu}=\phi_{,\mu}\phi_{,\nu}-(1/2)g_{\mu\nu}g^{\alpha\beta}\phi_{,\alpha}\phi_{,\beta} $. Some components of the Einstein and energy-momentum tensors for ϕ are shown below.

Here, $ (_{,r}) $ and $ (_{,t}) $ denote partial derivatives with respect to the coordinates r and t, respectively.

We define

and set $ G=1 $. Then, the equations $ G^{t}_{t}=8\pi T^{t}_{t} $ and $ G^{r}_{t}= 8\pi T^{r}_{t} $ respectively generate

The combination of $ G^{r}_{r}=8\pi T^{r}_{r} $ and Eq. (3) yields

With Eq. (2), we obtain

From the conservation of the energy-momentum tensor, $ T^{\mu\nu}_{;\mu}=0 $, we have

The Misner-Sharp mass is defined as [18]

The initial conditions for ϕ are set as $\phi|_{t_i}= a\exp[-(r/\sigma)^{2}]$ and $ \phi_{,t}|_{t_i}=0 $. The regularity of Eq. (3) at the center requires that $ A(r=0,t)=1 $. We choose $\delta(r=0, t)=0$, which implies that the coordinate time is equal to its proper time at the center. In the simulation, we integrate Eqs. (3) and (5)−(7) using the fourth-order Runge-Kutta method. A mesh refinement algorithm is implemented. We approach the critical solution by fine-tuning the initial profile of the scalar field via the bisection method. For details on the numerics, see Ref. [19].

We rewrite the metric (1) as

For convenience, we adjust the time coordinate such that $ t=0 $ when the naked singularity forms and define $X(r,t)\equiv \sqrt{2\pi}(r/\beta)\phi_{,r}$, $ Y(r,t)\equiv\sqrt{2\pi}(r/\alpha)\phi_{,t} $, $ \rho\equiv \ln r $, $ T\equiv\ln(-t) $, and $ u\equiv t/r $. Then, the equations for ϕ (6) and (7) can be respectively rewritten as

In critical collapse, the period in terms of the coordinate time t exponentially decreases. Consequently, the metric functions and matter field in the late stage of collapse and large-radius region for which $ |t/r|\ll 1 $ appear to be "frozen" rather than propagating [20, 21]. In Ref. [17], the authors offered one ansatz: in this region, the last terms in Eqs. (10) and (11) are negligible in comparison with the first terms. Moreover, treating α and β as constant, the authors obtained the following solutions:

where

The expressions (12) are consistent with the numerical results. However, some treatments in the above have not been fully justified. In addition, although approximate expressions for X and Y were obtained, the results for the metric functions and scalar field remain absent. We address such issues below.

In Ref. [17], some terms in Eqs. (10) and (11), that is, $ -u(\alpha Y)_{,u} $, $ -u(\alpha X)_{,u} $, $ \beta_{,u}X $, $ \alpha_{,\rho}Y $, $ \beta_{,u}Y $ , and $ \alpha_{,\rho}X $, were dropped. In fact, as shown in Figs. 1 and 2, in the large-radius region $ (r>10^{-3}) $, the absolute values of the terms $ -u\alpha_{,u}Y $, $ -u\alpha_{,u}X $, $ \alpha_{,\rho}Y $ , and $ \alpha_{,\rho}X $ can sometimes be greater than those of other terms. However, the terms dropped approximately cancel. Consequently, the equations constructed by the remaining terms roughly hold:

Figure 1.  (color online) Numerical results for Eqs. (10) and (14). (a) Results for Eq. (10). (b) Results for Eq. (14). In the large-radius region $ (r>10^{-3}) $, the absolute values of $ \alpha_{,\rho}Y $ and $ u\alpha_{,u}Y $ are sometimes greater than those of other terms in Eq. (10). However, the three terms $ \beta_{,u}X $, $ -\alpha_{,\rho}Y $ and $ u(\alpha Y)_{,u} $ roughly cancel. Then, Eq. (14) approximately holds.

Figure 2.  (color online) Numerical results for Eqs. (11) and (15). (a) Results for Eq. (11). (b) Results for Eq. (15).

From this perspective, treating α and β as constant in Ref. [17] is effectively valid.

In the analytic results (12) obtained in Ref. [17], only the phases in the sine functions are functions of r and t, and the amplitudes are constant. We examine the numerical results for ϕ and find that besides the phase, the amplitude for ϕ is also a function of r and t. We find that the field ϕ admits the following approximate expression:

The quantity $ [H(r,t)] $ has the following features:

(i) We first define $ H(r,t) $ in the usual manner:

Then, we have

(ii) $ [H(r,t)] $ is defined as

Here, $ \bar{\omega} $ and $ \bar{\alpha} $ are constant values of $ \omega(r,t) $ and $ \alpha(r,t) $. This is related to the treatment of α and β as constant in the discussions of Eqs. (14) and (15).

(iii) From Eq. (16), we obtain the expression for $ \phi_{,t} $,

where $ \tan\varphi_{1}\equiv C_{3}/C_{2} $. Regarding the quantity $ H_{,t}(=\omega\alpha/r+ \omega\alpha_{,t}t/r+\omega_{,t}\alpha t/r) $, the numerical results show that $ |{\omega\alpha_{,t}t}/{r}| $ is sometimes greater than $ {\omega\alpha}/{r} $. However, comparing expression (20) with the numerical results for $ \phi_{,t} $, we always obtain

This implies that in $ [H]_{,t} $ , the contributions from $ \omega\alpha_{,t}t/r $ and $ \omega_{,t}\alpha t/r $ are negligible. This should be related to the fact that the respective reductions of Eqs. (10) and (11) to Eqs. (14) and (15) are equivalent to treating α and $ \beta(\equiv \sqrt{1+\omega^{-2}}) $ as constant.

(iv) The numerical results in Fig. 3(a) show that in the large-radius region, the equation of motion for ϕ (7) is reduced to

Figure 3.  (color online) Numerical results for the equation of motion for ϕ (7) and the transition region. (a) Results for Eq. (7). In the large-radius region $ (r>10^{-3}) $, $ A^{-1}{\rm e}^{\delta}\phi_{,tt}\approx-(A^{-1}{\rm e}^{\delta})_{,t}\phi_{,t} $. This is very different from the flat-spacetime form $ \phi_{,tt}=r^{-2}(r^{2}\phi_{,r})_{,r} $. In the central region $ (r<10^{-5}) $, $A^{-1}e^{\delta}\phi_{,tt}\approx A {\rm e}^{-\delta}r^{-2}(r^{2}\phi_{,r})_{,r}$. Considering that in this region $ A\approx1 $ and $ \delta\approx 0 $, $ \phi_{,tt}\approx r^{-2}(r^{2}\phi_{,r})_{,r} $. (b) Transition region between the central and large-radius regions:$ [r_1,\; r_2] $. At $ r=r_{1} $, $ |C_{3}H|\sim|\omega\ln r| $. At $ r=r_2 $, $ |C_{3}H_{,r}|\sim\omega/r $.

Using Eq. (22) and the numerical results of $ |\delta_{,t}|\gg |A_{,t}| $, we have

Combining Eqs. (20), (21), and (23) and the numerical results of $ |\delta_{,t}|\gg H_{,t} $ generates

Namely, we treat α effectively as constant at the first-order accuracy, and the dynamical feature of α begins to take effect since the second-order temporal derivative of $ [H(r,t)] $.

(v) As shown in Fig. 3(b), at the late stage of critical collapse, in the large-radius region at which $ |t/r|\ll 1 $, $ |H|\ll 1 $, $ |H|\ll|\omega\ln r| $, and $ |H_{,r}|\ll 1/r $. Therefore, with Eq. (16), $ [H] $ mainly contributes to the temporal derivatives of ϕ, rather than to the field ϕ and its spatial derivatives.

The numerical results show that $ C_{1}\approx0.058 $, $ C_{2}^2+C_{3}^2\approx 1 $, and $ \varphi_{1}\approx 1.08 $. As shown in Figs. 4(a) and (b), the expressions for ϕ (16), $ \phi_{,t} $ (20), and $ \phi_{,tt} $ (23) agree well with the numerical results. Note that when we compare the analytic expressions with the numerical results in the large-radius region, the tail part, e.g., the region for $ r>1 $ in Fig. 4, is excluded.

Figure 4.  (color online) Comparison of the analytic and numerical results for the scalar field and metric functions in the large-radius region. (a) Results for ϕ (16) and $ \phi_{,t} $ (20). (b) $ \phi_{,tt}\approx-\delta_{,t}\phi_{,t} $. (c) Results for $ m/r $ (27). (d) Results for δ (30).

With Eqs. (16) and (20), we can rewrite Eq. (4) as

where $ C_4=4\pi\omega{C_{1}^2}\sqrt{C_{2}^2+C_{3}^2} $. Via integration, we have

Then, using Eq. (8) and the fact that $ |H|\ll 1 $, we obtain

where $C_6\approx {\rm e}^{C_5}C_{4}/(4C_3)\approx {\rm e}^{C_5}\pi\omega{C_{1}^2}\sqrt{C_{2}^2+C_{3}^2}/C_{3}$ and $C_7= (1/2)(1-{\rm e}^{C_{5}})$. As shown in Fig. 4(c), the expression for $ m/r $ (27) is consistent with the numerical results. The fitting results are $ C_6=0.013360\pm 0.000009\approx 1/75 $ and $ C_7\approx 0.065480\pm 0.000007\approx 1/15 $.

With Eq. (8), we can rewrite Eqs. (3) and (5) as

Then, the solution for δ can be expressed as

where $ C_8\approx-{2C_7}/(1-2C_7)-2C_{6}^2 $ and $C_9\approx-(C_6+ 8C_{6}C_{7})/\omega$. As shown in Fig. 4(d), the expression for δ (30) is consistent with the numerical results.

In Ref. [17], the quantities α and β were treated as constant. The approximate expressions for X and Y obtained in this manner agree well with the numerical results. Then, it was stated that in this circumstance, the spacetime is effectively flat. In fact, X and Y are combinations of the metric functions and derivatives of the scalar field, rather than the scalar field. To check whether the spacetime is effectively flat, it may be more appropriate to directly investigate the behavior of the equation of motion for the scalar field (7). As shown in Fig. 3(a), in the large-radius region, Eq. (7) is reduced to Eq. (22), which is clearly different from the flat-spacetime form $ \phi_{,tt}= r^{-2}(r^{2}\phi_{,r})_{,r} $. Therefore, the spacetime in this region is not effectively flat.

As shown in Fig. 3(a), in the central region, the absolute values of the terms $(A^{-1}{\rm e}^{\delta})_{,t}\phi_{,t}$ and $(A^{-1}{\rm e}^{\delta})_{,r}\phi_{,r}$ in Eq. (7) are considerably lower than those of $A^{-1}{\rm e}^{\delta}\phi_{,tt}$, $A{\rm e}^{-\delta}\phi_{,rr}$, and $(2/r)A{\rm e}^{-\delta}\phi_{,r}$. Moreover, in this region, $ A\approx1 $ and $ \delta\approx 0 $. Consequently, Eq. (7) is reduced to the flat-spacetime form

Regarding Eq. (31), we discuss the following:

(i) Equation (31) implies that in the central region, the scalar field ϕ evolves almost as in flat spacetime, not experiencing gravitational effects.

(ii) Besides critical collapse, we check the evolution of the scalar field in two other types of collapse (dispersion and black hole formation) and obtain similar results as (31).

(iii) The result (31) is closely related to the asymptotic behaviors of the metric functions and scalar field near the center. Under the smoothness requirement at the center, the metric functions and scalar field have the following power series expansions near the center [19]:

With Eqs. (3), (4), (5), (8), and (32), we obtain the following asymptotic expressions:

which are also shown in Fig. 5. With Eqs. (32) and (33), we can straightforwardly simplify Eq. (7) to (31).

Figure 5.  (color online) First-order temporal and spacial derivatives for A, δ, and ϕ. As discussed in Eq. (33), near the center, $ A_{,t}\propto r^2 $, $ \delta_{,t}\propto r^2 $, $ \phi_{,t}\approx \phi_{0}'(t)+\phi_{2}'(t)r^2 $, $ A_{,r}\propto r $, $ \delta_{,r}\propto r $, and $ \phi_{,r}\propto r $.

(iv) It is known that in critical collapse, the Ricci curvature scalar R in the central region is very high and will eventually diverge. This does not contradict the result in (31) because Eq. (31) is caused by the fact that the first-order derivatives of the metric components asymptote to zero. Conversely, as discussed below, the major terms constructing R include the second-order derivatives of the metric components, first-order derivatives of the metric components divided by r, and $ 2(1-A) $ divided by $ r^2 $, which are very large and will eventually diverge.

For the metric (1), the Ricci curvature scalar can be written as

With Eqs. (3), (4), (5), (8), (32), and (33), we obtain asymptotic expressions for all the terms on the right-hand side of Eq. (34):

The first five terms are dominant and have the same order of magnitude as $ 8\pi(\phi_{,t})^2 $, and the remaining terms are proportional to $ r^2 $ or $ r^4 $ and are negligible.

As shown in Fig. 3(b), the transition region between the central and large-radius regions can be expressed as $ r\in [r_1,\; r_2] $. At $ r=r_{1} $, $ |C_{3}H|\sim|\omega\ln r| $; at $ r=r_2 $, $ |C_{3}H_{,r}|\sim \omega/r $.

Analytic solutions are important for understanding the dynamics of gravitational collapse. Owing to the complexity of the Einstein equations, seeking the analytic solutions to the equations has been a very challenging issue. In successful circumstances, the equations are usually reduced to ODEs. In critical collapse, the equations remain as PDEs, whereas in the large-radius region and late stage of evolution, the spatial and temporal contributions are separate to some extent. This enables us to obtain approximate analytic expressions for the metric functions and matter field.

The boundary conditions at the center play a key role in the dynamics in the central region. In this region, owing to the boundary conditions, the terms related to gravitational effects in the equation of motion for the scalar field are negligible, such that the equation is reduced to the flat-spacetime form.

The authors are very thankful to Xiao-Kai He, Junbin Li, and Cheng-Yong Zhang for helpful discussions.