Critical collapse of an axisymmetric ultrarelativistic fluid in $2+1$ dimensions

In the spherically symmetric case, we showed in Ref. [6] that the critical phenomena depend on the value of $\kappa$. In particular, for $\kappa\lesssim 0.42$, we find typical type I behavior, where the mass and curvature do not scale and the critical solution near the center is static. How is this picture modified in the presence of angular momentum?

We start answering this question by showing, in Fig. 1, a log-log plot of the apparent horizon mass $M_{\text{AH}}$ (green group of curves), the adimensionalized apparent horizon angular momentum $J_{\text{AH}}/\ell$ (blue group of curves), and the maximum of curvature $-\Lambda\rho_{\text{max}}^{-1}$ (orange group of curves), against $p-p_{\star}$, for different values of $p_{\omega}$. For all these plots, we consider centered initial data with $\kappa=0.4$. The range of $p_{\omega}$ is $p_{\omega}=0.01$, $0.1$, $0.2$, $0.4$, $0.6$, $0.8$ and $1.0$. This corresponds to a range from “small” angular momentum (in the sense that $J_{\text{AH}}/(M_{\text{AH}}\ell)\ll 1$) to “large” angular momentum (for example, $p_{\omega}=1.0$ corresponds to $J_{\text{AH}}/(M_{\text{AH}}\ell)\simeq 0.64$).

This behavior, as in spherical symmetry, persists up to some critical value of $\kappa$ between $\kappa=0.42$ and $\kappa=0.43$. For comparison, we have also added an evolution with $\kappa=0.43$ (dotted curves) and $p_{\omega}=0.01$. The type I behavior also holds for the off-centered and ingoing initial data.

As in the spherically symmetric case, we consider the time-dependent quantity $R_{M}$, defined by

as a measure of the length scale of the solution. Similarly, we define the central density $\rho_{0}$, mass $M_{\text{OB}}$, and angular momentum $J_{\text{OB}}$ at the numerical outer boundary,

In Fig. 2, we plot in a linear-log plot, $\sqrt{-\Lambda\rho_{0}^{-1}}$, $J_{\text{OB}}/\ell$, $\sqrt{M_{\text{OB}}}$, and $R_{M}/\ell$ against central proper time $t$, for sub6 to sub15 initial data. We consider here centered initial data with $p_{\omega}=1.0$. Note that we have shifted $\log_{10}J_{\text{OB}}/\ell$ by a constant $c=1.2$ for clarity. As is typical for critical phenomena, the more fine-tuned the initial data to the black-hole threshold, the longer the critical regime, before the growing mode of the critical solution becomes dominant. In order to avoid cluttering, we truncate the plots after it is clear that the growing mode causes the curves to “peel off” from the critical regime.

For type I phenomena, it is not the mass and curvature that scale. Instead, it is the lifetime of the intermediate regime where the solution is approximated by the critical solution. Assuming that the critical solution is stationary, we can make the following ansatz for its linear perturbations:

where $Z$ stands for any dimensionless metric or matter variable, such as $R^{2}\rho$, $M$, or $\alpha$.

We assume the critical solution has a single growing mode, ${\rm Re}\,\sigma_{0}>0$. Since the solution is exactly critical at $p=p_{\star}$, this implies that $C_{0}(p)\sim p-p_{\star}$. We define the time $t=t_{p}$ to be the time where the growing perturbation becomes nonlinear. This occurs when

and so

The exponent $\sigma_{0}$ can be read off from Fig. 2. We find $\sigma_{0}\simeq 8.29$, close to its value in spherical symmetry (which was $\simeq 8.84$) [6].

In [11], we showed the existence of a two-parameter family of rigidly rotating, stationary star solutions for any causal equation of state $P=P(\rho)$. These solutions are analytic everywhere including at the center, and have finite total mass $M$ and angular momentum $J$. The two free parameters can be taken to be $s$, giving the overall length scale of the star, and $\tilde{\Omega}$, parametrizing the rotation of the star. The latter is defined so that $\tilde{\Omega}=0$ corresponds to nonrotating stars. We can write these exact solutions as

