On Critical Exponents for Self-Similar Collapse

We explore systematically perturbations of self-similar solutions to the Einstein-axion-dilaton system, whose dynamics are invariant under spacetime dilations combined with internal $SL(2,R)$ transformations. The self-similar solutions capture the enticing behavior critical systems on the verge of gravitational collapse, in arbitrary spacetime dimensions. Our methods rest on a combination of analytical and numerical tools, apply to all three conjugacy classes of $SL(2,R)$ transformations and allow accurate estimates of the corresponding Choptuik exponents. It is well known that these exponents depend on the spacetime dimension and on the matter content. Our main result is that they also attain different values, even within a given conjugacy class, for the distinct types of critical solutions that we recently identified in the Einstein-axion-dilaton system.


Introduction
The study of spherically-symmetric collapse of scalar fields started in [1] and the emergence of a critical behavior in gravitational collapse was first considered by 1 arXiv:1912.06103v3 [hep-th] 13 Feb 2022 Choptuik in [2], who showed that solutions on the verge of gravitational collapse feature some kind of spacetime self-similarity under special dilations. In addition, near-critical solutions display suggestive scaling laws. The example of [2] dealt with a minimally-coupled real scalar with initial conditions modulated by a parameter p related to the field amplitude. If the choice p = p crit results in the critical solution, super-critical initial conditions, which obtain for p > p crit , exhibit scaling laws for the Schwarzschild radius and hence for the mass of the final black hole, of the form where numerically one finds for the Choptuik exponent the value γ 0.37 [2].
Let us stress, however, that in a generic dimension d ≥ 4 the definition must be generalized as follows [5,6] Various numerical experiments with different choices of matter content can be performed along the same lines. For instance [7,8,9,10] found more types of critical solutions for a massless scalar field. The case of a complex scalar field was investigated in detail in [11], and the critical collapse of a radiation fluid was considered in the interesting works [12,13,5,14]. For spherically symmetric radiation fluids [13] obtained a very similar result, γ 0.36, which led to the natural conjecture [15] that γ ought to be universal for four-dimensional gravity coupled to any matter content.
It was then understood [5,14,16] that the Choptuik exponent can be extracted more conveniently relying on perturbations of the self-similar solutions. To this end, one perturbs a self-similar field h 0 according to where the perturbation h −κ has scaling dimension −κ ∈ C, to then determine the values of κ corresponding to the allowed modes. The most relevant mode κ * , for which Re(κ) attains the maximum value, is related to the Choptuik exponent according to This approach allowed for more precise determinations of the Choptuik exponents. For the radiation fluid, one thus arrived [5] at a better estimate, γ 0.3558, which clearly differs from the value that was originally obtained for the real scalar field. The authors of [16] then determined the Choptuik exponent for a complex scalar field, obtaining the value γ 0.3871, which is even more difficult to reconcile with the original result for a real scalar field. Therefore, the original hope that the value of γ reflect some sort of universality, independently of the matter content, was not fulfilled.
The extension to the case of axial symmetry was then considered in [17], while the authors of [18] examined critical collapse in the presence of gravitational shock waves. Non-linear sigma models provide an interesting class of slightly more complex systems. Specifically, [19] investigated critical collapse in four dimensions in the presence of sigma models related to the two-sphere and the hyperbolic plane H 2 , whose critical solutions are self-similar up to compensating U (1) internal rotations. The H 2 sigma model is particularly interesting, since it describes the coupling of an axion-dilaton system to gravity. The corresponding Choptuik exponent, γ 0.2641, is significantly lower.
This paper is devoted to the study of continuously self-similar solutions of the Einstein-axion-dilaton system in arbitrary dimensions. The internal SL(2, R) symmetry of the scalar sector has an interesting structure, due to the simultaneous presence of compact and non-compact one-parameter subgroups, and we thus generalize previous discussions considering all possible classes of continuously self-similar collapses. In this fashion, we discover several families of critical spacetimes in various dimensions. In a previous work [20] we have studied selfsimilar solutions to the Einstein-axion-dilaton system in several dimensions, and with the compensating internal transformation in the three conjugacy classes of sl(2, R), generalizing [21,22]. Surprisingly, we discovered that even for a given choice of dimension, matter content and conjugacy class of compensating transformation, entire families of physically-distinct self-similar solutions exist.
For example, the axion-dilaton with U (1)-compensated dilations (i.e., "elliptic" class), which was first considered in [19] in four dimensions, exhibits three distinct collapse solutions in five dimensions.
In this work we perturb the different classes of solutions discovered in [20], providing accurate estimates for the Choptuik exponents. Our main conclusion is that the Choptuik exponents depend not only on the choices of matter Lagrangian and spacetime dimension, but they also attain different values for the different critical points that may arise within a given system. We could reach these conclusions for the three conjugacy classes relying on a robust analytic and numerical formalism for linear perturbations of self-similar solutions in arbitrary dimensions, which we developed here and in [20]. We are thus able to recover the known value [19] of γ for the unique four-dimensional elliptic critical collapse, within numerical accuracy, which supports our confidence in the current setup.

