Nonaxisymmetric instability and fragmentation of general relativistic quasitoroidal stars

In a recent publication, we have demonstrated that differentially rotating stars admit new channels of black hole formation via fragmentation instabilities. Since a higher order instability of this kind could potentially transform a differentially rotating supermassive star into a multiple black hole system embedded in a massive accretion disk, we investigate the dependence of the instability on parameters of the equilibrium model. We find that many of the models constructed exhibit nonaxisymmetric instabilities with corotation points, even for low values of $T/|W|$, which lead to a fission of the stars into one, two, or three fragments, depending on the initial perturbation. At least in the models selected here, an $m=1$ mode becomes unstable at lower values of $T/|W|$, which would seem to favor a scenario where one black hole with a massive accretion disk forms. In this case, we have gained evidence that low values of compactness of the initial model can lead to a stabilization of the resulting fragment, thus preventing black hole formation in this scenario.

Figure 1

Time evolution of the $L^2$ norm of the Hamiltonian constraint for different resolutions (the numbers refer to a grid point on a single patch of our mesh-refinement grid). The time is normalized to the dynamical time scale $t_{\mathrm{dyn}}=R_e\sqrt{R_e/M}$.

Figure 2

Decadic logarithm of the density in a meridional plane of the model constructed with the parameters in Table . The model is of quasitoroidal nature.

Figure 3

${\rho }_{\max }/{\rho }_c$ (top) and $T/|W|$ (bottom) for the $R$ sequence. Beyond an axis ratio of $r_p/r_e\approx 0.3423$, the initial data solver converges to spheroidal models.

Figure 4

Polytropic models constructed with the RNS code. The plots show the parameter plane spanned by the axis ratio $r_p/r_e$ and the differential rotation parameter $A$, constrained by ${\rho }_c={10}^{-7}$ and $\Gamma =4/3$, and resolved by construction of 6400 ($80\times 80$) models. Top left: contours of the equatorial stellar over Keplerian angular velocity ${\Omega }_e/{\Omega }_K$ (the mass-shedding limit is the isosurface ${\Omega }_e/{\Omega }_K=1$). Top right: decadic logarithm of the ratio of maximal over central density ${log}_{10}({\rho }_{\max }/{\rho }_c)$. Middle left: rotation parameter $T/|W|$. Middle right: decadic logarithm of the inverse equatorial compactness $R_e/M$. Bottom left: normalized angular momentum $J/M^2$. Bottom right: Selection of evolved initial models. The thick continuous line marks the jump between quasitoroidal and spheroidal models apparent in the top left to bottom left plots, the thick dashed line is the mass-shedding limit, and the thin dashed line indicates an approximate division between spheroidal and quasitoroidal models. The models selected for numerical evolution are marked by circles.

Figure 5

Development of the fragmentation instability in the reference polytrope. The left sequence of plots shows the decadic logarithm of rest-mass density in the equatorial plane, the right sequence of plots shows the lapse function. The snapshots correspond to times $t/t_{\mathrm{dyn}}=0$ (top), 6.28 (middle), and 7.48 (bottom). The initial model, perturbed by Eq. (17) with ${\lambda }_m=1$ for $m=1,\dots ,4$, develops a spiral-arm instability and a collapsing fragment. The domain outside the star initially has the density of the artificial atmosphere (${\rho }_{\mathrm{atm}}\approx 5\times {10}^{-10}$), and the artifacts outside the star are caused by interactions of the atmosphere with the stellar surface and the outer boundaries of the computational domain. Here, and in all evolution sequences which follow, the extent of the spatial coordinate domain plotted is the same in all snapshots. Also, note that the color map is adapted to the range of function values in each plot.

Figure 6

Mode amplitudes versus time, extracted at $r=0.25R_e$, the radius of highest rest-mass density in the initial model, for different initial perturbations. The amplitude $A_m$ is the $m$th harmonic Fourier projection of the density, normalized to the average value. The perturbations corresponding to each plot are (cf. also Eq. (17)): ${\lambda }_m=1$ for $m=1,\dots ,4$ (top left), ${\lambda }_m={\delta }_{m1}$ (top right), ${\lambda }_m={\delta }_{m2}$ (middle left), ${\lambda }_m={\delta }_{m3}$ (middle right), ${\lambda }_m={\delta }_{m4}$ (bottom left), and ${\lambda }_m=0$ for $m=1,\dots ,4$. For details see text.

