Type II critical phenomena of neutron star collapse

We investigate spherically symmetric, general relativistic systems of collapsing perfect fluid distributions. We consider neutron star models that are driven to collapse by the addition of an initially “ingoing” velocity profile to the nominally static star solution. The neutron star models we use are Tolman-Oppenheimer-Volkoff solutions with an initially isentropic, gamma law equation of state. The initial values of (1) the amplitude of the velocity profile, and (2) the central density of the star, span a parameter space, and we focus only on that region that gives rise to type II critical behavior, wherein black holes of arbitrarily small mass can be formed. In contrast to previously published work, we find that—for a specific value of the adiabatic index ($\Gamma =2$)—the observed type II critical solution has approximately the same scaling exponent as that calculated for an ultrarelativistic fluid of the same index. Further, we find that the critical solution computed using the ideal-gas equations of state asymptotes to the ultrarelativistic critical solution.

Figure 1

Scaling behavior for supercritical—or black hole-forming—solutions. The top plot illustrates how the points from a series of supercritical runs follow the scaling law for the black hole mass (1), while the bottom plot shows how the data deviate from our best fit to this scaling law. The two dotted lines delineate the data used in making the best fit; this data is plotted separately in Fig. 2. Black holes were assumed to have formed when $\underset{r}{max}(2m(r,t)/r)\ge 0.995$. The gaps between some of the points represent those runs that crashed before $\underset{r}{max}(2m/r)$ reached this value. Smoothing was used for $\ln |p-p^{\star }|<-19.3$, which is also where we start our fit. These runs used ${\rho }_c=0.05$, $U=U_1$ and an initial grid defined by $\{N_a,N_b,N_c,\Delta r_a,\mathrm{\text{level}}\}=\{300,500,20,0.005,0\}$.

Figure 2

Best-fit for the scaling behavior of black hole masses near the threshold. The top plot shows calculated masses and the fitting line, while the bottom plot shows the deviation between the two. The scaling exponent for this fit, which is simply the slope of the line, is $\gamma =0.94$.

Figure 3

Scaling behavior in $T_{\max }$ for subcritical solutions, i.e. those not forming black holes. The line shown here is the best fit for the expected scaling law (28), using data only from the solutions closest to criticality. These runs used ${\rho }_c=0.05$, $U=U_1$ and an initial grid defined by $\{N_a,N_b,N_c,\Delta r_a,\mathrm{\text{level}}\}=\{300,500,20,0.005,0\}$.

Figure 4

Three scale-free quantities of near-critical solutions plotted versus ${\mathcal{X}}_a$ for the ideal-gas system (dashed line) and the ultrarelativistic system (solid line).

Figure 5

Deviation over time of those quantities displayed in Fig. 4. Here, $||f||$ denotes the ${\ell }_2$-norm of the function $f$. The deviations for $a$ (solid line), $\omega$ (dotted line), and $v$ (dashed line) are shown. The ${\ell }_2$-norms of these differences are computed at every time satisfying ${\mathcal{X}}_a<2$, and then logarithms of those norms are plotted as a function of self-similar time $\mathcal{T}$. Note that the sense of physical time is opposite to that of $\mathcal{T}$; that is, $\mathcal{T}\to -\infty$ as the solution approaches the accumulation time. As the evolution proceeds from the initial time, the two solutions asymptote toward each other. For $\mathcal{T}<-13$, the deviation between the two solutions increases as the ideal-gas near-critical solution departs from the asymptotic critical solution and eventually disperses from the origin.

Figure 6

Time sequences of $\omega$ for the most nearly critical solutions obtained with the ideal-gas EOS (dashed line) and the ultrarelativistic EOS (solid line). Both functions have been transformed into self-similar coordinates, based upon their respective accumulation times and respective values of $r_{a_{\max }}$. Approximate values of $\mathcal{T}$ are shown in the upper-left corners of the frames. Note that the ultrarelativistic $\omega$ is varying slightly frame-to-frame, contrary to appearances. Relative to the intrinsically ultrarelativistic solution, it takes more time for the ideal-gas solution to become self-similar since the length scale set by ${\rho }_{\circ }$ in the latter case only becomes insignificant for $P/{\rho }_{\circ }\gg 1$.

Figure 7

Scaling behavior in $T_{\max }$ near the critical solution for runs using different values of $\delta$. The scaling behaviors are shown for the original configuration (circles, solid line), a configuration with ${10}^2$ times the original floor value (squares, dotted line), and one with ${10}^4$ times the original floor value (triangles, dashed line). The scaling exponents $\gamma$ for these runs are listed in Table .