Self-similar Einstein-axion-dilaton solutions
We begin by summarizing, in this section, the main results of [20]. Our discussion starts with a description of the Einstein-axion-dilaton system, to then review the concept of continuously self-similar (CSS) configurations, the procedure for determining CSS solutions and the detailed solutions found in four and five dimensions in [20], which provide the backgrounds for the perturbations examined here.

Einstein-axion-dilaton system and self-similar fields
The Einstein-axion-dilaton system in d spacetime dimensions (d ≥ 4) consists of gravity coupled to an axion a and a dilaton φ, and is described by the action which is the same as the low energy effective action of type II string theory [23,24] where the axion-dilaton complex field τ ≡ a + ie −φ . Note that the SL(2, R) symmetry reduces to the SL(2, Z) subgroup once quantum effects are considered, and this S-duality is conjectured to be a non-perturbative symmetry of IIB strings [25,26,27]. From the action (6) one obtains the equations of In [20] and in the present paper we will only consider spherically symmetric configurations and perturbations 1 . The general spherically symmetric metric in d dimensions can be cast in the form where q ≡ d − 2.
The assumption of continuous scale invariance for the metric translates into a homogeneous scaling of the line element under dilations, as follows. If (t, r) → (Λt, Λr) , Λ > 0 then which implies the metric functions are themselves scale-invariant, so that where z ≡ −r/t is a scale-invariant coordinate. The scalar τ , rather than being directly scale-invariant, can be invariant, more generally, up to an SL(2, R) transformation, so that We term a configuration (g, τ ) satisfying eqts. (13), (14) to be continuously 1 We suspect it should be possible to prove that non-spherical perturbation modes (either with angular dependence, or with polarisation indices on the sphere) never include the most relevant modes. Intuitively, just as in time-translation invariant systems, the modes with angular momentum have increased energy, we expect that in a scale-invariant system angular momentum necessarily increases − Re κ. If this is indeed the case, then for the sake of determining the critical exponent spherically-symmetric perturbations are sufficient.

