Black hole formation through fragmentation of toroidal polytropes

We investigate new paths to black hole formation by considering the general relativistic evolution of a differentially rotating polytrope with toroidal shape. We find that this polytrope is unstable to nonaxisymmetric modes, which leads to a fragmentation into self-gravitating, collapsing components. In the case of one such fragment, we apply a simplified adaptive mesh refinement technique to follow the evolution to the formation of an apparent horizon centered on the fragment. This is the first study of the one-armed instability in full general relativity.

The formation of black holes from neutron stars, iron cores or supermassive stars is expected to be associated with a characteristic gravitational wave signal which may give information about the collapse dynamics and the physical environment of such objects.Therefore, and given that gravitational wave detectors are already taking data or are coming online, it is of prime importance to understand the dynamical features of the gravitational collapse of hydrodynamical systems.
The prototypical model of stellar collapse is an equilibrium polytrope subject to a radial or quasi-radial perturbation growing on a dynamical timescale.In spherical symmetry, every general relativistic polytrope with index N = 3 is unstable to radial oscillations [1] -in turn, there exists a critical N c < 3 for which the star is marginally stable.Without spherical symmetry, rotation can increase this critical value again [2].The black hole formation from the collapse of uniformly and differentially rotating polytropes induced by this instability is a well-investigated phenomenon, either with restriction to axisymmetry [3,4,5,6,7,8,9] or without [10,11,12,13].In the gauge choices usually employed, the dynamical behaviour of the system shows a radial contraction of the star, accompanied by the formation of an apparent horizon at late times.
Black hole formation from a dynamically unstable nonaxisymmetric mode, however, has not been modelled so far.Scenarios range from the development of a bar mode, subsequent transport of angular momentum into the shell and collapse of the central object, to fragmentation and off-center production of one or several black holes.In Newtonian theory, instabilities and fragmentation have received considerable attention, specifically in the context of binary formation from protostellar disks (e.g.[14,15,16,17,18] and references therein) and compact object production in stellar core collapse (e.g.[19,20,21] and references therein).In [11], the authors also report signatures of an m = 4 fragmentation behaviour in the collapse and centrifugal bounce of an N = 1 polytrope, but could not determine the final state due to resolution issues.
The cooling evolution of supermassive stars can be approximately described by the N = 3 mass-shedding sequence when the angular momentum transport timescales are short compared to the cooling timescale [22], so that uniform rotation is enforced.This sequence has a turning point for the onset of a quasi-radial instability, and numerical experiments confirm that the collapse remains axisymmetric [23].If the star is differentially rotating, the cooling sequence is less constrained and might end in a transition to nonaxisymmetric instability [24,25,26].The canonical expectation that a supermassive star produces one central black hole with a low-mass accretion disk might thus not be appropriate for differentially rotating configurations.
In this letter, we consider the production of a black hole through the fragmentation of a general relativistic polytrope.We focus on N = 3 polytropes, which are associated with precollapse cores of massive stars or supermassive stars.The soft equation of state also enhances the instability of the fragments compared to the common choice N = 1 for neutron stars.To represent this process accurately on a grid, we make use of an adaptive mesh refinement technique, since a possibly highly deformed apparent horizon needs to be located in some region of the domain which is unknown in advance.
The recent investigation of the collapse of differentially rotating supermassive stars by Saijo [27] was based on a sequence of relativistic N = 3 polytropes with a parameterized rotation law of the commonly used form j(Ω) = A 2 (Ω c − Ω), where Ω c is the angular velocity at the center, and the parameter A specifies the degree of differential rotation (A → ∞ is uniform rotation).The sequence selected was constrained by a constant central density ρ c = 3.38 × 10 −6 in units K = G = c = 1, and the choice A/r e = 1/3, where r e denotes the equatorial coordinate radius.
To examine the indirect collapse by fragmentation of a polytrope with toroidal shape, we choose a model with the same central density as in Saijo's [27] models, but with a ratio of polar to equatorial coordinate radius r p /r e = 0.24.The ratio of rotational kinetic energy to gravitational binding energy is T /|W | = 0.227.While the critical limit for the dynamical f-mode instability in uniform density, uniformly rotating Maclaurin spheroids is (T /|W |) dyn = 0.2738 (e.g.[28]), recent investigations of the stability of soft (N ∼ 3) differentially rotating polytropes in Newtonian gravity [29,30,31,32] have shown that the Maclaurin approximation is inappropriate for such systems, and generally find the critical (T /|W |) dyn to be below the Maclaurin value.Consequently, the toroid-like star considered in this study might be unstable to nonaxisymmetric perturbations.Here we present the first investigation of this instability in full general relativity, showing that relativistic effects are significant for the final outcome, as we observe that black holes can be produced.
All simulations have been performed in full general relativity.The only assumption on symmetry is a reflection invariance with respect to the equatorial plane of the star.The gauge freedom is fixed by the generalized 1+log slicing condition for the lapse function [33] with f (α) = 2/α, and by the hyperbolic-type condition suggested in [8] for the shift vector.
The computational framework is the Cactus code (www.cactuscode.org),which also provides a module to solve the geometric part of the field equations in the well-known BSSN form [34,35,36].In addition, the Carpet driver [37] is used for mesh refinement in Cactus.The hydrodynamics part of the field equations is evolved using the high-resolution shock-capturing PPM-Marquina implementation in the Whisky module [12,38], and a gamma law equation of state (P = ρǫ/N ).We are thus using a set of well-tested tools to evolve the general relativistic hydrodynamics equations.
To numerically construct the axisymmetric initial model described in the introduction, we use the RNS initial data solver [39] with a radial resolution of 601 and an angular resolution of 301 points.With the parameters described above, the model has toroid-like structure, with an off-center density maximum, but a non-zero central density.After mapping the model to the hierarchy of Cartesian grids provided by Carpet, a small perturbation of the form is applied with λ m = 0, 1 and λ = i λ i .In addition, the polytropic constant K is reduced by 0.1% to induce collapse if the model is radially unstable.After perturbing the model, the constraint equations are not solved again, since the amplitude B is chosen such that the violation of the constraints by the initial perturbation is about an order of magnitude smaller than that caused by the systematic error induced by the m = 4 symmetry of the Cartesian grid.
For most simulations, a fixed box-in-box mesh refinement with 5 levels is used to accurately resolve the central highdensity ring.The three innermost grids cover the star, while the two outermost ones push the outer boundaries to 6.4r e .The typical resolution used was 65 × 65 × 33 per grid patch, leading to a central resolution of δ x ≈ 10 −2 r e .However, runs with 89 × 89 × 45, 97 × 97 × 52 and 131 × 131 × 65 points per grid patch were also performed to test the resolution dependence of the solution; for the last setup, a simulation with a uniform grid setup would need to cover the equatorial plane of the star alone with 320 grid points to achieve the same central resolution.
To determine the amplitude of a specific mode in the equatorial plane, we perform a projection onto Fourier modes at certain coordinate radii [40,41].Care must be taken in interpreting the results as soon as the system starts to deviate significantly from axisymmetry, since the interpretation of the projection curve as a circle assumes ∂ φ to be a Killing vector.In addition to the Fourier extraction, we monitor the evolution of the rest mass and the ADM constraints.
In Table I, we have listed the parameters for the different simulation runs.For model N, the evolution of the moduli of the equatorial Fourier components at the initial radius of highest density is shown in Fig 1 .It is evident that, initially, the m = 4 component (thin dotted line) induced by our Cartesian grid is dominant.However, the star is unstable to m = 1 (thick solid line) and m = 2 (thick dashed line), and these modes consequently grow into the nonlinear regime, their efolding times being rather close.
To test the dependence of the results on the grid resolution, we have evolved the same initial data with different grid parameters as listed in Table I.Runs H1, H2 and H3 are high resolution versions of N. The setups I1 and I2 test the dependence on the amplitude B of the initial data perturbation, with I1 using an amplitude 10 times smaller than in setup N, and I2 using B = 0. M1 and M2 are variants of N, where only a specific unstable mode is imposed.
The rest mass is conserved numerically within 1.8% for setup N and to within 0.2% for setup H3.An approximate measurement of the e-folding times and mode frequencies can be obtained within an error of 5 − 10% related to ambiguities  To establish whether a black hole is formed by a fragment it is necessary to cover the fragment with significantly more resolution than affordable by fixed mesh refinement.Hence we have implemented a simplified adaptive mesh refinement scheme to follow the system to black hole formation: In this scheme, a tracking function, here provided by the location of a density maximum, is used to construct a locally fixed hierarchy of grids moving with the fragment.Additional refinement levels are switched on during contraction, until an apparent horizon is found.
Since the e-folding times for m = 1 and m = 2 turn out to be close, the number and interaction behaviour of the fragments in the non-linear regime depend sensitively on the initial perturbation.Thus setups M1 and M2 are used the follow the formation and evolution of a specific number of fragments.
The time evolution of the equatorial plane density for setup M1 is shown in Fig 2 .While the initial model is axisymmetric, it has already developed a strong m = 1 type deviation from axisymmetry at t = 6.43tD , which consequently evolves into a collapsing off-center fragment.At t = 7.45t D , we find an apparent horizon, using the numerical code described in [42].The horizon is centered on the collapsing fragment at a coordinate radius of r AH ≈ 0.16r e , and has an irreducible mass of M AH ≈ 0.24M star .Its coordinate representation is significantly deformed: its shape is close to ellipsoidal, with an axes ratio of ∼ 2 : 1.1 : 1.The apparent horizon is covered by three refinement levels and 50 to 100 grid points along each axis.
The evolution of the setup M2 is shown in Fig 3 .Two orbiting and collapsing fragments are forming.However, even with the adaptive mesh refinement method we use, constraint violations prevent us from continuing the simulation to the FIG.2: Time evolution of the equatorial plane density for setup M1.Shown are isocontours of the logarithm of the rest-mass density.The four snapshots extend to 0.37re and are taken at t/tD = 0, 6.43, 7.14, and 7.45, respectively.They show the formation and collapse of the fragment produced by the m = 1 instability.The last slice contains an apparent horizon demarked by the thick white line.Note that the shades of grey used for illustration are adapted to the current maximal density at each time, and that darker shades denote higher densities.formation of apparent horizons.Cell-based adaptive mesh refinement, a better choice of gauge, or methods based on numerical analysis (e.g.[43]) might be required in this case.
To summarize, we have studied fragmentation and black hole formation of a general relativistic equilibrium polytrope of index N = 3.The polytrope has been shown to be unstable to cylindrical perturbations of the form r sin(mφ), with m = 1, 2, which consequently grow to one or more selfgravitating fragments.We have applied an adaptive mesh refinement method to resolve the system accurately.In the m = 1 case, we have found an apparent horizon in the spacetime, indicating that a black hole has formed.
The dynamics of a nonaxisymmetric single-star collapse of this type is significantly different from that of the quasi-radial cases usually investigated.From the often considered case of quasi-radial collapse, to bar formation and subsequent collapse, to fragmentation and fragment inspiral we have a range of possible dynamical scenarios, which may be connected to discernable observable features in their gravitational wave signature.In that sense, the evolution presented here can be considered as an example of such processes.
Extensive discussion of this work as well as further results for the binary case will be presented elsewhere.In this case, the fragmentation process could transform a star into a binary black hole merger system with a massive accretion disk around it.

FIG. 1 :
FIG. 1: Time evolution of the mode amplitudes in the standard grid setup N. The amplitudes are obtained from a Fourier decomposition of the density profile on the equatorial plane circle at ̟ = 0.25re, the initial radius of highest density.Shown are the m = 1 (thick solid line), m = 2 (thick dashed line), m = 3 (thin dash-dotted line) and m = 4 (thin dotted line) mode amplitudes.The time is normalized to the dynamical timescale tD = re re/M .

FIG. 3 :
FIG. 3: Time evolution of the equatorial plane density for setup M2.The snapshot times are the same as in Fig 2. In this case, two fragments are forming.Constraint violations have forced us to terminate the simulation before apparent horizons could be located, as explained in the text.

TABLE I :
Parameters for the simulation runs.Different setups were chosen to confirm the results, including resolution tests, changes in the initial perturbation, and two setups with adaptive mesh refinement (AMR) to investigate black hole formation.