Three-dimensional relativistic simulations of rotating neutron-star collapse to a Kerr black hole

We present a new three-dimensional fully general-relativistic hydrodynamics code using high-resolution shock-capturing techniques and a conformal traceless formulation of the Einstein equations. Besides presenting a thorough set of tests which the code passes with very high accuracy, we discuss its application to the study of the gravitational collapse of uniformly rotating neutron stars to Kerr black holes. The initial stellar models are modeled as relativistic polytropes which are either secularly or dynamically unstable and with angular velocities which range from slow rotation to the mass-shedding limit. We investigate the gravitational collapse by carefully studying not only the dynamics of the matter, but also that of the trapped surfaces, i.e., of both the apparent and event horizons formed during the collapse. The use of these surfaces, together with the dynamical horizon framework, allows for a precise measurement of the black-hole mass and spin. The ability to successfully perform these simulations for sufficiently long times relies on excising a region of the computational domain which includes the singularity and is within the apparent horizon. The dynamics of the collapsing matter is strongly influenced by the initial amount of angular momentum in the progenitor star and, for initial models with sufficiently high angular velocities, the collapse can lead to the formation of an unstable disc in differential rotation. All of the simulations performed with uniformly rotating initial data and a polytropic or ideal-fluid equation of state show no evidence of shocks or of the presence of matter on stable orbits outside the black hole.

Figure 1

Gravitational mass shown as a function of the central energy density for equilibrium models constructed with the polytropic EOS, for $\Gamma =2$ and polytropic constant $K_{\mathrm{I}\mathrm{D}}=100$. The solid, dashed and dotted lines correspond to the sequence of nonrotating models, the sequence of models rotating at the mass-shedding limit and the sequence of models that are at the onset of the secular instability to axisymmetric perturbations. Also shown are the secularly (open circles) and dynamically unstable (filled circles) initial models used in the collapse simulations.

Figure 2

$L_2$ norm of the Hamiltonian constraint violation for the initial model D1 shown as a function of time. The different lines refer to different grid resolutions, but in all cases the IVP was solved after the pressure was uniformly decreased to trigger the collapse.

Figure 3

Collapse sequence for the slowly rotating model D1. Different panels refer to different snapshots during the collapse and show the isocontours of the rest-mass density and velocity field in the $(x,y)$ plane (left column) and in the $(x,z)$ plane (right column), respectively. The isobaric surfaces are logarithmically spaced and a reference length for the vector field is shown in the lower-right panel for a velocity of $0.2c$. The time of the different snapshots appears in the top right corner of each panel and it is given in ms, while the units on both axes are expressed in km. See main text for a discussion and Fig. 5 for a comparison with the collapse of a rapidly rotating model.

Figure 4

Magnified view of the final stages of the collapse of model D1. A region around the singularity that has formed is excised from the computational domain and is shown as an area filled with squares. Also shown with a thick dashed line is the coordinate location of the apparent horizon. Note that because of rotation this surface is not a two-sphere, although the departures are not significant and cannot be appreciated from the Figure (cf. Figure 16 for a clearer view).

Figure 5

Collapse sequence for the rapidly rotating model D4. The conventions used in these panels are the same as in Fig. 3, which can be used for a comparison with the collapse of a slowly rotating model. Note that a region around the singularity that has formed is excised from the computational domain and is indicated as an area filled with squares. Also shown with a thick dashed line is the coordinate location of the apparent horizon.

Figure 6

Magnified view of the final stages of the collapse of model D4. Note that a region around the singularity that has formed is excised from the computational domain and this is indicated as an area filled with squares. Also shown with a thick dashed line is the coordinate location of the apparent horizon. Note that because of the rapid rotation, this surface has significant departures from a two-sphere (cf. Figure 16 for a clearer view).

Figure 7

Ratio of the proper polar radius to the proper equatorial radius for all the initial models. Each curve ends at the time when, for each simulation, all the matter along the $z$-axis has fallen below the apparent horizon.

Figure 8

Rest-mass density of model D4 normalized to the initial value at the stellar center. The profiles are measured along the $x$-axis on the equatorial plane and refer to different times (see main text for details). Line 5, shown as dotted, corresponds to the time when the apparent horizon is first found. The inset shows a magnified view of the final stages of the evolution using a logarithmic scale and also the location of the excised region as it grows in time.