self-similar (CSS).
Physically distinct cases correspond to the different conjugacy classes of dM dΛ Λ=1 . For the three cases of elliptic, parabolic and hyperbolic elements one can parametrize the CSS axion-dilaton field as , elliptic in terms of an unknown real parameter ω and a scale-invariant field f (z), which satisfies |f (z)| < 1 in the elliptic case and Im f (z) > 0 in the parabolic and hyperbolic cases. Note that ω can be assumed to be positive as the SL(2, R) symmetry τ → − 1 τ leaves the form of the ansätze (15) invariant, but changes the sign of ω.
At this point, making use of the CSS ansätze in the equations of motion (7), (8), one can recover the CSS equations u(z), b(z), f (z). In fact, u(z) can be eliminated in favour of b(z) and f (z), and after some simplifications one can construct a final system of ordinary differential equations (ODEs) These equations have several singular points; the relevant range for z, which contains the infinite past, lies between the two singularities In particular, z = z + is a null surface and is in fact a type of event horizon (homothetic horizon). This surface should be merely a coordinate singularity, and thus the axion-dilaton field τ should be regular across it. This condition is equivalent, in fact, to demanding that f (z) remain finite as z → z + ; the vanishing of the divergent part of f (z) is a complex-valued constraint at z + For what concerns z = 0, using the regularity of τ and residual symmetries in the equations of motions, one can determine the initial conditions where x 0 is a real parameter. At this point, one is left with two unknown parameters (ω, x 0 ) and two constraints (the real and imaginary parts of G). The system is thus completely determined, and will generically feature a discrete solution set. In [20] we investigated this system numerically in four and five dimensions, finding several CSS solutions that we now review briefly.

Solutions in four and five dimensions
In this section we review the CSS solutions in four and five dimensions for the three conjugacy classes of sl(2, R) found 2 in [20]. They are determined by the values of the parameters ω and x 0 , as explained in Section 2.1, and we also note the position of the z + singularity. The actual profiles of these solutions can be readily reproduced integrating numerically the CSS equations of motion.
Notice that multiple distinct solutions were found in [20] for some specific choices of dimensions and conjugacy classes. The set of 11 CSS solutions in four and five dimensions discovered in [20] corresponds to the following parameters: 2 No solutions were identified in the parabolic class, within our numerical precision. Let us also note that in [20] we tentatively identified at least a fourth solution δ in the four-dimensional hyperbolic class. Since their numerical accuracy was very limited already in the unperturbed solutions, we exclude these cases from our current investigation. Moreover, we were unable to ascertain whether additional hyperbolic solutions were present in a region essentially inaccessible due to numerical errors. [20] contains a more precise description of the solution sets in parameter space, along with graphical representations.

Perturbation theory
In this section we would like to outline a comprehensive method that can allow one to investigate perturbations in general numbers of dimensions and for any matter content. We apply this formalism to the axion-dilaton system and identify linearized perturbations of the equations of motions for elliptic and hyperbolic ansätze. We shall follow on the lines of previous discussion of the perturbation theory of self-similar backgrounds [28] 3 . We will, however, improve upon the methodology by implementing several algebraic manipulations for the purpose of simplifying an otherwise extremely cumbersome and opaque computation.
We emphasize that by means of the following method we are able to reconstruct the well-known [30] result for the Choptuik exponent for the four-dimensional elliptic solution, namely We shall consider this as a strong indication that our algebraic approach, numerical procedure and implementation of the latter are all sound, and that additional results are to be trusted.

Perturbation ansätze
Given an exactly critical analytical solution h, one can perturb it to discover the critical exponent γ letting where ε is a small number. By expanding the equations of motion in powers of ε, variable. This suggests the Fourier-Laplace decomposition where κ spans with generality over a discrete set of modes in C. In κ-space, the equation one thus obtains is algebraic. We can therefore from now on omit the with the understanding that once one identifies the spectrum of κ solving the corresponding equations for h 1 (z), the general solution to the first-order equations is given by a linear combination of these modes.
The spectrum of κ will be generally not real. Modes with Re κ < 0 represent instabilities, and their existence would imply that the solution is not truly critical.
Unstable CSS collapses were identified for example in [19] in four-dimensional σ-models. Given that the axion-dilaton system in four dimensions in the elliptic class is known to be stable, we have assumed that the same property holds for all cases studied in this work, and we have only confined our attention to Re κ > 0.
We are actually interested specifically in the mode κ * with largest real part, or in other words the fastest decaying mode, which is related to the Choptuik exponent via [5,14,16] We have assumed that, as for the four-dimensional elliptic critical collapse, κ * is real. This restriction is required due to the demanding computational workload, although some of the self-similar backgrounds discussed here might have complex fastest decaying modes. If so, our determination of the Choptuik exponent would be an overestimate.

