Dynamics of critical collapse

Critical collapse of a massless scalar field in spherical symmetry is systematically studied. We combine numerical simulations and asymptotic analysis, and synthesize critical collapse, spacetime singularities, and complex science. First set of approximate analytic expressions near the center are obtained. We observe that, near the center, the spacetime is nearly conformally flat, the dynamics is not described by the Kasner solution, and the Kreschmann scalar is proportional to r-5.30\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{-5.30}$$\end{document}, where r is the areal radius. These features are significantly different from those in black hole singularities. It is speculated that the scalar field in critical collapse may be a special standing wave.


Introduction
Complex systems widely exist in Nature. The interactions between the parts, alternatively the nonlinearities, in a complex system make the system more than the sum of its parts, causing the emergence of the intriguing critical phenomena and discrete scale invariance in some circumstances [1,2]. The critical phenomena in gravitational collapse caused by the nonlinearities of Einstein equations were originally discovered by Choptuik via numerical simulations [3]. A naked singularity is formed in critical collapse [4]. So critical collapse connects two basic fields in science: critical phenomena and spacetime singularities.
Choptuik simulated critical collapse of a massless scalar field in spherical symmetry in general relativity. When the scalar field is weak, the scalar field will implode and then disperse, leaving a flat spacetime behind; while when the scalar field is strong enough, it will implode and form a black hole (BH). By fine tuning the strength of the scalar field, a critical solution which is on the threshold of BH formation can be obtained. The results for critical collapse include discrete a e-mail: sps_guojq@ujn.edu.cn b e-mail: sps_zhanghs@ujn.edu.cn self-similarity (DSS), universality, and mass scaling law. The DSS states that the scalar field oscillates periodically at everdecreasing time and length scales related by a factor of e with ≈ 3.44. Regarding the universality, all the families of the near-critical evolutions approach the same solution. The mass scaling law describes that, in the super-critical case, the mass of the tiny BH has a scaling relation with the parameter of the scalar field, M BH B| p − p * | γ , where B is a familydependent parameter, p is one parameter for the scalar field, describing the strength of the scalar field, p * is the critical value for p, and γ is a universal scaling exponent γ ≈ 0.37.
Considering that the results on critical collapse by Choptuik are numerical, it is natural to study this issue analytically. The existence of a real analytic solution to critical collapse was proved in Ref. [20]. In Refs. [21][22][23], under the requirements of discrete scale invariance, analyticity, and an additional reflection-type symmetry, the critical collapse is reduced to an eigenvalue problem. The rescaling factor becomes an eigenvalue and is solved numerically by a relaxation algorithm with high precision. Motivated by understanding the nature of the echoing of the scalar field, Price and Pullin found that although the oscillatory behavior of the scalar field seems to come from the nonlinearities of general relativity, it can be approximated by a scalar field solution in flat spacetime [24]. In addition, analytic models for the continuous self-similar collapse in spherical symmetry (Roberts solution) [25][26][27] and cylindrical symmetry [28] were presented. In Ref. [29], with analytic perturbation methods, it was found that a generic perturbation departs from the Roberts solution in a universal way. For reviews on critical collapse, see Refs. [30][31][32].
Despite the above efforts on analytic studies, due to the complexity of Einstein equations, ever since the publication of Choptuik's numerical results in 1993, the analytic solution remains unknown, and the nature of critical collapse is far from being well understood. In this paper, we explore the dynamics of critical collapse by combining numerical simulations and asymptotic analysis. Considering that critical collapse ends up with a naked singularity, we connect critical collapse to the results on the dynamics near spacetime singularities that have been obtained before, including the Belinskii, Khalatnikov, and Lifshitz (BKL) conjecture and BH formation. With such efforts, some approximate analytic solutions near the center are obtained. It is found that, near the center, the spacetime is nearly conformally flat, and the dynamics is not described by the Kasner solution.
The paper is organized as follows. In Sect. 2, we describe the framework. Section 3 presents the nearly conformal flatness of the spacetime in critical collapse. In Sect. 4, approximate analytic information on critical collapse is reported. We compare critical collapse with the BKL conjecture and BH formation in Sect. 5. In Sect. 6, the connection between the scalar field in critical collapse and standing waves is argued. In Sect. 7, the results are discussed. Throughout the paper, we set G = c =h = 1.