Figure 8

Comparison of the scaling behavior in $T_{\max }$ obtained with two different Riemann solvers. The “smoothed Roe” line corresponds to the original calculations. The other (dotted) line was generated using the Marquina method, with other computational methods and parameters identical to the original calculations. The scaling exponents, $\gamma$, for these runs are listed in Table .

Figure 9

Scaling behavior in $T_{\max }$ near the critical solution for runs using different levels of resolution. The runs were made with ${\rho }_c=0.05$, $U=U_1$, and the solid line with circles was generated from runs using the original configuration. The $\mathrm{\text{level}}=1$ (squares, dotted line) and $\mathrm{\text{level}}=2$ (triangles, dashed line) runs, respectively, used computational grids that were locally 2 times and 4 times as refined. The scaling exponents, $\gamma$, for these runs are listed in Table .

Figure 10

Logarithm of scaled, independent residuals of the Hamiltonian constraint (21) and slicing condition (22) for three levels of resolutions calculated from solutions in the self-similar regime plotted versus ${\mathcal{X}}_a$. The dotted (dashed) lines are from a run which used 2 (4) times the local spatial and temporal resolutions of the original run, which is represented by the solid lines; the dotted (dashed) residual was scaled by a factor of 4 (16) in order to make the $O(\Delta r^2)$ convergence of the solution more apparent. Each distribution is from a solution that has been tuned to $\ln |p-p^{\star }|\simeq -19$ with respect to the value of $p^{\star }$ for each resolution, and function values at every tenth grid point are shown. The physical velocity of the fluid for the $l=0$ run is shown in the bottom frame in order to facilitate comparison of features in the independent residuals to those in the solution.

Figure 11

Scaling behavior in $T_{\max }$ for several families of initial data, which differ only in one respect from the original configuration (filled circles). Open triangles represent the scaling law found from using a different initial velocity profile, $U=U_2$. Filled triangles, $\times$’s, open circles, open squares, and filled squares represent—respectively—scaling laws using ${\rho }_c=0.01$, 0.02, 0.03, 0.04, and 0.0531. The scaling exponents, $\gamma$, for these runs are listed in Table .

Figure 12

Three scale-free quantities of near-critical solutions plotted versus ${\mathcal{X}}_a$ for various initial data sets. The ultrarelativistic solution is the solid curve. The solutions that started with $U=U_1$ and ${\rho }_c=0.01$, 0.02, 0.03, 0.04, 0.05, 0.0 531 are represented, respectively, as dots, short dashes, long dashes, dotted-short dashes, and short-long dashes. The functions from the evolution that started with $U=U_2$ are shown in dotted-long dashes.

Figure 13

Logarithm of the difference between each of the scale-free functions of an ideal-gas EOS near-threshold solution and those from the ultrarelativistic solution. The differences are shown for all initial data sets, and are plotted versus ${\mathcal{X}}_a$. The solutions that started with $U=U_1$ and ${\rho }_c=0.01$, 0.02, 0.03, 0.04, 0.05, 0.0 531 are represented, respectively, as a solid curve, dots, short dashes, long dashes, dotted-short dashes, and short-long dashes. The functions from the evolution that started with $U=U_2$ are shown in dotted-long dashes.

Figure 14

Scaling behavior in $T_{\max }$ for different values of $\Gamma$. The original (solid line with circles) was made from runs with ${\rho }_c=0.05$ and $\Gamma =2$ as before. Both the $\Gamma =1.88$ (dashed line with triangles) and $\Gamma =1.90$ (dotted line with squares) runs used ${\rho }_c=0.013$ with $\{N_a,N_b,N_c,\Delta r_a\}=\{420,700,20,0.005\}$. The scaling exponents are given in the text.

Figure 15

Contour plot of ${log}_{10}$ of the relative error between the exact value of $v$ and the value obtained using our new primitive variable solver. Here, $W=1151$. The color map is shown at the bottom. Note that black regions represent portions of $(P,{\rho }_{\circ })$ space where the new solver results in errors of order 100%, while white regions signal errors that are at the level of machine precision or below.

Figure 16

Comparison of the accuracy achieved with our primitive variable solver and that used in 25. Shown are frames like the one shown in Fig. 15, but for different values of $W$ (rows) and for different variables and methods (columns). The label “Old” indicates that the frames in that column were calculated using the method described in 25, while the label “New” indicates that those frames were calculated using our new primitive solver. The left two frames show ${log}_{10}$ of the relative error in $v$, while the right two frames show ${log}_{10}$ of the relative error in ${\rho }_{\circ }$. The color coding is the same as in Fig. 15.