where $\check{Z}$ is the corresponding exact stationary solutions. The cosmological constant enters this picture through a specific combination, parametrized by the dimensionless quantity $\mu$,

In what follows, all quantities pertaining to the exact stationary solutions have a check symbol, as in (37).

The parameters $\mu$ and $\tilde{\Omega}$ correspond to $\mu$ and $\Omega$ in [11]. (The tilde is used to distinguish $\tilde{\Omega}$ from the unrelated conserved variable $\Omega$). Regular stationary solutions exist only for $0\leq\mu\leq 1$, and $\mu$ can be interpreted as parametrizing the competition of the attractive acceleration induced by a negative cosmological constant and the centrifugal acceleration.

In Fig. 3, we compare our best subcritical numerical solution at $t=2.1$, where the solution is approximately stationary, to the family of stationary solutions. Since the stationary solutions form a two-parameter family, we need to fit those two parameters. For this, we match the central density and total angular momentum, using the angular momentum at the numerical outer boundary as a proxy for the total angular momentum. The matching conditions are

Note that there is a slight abuse of language here; the total angular momentum (and mass) of the system is a conserved quantity. The angular momentum (and mass) at the numerical outer boundary (at finite radius) differ from it, since a bit of spin (and mass) are radiated away through this boundary at the start of the evolution, when the solution has yet to enter the critical regime. At late times, the approximately time-independent value of the angular momentum (and mass) at the outer boundary during the critical regime is a fixed fraction of the constant total spin (and mass) of the system.

Our matching of the angular momentum therefore introduces a small systematic error. At first, we attempted to ask for the value of $R_{M}$ between the numerical and exact solutions to agree. This is also possible, but $R_{M}$ depends very weakly on rotation. Such a matching condition is therefore very sensitive to numerical error.

In the example of Fig. 3, our matching procedure gives $s\simeq 0.0882$ and $\tilde{\Omega}\simeq 0.0498$. We then find very good agreement for $\rho$, $M$ and $a/\alpha$ at all $R$. We find slightly less good agreement for the angular momentum. Since $J$ is relatively small, it is likely that this is a numerical error, although its precise nature is unknown.

We remark here that, at first, we attempted to push the numerical outer boundary further out. We noticed however, for $R_{\text{max}}\gtrsim 4$, there is a clear numerical error that builds up and causes mass and spin to slowly radiate away. For example, for $R_{\text{max}}\simeq 1.25$, the total mass and angular momentum stay constant in the critical regime within $\ll 1\%$ (see Fig. 7). On the other hand, for $R_{\text{max}}\simeq 4$, they decrease slowly and approximately linearly. As compared to the near constant value that they take for $R_{\text{max}}\simeq 1.25$, they further decrease by about $29\%$ and $16\%$ for $R_{\text{max}}\simeq 4$, respectively. This is measured up until the time when the “star” disperses. This decay rate is also different for larger $R_{\text{max}}$. We do not have an explanation for this numerical error.

In spherical symmetry, we have given strong evidence that the critical phenomena are type II for $\kappa\geq 0.43$ (the apparent-horizon mass and curvature scale as some power law and the critical solution shrinks quasistatically to arbitrarily small size); see [6]. As we did for $\kappa\leq 0.42$, we now consider families of initial data with different initial angular momentum $p_{\omega}$. In the following, we focus our attention on the case $\kappa=0.5$.

In Fig. 4, we plot $-\Lambda\rho_{\text{max}}^{-1}$ (orange group of curves), $M_{AH}$ (green group of curves) and $J_{AH}/\ell$ (blue group of curves) as a function of $p-p_{\star}$ for centered initial data and with different $p_{w}$, namely, $p_{w}=0.02$, $0.04$, $0.06$, $0.08$, $0.1$, $0.16$, $0.2$, $0.3$, $0.4$ and $0.5$ (solid lines). The dashed-dotted lines were obtained with $p_{w}=0.2$, but using the monotonized central (MC) limiter.