Framework
In this section, we present the framework for numerical simulations of critical collapse. Critical collapse of a massless scalar field in spherical symmetry is simulated in double-null coordinates, where u = (t − x)/2 and v = (t + x)/2. Consider a massless scalar field ψ with the energy-momentum tensor T μν = ψ ,μ ψ ,ν − (1/2)g μν g αβ ψ ,α ψ ,β . Then the equations of motion can be written as [33,34] The constraint equations are [33,34] r ,t x + r ,t σ ,x + r ,x σ ,t + 4πr ψ ,t ψ ,x = 0, r ,tt + r ,x x + 2r ,t σ ,t + 2r ,x σ ,x + 4πr (ψ 2 ,t + ψ 2 ,x ) = 0. (6) In this paper, r ,t is defined as r ,t ≡ ∂r (t, x)/∂t, and other quantities, e.g., r ,x , r ,tt , etc, are defined analogously. For numerical stability concern, the Misner-Sharp mass m is used as an auxiliary variable as has been successfully implemented in Ref. [35], Then Eqs. (2) and (3) can be rewritten as The dynamics of m is described by The quantity m ,x is The derivations of Eqs. (10) and (11) are given in the Appendix. In the simulations, Eqs. (4), (8), (9), and (10) are integrated numerically with finite difference method and leapfrog integration scheme. We impose r ,tt = r ,t = σ ,t = ψ ,t = 0 at t = 0. The initial profile of ψ is set as with a being tuned as a = 0.0908379681, b = 0.01, and x 0 = 0.25. We set r = σ = m = 0 and r ,x = 1 at the origin (x = 0, t = 0). The expressions for σ ,x , r ,x x , and m ,x can be obtained from Eqs. (6), (8), and (11), respectively. We obtain the quantities r , σ , and m on the initial slice of t = 0 by numerically integrating such expressions from x = 0 to x = 2 with the fourth-order Runge-Kutta method. From Eqs. (4), (8), and (9), the values of ψ, r , and σ at t = t can be determined via a second-order Taylor expansion. Take the first-order time derivative of Eq. (10), one can obtain m ,tt . With m ,t and m ,tt , m at t = t is obtained via a second-order Taylor series expansion.
Regarding the boundary conditions, we always set r = m ≡ 0 at x = 0. Then there are always r ,t = r ,tt ≡ 0 at x = 0. From Eqs. (4) and (8), one obtains respectively ψ ,x ≡ 0 and r ,x x ≡ 0 at x = 0. Then with Eq. (6), there is The simulations of critical collapse need to be highly accurate. In order to achieve this objective, adaptive or fixed mesh refinement techniques are usually used. While in this paper, we take an even simpler approach. From the beginning, we use very small spatial and temporal grid spacings , not making any mesh refinement throughout the whole simulations. It turns out that the numerical results by this approach can show the basic features of critical collapse and are adequate for us to study the dynamics of critical collapse. The code is second-order convergent.