Figure 17

Conserved variable $D(\mathcal{X},\mathcal{T})$ from the most nearly critical evolution obtained with the use of the approximate Roe solver without smoothing. $\mathcal{X}$ and $\mathcal{T}$ are defined by Eqs. (30, 31), respectively. The dashed line indicates the location of the sonic point, $\mathcal{X}=0$. No refinement takes place during the period shown here, and $\Delta r_a\simeq 1.55\times {10}^{-7}$. From left to right and top to bottom, the $\mathcal{T}$ values of the frames are $-10.4109$, $-10.4977$, $-10.5916$, $-10.6938$, $-10.7822$, $-10.7823$, $-10.7824$, $-10.7825$, $-10.7826$. The evolution started with a TOV star of central density ${\rho }_c=0.05$ that was perturbed using profile $U_1$ [see (26)] with the overall amplitude factor, $U_{\mathrm{amp}}$, tuned to produce near-critical evolution.

Figure 18

Characteristic speeds of the fluid for the most nearly critical solution obtained with the approximate Roe solver without smoothing. The wave speeds are plotted here as functions of the self-similar coordinate $\mathcal{X}$, and are shown at $\mathcal{T}=-10.6938$. A closer view of the characteristic speeds near the sonic point is shown as an inset in the lower right of the plot, revealing the severity of the discontinuity in ${\lambda }_{-}$ as discussed in the text.

Figure 19

Pressure and rest-mass density of the most nearly critical solution obtained with the approximate Roe solver without smoothing. Both quantities are plotted as a function of the self-similar coordinate $\mathcal{X}$, and are shown at $\mathcal{T}=-10.6938$. Interior to the sonic point, $\mathcal{X}=0$, the fluid is clearly in the ultrarelativistic limit, with $P/{\rho }_{\circ }\approx {10}^4$. However, beyond the sonic point, the flow is not ultrarelativistic—in fact, $P<{\rho }_{\circ }$ in most of the domain exterior to $\mathcal{X}=0$. A closer view of the distributions near the sonic point is shown in the inset and more clearly illustrates the formation of an expansion shock as discussed in the text.

Figure 20

Riemann solution using the approximate Roe method with initial data $\{P,v,{\rho }_{\circ }\}(x<0.5)=\{100,0,1\}$, $\{P,v,{\rho }_{\circ }\}(x>0.5)=\{1,0,1\}$, $\Gamma =5/3$, using 200 cells. $P(x)/120$ (circles), ${\rho }_{\circ }(x)/6$ (triangles), and $v(x)$ ($\times$’s) are plotted at $t=0.4$. The lines correspond to the exact solution.

Figure 21

Same as in Fig. 20 except using the Marquina method.

Figure 22

Convergence test for the fluid variable $D$. The top panel shows $\ln |D^{8h}-D^{4h}|$ (solid line), $\ln |4(D^{4h}-D^{2h})|$ (dotted line), $\ln |16(D^{2h}-D^h)|$ (dashed line) which have been scaled so that they will coincide if they are well described by Richardson expansions of the form (e1). Quantities with superscript “$lh$” have been calculated using grid point separations $l$ times larger than the fiducial $l=1$ quantities. The bottom panel shows $D(r,0)$ (dashed line) and $D(r,t)$ (solid line), where $t$ is the time at which we performed the convergence test. The initial data consisted of a self-gravitating fluid specified by a Gaussian function for ${\rho }_{\circ }$ centered at $r=0.1$ with an initial linear velocity profile. The initial grid used for the coarsest solution shown is defined by the parameters $\{N_a,N_b,N_c,\Delta r_a\}=\{200,300,20,0.005\}$.

Figure 23

Convergence test for the fluid variable $\Pi$. The top panel shows the scaled error estimates described in the caption of Fig. 22. The bottom panel shows $\Pi (r,0)$ (dashed line) and $\Pi (r,t)$ (solid line), where $t$ is the time at which convergence is tested.

Figure 24

Convergence test for the fluid variable $\Phi$. The top panel shows the scaled error estimates described in the caption of Fig. 22. The bottom panel shows $\Phi (r,0)$ (dashed line) and $\Phi (r,t)$ (solid line), where $t$ is the time at which convergence is tested.