# Critical collapse of ultra-relativistic fluids: damping or growth of aspherical deformations

###### Abstract

We perform fully nonlinear numerical simulations to study aspherical deformations of the critical self-similar solution in the gravitational collapse of ultra-relativistic fluids. Adopting a perturbative calculation, Gundlach predicted that these perturbations behave like damped or growing oscillations, with the frequency and damping (or growth) rates depending on the equation of state. We consider a number of different equations of state and degrees of asphericity and find very good agreement with the findings of Gundlach for polar modes. For sufficiently soft equations of state, the modes are damped, meaning that, in the limit of perfect fine-tuning, the spherically symmetric critical solution is recovered. We find that the degree of asphericity has at most a small effect on the frequency and damping parameter, or on the critical exponents in the power-law scalings. Our findings also confirm, for the first time, Gundlach’s prediction that the modes become unstable for sufficiently stiff equations of state. In this regime the spherically symmetric self-similar solution can no longer be recovered by fine-tuning to the black-hole threshold, and one can no longer expect power-law scaling to hold to arbitrarily small scales.

###### pacs:

04.20.Jb, 04.70.Bw, 98.80.Jk, 04.25.dg## I Introduction

Critical phenomena in gravitational collapse were first reported in the seminal work of Choptuik Choptuik (1993). Consider initial data that are parametrized by some parameter , and assume that the data will either collapse to form a black hole – for sufficiently large values of , say – or will disperse and leave behind flat space – for sufficiently small values of . Then there must exist a critical value of the parameter, , which marks the onset of black-hole formation and separates supercritical from subcritical data.

Critical phenomena refer to properties of the solution in the vicinity of . Specifically, Choptuik found that, for close to , the evolution will pass through a phase of self-similar contraction. Ultimately the evolution will leave this self-similar phase to either form a black hole, or to disperse to infinity. The closer is chosen to , the longer the evolution will follow the self-similar contraction, and hence the smaller the length scale at which it will diverge from the self-similar solution. This length scale sets the length scale for any dimensional quantity characterizing the solution. For supercritical data, in particular, this results in the famous power-law scaling laws for the black hole mass

(1) |

where is the critical exponent. The value of the critical exponent can be found from the growth rate of perturbations of the self-similar critical solution (see, e.g., Maison (1996)).

In his original work, Choptuik Choptuik (1993) adopted a mass-less scalar field as the matter model, for which the critical solution displays a discrete self-similarity. Choptuik’s discovery triggered numerous follow-up studies, creating a large body of literature on an entire new field of research (see, e.g., Gundlach (2003); Gundlach and Martín-García (2007) for reviews). Particularly important for our purposes is the discovery of critical phenomena in the collapse of a radiation fluid by Evans and Coleman Evans and Coleman (1994). A radiation fluid is a special case of an ultra-relativistic fluid, i.e. a fluid whose equation of state is

(2) |

where is the pressure, the energy density, and is a dimensionless constant which, for a radiation fluid, takes the value . For ultra-relativistic fluids the critical solution is continuously self-similar – as opposed to the discrete self-similarity found for the massless scalar wave – which makes it easier in some ways to analyze this collapse. The critical exponent for a radiation fluid was determined analytically by Koike et al. (1995) as well as Maison (1996), who considered a number of different values of . In particular, these studies showed that the critical exponent is not universal, but depends on the matter model. Neilsen and Choptuik Neilsen and Choptuik (2000) generalized the studies of Evans and Coleman (1994) by performing numerical simulations of critical collapse of ultra-relativistic fluids with different values of , finding very good agreement in the critical exponents with Maison (1996).

Until recently, most numerical studies of critical collapse assumed spherical symmetry (but see Abrahams and Evans (1993); Choptuik et al. (2003, 2004); Sorkin (2011) for some notable exceptions). This is not surprising, since spherical symmetry makes it easiest to resolve the small structures and tiny black holes that form in critical collapse close to criticality. At the same time, some interesting questions in the context of critical collapse cannot even be addressed in spherical symmetry – relating, for example, to the effects of angular momentum Gundlach (1998, 2002a) or aspherical deformations Martín-García and Gundlach (1999); Gundlach (2002b) on the critical solution.