The spacetime near the center is nearly conformally flat
The numerical results for the evolution of r , σ , m, and ψ are shown in Fig. 1. Near the center, σ goes to +∞. The scalar field ψ oscillates, and at the same time, moves toward the center under gravity. It is noticeable that there are at least two basic parameters for a wave: period and amplitude. However, the amplitude of the scalar field has not been paid as much attention as the period has in the literature in the past. As shown in Figs. 1e and 4b, the amplitude is about 0.61. We note that, in Choptuik's original work, where the polarslicing coordinates were used, the amplitude is about 0.45 [9]. In the work by Akbarian and Choptuik [36] and the work by Baumgarte [37], where the Baumgarte-Shapiro-Shibata-Nakamura formulation of Einstein's equations in spherical polar coordinates, the 1 + log slicing condition for the lapse and the Gamma-driver condition for the shift were implemented, the amplitude is about 0.61. One may obtain additional hints on the nature of and possible analytic solution Noting that the dynamics is simply determined by the equations, we examine the behavior of each term in each equation by plotting the numerical results. We investigate the contribution of each term in the equation of motion for r (2) on the slice of x = 4.75 × 10 −4 . Remarkably, as shown in Fig. 2a, near the center, ever since the beginning of the collapse, there is always e −σ ≈ r ,x . Besides, analytic results means relations between different quantities, so we check the relations between quantities by plotting the corresponding numerical results, e.g., r ,x vs. r . We plot r ,x vs. r on the same slice, and find that Substitution of Eq. (13) into Eq. (1) yields So, near the center, the spacetime is nearly conformally flat (NCF). Some comments on the nearly conformal flatness feature are made as follows: (i) We also simulate collapse toward BH formation and flat spacetime (dispersion), and find that Eq. (13) is also valid in these two cases. So the NCF region also exists near the center in BH formation and dispersion cases. (ii) In Ref. [38], the dynamics near the center in axisymmetric system was studied in detail. The local flatness near the center yields the conformal flatness for the spatial part of the metric. This feature and the homogeneity property [38] will generate the nearly conformal flatness for the spacetime near the center. As shown in Fig. 1b, near the center, for slice t = Constant, there are σ ≈ Constant and r (t, x) ∝ x, which should be related to the boundary conditions of r ,x x = σ ,x ≡ 0 at x = 0, see Sect. 2 and Fig. 3. In the Oppenheimer-Snyder spherical dust collapse, the density ρ is uniform, the pressure is zero, and the energymomentum tensor is T μ ν = ρ(−1, 0, 0, 0). The interior spacetime is described by the closed Friedmann-Lemaître-Robertson-Walker metric [39,40]. Obviously, the spacetime is conformally flat. In spherical scalar collapse, as discussed in Sect. 2, near the center, ψ ,x = 0. So near the center, the energy-momentum tensor becomes T μ ν ≈ (1/2)e 2σ ψ 2 ,t (−1, 1, 1, 1) and is almost uniform. Therefore, spherical scalar collapse near the center is similar to the Oppenheimer-Snyder spherical dust collapse, which naturally leads to the nearly conformal flatness feature near the center. (iii) The nearly conformal flatness means that the Weyl tensor C αβμν , the scalar C W ≡ C αβμν C αβμν , and the tidal force are nearly zero. Then for the metric (1), using Eq. (2), there is (iv) Using σ ≈ − ln (r/x) and r (t, x) ≈ D(t)x, one obtains Then with e −σ ≈ r ,x , the first line of Eq. (15) leads to C W ≈ 0, verifying that the spacetime in which e −σ ≈ r/x is indeed nearly conformally flat. Moreover, Eqs. (17) and (18) imply that, for the dynamics near the center, the coordinates t and x are not symmetric. (v) We also compute the Weyl tensor in the Roberts solution, which is one type of continuously self-similar critical collapse of a scalar field [25][26][27], and obtain nonzero results. The solution can be written as below [26]: where k is a parameter. For this solution, which is nonzero for non-Minkowskian spacetime.

Approximate analytic expressions
In this section, we present partial analytic information obtained via combination of numerical simulations and asymptotic analysis. In Fig. 4a, we plot −σ vs. − ln τ on the slice of x = 4.75 × 10 −4 . τ is the proper time, defined as τ ≡ ξ 0 e −σ dξ , where ξ ≡ t − t 0 and t 0 is the coordinate time when r goes to zero. As shown in Fig. 4a, −σ is a superposition of a linear function and a periodic function of − ln τ . In fact, in many discrete scale invariance systems, the log-periodic oscillations exist in the time dependence of the energy release as the impending rupture is approached. The typical time-to-failure formula is [1,2] E ∼ (t r − t) n 1 + β cos 2π where E is the energy released or some other variable describing the on-going damage, t r is the rupture time, n is a critical exponent, λ is a preferred scaling factor, and ϕ is a phase. The critical collapse system is also a discrete scale invariance system. So we fit the numerical results of −σ vs. − ln τ according to the logarithm of Eq. (23), Noting that the period for σ in Fig. 4a is T σ ≈ 1.795, we fix the angular frequency ω in Eq. (24) as ω = 2π/T σ = 3.50.
As shown in Fig. 4a, the numerical results can be well fit by this formula. The close agreement between the numerical results and the fitting formula strongly implies that critical collapse is indeed one more subfield of complex science. With the fitting result (24) for σ = σ (τ ) and ξ = τ 0 e σ dτ , one can obtain the approximate analytic expression for ξ = ξ(τ ). The fitting result for a in Eq. (24) is −0.5028 ± 0.0015. For simplicity, we set a = −1/2. In addition, let y = ω ln τ − c and use Taylor series expansion, there is where H = 2b/[(2 + b 2 ) √ 1 + 4ω 2 ] ≈ 0.03, and δ = arccos[2ω/ √ 1 + 4ω 2 ] ≈ 0.14 radians. Noting that H ≈ 0.03 1, the major part for the inverse function τ = τ (ξ) is (26) Substitution of Eq. (26) into the second term in the square brackets in Eq. (25) yields an approximate analytic expression for τ = τ (ξ) As shown in Fig. 6, Eq. (28) matches well with the numerical results for ξ < 0.15, which corresponds to formal critical collapse. There is a large deviation between the analytic expression (28) and the numerical results for ξ > 0.15, which should be because that this region is not a formal critical collapse region, but an entrance to critical collapse. Also see Figs. 2 and 4. The numerical results on ψ vs. − ln τ are plotted in Fig. 4b. Roughly speaking, ψ is a linear function of ln τ and switches between α ln τ and −α ln τ after every half-period with α ≈: 0.6 ∼ 0.8. The numerical results on ψ vs. − ln τ are fit according to fourth-order Fourier sine and cosine series. The fitting results are shown in the caption of Fig. 4. In Sect. 5, we will give one preliminary analytic explanation to the slopes of (−σ, ψ) vs. − ln τ .