Pure gauge modes
It is important to identify which values of κ correspond to unphysical gauge modes. These will be solutions of the perturbed equations that however are not physically distinguishable from the background. One can generate these modes, and thus identify them, performing infinitesimal gauge transformations on the self-similar solutions. We remark that the following discussion is valid in any dimension.
Consider first the symmetry transformations f (z) → e i f (z) for the elliptic class, and f (z) → e f (z) for the hyperbolic class. If 1, they reduce to while b 0 (z) is unaffected. Comparing with the perturbation ansätze (26), one can deduce that the following mode has been generated: One can thus see that this unphysical mode corresponds to κ = 0, and consequently this root must be excluded from the analysis.
In fact, a similar mode would also be present in the parabolic class, if one were to identify self-similar solutions there, but as we have anticipated we have found no solutions belonging to this class. At any rate, under f (z) → f (z) + one would generate the mode b 1 (z) = 0, f 1 (z) = 1, which has also κ = 0. Even more generally, in any model with a matter sector possessing internal symmetries, these will always give rise to pure gauge perturbations that are scale-invariant, and thus again with κ = 0.
Time translations at constant r, t → t + , are another class of gauge transfor-mations. In this case r → r, and therefore z → z − z t , so that Comparing with the ansatz, one can read off again the value κ = 1. Therefore, this represents an additional unphysical gauge mode.

Linearized equations of motion in any dimension
We shall now set up the perturbation theory for the self-similar solutions discussed above. We first perform this step for the elliptic class of solutions in an arbitrary dimension d = q + 2 ≥ 4

Linearized equations for the elliptic class
The perturbation ansatz (26), specified to the u, b, τ functions, reads Using the ansätze (31), (32) for the metric functions, one can compute the Ricci tensor for the metric as a function of ε. The dependence on the spacetime dimension d = q + 2 can be determined following the procedure described in [20]. One can thus find the zeroth-order and first-order parts from the limiting behaviors The same can be done with the right-hand sideT ab of the field equations, introducing the axion-dilaton perturbations (33), (34), and therefore, in a similar fashion,T The Einstein Field Equations (EFEs) must hold order by order, so that It is now possible to use part of these equations to eliminate u(t, r) from the others. Specifically, one can use R We omit the explicit form of u 1 (z), since it is rather cumbersome, and henceforth we shall also assume that u(t, r) and its first derivatives have been expressed, in all equations, in terms of the other variables. One is thus left with b(t, r) as the only remaining metric function, and it is now interesting to combine the field equations in such a way as to eliminate the second-derivative terms in b(t, r). This procedure identifies the so-called Hamiltonian constraint, which here results from the combination The lowest-order contribution is a first-order equation linking b 0 to b 0 f 0 , f 0 (in fact, one thus recovers (16)), Note that the right-hand side of (45) is real, as expected. In a similar fashion the first correction, where At this point we can make an important remark: eq. (47) is evidently invariant under dilations (t, r) → (e λ r, e λ t), so that, changing coordinates to (t, z), the factors of t will cancel and the result will only depends on z. All in all, one is actually dealing with a real and linear ordinary differential equation.
One can now introduce the perturbation ansätze in the τ equation of motion (8).
After extracting the zeroth and first-order parts, one can finally replace all instances of u 0 , u 0 , u 1 , u 1 as above. The resulting zeroth-order part is an  17): (49) Here and again the first-order part will involve b 1 , b 1 , f 1 , f 1 , f 1 linearly, and the zeroth-order functions and their derivatives non-linearly. In detail where Note, again, that the equation is explicitly scale-invariant, and therefore it turns into an ordinary differential equation after a change of coordinates to (z, t).
Once the zeroth-order solutions are obtained integrating numerically the unperturbed equations, and after replacing b 1 via eq. (47), one is left with a system of ordinary linear differential equations that we write schematically in the form (54) B 1 and F 1 are linear functions that still depend non-linearly on the unperturbed solution. In addition, there is still a quadratic dependence on κ. Note that these equations inherit the singularities of the unperturbed system of eqs. (16) and (17). In particular, they are also singular for z = 0 and b 2 (z) = z 2 .