After numerical simulations of binary black holes in three spatial dimensions became possible (see Pretorius (2005); Campanelli et al. (2006); Baker et al. (2006)), most code development efforts in the numerical relativity community focused on the simulation of binaries. These codes have only rarely been used to study critical collapse (but see Alcubierre et al. (2000); Hilditch et al. (2013); Healy and Laguna (2014); Deppe et al. (2018)). A separate, more recent code development effort has resulted in methods for numerical relativity in curvilinear coordinates Montero and Cordero-Carrión (2012); Baumgarte et al. (2013); Montero et al. (2014); Baumgarte et al. (2015) (see also Ruchlin et al. (2018); Mewes et al. (2018) for more recent implementations of these techniques). Spherical coordinates, in particular, are well suited for simulations of critical collapse, and in fact, the code of Baumgarte et al. (2013); Montero et al. (2014); Baumgarte et al. (2015) has already been used to study critical phenomena in the gravitational collapse of both aspherical radiation fluids Baumgarte and Montero (2015) and rotating ultra-relativistic fluids Baumgarte and Gundlach (2016); Gundlach and Baumgarte (2016, 2018).

In this paper we expand the calculations of Baumgarte and Montero (2015), and perform numerical simulations of the critical collapse of ultrarelativistic fluids in the absence of spherical symmetry. Gundlach Gundlach (2002b) predicted from perturbative calculations that deviations from sphericity should behave as either damped or growing oscillations (see Eq. (6) below), where the oscillation frequency and the damping (or growth) rate depend on the value of in the equation of state (2). In Baumgarte and Montero (2015) this behavior was confirmed for a radiation fluid with , albeit only with a modest accuracy and only for two different values of the asphericity. Here we expand these calculations in several ways.

We first redo the calculations for Baumgarte and Montero (2015) with better grid resolution (see Section III below), which allows us to track the self-similar solution for longer and measure both the frequency and the damping (or growth) rate of the deformation more accurately. We also consider more and larger values of the deviation from spherical symmetry, so that we can examine its effect on the above parameters more systematically. Perhaps most importantly we also consider more general ultra-relativistic fluids, i.e. different values of in the equation of state (2), and find good agreement with the dependence of and on as predicted perturbatively by Gundlach (2002b) (see Fig. 8 below). In particular we confirm Gundlach’s result that for the modes become unstable and grow.

Our paper is organized as follows. In Section II we review some properties of the continuous self-similar solution encountered in the critical collapse of an ultra-relativistic fluid. In Section III we present some basic equations and our choice of initial data, we describe our numerical code, and discuss the diagnostics used to analyze our results. In Section IV we present these results, both for radiation fluids with and for more general ultra-relativistic fluids. We summarize and discuss our findings in Section V. Throughout this paper we adopt geometrized units with .

## Ii The critical solution

In Section III.2 below we will introduce a two-parameter family of initial data for ultra-relativistic fluids, with one parameter, , describing the overall density and a second parameter, , governing the deviation from spherical symmetry. We assume that for some part of this parameter space the evolution will pass through a self-similar phase. As discussed in more detail in Gundlach and Baumgarte (2018) (where the second parameter was the spin rate rather than ), the evolution can then be described as passing through three distinct phases (see also Fig. 3 below for an example).

In Phase 1, the evolution evolves from the initial data to the self-similar solution. During Phase 2, the evolution can be described as the critical, self-similar solution, plus a perturbation. The critical solution is unstable to at least one growing perturbative mode, which is spherically symmetric. Such a growing mode will cause the evolution to ultimately deviate from the critical solution, marking the transition to Phase 3. During Phase 3, the solution either disperses to infinity or collapses to a black hole.

For ultra-relativistic fluids with the equation of state (2) the (spherically symmetric) critical solution displays a continuous self-similarity Evans and Coleman (1994), describing a continuous contraction of the solution to an accumulation event (see, e.g., Fig. 1 in Neilsen and Choptuik (2000) for an illustration). We denote the proper time of the accumulation event, as measured by an observer at the center, as . The length scale of the critical solution at a proper time is then given by . A dimensionless quantity describing the critical solution can then depend on the ratio

(3) |

only, where is the areal radius (see Fig. 1 below for a demonstration). A dimensional quantity with units , on the other hand, where carries units of length, has to scale with . For the energy density at the center, for example, we expect

(4) |

The length scale of global properties of the solution, e.g. the black hole mass in supercritical evolutions or the maximum density in subcritical evolutions, is determined by the length scale of the critical solution at the transition from Phase 2 to Phase 3.

Spherically symmetric perturbations of the critical solution were considered by Maison (1996). There exists exactly one unstable, spherically symmetric mode, and the growth rate of this unstable mode, , determines the critical exponent in the mass scaling law (1) (see also Gundlach and Martín-García (2007); Gundlach and Baumgarte (2018) for a review of these arguments). In fact, as first realized by Garfinkle and Duncan (1998), the same arguments result in a power-law scaling for other quantities as well. For example, in subcritical evolutions the maximum energy density encountered during the evolution, , satisfies

(5) |

where, on dimensional grounds, .