Critical collapse vs. dynamics near spacetime singularities
A naked singularity is formed in critical collapse. Therefore, it is meaningful to connect critical collapse to the results on the dynamics near spacetime singularities that have been obtained before, including the BKL conjecture and BH formation. We will do so in this section.

Critical collapse vs. BH formation
In Table 1, we compare the dynamics of critical collapse near the center and that at x = 0 along the singularity curve in BH formation. Some discussions are given below. In BH formation, the scalar field and gravity are strong. Then, the scalar field is simply absorbed into the central singularity without reflections, ψ = α ln τ or ψ = −α ln τ with α ≈ 0.16, and σ = −(1/3) ln τ . In critical collapse, the scalar field and gravity are weak, and there is a balance between gravity and reflections at the center. Consequently, ψ switches between α ln τ and −α ln τ with α ≈: 0.6 ∼ 0.8, and σ is a superposition of a linear function and a periodic function of ln τ with the slope for the linear part being about −1/2.
The slopes of (σ, ψ) vs. ln τ in critical collapse are at the same order of magnitude as those in BH formation, respectively. However, as shown in Eqs. (34) and (35), the values of the slopes in BH formation can be obtained analytically. In this sense, for the first time, we give one rough analytic explanation to the slopes of (σ, ψ) vs. ln τ in critical collapse.
The presence of the cosine function in the second line of Eq. (28) implies that the dynamics in critical collapse is not described by the Kasner solution. Therefore, the BKL conjecture is not valid for the naked singularity formed in critical collapse, which is different from BH singularities.
As shown in Eq. (36), in BH formation, in the case of C = ± √ 3/2, the Kreschmann scalar K ∝ r −12 . We plot the Kreschmann scalar in critical collapse in Fig. 8, and observe that K ∝ r −5.30 , much weaker than the BH formation case, and also weaker than the Schwarzschild BH case where K ∝ r −6 .
In summary, critical collapse and BH formation share some common features: after all they both are on the dynamics in the vicinities of singularities. On the other hand, there are some remarkable differences between the two: the scalar field and gravity are weaker in critical collapse than in BH formation.

Scalar field in critical collapse vs. standing waves
In critical collapse, the scalar field moves toward the center under gravity. At the same time, the scalar field is reflected at the center. The resulting scalar field is discretely self-similar. This picture makes us to speculate that the scalar field in critical collapse may be a special standing wave. The word 'special' is due to the fact that, in critical collapse, the scalar wave keeps shrinking toward the center. We think that this  ψ is reflected at x = 0, and switches between α ln τ and −α ln τ with α ≈: 0.6 ∼ 0.8 after each half-period.
The spacetime is nearly conformally flat.
The spacetime is nearly conformally flat.
The BKL conjecture applies.
The BKL conjecture does not apply.
Kreschmann scalar ∝ r −12 Kreschmann scalar ∝ r −5.30 special standing wave picture helps to interpret the DSS feature and to seek the full analytic solution to critical collapse. In fact, standing waves are quite common in astrophysics: e.g., (i) The 5-minute solar oscillation in local surface velocities discovered by Leighton [53,54]. (ii) The density wave theory on the spiral structure of disk galaxies by Lin and Shu [55][56][57][58][59]. (iii) The microwave cavity hypothesis on ball lighting by Kapitsa [60,61].

Discussions
We obtained the first approximate analytic expression for the metric near the center in critical collapse. We observed that in spherical scalar collapse, the spacetime near the center is nearly conformally flat. The interplay between critical collapse and complex science deserves to be explored fur-ther, which is meaningful for both sides. In critical collapse, gravity and reflections at the center compete and compromise. Similar mechanisms widely exist in other branches of complex science [62]. The BKL conjecture may still be a guiding principle to achieve deeper understanding on critical collapse.
Critical collapse is an interdisciplinary subject, connecting gravitation (including BH physics and spacetime singularities), complex science, and partial differential equations, etc. Studies in critical collapse can enrich each of these, and deserve further efforts.