This plot illustrates multiple points. First, for “small” initial spin, we find typical type II phenomena where the density, mass and spin scale like power laws. Second, the angular momentum $J_{\text{AH}}$ scales more slowly than $M_{\text{AH}}$. In particular, the spin-to-mass ratio of the resulting black hole increases as we fine-tune to the black-hole threshold; see also Fig. 6. For larger initial spin such as $p_{w}\geq 0.2$ (or respectively more fine-tuning for smaller $p_{\omega}$), the spin-to-mass ratio approaches extremality. It does not go beyond, however, because the critical solution becomes stationary and, as a result, the power-law scaling with respect to $p-p_{*}$ also smoothly levels off.

The phenomena highlighted above have also been checked to hold for the off-centered and initially ingoing data, although we chose not to include them in the plot to avoid cluttering.

Note that for the evolution using the MC limiter (dashed-dotted), the scaling stops abruptly and the spin-to-mass ratio seems to remain constant with further fine-tuning. For reasons that will be made clearer later, this is expected to be a numerical artifact: for relatively large spin-to-mass ratio, numerical instabilities are generated at the numerical outer boundary, similar to those we faced when simulating type I phenomena in the spherically symmetric case. These instabilities travel inward and produce shocks near the center, which are absent using a more dissipative limiter such as the first-order Godunov limiter.

In Fig. 5, we explore the relationship between $J_{\text{AH}}$ and $M_{\text{AH}}$ more explicitly, for the critical solution and the final black hole. On the left, we show a parametric plot of $M_{\text{AH}}$ against $J_{\text{AH}}/\ell$ for the same initial data as in Fig. 4 (and using the same convention for the lines). We have added the bisections from the off-centered and ingoing initial data with $p_{w}=0.02$, $0.04$ and $0.06$ (dotted and dashed lines, respectively). The red line corresponds to $J_{\text{AH}}=M_{\text{AH}}\ell$. The trajectories in the $(M_{\text{AH}},J_{\text{AH}})$ plane are clearly family dependent, although as we increase the fine-tuning and the black holes become smaller (the bottom left corner) all the curves become approximately parallel to each other. This provides some evidence that the critical phenomena are controlled by a unique one-parameter family of critical solutions (universality).

On the right, we show a parametric plot of the mass $M_{\text{OB}}$ against the angular momentum $J_{\text{OB}}/\ell$, for our best supercritical data. As expected, both plots are qualitatively very similar, since one expects the mass and spin at the apparent horizon to be in some fixed ratio to their values at the numerical outer boundary, at the time when the solution veers off from the critical solution. It turns out that, for close to critical data, $M_{\text{AH}}\simeq M_{\text{OB}}$ and $J_{\text{AH}}\simeq J_{\text{OB}}$ at this time, and so this fixed fraction is almost one.

In Fig. 6, we use again the same data as in Fig. 5. On the left plot, we show the spin-to-mass ratio $J_{\text{AH}}/(M_{\text{AH}}\ell)$ as a function of $p-p_{\star}$ for different levels of fine-tuning. The only addition is a second bisection using a MC limiter (a second dashed-dotted curve), but where the numerical outer boundary is at $R_{\text{max}}\simeq 2.65\simeq 4.17\ell$. The fact that both dashed-dotted lines do not level off at the same spin-to-mass ratio gives evidence that the unphysical numerical outer boundary spoils the numerical results in that highly rotating regime. Instead, the results from the Godunov limiter, which removes the aforementioned instabilities, are more plausible; in the limit of fine-tuning, the black hole is extremal.