Aspherical perturbations were considered by Gundlach Gundlach (2002b). In particular, Gundlach, considered modes of order , and showed that these modes display an oscillatory behavior described by

(6) |

where the time is given by^{1}^{1}1Strictly speaking, this should be , where we have yet to determine an overall length scale .

(7) |

and where the damping coefficient and the angular frequency depend on the constant in the ultra-relativistic equation of state (2) as well as the mode of the perturbation. In this paper we will focus on polar modes. Gundlach found that, for these modes, is negative for , resulting in a damping of these modes, but that becomes positive for , in which case the mode grows.

## Iii Basic equations and numerical method

### iii.1 Formalism

We solve Einstein’s equations, expressed in a reference-metric formulation Bonazzola et al. (2004); Shibata et al. (2004); Brown (2009); Gourgoulhon (2012) of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998), in spherical coordinates . In particular, this formalism adopts a “3+1” decomposition of the spacetime in which the line element is written as

(8) |

where is the spacetime metric, the lapse function, the shift vector, and the spatial metric (see, e.g., Alcubierre (2008); Baumgarte and Shapiro (2010); Gourgoulhon (2012) for textbook introductions.) We also adopt a conformal decomposition of the spatial metric,

(9) |

where is the conformal factor and the conformally related metric. In a reference-metric approach we also introduce the flat metric in whatever coordinate system is used; for our applications the reference metric is the flat metric expressed in spherical coordinates. Details of this formalism, and its implementation in our code, can be found in Baumgarte et al. (2013, 2015).

We assume that the matter is described by a perfect fluid with stress-energy tensor

(10) |

and that the equation of state is that of an ultra-relativistic fluid, i.e. Eq. (2).

### iii.2 Initial data

We adopt the same initial data as in Baumgarte and Montero (2015). Specifically, we choose a moment of time symmetry and assume that the metric is conformally flat initially, . We then choose the initial density distribution as

(11) |

where we have abbreviated

(12) |

and where is the second Legendre polynomial.

The initial data (11) form a two-parameter family, parametrized by , which governs the overall density amplitude, and , which determines the deviation from spherical symmetry. Note that in the above expressions is the (isotropic) coordinate radius. In the functions the product then becomes the areal radius in spherical symmetry. In this limit, the density distribution is centered on an areal radius of , and drops off on a length scale of . We set to unity in our code, , and hence report all dimensional quantities in units of . We introduced the functions in order to ensure that the density and its derivatives are well-behaved at the origin.

The conformal factor in (11) satisfies the Hamiltonian constraint

(13) |

where is the flat Laplace operator associated with the flat metric . We construct solutions to (11) and (13) using an iterative process.

In the limit of spherical symmetry () and for the density distribution (11) reduces to that adopted by Evans and Coleman (1994). In this case the total gravitational mass can be found analytically, , so that becomes a non-dimensional measure of the initial strength of the gravitational fields.

We complete the specification of the initial data by choosing and .

### iii.3 Numerical code

We evolve the gravitational field and fluid variables using the code described in Baumgarte et al. (2013, 2015). The code implements the BSSN equations Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998) in spherical coordinates, adopting a reference-metric formalism Bonazzola et al. (2004); Shibata et al. (2004); Brown (2009); Gourgoulhon (2012) together with a rescaling of all tensorial quantities in order to handle all coordinate singularities analytically and thereby allow for a stable evolution. We similarly adopt a reference-metric approach in solving the equations of relativistic hydrodynamics Montero et al. (2014); Baumgarte et al. (2015).

The code does not make any symmetry assumptions, but since we evolve axisymmetric initial data we set to zero all derivatives in the direction and use only one grid point in the azimuthal direction. For spherically symmetric evolutions (i.e. ) we adopt the minimum number of grid points in the direction, . In the absence of spherical symmetry () we adopt equatorial symmetry, and, unless noted otherwise, resolve the remaining hemisphere with a very modest number of grid points (see Fig. 5 below and the related text for a discussion). We use radial grid points in our simulations. As in Baumgarte and Montero (2015), these grid-points are allocated logarithmically (see Appendix A in Baumgarte and Montero (2015)), with the ratio in size between neighboring grid cells chosen to be . For 312 grid points, this means that the ratio between the size of the innermost and the outermost grid cells is .

Unlike in Baumgarte and Montero (2015), however, we also use re-gridding to improve the resolution during the evolution. This allows us to track the critical solution for longer, and to follow the oscillations (6) of aspherical modes for multiple periods. Specifically, we estimate the typical length scale of the solution at the origin from , and compare with the size of the innermost grid cell. Whenever exceeds a certain cut-off, chosen to be 0.05 in our simulations here, we shrink the entire grid by moving the outer boundary to a smaller value, and interpolating all data to the new grid. Unless noted otherwise we start with the outer boundary at (in units of ) and allow the grid to shrink down to in 10 steps. In some cases (see Section IV.2) we also performed higher-resolution runs and allowed the outer boundary to contract to in 20 steps. In either case we ended all simulations before the center comes into causal contact with the outer boundary, and becomes affected by numerical error originating there.