Figure 9

Evolution of the mass fraction versus time during the collapse of D1 (left panel) and D4 (right panel). The rest mass is measured within two-spheres of coordinate radii $r_i=0.2,0.4$ and 0.6 times the initial stellar equatorial circumferential radius $R_e$ (cf. Table ). Marked with filled dots on the different lines are the times at which the apparent horizon is first found (the data refers to a simulation with ${96}^2\times 48$ gridzones). The insets in both panels show, on a logarithmic scale, the evolution of the normalized rest mass outside the apparent horizon. Note that this is appreciably nonzero for a rather long time in case of model D4.

Figure 10

Evolution of the averaged angular velocity $⟨\Omega ⟩$ (left panel) and of the averaged angular momentum per unit mass $⟨h{u}_{\varphi }⟩$ (right panel). Both quantities are measured at the stellar equator, are normalized to the initial value at the stellar surface and refer to both models D1 (upper parts) and D4 (lower parts).

Figure 11

Convergence of the measure of the black-hole mass as the resolution is increased. The curves refer to estimates using the event-horizon equatorial circumference [i.e. Eq. (5.1)] and have been obtained using ${288}^2\times 144$ and ${192}^2\times 96$ zones, respectively. Shown in the small inset are the results for model D1, while those for model D4 are in the main panel.

Figure 12

Evolution of the event-horizon mass $M=C_{\mathrm{e}\mathrm{q}}/4\pi$ for models D1 and D4. Different lines refer to the different initial guesses for the null surface and are numbered 0, 1 and 2. Note that they converge exponentially to the correct event-horizon surface for decreasing times.

Figure 13

Fitting the oblateness of the event horizon to QNMs of a Kerr black hole. The fit is shown with the solid line, while the open circles represent the computed values of $C_r$. The estimate for $C_r$ of a Kerr black hole having the fitted value of $a/M_{hor}$ is shown with a dashed line.

Figure 14

Comparison of the different measures of the angular momentum for the cases D1 (left panel) and D4 (right panel). The estimate using the fit to the circumference ratio (see left panel of Fig. 13) is also shown. The dynamical-horizon spin measurement is considerably more accurate at late times as the event-horizon surfaces will diverge exponentially at this point. Shown with the horizontal short-long dashed lines are the values of $(J/M^2{)}_{\mathrm{A}\mathrm{D}\mathrm{M}}$ in the two cases as measured from the initial data (see main text for details).

Figure 15

Comparison of the five different approaches used in the measure of the total black-hole mass for the collapse of models D1 and D4. Different lines refer to the different methods discussed in the main text. The left panel (model D1), shows that the different methods are overall comparable when the rotation is slow, but that differences are already present (this is as shown in the inset). The right panel (model D4) shows that the different measures can be considerably different when the rotation is large.

Figure 16

Evolution of the most relevant surfaces during the collapse for the cases D1 and D4. Solid, dashed and dotted lines represent the worldlines of the circumferential radii of the event horizon, of the apparent horizon and of the stellar surface, respectively. Note that for the horizons we plot both the equatorial and the polar circumferential radii, while only the equatorial circumferential radius is shown for the stellar surface. Shown in the insets are the magnified views of the worldlines during the final stages of the collapse.

Figure 17

Solution of a Riemann problem set on the main diagonal of the cubic grid. The figure shows the comparison of the hydrodynamical variables evolved by Whisky, indicated with symbols, with the exact solution. The numerical simulation was obtained with the van Leer reconstruction method and the Roe solver, on a grid with ${140}^3$ points.

Figure 18

Central mass-density, normalized to the initial value, in a stable TOV star ($M=1.4M_{\odot }$ and polytropic exponent $\Gamma =2$) evolution at different resolutions. PPM and Marquina were used for all runs.

Figure 19

Power spectrum of the central mass-density evolution of an $M=1.4M_{\odot }$, $\Gamma =2$ stable TOV star performed with ${128}^3$ grid points. The units of the vertical axis are arbitrary.

Figure 20

Normalized central mass-density evolution of an $M=1.4M_{\odot }$, $\Gamma =2$ unstable TOV star performed with ${96}^3$ grid points.