Similarly, in the right plot, we show the spin-to-mass ratio evaluated at the numerical outer boundary, for our best supercritical data, against proper time $t$. Both left and right plots are qualitatively similar. This is expected since the more fine-tuned the initial data are to the black-hole threshold, the longer it takes for the growing mode to dominate the perturbation, and so the more the spin-to-mass ratio can grow, before the solution collapses.

If we define the time where the critical regime starts to be the time from which the central density shows critical scaling (at least for small rotation), then this occurs at $t\simeq 0.8$ for centered initial data (see, for example, Fig. 7), while it occurs at $t\simeq 1.0$ for ingoing and off-centered initial data. In the right plot of Fig. 6, the different trajectories do not cross during this critical regime. This again provides evidence that, in a suitable adiabatic limit, a given pair $(M_{\text{AH}},J_{\text{AH}})$ in a quasistationary sequence uniquely determines its evolution.

We have seen in Fig. 5 that for sufficient fine-tuning (on the left) or at sufficiently late times (on the right), the trajectories in the $J$-$M$ plane do not cross, providing evidence that there is a universal one-parameter family of critical solutions that fibrates the $J$-$M$ plane. Independently, the fact that we can make arbitrarily small black holes or arbitrarily large central densities by fine-tuning generic one-parameter families of initial data suggests that each member of this family of critical solutions has precisely one unstable perturbation mode.

As in the spherically symmetric case, the profile of the solution can be roughly split into two regions: one central region where most of the density lies and whose size shrinks in time and an outer region (“atmosphere”) where the mass and spin are approximately constant in space and decrease to zero (exponentially) in time.

At a stage where the effect of spin can still be regarded as perturbative, one would expect the critical solution to be approximated by the quasistatic solution in spherical symmetry, plus a (unique, growing) nonspherical perturbation proportional to $p_{w}$ that carries the angular momentum. In particular, one would expect that the solution shrinks exponentially quickly in time.

To confirm this, we plot, on the left of Fig. 7, the logarithms of $R_{M}(t)/\ell$, $\sqrt{-\Lambda\rho_{0}^{-1}(t)}$, $\sqrt{M_{\text{OB}}(t)}$ and $J_{\text{OB}}/\ell$ for sub6, sub9, sub12, and sub15 centered data with $p_{w}=0.01$. For those plots, we used the MC limiter, as even at sub15, the spin-to-mass ratio still remains relatively small, and the aforementioned instabilities are not present.

We find, as anticipated from our study in spherical symmetry, that these quantities are exponential functions of $t$. $J_{\text{OB}}$ scales slower then $M_{\text{OB}}$ as it is the case for the apparent-horizon mass and spin, so that the spin-to-mass ratio $J_{\text{OB}}/(M_{\text{OB}}\ell)$ increases. Note that right after the critical regime, the trajectories for each variable at different levels of fine-tuning align, up to a rescaling and a shift in time. This kind of behavior is expected in a spacetime where the cosmological constant is dynamically irrelevant and where the field equations become approximately scale invariant.

Note that after the critical regime, $J_{\text{OB}}$ enters what seems to be a second regime, still decreasing exponentially, but noticeably more slowly than during the critical regime.

On the right of Fig. 7, we plot the logarithm of $R_{M}(t)/\ell$ and $\sqrt{-\Lambda\rho_{0}^{-1}(t)}$ for our best subcritical evolution (sub15) with $p_{w}=0.01$, $0.05$, $0.1$, $0.16$, $0.2$ using again centered initial data, but with the Godunov limiter. The solution first shrinks and displays typical type II phenomena for “small” $p_{w}$. During this phase, $R_{M}(t)/\ell$ and $\sqrt{-\Lambda\rho_{0}^{-1}(t)}$ are completely independent of $p_{\omega}$, as we would expect while angular momentum is a perturbation.