As in Baumgarte and Montero (2015) we carry out our simulations with “moving puncture” coordinate conditions. Specifically, we adopt the 1+log slicing condition

(14) |

(see Bona et al. (1995)) for the lapse , where is the mean curvature, and a Gamma-driver condition for the shift vector (see Alcubierre et al. (2003)) as presented in Thierfelder et al. (2011).

### iii.4 Diagnostics

We monitor several different quantities in order to analyze our numerical results.

For supercritical data we locate outer-most trapped surfaces. Once the newly formed black holes have settled down into an equilibrium state, the location of these surfaces coincides with that of an event horizon. We then determine the black hole mass from the proper area of the horizon. In the vicinity of the black hole threshold we can then make fits to (1) in order to determine the critical parameter as well as the critical exponent .

We also track the central energy density as a function of the proper time measured by an observer at the center. We fit to the expected behavior (4), both as an indirect verification of the self-similar contraction during Phase 2, and to find the proper times during which the solution passes through this phase.

For subcritical evolutions we also record , and fit these data to (5) to determine the critical exponent .

Measuring deviations from sphericity is more subtle (see also Choptuik et al. (2003) for a discussion). For simplicity we adopt the following approach to measure the degree of asphericity of our spatial slices. We start by computing, for every grid point , the proper distance from the origin along the coordinate line of constant and . In general, this definition of will depend on the choice of spatial coordinates. The polar () and equatorial () directions, however, take on an invariant meaning in our axially and equatorially symmetric simulations, so that measured in these directions is independent of spatial coordinates. We next define the dimensionless density variable

(15) |

We chose to use the proper distance in this definition, rather than the areal radius as in Evans and Coleman (1994), since it is easier to define the former in the absence of spherical symmetry.^{2}^{2}2We also considered the approach of Baumgarte and Montero (2015), who defined with . The radial variable is gauge-dependent, but does reduce to in the limit of spherical symmetry. As expected, differs from in our simulations, but results for the coefficients and as computed from or agree to within our estimated errors. Even though our definition (15) is slicing-dependent, displays self-similar behavior during Phase 2 in our simulations with the 1+log slicing condition (14), as demonstrated in Fig. 1.

At each instance of time we measure the maximum values of in the axial direction, , and the maximum value in the equatorial plane, (see Fig. 4 below for an illustration), and then compute the difference

(16) |

as a measure of the departure from spherical symmetry.

Given our assumption of equatorial symmetry, is affected by all even modes with . To linear order in the deformation parameter our initial density distribution (11) will produce an mode with an amplitude proportional to . Higher-order modes enter through non-linear coupling, and therefore become important only for large values of . However, we expect these higher-order modes to decay more rapidly than the modes (see Figs. 12 and 13 in Gundlach (2002b)), so that, at sufficiently late times, becomes a measure of the modes.

## Iv Results

### iv.1 Radiation Fluids

We first consider a radiation fluid with . We choose different values of the deformation parameter in the initial data (11), and then vary the amplitude parameter in order to bracket its critical value (see Table 1). We find that the difference between and its value in spherical symmetry, , is approximately quadratic in ,

(17) |

with . As discussed in Section III.4 we measure the black hole mass for supercritical data, and the maximum value of the central density for subcritical data. We then fit these two quantities to the scaling relations (1) and (5) in order to obtain the critical parameter as well as the critical exponents and . In addition to numerical error, these quantities also depend on which data points are included in the fits (see also Appendix A in Deppe et al. (2018) for a discussion). Data too far away from the critical parameter no longer satisfy the scaling relations (1) and (5), while data points too close to the critical parameter are more strongly affected by numerical error, because the increasingly small structures can no longer be resolved. In Fig. 2 we show an example for , where we do find excellent agreement of our data with the fits over multiple orders of magnitude. We also list the results of our fits, for all considered values of , in Table 1. The critical exponents and do not appear to depend on , and agree very well with the perturbative value of 0.3558 (see Maison (1996)).

perturbative | 0.3558 | -0.3846 | 3.6158 | |||
---|---|---|---|---|---|---|

0 | 0.12409 | 6.449 | 0.357 | 0.357 | – | – |

0.01 | 0.12409 | 6.449 | 0.355 | 0.356 | -0.36 | 3.64 |

0.1 | 0.12410 | 6.450 | 0.357 | 0.356 | -0.36 | 3.64 |