Figure 7

Comparison of mode amplitude versus time for different initial perturbations. The upper panel shows the mode amplitude $A_1$, and the lower one shows $A_2$. Note that, in the case ${\lambda }_m={\delta }_{m1}$, the $m=2$ mode is dominated by nonlinear effects from the fragmentation.

Figure 8

Similar to Fig. 5, but now for a perturbation ${\lambda }_m={\delta }_{m1}$. Shown is the decadic logarithm of the density in the equatorial plane. The snapshots correspond to times $t/t_{\mathrm{dyn}}=6.28$ (left), 7.11 (middle), and 7.48 (right). While the $m=2$ mode is now suppressed, the qualitative evolution is similar to that displayed in Fig. 5.

Figure 9

Similar to Fig. 5, but for a perturbation ${\lambda }_m={\delta }_{m2}$. Shown is the decadic logarithm of the density in the equatorial plane. The snapshots correspond to times $t/t_{\mathrm{dyn}}=6.28$ (left), 7.66 (middle), and 8.85 (right). Two fragments develop and encounter a runaway instability while orbiting each other.

Figure 10

Local growth times ${\tau }_m(t)≔(d\ln A/dt{)}^{-1}$ (top) and frequencies $d{\varphi }_m/dt$ (bottom) for modes $m=1$ and $m=2$, both taken from simulations with the corresponding initial perturbation. The plots denoted by “time averaged” are running averages over $t_{\mathrm{dyn}}$.

Figure 11

Contour plot of the function ${log}_{10}(|\rho (x,y,z,t)-{\rho }_0(x,y,z)|+ϵ)$ in the equatorial plane, where $\rho$ denotes the rest-mass density of the $m=1$ evolution shown in Fig. 8, ${\rho }_0$ denotes the density function of the equilibrium polytrope, and $ϵ$ is a small number ($ϵ={10}^{-10}$ in this plot).

Figure 12

Evolution of the mode amplitudes $A_1$ and $A_2$ in the reference polytrope perturbed by ${\lambda }_m={\delta }_{m1}$, for different mode extraction radii. The radius of highest initial density is indicated.

Figure 13

Same as Fig. 12, but for the mode amplitudes $A_3$ and $A_4$.

Figure 14

Local mode frequency $t_{\mathrm{dyn}}d\varphi /dt$ for the modes $m=1$ (top) and $m=2$ (bottom) in the reference polytrope, extracted at different radii in the equatorial plane.

Figure 15

Angular velocity of the reference polytrope over the $x$ axis (black line), and approximate location (with error bar) of the pattern speed of the $m=1$ mode (red, upper rectangle), and the $m=2$ mode (blue, lower rectangle). Both modes are inside the corotation band.

Figure 16

Evolution of the mode amplitudes $A_1$ (top) and $A_2$ (bottom) for different grid resolutions. The grid sizes in the legend refer to the four outer grid patches; the innermost patch covering the high-density central toroidal region of the star has a resolution of $65\times 65\times 33$, $97\times 97\times 49$, $129\times 129\times 65$, or $193\times 193\times 97$, correspondingly.

Figure 17

Same as Fig. 16, but for the evolution of the maximum of the rest-mass density ${\rho }_{\max }$.

Figure 18

Same as Fig. 16, but for the evolution of the total rest mass $M_0$. The values are normalized to the initial rest mass $M_0^{\mathrm{ID}}$ obtained from the initial data solver. The last graph is from a simulation with a less dense artificial atmosphere.

Figure 19

Evolution of the mode amplitudes $A_1$ and $A_2$ in the reference polytrope, for artificial atmospheres of different density ${\rho }_{\mathrm{atm}}/{\rho }_c$.

Figure 20

Evolution of the mode amplitudes $A_1$ and $A_2$ for different members of the $R$ sequence (cf. Table ).

Figure 21

Evolution of the mode amplitudes $A_1$ and $A_2$ for different members of the $G$ sequence (cf. Table ).

Figure 22