Linearized equations for the Hyperbolic class
The same procedure can be followed for the hyperbolic class, again in an arbitrary number of dimensions. In this case eq. (33) has to be replaced by while all the other assumptions are unchanged. The form of u 0 and u 1 is the same as in the elliptic case, since these quantities are independent of the matter content. One thus obtains for the hyperbolic class of solutions the following form for u 0 : As before, for brevity we omit the explicit expression of u 1 (z).
Following similar steps as before, the unperturbed equations of motion in the hyperbolic class are found to be Retracing the steps of the previous section yields the linearized first-order equa- where Finally, the linearized equation for f 1 reads where L 6 = −r 2 b 0 + t 2 b 3 0 . Once more, the equations for b 1 and f 1 are scaleinvariant, and thus can be turned into ordinary differential equations in z.
In both cases, the modes are identified determining the values of κ that lead to smooth solutions of the perturbed equations (53), (54) that satisfy certain boundary conditions, which we now turn to describe.

Boundary conditions for perturbations
Eqs. (53), (54) are to be supplemented by suitable boundary conditions. At z = 0, thanks to the freedom to rescale the time coordinate, one can set In addition, regularity of the axion-dilaton implies So that one is left only with an unknown complex parameter, f 1 (0).
At z + one can also demand that the axion-dilaton (and, thus, its perturbation) be regular across the homothetic horizon, which we emphasize is only a coordinate singularity. We argue that regularity is encoded in the finiteness of the second derivatives ∂ 2 r f (t, r), ∂ r ∂ t f (t, r), ∂ 2 t f (t, r) as z → z + , since t and r are themselves regular there. This implies that f 0 (z) and f 1 (z) must be finite as z → z + . In practice, we define a parameter β = b 0 (z) − z and then expand f 0 and f 1 near the singularity, as where we have introduced the schematic notation h 0 = (b 0 (z + ), f 0 (z + ), f 0 (z + )), ). The vanishing of G(h 0 ) is simply the unperturbed complex value constraint (20) at z + , and we shall assume that background solutions satisfy this condition. We also verify that which is a consistency check. One is thus only left with the complex-valued which is linear in h 1 . We do not display the specific forms of this constraint for the elliptic and hyperbolic classes, since they are particularly unwieldy, but we remark that they are in principle easily obtained from the zeroth-order and first-order equations of motion according to eqs. (63) and (64). We solve this constraint for f 1 (z + ) in terms of f 1 (z + ), b 1 (z + ), κ and h 0 . This condition has thus reduced the number of free parameters in the boundary conditions at z + to a real number b 1 (z + ) and a complex f 1 (z + ). Summarizing, there are six real unknowns, which are κ and the five-component vector: and the system of linear ODE's eqs. (53), (54), whose total real order is five.

Numerical procedure
We define a "difference function" D(κ; X) (X is as defined in (67)) through the following procedure: • First solve the linear equation (66) to determine f 1 (z + ).
• With complete boundary conditions at z = 0 and z = z + , integrate forward from z = 0 to an intermediate point z mid , and also backward from z + to z mid .
• Collect in a vector (b 1 , Re f 1 , Im f 1 , Re f 1 , Im f 1 ) computed at z mid the difference between the two results.
We note that D(κ; X) is linear in X. Indeed, since (66) is linear in X, so are the boundary conditions for the ODE system, and since the solutions to a linear ODE system are linear in the boundary conditions themselves, so will be the differences of the values at z mid as functions of X. Since D is linear, it admits a representation as a matrix product: with A(κ) a 5 × 5 real matrix depending non-linearly on κ.
The search for solutions is thus reduced to the task of finding zeroes of D(κ; X). We compute the columns of A(κ) numerically through evaluations on the basis vectors of R 5 : which amounts to five evaluations of D(κ; X). We then perform an easy root search for the determinant as a function of κ. The root with the largest value is related to the Choptuik exponent through eq. (27).