0.25 | 0.12417 | 6.451 | 0.356 | 0.357 | -0.36 | 3.64 |

0.5 | 0.12444 | 6.460 | 0.356 | 0.357 | -0.36 | 3.64 |

1.0 | 0.12554 | 6.496 | 0.356 | 0.357 | -0.37 | 3.65 |

We next fit the central density to the scaling (4) in order to obtain a value for , the proper time of the accumulation event (see Sect. II). Results for are again listed in Table 1. An example for is shown in Fig. 3, where we plot as a function of for pairs of subcritical and supercritical evolutions that bracket the critical value of with different accuracy. Note that time advances from right to left in this figure. We can clearly distinguish the three phases of the evolution that we described in Section II. Phase 1 starts with the initial data on the bottom right of the figure. During Phase 2, the evolution follows the self-similar critical solution; during this part of the evolution the central density is well approximated by the fit (4). In Fig. 3, Phase 2 starts around , when the central density approaches the fit marked by the solid line. Initial data that are better fine-tuned to the critical solution will follow the critical solution longer and to higher density; this is clearly visible in the figure. Ultimately, the evolution starts to deviate from the critical solution, with the central density either increasing more rapidly for supercritical evolutions, or dropping to smaller values for subcritical evolutions. The departure from the critical solution marks the transition from Phase 2 to Phase 3. Graphs like the one shown in Fig. 3 allow us to determine the time brackets during which our evolution follows the self-similar critical solution; this will be important in our analysis of aspherical deformations.

An example of the evolution of aspherical deformations is shown in Fig. 4 (see also ani () for an animation). Specifically, we show plots of the dimensionless density function , defined in (15), at six different instances of time during Phase 2. The wire-frame shows a subcritical, spherically symmetric evolution close to the critical solution (with ) as a function of and . We adjust the scale of the and -axes in Fig. 4 so that this function does not appear to change its shape at all, reflecting the self-similar contraction. Even though self-similarity is defined in terms of the gauge-invariant areal radius (see Eq. (3)), the coordinate radius used here apparently serves as an excellent proxy.

The colored surface in Fig. 4 represents an subcritical evolution with , with again within about of the critical value (values of for different values of are listed in Table 1). During Phase 2, this aspherical solution appears to perform damped oscillations around the spherically symmetric critical solution. We show snapshots of at six different subsequent times at which the deviation from the spherical solution are approximately largest; these snapshots show that the fluid appears to “slosh” back and forth between the poles and equator.

At any instance of time we measure the maximum values of in both the axial and equatorial directions, and . These values are marked by the green dots in Fig. 4. Computing then yields a measure of the aspherical deformation, as discussed in Section III.4. In Fig. 5 we plot as a function of (central) proper time for , for four different angular resolutions . The results for and can hardly be distinguished in this graph, except for very late times (see below). This demonstrates that for the numerical errors resulting from our low angular resolution are small, and highlights the advantages of spherical coordinates for these simulations.

We note that for different , the critical parameter takes slightly different numerical values. While we bracketed to for each resolution, the evolutions shown in Fig. 5 represent slightly different values of for each , meaning that they will depart from the critical solution at slightly different times. This explains the apparent non-convergent behavior at very late times.

Fig. 5 shows that performs a damped oscillation with decreasing period. The latter is not surprising, because the period of the oscillation should be related to the length scale of the unperturbed critical solution, which decreases with (see Section II). It is therefore more natural to display as a function of the time , see Eq. (7), where can be determined from the fit to (see Fig. 3). This graph, shown in Fig. 6, visually appears like the exponentially damped oscillation (6) predicted by Gundlach (2002b). We can now make fits^{3}^{3}3Strictly speaking, we include an offset in these fits, i.e. we fit to , but always find to be small. to (6) in order to find estimates for the damping coefficient and the angular frequency .

Since (6) describes deformations of the self-similar critical solution, fits to this behavior should be performed only during Phase 2, as identified, for example, in Fig. 3. In Section IV.2 below we will also see that, for stiffer equations of state, is still dominated by transients at early times during Phase 2, and that the behavior (6) dominates only at later times. It is therefore not surprising that the resulting fitted parameters and depend on the time window over which the fits were carried out. In Fig. 6, for example, the fit depends on whether the first peak, around , is included in the fit or not, which leads to changes in the damping parameter of several percent, and slightly smaller changes in the angular frequency .

In practice we perform these fits in two different ways. One fit assumes that (in ) is given from the fit to the central density (which, of course, also depends on the time window used in that fit), while the other fit simultaneously varies in the fit to (6). As a self-consistence test we than vary the time window until both fits result in parameters that are close to each other (typically less than 1% difference for a radiation fluid). The values reported in Table 1 are average values between the two fits, and carry an error of at least several percent, larger than those for the critical exponents and .