As the spin-to-mass ratio has become sufficiently large, both $R_{M}$ and $\rho_{0}^{-1}$ decrease more slowly and start to level off. For the same fine-tuning, this corresponds to larger values of $p_{w}$; see the inset on the right of Fig. 7. We do not have a satisfactory dynamical explanation for this.

Together with Fig. 4, one expects that in the limit of perfect fine-tuning, the critical solution becomes stationary at late times, and correspondingly, on the black-hole side of the threshold of collapse, the black hole becomes extremal.

In spherical symmetry, the critical solution is quasistatic, meaning that it adiabatically goes through the sequence of static solutions, with $s$ now a function of $t$, and $|\dot{s}(t)|\ll 1$. Furthermore, $\mu\sim s^{2}$, and it follows that in the quasistatic ansatz, $\mu(t)\sim s(t)^{2}$. Since the solution shrinks exponentially, we have the following ansatz for $s(t)$:

for constants $s_{0}$ and $\nu$. To leading order in $\dot{s}$, the quasistatic ansatz then takes the form

and similarly for other suitably rescaled variables.

As we hinted before, one can expect, at least in the regime where the effects of angular momentum are still perturbative, that the solution can be thought to shrink adiabatically to zero size, going through a sequence of stationary solutions. The picture in spherical symmetry can then be straightforwardly generalized to axisymmetric initial data. Specifically, we assume that, to leading order in $\dot{s}$, the critical solution can be well approximated by

and so on. From the exponential form of $s(t)$, we further make a relatively agnostic ansatz for $\tilde{\Omega}$ of the form

where $\tilde{\Omega}_{0}$ and $\phi$ are constants. This exponential ansatz is justified from the fact that, for the family of stationary solutions, the angular momentum at infinity satisfies (39). In particular $J\sim\tilde{\Omega}$. From Fig. 7, $J_{\text{OB}}$, seen as a proxy for the corresponding angular momentum at infinity, decays exponentially, thus suggesting the exponential form for $\tilde{\Omega}$. Together with the ansatz for $s(t)$ (41), we have now four parameters to fit: $s_{0}$ and $\nu$ (as in spherical symmetry) and $\tilde{\Omega}_{0}$ and $\phi$.

For $p_{w}=0.01$, where we only observe type II behavior to our level of fine-tuning, we choose to fix them by requiring the central density and spin at infinity to match those of the stationary solutions at times $t\simeq 1.0$ and $1.4$, where we believe that the solution is in its critical regime; see Fig. 7. For the spin at infinity, we take the spin at the numerical outer boundary as a proxy. This is justified because in the atmosphere, where the density is small, the mass and spin are approximately constant in space. We find the following values:

As expected, the above values of $s_{0}$ and $\nu$ are consistent with their values in the spherically symmetric case $p_{w}=0$. We have repeated this procedure for the ingoing and off-centered initial data as well, and summarize the values of $\nu$ and $\phi$ in Table 1. Due to the smallness of $\phi$, the relative variation of $\phi$ is important. In particular, it is difficult to confidently say if those variations are purely a numerical error. However, because Figs. 5 and 6 suggest some universality for the critical solution, we are inclined to believe it is so. That is, we believe that $\phi$, just as $\nu$, is family independent. We find good agreement between the numerical and exact solutions.

In Fig. 8, we compare the leading-order quasistationary solution for $J$, of the same form as in Eq. (43) (dotted lines), with the numerical data (solid lines) at constant time intervals $t\simeq 1$, $1.1$, $1.2$, $1.3$ and $1.4$.

Returning back to the left plot of Fig. 7, the quasistationarity of the critical solution can be used to explain the slope of $J_{\text{OB}}$. Specifically, from (39), we have

In the spherically symmetric case, we showed that $M_{\text{OB}}\sim s^{2}$. At the level where angular momentum is still viewed as a perturbation, this is still expected to hold, so that

or

where

is small. For the centered initial data for example, (45) gives $\Delta\simeq 0.144$.