Momentum-space singularity
One could notice that the perturbed equations of motion (47), (54) are singular whenever the factor present in the denominators, vanishes. For fixed κ, a zero can be encountered in the integration, and the numerical procedure fails at that point. We actually find that there are entire regions for κ, in the different solutions, which lead to failure. We can roughly pinpoint the location of these intervals determining the value of κ for which a singularity appears directly at z + : In these case the divergence affects directly the boundary condition (66).
We find evidence that there is generically a whole range of values for κ near κ sing , within which the singularity is reached somewhere between z = 0 and z = z + .
κ sing appears to usually lie either within or at the right edge of this interval. The numerical procedure more specifically either fails as the integration cannot be brought to completion or because numerical instabilities are not properly taken into account, and the output of the difference function is then incorrect.
One could possibly bypass this issue performing in advance suitable algebraic reparametrizations. However, the most relevant mode κ * lies, in almost all cases, outside the failure region (either before or after), so that this problem does not affect our estimates of the Choptuik exponents (see Table 1).

Results
We can now present the numerical results for the Choptuik exponents of the critical solutions, in four and five dimensions, in the elliptic and hyperbolic classes introduced in Section 2.2.
We first test the technique on the unique four-dimensional elliptic solution (4E).
In to be which implies a Choptuik exponent although, strictly speaking, our determination is guaranteed to be precise only up and including the first two decimal digits 4 .
Remarkably, this result is in perfect agreement with the existing literature, see [30]. We take this non-trivial match as evidence that our methodology and our numerical set up are reliable. For this solution, κ sing = 1.224, and the integration fails around 1 κ 1.4, well below the position of the most relevant mode. This is visible as a range of κ values in Fig. 1 for which Mathematica could not complete the computation of det A or returned a grossly incorrect output. In any case, because of the considerable distance of this region from the final crossing, we have reasons to consider our numerical determination of κ * reliable.
Our results for the other Choptuik exponents can be found in Table 1, together with some additional details on the failure regions.
Note that our results are in contrast with those of [21,22]. In fact, as there exists a mismatch already at the level of the self-similar solutions between [20] and [21,22], there is no expectation of overlap of the results of perturbation

Conclusions
In this paper, we have outlined a numerical set up for perturbing self-similar collapse solutions to the Einstein-axion-dilaton system in arbitrary spacetime dimensions. In all cases we focus on the most relevant mode, which determines the Choptuik critical exponent. Not only does the method work for different conjugacy classes of the SL(2, R) transformation, but it can potentially encompass more general choices of matter content. In terms of results, what is perhaps most striking is that not only do the Choptuik exponents depend on the matter content (the axion-dilaton sector, in this case) and on the choice of spacetime dimension, but they also differ sizably for the various solutions of self-critical collapse that were recently found in [20].
Therefore, our results provide some clear evidence that the original expectation that Choptuik phenomena exhibit some form of universality [14] is not fulfilled.
The current analysis does not exclude that some universal behavior lies hidden in combinations of the critical exponents and other parameters of the theory.
Nevertheless, this analysis and previous ones provide some definite evidence against transferring standard expectations of Statistical Mechanics to the case of critical gravitational collapse. Some uncertainties, however, still remain. In particular, the last value of γ in the table 1 is sizably larger than the others and emerges from a case where the value of k * lies below the κ = 1 gauge mode. We hope to return to these important issues in the future.