### iv.2 Other ultra-relativistic fluids

We next analyze the dependence of our results on the stiffness of the equation of state, i.e. on the constant in Eq. (2). Specifically, we consider , 0.4, 0.5 and 0.6 in addition to for the radiation fluid of Sect. IV.1. For each value of we choose different values of the deformation , and then bracket the critical parameter . We again perform fits to (1) and (5) to find the critical exponents and , and to (6) to find and . Numerical values for our fits are provided in Tables 2 through 5.

perturbative | 0.2614 | -1.296 | 5.1884 | |||
---|---|---|---|---|---|---|

0 | 0.10772 | 9.85 | 0.256 | 0.263 | – | – |

0.01 | 0.10772 | 9.85 | 0.261 | 0.263 | -1.2 | 5.2 |

0.1 | 0.10773 | 9.86 | 0.257 | 0.263 | -1.2 | 5.2 |

0.5 | 0.10806 | 9.89 | 0.262 | 0.265 | -1.2 | 5.3 |

As for the radiation fluid, we find that the critical exponents and agree well with the perturbative values of Maison (1996), and show very little dependence on , certainly within what we estimate to be our numerical and fitting errors. Also as for the radiation fluid, it is again significantly more challenging to determine the coefficients and for the deviation from the critical solution, but our values nevertheless agree quite well with the perturbative values provided by Gundlach (2002b).

perturbative | 0.4035 | -0.1715 | 3.07312 | |||
---|---|---|---|---|---|---|

0 | 0.12795 | 5.59 | 0.405 | 0.403 | – | – |

0.01 | 0.12795 | 5.58 | 0.403 | 0.403 | -0.16 | 3.10 |

0.1 | 0.12796 | 5.58 | 0.406 | 0.403 | -0.17 | 3.12 |

0.5 | 0.12832 | 5.59 | 0.404 | 0.403 | -0.17 | 3.11 |

1.0 | 0.12946 | 5.62 | 0.406 | 0.403 | -0.18 | 3.15 |

In Fig. 7 we show an example of these fits for all five different values of , all for . Several general trends are clearly visible in this Figure. We first notice that, for larger , the deformations are damped more slowly, i.e. increases, and the angular frequency of the deformations decreases. Moreoever, we find that changes sign around , so that the modes become unstable and grow for . All of these observations are consistent with the perturbative results of Gundlach (2002b), as demonstrated by the direct comparison our numerical values for and with those of Gundlach (2002b) in Fig. 8.

perturbative | 0.4774 | 0.0135 | 2.44 | |||
---|---|---|---|---|---|---|

0 | 0.13087 | 4.726 | 0.479 | 0.476 | – | – |

0.5 | 0.13127 | 4.733 | – | 0.475 | 0.003 | 2.48 |

The graphs in Fig. 7 also show that, for larger , the behavior (6) emerges only for later times , while earlier times appear to be dominated by transients (and possibly modes of higher order which, for larger , also decay more slowly; see Figs. 12 and 13 in Gundlach (2002b)). Since the period of the oscillations (6) also increases with increasing , and an accurate determination of and requires tracking the oscillation for multiple oscillation periods, these simulations become increasingly challenging for larger . For our usual grid setup with 10 regrids, shrinking the outer boundary to (see Section III), did not provide sufficiently accurate results. We therefore allowed up to 20 regrids down to an outer boundary of in these cases, but performed these simulations for only. Stopping the simulations before the center comes into causal contact with the outer boundary also meant that the simulations did not allow for enough time for the horizons to settle down – we therefore did not compute from these simulations. We note, however, that lower-resolution simulations with only 10 regrids showed very good agreement of with the perturbative values.

As shown in Fig. 7 the amplitude of the oscillations changes very little for (but appears to grow very slowly), while for the oscillations grow very clearly. This behavior is consistent with the perturbative findings of Gundlach (2002b), who found that changes sign at about . As we will discuss in the following section, this sign change has profound consequences for the behavior of solutions close to the black hole threshold.

perturbative | 0.5556 | 0.112 | 1.968 | |||
---|---|---|---|---|---|---|

0 | 0.13157 | 4.171 | 0.560 | 0.555 | – | – |

0.5 | 0.13201 | 4.176 | – | 0.554 | 0.102 | 2.03 |

## V Summary and Discussion