Evolution of the mode amplitudes $A_1$ and $A_2$ for different members of the $C$ sequence (cf. Table ).

Figure 23

Evolution of the maximum of the rest-mass density $\rho$ and the minimum of the lapse function $\alpha$ for different members of the $C$ sequence.

Figure 24

Evolution of density in the equatorial plane of the model $C8$. Shown is the decadic logarithm of the rest-mass density. The snapshots were taken at times $t/t_{\mathrm{dyn}}=0$ (top), 6.28 (middle), and 8.28 (bottom). In contrast to the more compact model $C1$, which is the reference polytrope investigated earlier (cf. Fig. 5), the fragment reexpands after a maximal compression.

Figure 25

Stability of quasitoroidal models with ${\rho }_c={10}^{-7}$ (cf. Fig. 4). A Latin number denotes the highest azimuthal order of the unstable modes, i.e. I for $m=1$ unstable, II for $m=1$, 2 unstable, and III for $m=1$, 2, 3 unstable. Models denoted by “(I)” are either long-term unstable with growth times $\tau \gg t_{\mathrm{dyn}}$, or stable (see text), and models denoted by A exhibit an axisymmetric instability. The (red) line in the lower left is the approximate location of the sequence $J/M^2=1$ (cf. Fig. 4), and the three (blue) lines inside the quasitoroidal region are the approximate locations of sequences with $T/|W|=0.14$ (right), $T/|W|=0.18$ (middle), and $T/|W|=0.26$ (left).

Figure 26

Remnants of the models from Fig. 25, which are unstable with respect to nonaxisymmetric modes. The nonlinear behavior has been analyzed by observing the evolution of the function $\min \alpha$ (see also Fig. 23). Models which show a minimum in this function are marked by B for “bounce,” while models exhibiting an exponential collapse of the lapse are marked by C for “collapse.”

Figure 27

Mode amplitudes in the model $A0.1R0.15$ (cf. Table ), extracted at the radius of highest initial rest-mass density.

Figure 28

Mode amplitudes in the model $A0.1R0.50$ (cf. Table ), extracted at the radius of highest initial rest-mass density.

Figure 29

Mode amplitudes in the model $A0.3R0.15$ (cf. Table ), extracted at the radius of highest initial rest-mass density.

Figure 30

Mode amplitudes in the model $A0.3R0.50$ (cf. Table ), extracted at the radius of highest initial rest-mass density.

Figure 31

Mode amplitudes in the model $A0.6R0.15$ (cf. Table ), extracted at the radius of highest initial rest-mass density.

Figure 32

Equatorial density evolution of the model $A0.6R0.15$, with an $m=3$ perturbation. (Note that, in Fig. 31, a perturbation with ${\lambda }_m=1$ has been used instead). Shown is the decadic logarithm of the rest-mass density. The snapshots were taken at times $t/t_{\mathrm{dyn}}=0$ (top), 6.28 (middle), and 7.60 (bottom). Three fragments develop and subsequently encounter collapse similar to the two-fragment case (cf. Fig. 9). The evolution of the model perturbed with $m=1$ and $m=2$ is similar to the corresponding one in the reference polytrope.

Figure 33

Equatorial density evolution of the model $A0.1R0.15$. Shown is the decadic logarithm of the rest-mass density. The snapshots were taken at times $t/t_{\mathrm{dyn}}=0$ (top) and 6.28 (bottom). The model exhibits an axisymmetric instability.

Figure 34

Angular velocity of the polytropes $A0.3R0.15$ (top), $A0.3R0.30$ (middle), and $A0.3R0.40$ (bottom) over the $x$ axis (black line), and approximate location of the pattern speed of the $m=1$ mode (red, upper rectangle) and the $m=2$ mode (blue, lower rectangle), cf. also Fig. 15.

Figure 35

Long-term evolution of the mode amplitudes for the model $A0.2R0.45$, which is unstable to an $m=1$ perturbation. The mode, however, grows rather slowly over a time of $20t_{\mathrm{dyn}}$.

Figure 36

Amplitude of the $m=1$ mode for different models from a sequence with constant $T/|W|$ limiting in the model $A0.2R0.40$, which has ${\rho }_c={10}^{-7}$; cf. also Table .