We perform numerical simulations of the gravitational collapse of ultra-relativistic fluids to study critical phenomena in the absence of spherical symmetry. Specifically, we consider initial data that, to lowest order, describe polar deformations of otherwise spherically symmetric fluid distributions. We evolve these fluids dynamically – using a numerical code that adopts spherical coordinates – and measure the deviations from spherical symmetry as a function of time. We vary the stiffness of the equation of state, parametrized by in (2), and consider different degrees of asphericity, parametrized by in (11). We find that the deviations are well described by oscillations that are exponentially damped or growing, see (6), in very good agreement with the perturbative results of Gundlach (2002b). In particular, we confirm that for stiffer equations of state, i.e. for larger , the growth rate in (6) increases, and that the oscillation frequency decreases (see Fig. 8). We also find that and depend at most very weakly on the degree of deformation .

For sufficiently soft equations of state, with , is negative, so that the oscillations are damped. We find (see Fig. 8); the value reported in Gundlach (2002b) is . In the limit of perfect fine-tuning, damped oscillations have an infinite amount of time to decay (as measured by the logarithmic time defined in (7)), so that the spherically symmetric critical solution is recovered, and the evolution is dominated by the unstable spherically symmetric mode (see Section II). For we therefore expect the characteristic power-law scaling to hold to arbitrarily small scales.

To the best of our knowledge, our simulations are also the first numerical confirmation that, for , the coefficient in (6) is positive, so that deformations grow rather than decay. This has profound implications for the nature of the black hole threshold, since, in the limit of perfect fine-tuning, the spherically symmetric critical solution can then be recovered only when exactly. For all non-zero , aspherical deformations will grow, and, for sufficient fine-tuning to the black hole threshold, will ultimately dominate the evolution. Therefore, power-law scaling can no longer be expected to hold to arbitrarily small scales. However, the aspherical mode grows only slowly in comparison with the spherical mode, so that, for a given value of , the effects of the aspherical mode on the power-law scalings can be observed for exquisite fine-tuning only. For and we observe some small deviations from power-laws for , but we are not confident that these are caused by true departures from power-laws rather than by numerical error.

A related effect was reported by Choptuik et al. (2003), who studied critical collapse of massless scalar fields in axisymmetry. They found that, for data with large aspherical deformations, and close to the black hole threshold, there exists a growing, aspherical mode. This growing mode ultimately leads to a bifurcation, so that two distinct regions on the symmetry axis collapse individually.

We also note that for , an mode describing rotation becomes unstable (see Gundlach (2002b); see also Gundlach and Baumgarte (2018) for numerical results). Combining these results for those for , we see that for ultra-relativistic fluids, all aspherical modes are stable for only (see also the discussion in Gundlach (2002b)), meaning that we can expect scaling laws for generic initial data to hold to arbitrarily small scales in this regime only. From a physics or astrophysics perspective, the most relevant example of an ultra-relativistic fluid is that of a radiation fluid with , which is almost exactly in the center of the above regime.

###### Acknowledgements.

It is a great pleasure to thank Carsten Gundlach for numerous elucidating conversations, for his comments on an earlier version of this manuscript, and for providing numerical values for the perturbative values of and . JC would like to thank Bowdoin College for hospitality. This work was supported in part by CAPES/Brazil through the Sandwich Doctorate Fellowship program, as well as by NSF grants PHY-1402780 and 1707526 to Bowdoin College. Numerical simulations were performed on the Bowdoin Computational Grid.## References

- Choptuik (1993) M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993).
- Maison (1996) D. Maison, Phys. Lett. B 366, 82 (1996).
- Gundlach (2003) C. Gundlach, Phys. Rept. 376, 339 (2003).
- Gundlach and Martín-García (2007) C. Gundlach and J. M. Martín-García, Living Rev. Rel. 10, 5 (2007).
- Evans and Coleman (1994) C. R. Evans and J. S. Coleman, Phys. Rev. Lett. 72, 1782 (1994).
- Koike et al. (1995) T. Koike, T. Hara, and S. Adachi, Phys. Rev. Lett. 74, 5170 (1995).
- Neilsen and Choptuik (2000) D. Neilsen and M. Choptuik, Class. Quantum Grav. 17, 761 (2000).
- Abrahams and Evans (1993) A. M. Abrahams and C. R. Evans, Phys. Rev. Lett. 70, 2980 (1993).
- Choptuik et al. (2003) M. W. Choptuik, E. W. Hirschmann, S. L. Liebling, and F. Pretorius, Phys. Rev. D 68, 044007/1 (2003).
- Choptuik et al. (2004) M. W. Choptuik, E. W. Hirschmann, S. L. Liebling, and F. Pretorius, Phys. Rev. Lett. 93, 131101 (2004).
- Sorkin (2011) E. Sorkin, Class. Quantum Grav. 28, 025011 (2011).
- Gundlach (1998) C. Gundlach, Phys. Rev. D 57, R7080 (1998).
- Gundlach (2002a) C. Gundlach, Phys. Rev. D 65, 064019/1 (2002a).
- Martín-García and Gundlach (1999) J. M. Martín-García and C. Gundlach, Phys. Rev. D 59, 064031/1 (1999).
- Gundlach (2002b) C. Gundlach, Phys. Rev. D 65, 084021/1 (2002b).
- Pretorius (2005) F. Pretorius, Phys. Rev. Lett. 95, 121101/1 (2005).
- Campanelli et al. (2006) M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96, 111101/1 (2006).
- Baker et al. (2006) J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Phys. Rev. Lett. 96, 111102/1 (2006).
- Alcubierre et al. (2000) M. Alcubierre, G. Allen, B. Brügmann, G. Lanfermann, E. Seidel, W.-M. Suen, and M. Tobias, Phys. Rev. D 61, 041501(R) (2000).
- Hilditch et al. (2013) D. Hilditch, T. W. Baumgarte, A. Weyhausen, T. Dietrich, B. Brügmann, P. J. Montero, and E. Müller, Phys. Rev. D 88, 103009/1 (2013).
- Healy and Laguna (2014) J. Healy and P. Laguna, Gen. Relativ. Gravit. 46, 1722/1 (2014).
- Deppe et al. (2018) N. Deppe, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, ArXiv e-prints (2018), eprint 1802.08682.
- Montero and Cordero-Carrión (2012) P. J. Montero and I. Cordero-Carrión, Phys. Rev. D 85, 124037/1 (2012).
- Baumgarte et al. (2013) T. W. Baumgarte, P. J. Montero, I. Cordero-Carrión, and E. Müller, Phys. Rev. D 87, 044026/1 (2013).
- Montero et al. (2014) P. J. Montero, T. W. Baumgarte, and E. Müller, Phys. Rev. D 89, 084043/1 (2014).
- Baumgarte et al. (2015) T. W. Baumgarte, P. J. Montero, and E. Müller, Phys. Rev. D 91, 064035/1 (2015).
- Ruchlin et al. (2018) I. Ruchlin, Z. B. Etienne, and T. W. Baumgarte, Phys. Rev. D 97, 064036 (2018).
- Mewes et al. (2018) V. Mewes, Y. Zlochower, M. Campanelli, I. Ruchlin, Z. B. Etienne, and T. W. Baumgarte, Phys. Rev. D 97, 084059 (2018).
- Baumgarte and Montero (2015) T. W. Baumgarte and P. J. Montero, Phys. Rev. D 92, 124065/1 (2015).
- Baumgarte and Gundlach (2016) T. W. Baumgarte and C. Gundlach, Phys. Rev. Lett. 116, 221103 (2016).
- Gundlach and Baumgarte (2016) C. Gundlach and T. W. Baumgarte, Phys. Rev. D 94, 084012 (2016).
- Gundlach and Baumgarte (2018) C. Gundlach and T. W. Baumgarte, Phys. Rev. D 97, 064006 (2018).
- Garfinkle and Duncan (1998) D. Garfinkle and G. C. Duncan, Phys. Rev. D 58, 064024/1 (1998).
- Bonazzola et al. (2004) S. Bonazzola, E. Gourgoulhon, P. Grandclément, and J. Novak, Phys. Rev. D 70, 104007/1 (2004).
- Shibata et al. (2004) M. Shibata, K. Uryū, and J. L. Friedman, Phys. Rev. D 70, 044044/1 (2004).
- Brown (2009) J. D. Brown, Phys. Rev. D 79, 104029/1 (2009).
- Gourgoulhon (2012) E. Gourgoulhon, 3+1 formalism in general relativity (Springer, New York, 2012).
- Nakamura et al. (1987) T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987).
- Shibata and Nakamura (1995) M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
- Baumgarte and Shapiro (1998) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 024007/1 (1998).
- Alcubierre (2008) M. Alcubierre, Introduction to 3+1 Numerical Relativity (Oxford University Press, New York, 2008).
- Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro, Numerical relativity: Solving Einstein’s equations on the computer (Cambridge University Press, Cambridge, 2010).
- Bona et al. (1995) C. Bona, J. Massó, E. Seidel, and J. Stela, Phys. Rev. Lett. 75, 600 (1995).
- Alcubierre et al. (2003) M. Alcubierre, B. Brügmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 67, 084023 (2003).
- Thierfelder et al. (2011) M. Thierfelder, S. Bernuzzi, and B. Brügmann, Phys. Rev. D 84, 044012 (2011).
- (46) See www.bowdoin.edu/tbaumgar/radfluidmovie.mp4 for an animation of these simulations.