A Finite Difference Equations and Approximate Riemann Solvers

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 ”in-going” 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 ()—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.


I Introduction

Critical phenomena in general relativity involves the study of solutions—called critical solutions—that lie at the boundary between black hole-forming and black hole-lacking spacetimes. (See (1); (2); (3) for reviews.) Published work in general relativistic critical phenomena began over a decade ago with a detailed numerical examination of the collapse dynamics of a massless scalar field, minimally coupled to the general relativistic gravitational field (4). This first study in critical phenomena touched upon the three fundamental aspects of black-hole-threshold critical behavior: 1) universality and 2) scale invariance of the critical solution with 3) power-law behavior in its vicinity. All three of these features have now been seen in a multitude of matter models, such as perfect fluids (6); (7); Brady et al. (2002), an Yang-Mills model (8); (9), and collisionless matter (11); (10) to name a few. It was eventually found that there are two related yet distinct types of critical phenomena: Type I and Type II, so named because of the similarities between critical phenomena in general relativity and those of statistical mechanics.

Type II behavior was the first to be discovered (4), and entails critical solutions that are either continuously self-similar (CSS) or discretely self-similar (DSS). Super-critical solutions—those that form black holes—give rise to black holes with masses that scale as a power-law,


implying that arbitrarily small black holes can be formed. Here, parameterizes a -parameter family of initial data with which one can tune toward the critical solution, located at , and is the scaling exponent of the critical behavior. Since as , this type of critical behavior was named “Type II” since it parallels Type II (continuous) phase transitions of statistical mechanics.

As in the statistical mechanical case, there is also a Type I behavior in gravitational collapse, where the black hole mass “turns on” at a finite value. As might be expected, Type I critical solutions are quite different from their Type II counterparts, tending to be meta-stable star-like solutions that are static or periodic. In this paper, attention is restricted to Type II behavior; results from our study of Type I transitions in our model are reported in a separate paper (12).

The accepted picture describing the scaling behavior seen in Type II critical collapse was suggested by Evans and Coleman (6), who computed the critical solution for a radiation fluid (fluid pressure, , and density, , related by ) in two distinct ways. First, using a code that solved the full set of partial differential equations (PDEs) for the fluid and gravitational field, and by tuning an initial data parameter as sketched above, Evans and Coleman were able to compute a strong field solution that sat at the threshold of black hole formation, as well as establish a scaling law of the form (1) . Furthermore, the results of this numerical experiment provided compelling evidence that the threshold solution was continuously self-similar. Second, by adopting the assumption of continuous self-similarity as an ansatz, Evans and Coleman reduced the set of PDEs governing their model to a set of ODEs, from which a precisely CSS solution was calculated. The solutions computed using these two completely different techniques were found to agree extremely well. Crucially, it was argued that the observed scaling behavior of the mass above threshold could be explained using linear perturbation theory about a background given by the CSS critical solution. Such an analysis was carried out by Koike et al. (13) (for the radiation fluid), who showed that the scaling exponent, , was the inverse of the Lyapunov exponent of the critical solution’s single, unstable eigenmode.

Subsequent work showed that was not a truly universal constant, but that its value could depend on the specifics of the matter model used. The first evidence for this non-universality in scaling behavior was given in concurrent works by Maison (14) and Hara et al. (15). Using similar methods to those of (13); (16), they found that for an “ultrarelativistic” fluid with equation of state (EOS)


is dependent on the adiabatic index, .

Most of these investigations, however, have involved ultrarelativistic fluids that are explicitly scale-free. The reason for the predominance of this type of fluid is due to the fact that Cahill and Taub (17) showed that only those perfect fluids which have state equations of the form of Eq. (2)—i.e. the so-called ultrarelativistic EOS—can give rise to spacetimes that admit a homothetic symmetry (i.e. which are self-similar). Hence, it is not completely unreasonable to expect that Type II, CSS critical solutions would only appear in such fluids, or at least in fluids that admit an ultrarelativistic limit. To study this conjecture, Neilsen and one of the current authors (7) considered the evolution of a typical perfect fluid with equation of state,


that introduces a length scale into the field equations. Here, is the pressure, is the rest-mass energy density, and is the specific internal energy density. It was argued in (7) that Type II critical collapse scenarios are typically kinetic energy dominated, entailing increasingly large central pressures that maintain the tenuous balance between the matter dispersing from the origin and collapsing to a black hole. Therefore, if was reasoned that should one be able to give the fluid sufficient kinetic energy, then it would naturally enter an ultrarelativistic phase. Specifically, if the fluid undergoes a collapse such that dynamically, then will effectively become negligible in the equations of motion (EOM) and the system will be able to follow a scale-free—hence self-similar—evolution. To see if this hypothesis was correct, compact distributions of perfect fluid, with () were collapsed, and the calculations were tuned to a threshold solution. The critical solution thus obtained by solving the full set of PDEs closely matched the precisely self-similar solution, which was calculated by assuming that a model governed by the ultrarelativistic EOM had an exact homothetic symmetry. Further, it was found that the scaling exponent, , defined by Eq. (1) matched that of the ultrarelativistic critical solution with . Since the ultrarelativistic fluid exhibited Type II phenomena for all considered values of the adiabatic index in the range , the results of (7) suggested that the Type II ideal-gas critical solution for any in that range should be the same as that for an ultrarelativistic fluid with the same .

This hypothesis is not without precedence, since several models have been found to exhibit DSS or CSS collapse, even when explicit length scales are present. For instance, one of us found Type II behavior for the case of a collapsing massive scalar field (18)—that is a scalar field with potential —even though the model has an explicit length scale set by . The heuristic argument presented in (18) is that the potential term is naturally bounded since itself is bounded in the critical regime, but that the kinetic term——diverges in the critical limit. Hence, the kinetic term overwhelms the potential term and essentially makes the critical evolution scale-free.

The single study exhibiting Type II behavior in perfect fluid collapse with an ideal gas EOS (7) remained unverified until work by Novak (19). To determine the possible range in masses of nascent black holes formed from stellar collapse, Novak performed a parameter space survey using Tolman-Oppenheimer-Volkoff (TOV) solutions with , varying the overall amplitude, , of an otherwise fixed initial coordinate velocity profile for the fluid in order to generate critical solutions. The Type II behavior observed was quantified by fitting to the typical black hole mass scaling relation (1), where was identified with . Significantly, Novak was able to observe scaling behavior even with a realistic equation of state formulated by Pons et al. (20). This was somewhat surprising, since there were some expectations that Type II phenomena would not be observed when using realistic equations of state (3). However, provided that the equation of state admits an ultrarelativistic limit, as is apparently the case for Novak’s calculations, the heuristic argument sketched above suggests that one should expect to see Type II transitions.

Although Novak observed Type II behavior, he did not find the same scaling exponent as had been observed for the ultrarelativistic fluid in the study described in (7). In addition, he claimed that was a function of 1) the central rest-mass density, , which (as described in the next section) parameterizes the initial star solution, and 2) the EOS used. He observed that the fit to Eq. (1) worsened as increased to that of the maximum mass solution, and that it eventually broke down completely. Specifically, he found for the ideal-gas EOS (3)


and when using the realistic EOS


These values are significantly different from the values most recently calculated with the ultrarelativistic fluid Brady et al. (2002) using a variety of methods:


Here we have taken the average of the three independent values calculated in Brady et al. (2002) and the quoted uncertainty is the standard deviation of those values, and does not include the systematic errors inherent in the distinct calculations.

However, Novak clearly states in (19) that his code was not designed to simulate the formation of very small black holes, and apparently was only able to tune to a precision of . In this paper, we reexamine the Type II behavior in this particular system in order to check the claims of (19), and obtain what we claim is an improved measurement of the scaling exponent. Different families of initial data are used to demonstrate the universality of our computed critical solution. Also, we compare the critical configuration calculated from a near-threshold neutron star collapse to the critical solution obtained using an explicitly ultrarelativistic fluid. To more accurately study CSS behavior as the black hole threshold is approached, we employ mesh refinement techniques and non-uniform discretization. Further, we implement methods that improve the accuracy of the transformation from so called conservative variables to primitive variables that is required in our numerical analysis.

The remainder of this paper is structured as follows. Section II introduces the model and equations used to describe our collapsing neutron stars. In Sec. III, we describe the numerical methods used to solve the coupled fluid and gravitational PDEs, and include a discussion of an instability we observe near the critical threshold. In Sec. IV, we analyze the observed Type II behavior and compare it to previously published work. We then conclude in Sec. V with some closing remarks and notes on anticipated future work.

Geometrized units, , are used throughout, and our tensor notation follows (21).

Ii Theoretical Model

As in many previous critical phenomena studies in spherical symmetry Brady et al. (2002); (4); (8); (7); (19), we employ the so-called polar-areal metric


where is often referred to as the lapse function. We use the perfect fluid approximation for the matter model of our neutron stars, so the stress-energy tensor takes the form


Here, is the 4-velocity of a given fluid element, is the isotropic pressure, is the energy density, is the rest-mass energy density, and is the specific internal energy.

Modern computational methods that cast the hyperbolic fluid equations of motion in conservation form, and that use information concerning the characteristics of the equations, have been used very successfully in the modeling of highly-relativistic flows with strong gravity (see Banyuls et al. (1997); (23); (24); (25); (26) for a small, but representative, selection of papers on this topic). Here we adopt the formulation used by Romero et al. (26), with a change of variables similar to that performed in (25). The formulation described in (26) has been used extensively for treating relativistic flows in the presence of strong gravitation fields, and appears to work quite well in such instances.

The EOM for the fluid are derived from the local conservation equations for energy and baryon number which are, respectively,


Rather than working directly with the components of the fluid 4-velocity, it is more useful to employ the radial component of the Eulerian—or physical—velocity of the fluid as measured by an Eulerian observer:


where . The associated “Lorentz gamma function” is defined by


and satisfies the usual relation


since the 4-velocity is time-like and unit-normalized, i.e. . Adopting a notation where bold face symbols denote state vectors, the fluid EOM can be cast in conservative form as


where is a state vector of conserved variables, and and are, respectively, the flux and source state vectors. Our choice for the conserved variables is the one used almost exclusively in the field Brady et al. (2002); (7); (25); (19); (26), namely with


and where is the specific enthalpy of the fluid. , , , and can be thought of as the rest-mass density, momentum density, total energy density, and internal energy density, respectively, as measured in a Eulerian-frame defined by the ADM slicing. We found that for extremely relativistic flows near the threshold of black hole formation, this formulation was not very stable. We therefore use a different formulation motivated by (25), where it was found that evolving allowed for a more precise calculation since in the ultrarelativistic regime. We thus define new variables


and the state vectors become


The elements of the vector are the set of primitive variables used. The source function takes the form




With this source, the governing equations with which we solve for our metric functions are the Hamiltonian constraint of the ADM Arnowitt et al. (1962) formulation


and the polar-areal slicing condition


In order to monitor how near the spacetime is to forming an apparent horizon, we employ the mass aspect function, ,


From Birkhoff’s theorem, we know that an apparent horizon forms in the limit . We note that Eqs. (20-21) were used to compute the fluid equation source terms (18-19), and that flat space equations are obtained by setting and .

With the conservation equations (9)-(10), the equation of state (EOS) closes the system of hydrodynamic equations. Largely due to the extensive nature of our parameter space survey, we restrict the current study to continuum state equations (i.e. we do not use tabulated equations of state). The polytropic EOS


is used only when calculating initial conditions for a star. Here, is taken to be constant (isentropic condition) while is the adiabatic index. After , we allow for the development of shocks and therefore only use the “ideal-gas” or “gamma-law” equation of state (3). Our initial neutron star models are approximated by solutions of the spherically-symmetric hydrostatic Einstein equations, called the Tolman-Oppenheimer-Volkoff (TOV) solutions (28); (30); (29). We simulate the stiffness of matter at super-nuclear densities by setting in all of the calculations discussed here. As pointed out by Cook et al. (31), the constant can be thought of as the fundamental length scale of the system, which one can use to scale any dynamical quantity with set values of to a system with different values . As with and , we set , thereby making our equations dimensionless. We note that for a specific value of , the TOV solutions generically constitute a one-parameter family, where the central value of the fluid density, , serves as a convenient parameter. The ADM masses of stars for a TOV family governed by the ideal-gas EOS typically depend on the value of , with for , and achieving a global maximum at some value . Stars with are stable against radial perturbations, while those with are dynamically unstable in radial perturbation theory. In the experiments described below, we generally work with stars that have central densities significantly less than .

After the initial, star-like solution is calculated, an in-going velocity profile is added to drive the star to collapse. In order to do this, we follow the prescription used in (32) and (19). The method entails specifying the coordinate velocity


of the star. In general, the profile takes the algebraic form:


The two profiles that were used in (19) are


where and is the radius of the TOV solution. Unless stated otherwise, the profile will be used for all the results herein. The velocity profile is added consistently to the TOV solution by recalculating and via Eqs. (20-21) once the profile has been assigned to a given star. Further details concerning the calculation of initial data can be found in (33); (12).

As mentioned above, previous critical phenomena studies of perfect fluids have focused on models governed by the so-called “ultrarelativistic” EOS (2). This can be thought of as an ultrarelativistic limit of (3), wherein or . In this limit, becomes insignificant and one is left with two equations of motion for the fluid, which can be easily derived from Eqs. (14)–(19) by ignoring the EOM for and setting . The full expressions are given in (33). In this paper, we only use the ultrarelativistic EOS in order to dynamically calculate ultrarelativistic Type II critical solutions. All other computations are performed using the ideal-gas EOS (3).

Iii Numerical Techniques & Computational Issues

In this section we discuss the numerical techniques we use to simulate the highly-relativistic flows encountered in the driven collapse of neutron stars. The simulations entail solution of a system of coupled, partial and ordinary differential equations that describe how the fluid and gravitational field evolve in time. In section III.2 we also discuss an instability that generically appears in our calculations that are very close to the black hole threshold—this instability ultimately limits how closely we can tune any given family parameter to criticality.

We use the Rapid Numerical Prototyping Language (RNPL) (34) to handle check-pointing, input/output, and memory management for all our simulations. Secondary routines are called to solve the fluid and geometric equations. We use second-order high-resolution shock-capturing (HRSC) methods to evolve the fluid. The discrete equations are derived using a finite volume approach and are detailed in Appendix A. We generally use a Roe-type method (35); (26) as our approximate Riemann solver. However, particularly in the investigation of the instability mentioned above, we have sometimes used the Marquina flux formula (36) and Harten and Hyman’s (37) entropy-fix for Roe’s method, to compare with the basic Roe solver. For accurate and stable resolution of shocks, we use the minmod slope-limiter (38) to reconstruct the primitive variables at cell interfaces. We have also implemented the the linear MC (39) and Superbee (40) limiters, as well as the high-order essentially non-oscillatory (ENO) scheme (41), but have found the minmod limiter to provide the most stable evolutions near the threshold of black hole formation while still resolving shocks adequately.

In order to track the continuously decreasing spatio-temporal scales typically seen in CSS phenomena, we use a nonuniform grid that refines as needed. The origin, the point upon which the matter collapses, is the natural setting for the smallest dynamical scales and we therefore only refine the innermost region. Our implementation was inspired by Neilsen (42) and is detailed in (33). The basic idea is to segment the discrete domain into three regions: an innermost uniform grid with the smallest grid spacing, , composed of cells, an adjacent intermediate grid with cells of sizes , and an outermost uniform grid of cells. Refinement occurs when the maximum of is attained within cell widths from the origin. Interpolation is performed with a -order ENO interpolation procedure written by Olabarrieta (43). We usually set , , and adopt an initial value for such that the outer boundary lies at about times the initial radius of the star.

Time integration is performed separately from the spatial discretization using the method of lines. Specifically, an explicit, two-step predictor-corrector technique, called Huen’s method, is used to time-advance the ODEs that result from the spatial discretization of our time dependent PDEs. Discrete timesteps, , for the ODE integration are constrained to magnitudes given by , ensuring that the Courant-Friedrichs-Levy condition for our scheme is not violated. Additional details concerning the time integration are given in App. A as well as in (33). Results from a shock tube test, which measures our code’s ability to evolve discontinuities, and a convergence test involving the evolution of a self-gravitating distribution of fluid are presented in App. C.

iii.1 Primitive Variable Calculation

Since only the conserved variables are evolved by the HRSC schemes discussed above, the primitive variables must be derived from the conserved variables after each predictor and corrector step in order to compute fluxes and source functions for the next evolution step. This involves inverting the three equations —given by the definitions of the conserved variables (15)—for the three unknown primitive variables, . While closed-form expressions for the inverted equations exist for the ideal-gas EOS, numerically solving the equations is far more efficient (44). At each grid point, we use a Newton-Raphson method to find the values of that minimize the residuals of the conserved variable definitions (15). Instead of solving the full -by- system at each point, an identity function —derived from Eq. (15)—is used as a residual, making the solution process one-dimensional. This makes the procedure much more efficient.

Our method for performing the inversion, which is discussed further in App. D, is based on one specially suited for spherically symmetric ultrarelativistic flows (42); (25), and uses a residual function based on the definition of :


Here, is the enthalpy. In order to increase the accuracy of our computation of , we use different methods for calculating the residual and its derivative in different regimes (including both the ultrarelativistic and nonrelativistic limits). The “nonrelativistic” and “intermediate” methods originated from (42); (25), where flows in the ultrarelativistic limit were also studied. However, we have found that in the ultrarelativistic limit, where , the intermediate method still gives imprecise results. This imprecision can be traced to a loss of precision in the calculation of the quantity


We thus expand all nonlinear expressions appearing in the conversion to primitive variables in powers of , which yields results with increased floating-point accuracy. In the other limit, , where the the system is nonrelativistic, we use expansions up to that similarly reduce the influence of round-off errors. In practice, the ultrarelativistic regime is defined by an adjustable parameter and the nonrelativistic regime by . For example, for all the results shown below, we used and ; these values ensure that the leading-order error terms in the ultrarelativistic and nonrelativistic expansions are below the intrinsic round-off error of the computations.

Comparisons of the accuracy of our improved method and the work presented in (42); (25) are shown in Fig. 1. In order to estimate the relative error in the primitive variable calculation for each method, we seed each solver with guesses for that are a fixed factor away from the true solution, , where is the exact solution, and denotes the -th component of the state vector . A constant set of seeds are used in order to put all calculations on equal footing. Although the convergence of both methods does depend on the size of the , our pre-specified initial values are generally further from the true solution values than they are in the context of an actual calculation. In addition, use of different seeds with magnitudes comparable to those given above yielded similar results, suggesting that these particular values of are appropriately representative. For each method, once a best estimate for is computed, the relative error for the solver is computed as . Fig. 1 shows the logarithm of the relative errors for the variables and (the trend of the error in is similar to that of ). The improved accuracy of our new method is most noticeable in the computation of , as demonstrated by the first two columns in the figure. We note that our method can accurately calculate for and for all tested, while the method described in (25) develops significant problems when and .

Figure 1: Comparison of the accuracy achieved with our primitive variable solver and that used in (25). The first (third) column shows of the relative error between the exact value of () and the value obtained using the method described in (25), while the second (fourth) column shows of the relative error between the exact value of () and the value obtained from the method described in Sec. III.1. In any given plot, the vertical (horizontal) axis is discretized into 50 uniformly-spaced values of (); the minimum (maximum) value used for both and is () . Each row contains plots of -space calculations that were performed with a given value of shown to the left of the row. The uniformly-spaced color map is shown at the bottom and indicates that the maximum (darkest shading) represents errors comparable to and exceeding , while the minimum (lightest shading) represents errors comparable to and below machine precision.

Even though the above methods improved the accuracy of the primitive variable calculation, significant errors still remain for highly-relativistic flows () where and differ by orders of magnitude—i.e. when or . In these regimes, machine precision limits the accuracy of the calculation since terms in and become numerically insignificant (relative to other terms) even though their presence is essential to the computation of a solution. The effect from round-off error in these regimes can easily be seen by plotting using different orders of numerical precision, which we have done using arbitrary precision arithmetic in Maple. In order to accurately calculate in these regimes, one would need a new algorithm that performs better in these regimes, or use higher-precision arithmetic in the simulation. Fortunately, the improvements we have made seem to be sufficient for our current purposes.

iii.2 Instability at the Sonic Point in the CSS Regime

In this section we provide a description of an instability that develops in our calculations in the vicinity of the sonic point in near-critical evolutions. When using the unmodified approximate Roe solver, this instability made it impossible for us to obtain consistent brackets about the critical parameter, , for . This significantly hindered our study, since we found that we needed to tune quite closely to the threshold solution in order to calculate an accurate value of the scaling exponent .

An example of the instability seen in evolutions using primitive variable reconstruction is shown in Fig. 2. The conserved variable is plotted in CSS coordinates and X defined later in the paper (Eqs. (33) and (34), respectively). The last five frames show data from the last 5 time steps before the code crashes, while the first four frames are more distributed in . Hence, we see that the feature at the sonic point exists for many discrete time steps before its growth produces a code crash.

Figure 2: Conserved variable from the most nearly critical evolution obtained with the use of the approximate Roe solver without smoothing. and are defined by Eqs. (34) and (33), respectively. The dashed line indicates the location of the sonic point, . No refinement takes place during the period shown here, and . From left-to-right and top-to-bottom, the values of the frames are , , , , , , , , . The evolution started with a TOV star of central density that was perturbed using profile (see (26)) with the overall amplitude factor, , tuned to produce near-critical evolution.

The instability manifests itself in different ways, depending on the type of cell reconstruction used. For example when using the conserved variables to reconstruct the solution at the cell borders, we find that the conserved variables themselves remain smooth, but that each of the primitive variables exhibits persistent oscillations near the sonic point that span of order 2-4 grid cells. On the other hand, reconstructing with the primitive variables leads to smooth but oscillations in . A third and final reconstruction method was tried with the so-called characteristic variables, which are the advected quantities in the equations found by diagonalizing the quasi-linear form of Eq. (14). This method was significantly more diffusive but less stable than reconstruction with primitive variables.

The instability is also sensitive to the slope limiter used, with Superbee and Monotonized Central-differenced (MC) limiters producing more spurious oscillations than the more diffusive minmod limiter. We also found no improvements by varying the order of the ENO reconstruction from through .

In terms of ruling out potential sources of the problem, we have ensured that the regridding procedure is not responsible for the instability. To accomplish this, we first evolved a system that was tuned near the critical solution. We extracted the grid functions at a specific time, , before the appearance of instability, and interpolated them onto a new grid fine enough so that no further refinement would be required in the subsequent evolution. The data was allowed to evolve from this time, and the instability developed in the same manner and at the same time, , as in the original run.

Moreover, we find that the instability does not “converge away.” We tuned the initial data towards criticality for three different levels of refinement, where refinement was done locally so that for all , and is the “level” of refinement. We find that as increases the oscillations associated with the instability do not significantly change in magnitude and remain confined to approximately the same number of grid cells. Also, the solutions eventually diverge at the sonic point in all cases.

In order to describe the likely source of the instability, we first need to provide a better description of the near-critical solution. When the initial data has been tuned close to the critical solution at the threshold of black hole formation, the behavior near the origin is self-similar up to the sonic point, , where the flow velocity equals the speed of sound, (50). For these solutions the fluid becomes ultrarelativistic—e.g. —for and we expect that in that region. Also, from previous ultrarelativistic studies using such as Brady et al. (2002); (7), we expect that for . Thus, about the sonic point, the characteristic speeds (44) should take the values given in Table 1.

Characteristic Speed
Table 1: Asymptotic values of the fluid’s characteristic speeds in the ultrarelativistic limit. The sonic point is located at .
Figure 3: 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 , and are shown at . 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 as discussed in the text.

In fact, this is exactly what we find when using the ideal-gas state equation, as seen in Fig. 3 and Fig. 4. In Fig. 4 we see that within the self-similar region, but that for .

From these plots we also find that the transition from the ultrarelativistic regime to the exterior solution—defined by , and characterized by an absence of self-similarity and decay to asymptotic flatness—is quite abrupt. For instance, the discontinuity in is resolved by only a few grid points, signifying the presence of a shock which can also be seen for in the plots of and shown in Fig. 4. Since and , the discontinuity represents a point of transonic rarefaction. Also, the shock appears to be an expansion shock, which is entropy-violating, since it travels into a region of higher pressure and density. The reason why we can have an entropy-violating shock develop is because it is coincident with a change in curvature: as the high-pressure matter leaves the confines of the potential, it freely expands and its internal energy is converted into bulk kinetic energy.

Figure 4: 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 , and are shown at . Interior to the sonic point, , the fluid is clearly in the ultrarelativistic limit, with . However, beyond the sonic point, the flow is not ultrarelativistic—in fact in most of the domain exterior to . 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.

LeVeque states in (45) that the Roe solver can lead to the wrong Riemann solution at transonic rarefactions (in flat spacetime) since the linearization that the Roe solver performs on the EOM leads to a Riemann solution with only discontinuities and no rarefaction waves. He illustrates this point in (46) using a boosted shock tube test that makes the rarefaction transonic. Other failures of Roe’s method that are attributed to its linearization have been shown by Quirk (47), and by Donat et al. (48) where an unphysical “carbuncle” forms in front of a relativistic, supersonic jet.

To see whether Roe’s method contributes to the instability, we implemented the Marquina method and the entropy-fix for the Roe scheme due to Harten and Hyman (37). A comparison between the three methods is shown in Fig. 5, where we have evolved a shock tube problem that emulates the fluid state about the sonic point of near-critical solutions. The initial conditions used for this test are and . These values are such that, initially, , , , and , closely approximating the values listed in Table 1. The Roe, Marquina and Roe with entropy-fix evolutions each used 400 points in the entire grid of domain (only a portion of the grid is shown here), and both used a Courant factor of . The exact solution was obtained from the Riemann solver provided in (49) with points, using the same range in and same initial conditions as the Roe and Marquina runs. The Marquina method and entropy-fixed Roe method produce similar—yet more diffused—solutions than the exact solver, while the basic Roe solver severely diverges from the exact solution near the transonic rarefaction during the first few time steps. Even though the initial divergence of the Roe solution eventually vanishes, a relic feature still exists and propagates away from the center. If we were to reverse the evolution of the Roe solution shown here, the sequence would be reminiscent of how the instability in grows near the sonic point of near-critical solutions (Fig. 2).

Figure 5: One-dimensional, slab-symmetric shock tube test to simulate the discontinuity observed near the sonic point of near-critical solutions. The rest-mass density computed using different Riemann solvers is plotted versus in each constant- frame. Solution time is shown in the upper-right part of each frame. The solid line without points is an exact Riemann solution, the dashed line with triangles corresponds to the solution obtained with the approximate Roe solver, squares represent the solution from Marquina’s method, and circles are from the Roe scheme with entropy fix.

This shock tube test suggests that our use of Roe’s method may underlie the observed instability. In order to address this possibility, we implemented both the Marquina solver and the entropy-fix in the general relativistic code and tuned towards the critical solution. We were able to tune to with Marquina’s method, which is approximately a factor closer to than we reached with Roe’s method. Also, the “bump” at developed at a later with Marquina’s method. However, the use of Marquina’s flux formula did not completely solve the problem since evolutions using it also eventually succumb to the instability, preventing us from tuning beyond . Moreover, Roe’s method with the entropy-fix, performed more poorly than when used without the entropy-fix, impeding any further tuning past .

A possible explanation for the failure of the entropy-fixing methods to stabilize near-critical evolutions may involve the fact that they do provide an entropy-fix. As mentioned before, these methods would evolve the critical solution’s expansion shock to a rarefaction fan if no source was present (see Fig. 5). However, the source in the EOM and/or the geometric factor in the flux must be what keeps the expansion shock from being dissipated away since the shock is a part of the critical solution. The competition between the approximate Riemann solutions of the entropy-fixing methods and the spacetime may exacerbate the instability. However, it is unclear why Marquina’s method is more successful than Roe’s method with an entropy-fix.

Even though Marquina’s method provides a marked improvement over Roe’s method, it does not eliminate the sonic-point instability. A first attempt to explicitly dissipate the instability involves applying artificial viscosity in the region. We follow Wilson’s (50) artificial viscosity method and set in alone, where


and is a user-specified parameter. However, the instability worsens as , and, since is not proportional to like other terms in , becomes irrelevant as the flow becomes more relativistic. As a result, moderate values of lead to insignificant changes in behavior near the sonic point, while extremely large values tend to amplify the instability and induce additional high-frequency modes near .

Since we find that the instability near becomes more severe as became more discontinuous, our second attempt to control the blow-up entails smoothing the conserved variables about at every predictor/corrector step of the fluid update. The smoothing is done once the discontinuity develops. Specifically, we start to use the smoothing procedure when , and at times when begins to be resolved over approximately or fewer zones. We can use the same time, , to begin smoothing for all runs since the evolution for is almost identical for all near-critical values of . The smoothing is performed over the first contiguous set of points, , that satisfy , where is some adjustable parameter which we set to . The smoothing operation replaces the quantity with .

We also find that the instability worsens as the number of points between the origin and the sonic point decreases, as occurs in those cases where the solution disperses from the origin instead of forming a black hole. Our ability to follow evolutions through to dispersal is necessary for calculation of the scaling exponent, , since, as described in the next section, one way of estimating involves measuring the scaling of the global maximum, , of the stress energy trace, , as a function of . We find that is generally attained at a time when the fluid is beginning to disperse. Consequently, we found it necessary to refine the grid whenever the discontinuity or reaches .

The diffusion introduced by the smoothing allows us to further tune toward the critical solution, eventually to . However, we are still unable to calculate the global maximum of , , for the most nearly critical runs even though we can identify them as being dispersal cases. The minimum value of for which we can calculate is about , as illustrated in Fig. 8. This is far smaller, however, than we can achieve without smoothing or with any other method we have tried. Surprisingly, smoothing about the sonic point did not make the Marquina evolutions any more stable.

Iv Simulation Results

In this section, we study the Type II, CSS critical solution found at the black hole-forming threshold of the parameter space described in (12). If not otherwise stated, the results in this section use (27) for the initial velocity profile, and the overall amplitude, of the velocity profile for the tuning parameter, . As previously mentioned, the stars that we able to drive to a Type II black hole threshold generally have central densities, , significantly smaller than the maximum value, , along the stable branch, which for is in our units. Although we were generally able to form black holes from stars with an initial rest-mass central density greater than , we have closely tuned towards critical solutions for only a handful of such initial states. (We were unable to form black holes from stars with using an initially-ingoing velocity profile.) In Table 2, we list the central densities of the stars for which Type II behavior was actually observed, and quantify how close to the critical value we were able to tune. The instability described in Sec. III.2 limited the tuning in all instances.

Table 2: Star solutions in which we observed Type II behavior, and the minimum black hole masses we were able to form from them. We denote the mass of the smallest black hole found for a given by , is the mass of the initial star solution, and is the relative precision reached in per star. The final columns lists the critical parameter values we obtained.

From Table 2 it is clear that the instability’s effect on our ability to find the critical parameter increases with decreasing . This is most likely due to the fact that sparser stars require greater in-going velocities in order to collapse, giving rise to more relativistic and, consequently, less stable evolutions. We note, however, that our results represent great improvement over the precision obtained in (19); the smallest black hole attained in that study was . The success of our code is most likely due to our use of adaptive/variable mesh procedures and the great lengths we went to combat the sonic point instability.

Unless otherwise stated, and for the remainder of the section, we focus on behavior seen with the star having central density .

Figure 6: 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. 7. Black holes were assumed to have formed when . The gaps between some of the points represent those runs that crashed before reached this value. Smoothing was used for , which is also where we start our fit. These runs used , and an initial grid defined by .
Figure 7: 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 .

To demonstrate the scaling behavior of , we show in Fig. 6 versus for a wide range of supercritical solutions. The slope of the trend is equal to the scaling exponent, . From the figure, we can clearly see that the scaling law provides a good fit only in the limit as expected (13). The jump seen at represents the point at which the fluid enters a dynamical phase where the center part of the star has enough kinetic energy to dominate the effective potential energy, whose magnitude is set by . In this regime, the fluid then follows a CSS-type evolution.

In addition, Fig. 6 is meant to illustrate problems in our calculations associated with the coordinate singularity that inevitably develops in our Schwarzschild-like coordinate system in super-critical evolutions. Computations were run for a set of parameter values, , distributed uniformly in —any gaps in the plotted data thus represent instances where our code crashed prematurely. We note that the flow velocity becomes discontinuous and nearly luminal when black hole formation is imminent, and this seems to amplify the instability mentioned in Sec. III.2. This results in the evolution halting before exceeds its nominal threshold of 0.995. However, for a set of parameter values delineated in the figure by dashed lines, we were able to find a good fit to a scaling law. For that subset of data, the fit, as well as the data’s deviation from the fit, are shown in Fig. 7. One measure of how well the black hole masses are described by a relation of the form (1) is that deviations from the fit are small and apparently random. The slope of the trend yields an estimated scaling exponent of .

Figure 8: Scaling behavior in for subcritical solutions, i.e. those not forming black holes. The line shown here is the best-fit for the expected scaling law (31), using data only from the solutions closest to criticality. These runs used , and an initial grid defined by .
Figure 9: Three scale-free quantities of near-critical solutions in self-similar coordinates for the ideal-gas system (dashed line) and the ultrarelativistic system (solid line).

As mentioned in the previous section, to obtain another estimate of the scaling exponent, we calculate how the global maximum, , of the stress-energy trace scales as , using subcritical computations. As Garfinkle and Duncan (51) pointed out for the case of spherically-symmetric massless scalar collapse, the global maximum of the Ricci scalar should be proportional to the inverse square of the fundamental length scale of the self-similar solution. Hence for near critical solutions below the threshold should follow the scaling law:


Using instead of the Ricci scalar is computationally more expedient since it does not require the calculation of second-order space and time derivatives of the metric functions.

By determining from a plot of versus , in addition to a fit to Eq. (1), we can get an estimate of the systematic errors in our estimation of for both methods. Estimation of from the scaling of also has the advantage that, in the limit , the code is more stable for subcritical, rather than supercritical, evolutions. The scaling behavior for can be seen in Fig. 8 where is plotted versus . The solutions far from criticality seem to smoothly asymptote toward the critical regime. The line shown in this plot uses only those points in the regime that provide the best linear fit; a closer view of the points used in the fit are shown, for instance, in Fig. 12. Since the slope of the line now represents (31), we find from this fit that , which agrees with the value found from the scaling of to within the estimated systematic error in our computations.

Although our calculated scaling exponents match well to results previously obtained for the ultrarelativistic fluid with , this does not necessarily say how well the ideal-gas critical solutions compare to the ultrarelativistic ones in detail. To obtain the ultrarelativistic critical solutions, we let an adjustable distribution of ultrarelativistic fluid free-fall and implode at the origin; specifically, the initial data for the fluid is set so that is a Gaussian distribution and , and the amplitude of the Gaussian is used as the tuning parameter. The scale-free functions from the critical solutions of the velocity-induced neutron star system and the ultrarelativistic system are shown in Fig. 9. Here, the quantity is another scale-free function determined from metric and fluid quantities via


In order to make the comparison between the two solutions, the grid functions were transformed into self-similar coordinates and :


where is the elapsed central proper time,


is the accumulation time of a given critical solution, and is the location of the sonic point. Since we find that computing a smooth near the critical point is difficult, we typically use ,


as our self-similar radial coordinate. Here, is the position of the local maximum of that lies closest to .

Figure 10: Deviation over time of those quantities displayed in Fig. 9. Here, denotes the -norm of the function . The deviations for (solid line), (dotted line), and (dashed line) are shown. The -norms of these differences are computed at every time satisfying , and then logarithms of those norms are plotted as a function of self-similar time . Note that the sense of physical time is opposite to that of ; that is, as the solution approaches the accumulation time. As the evolution proceeds from the initial time, the two solutions asymptote toward each other. For , 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.

Our results indicate that the ideal-gas system does asymptote to the ultrarelativistic self-similar solution in the critical limit. While the ultrarelativistic fluid enters a self-similar phase shortly after the initial time, the ideal-gas solution generally tends toward the critical solution relatively slowly, then eventually diverges from it. The agreement between the ideal-gas and ultrarelativistic solutions improves as , as expected, and Fig. 9 shows profiles at a time when the difference between the solutions was minimized. The -norms 2 of the deviations between the ideal-gas and ultrarelativistic scale-free functions plotted over time are shown in Fig. 10; it can be easily gleaned from this figure that the minimum of the average deviations occurs at approximately , which is the time at which we have displayed the profiles in Fig. 9. Also, Fig. 10 graphically illustrates how the ideal-gas solution asymptotes exponentially—in central proper time, —to the ultrarelativistic critical solution at early times. The deviations for the three functions seem to have the same qualitative trend, indicating that metric and fluid quantities asymptote to their ultrarelativistic counterparts.

Figure 11: Time sequences of 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 . Approximate values of are shown in the upper-left corners of the frames. Note that the ultrarelativistic 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 in the latter case only becomes insignificant for .

This exponential approach of the ideal-gas solution to the self-similar solution is better seen in the comparison of time sequences of extracted from the ideal-gas and ultrarelativistic computations, as shown in Fig. 11. Here, has already attained a self-similar form at the beginning of the displayed sequence, while becomes self-similar at later times, and remains self-similar only for a time interval .

iv.1 Universality and Consistency

This section describes several numerical experiments primarily designed to ensure that the results presented above are not artifacts of the computational techniques used. These computations also provide a measure of the systematic error in our calculation of . Moreover, in order to check previous claims that critical solutions generated from perfect fluid configurations having the same adiabatic index may not reside in the same universality class, we also measure for different initial conditions, while keeping constant. When making any comparisons, the methods, parameters, and initial data used to produce Figs. 68 will be referred to as the “original” configuration. A tabulation of the values of and calculated from the different simulation configurations discussed in this section is given in Table 3.

Figure 12: Scaling behavior in near the critical solution for runs using different values of . The scaling behaviors are shown for the original configuration (circles,solid line), a configuration with times the original floor value (squares, dotted line), and one with times the original floor value (triangles, dashed line). The scaling exponents for these runs are listed in Table 3

The effect on the scaling behavior due to the floor value, , used for the fluid (see App. A) is estimated first. Since the magnitude of the floor is set in an ad hoc fashion—i.e. without any physical basis—it is crucial to verify that any results are independent of it. To test this, we replicated the original runs using different values of the floor while keeping all other parameters fixed. The scaling behavior obtained using different floor values is illustrated in Fig. 12. The dotted and dashed lines correspond to floor values that are factors of and , respectively, larger than the original configuration, which itself used . The minimal influence of the floor on solutions in the critical regime is clearly seen by the fact that all points follow nearly the same best-fit line. In fact, Table 3 indicates that all estimated values of agree to within and that all estimates of coincide to within . The deviations of the calculated sets from their respective best-fit lines for the different floor values even follow the same functional form, suggesting that the observed “periodic” deviations from linearity are not due to the floor.

The absence of any dependence on the floor is not too surprising since the component of the fluid that undergoes self-similar collapse is never rarefied enough to trigger use of the floor. For instance, at a time when the central part of the star begins to resemble an ultrarelativistic critical solution, the minimum values of within the sonic point are, respectively, —far above the typical floor values used. Only for is the floor activated, and dynamics in this region cannot affect the interior solution once self-similar collapse begins due to the characteristic structure of near-critical solutions (see Table 1).

Figure 13: Comparison of the scaling behavior in 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, , for these runs are listed in Table 3

The effect from the Riemann solver used on the scaling behavior is seen in Fig. 13. We find that the scaling behavior of from the two methods is remarkably close. Even though the Roe method with smoothing allows us to determine for smaller values of , the deviations from the best-fit of the two data sets are of the same order of magnitude for common values of . From Table 3, we see that the respective values of agree to within and that values of agree to within . These differences are quite small—comparable to those found as a result of varying the floor. Hence, we conclude that the choice in Riemann solvers has little, if any, effect on the computed scaling behavior, indicating that the smoothed approximate Roe solver is adequate for our purposes.

Figure 14: Scaling behavior in near the critical solution for runs using different levels of resolution. The runs were made with , , and the solid line with circles was generated from runs using the original configuration. The (squares, dotted line) and (triangles, dashed line) runs, respectively, used computational grids that were locally times and times as refined. The scaling exponents, , for these runs are listed in Table 3

When using finite difference methods, it is vital to verify that the order of the solution error is the same as the order to which the derivatives are approximated by difference operators. For example, our HRSC scheme should be accurate in smooth regions and near shocks, so we should expect this scaling behavior of the error as is varied. First, we wish to see if our estimate for converges as the grid is refined. Figure 14 shows a plot of versus for the original configuration, along with others computed at higher resolutions. We first see that the three distributions follow lines of approximately the same slope (and which are thus shifted vertically relative to one another by constant amounts) while the deviation of the best-fits seems to increase slightly with resolution. Also, we can see that an increase in resolution permits us to follow the collapse through to dispersal for solutions closer to the critical threshold, allowing for the scaling law to be sampled at smaller . Even though the deviations from the best-fits for are quite small compared to the typical size of , it is a little worrisome that they are larger than those from the lowest resolution runs. However, this behavior can likely be attributed to the sonic point instability and the smoothing procedure used to dampen it. In particular, the “bump” at sharpens with increasing resolution spanning a roughly constant number of grid cells (see Sec. III.2 for more details). Consequently, the impact of the instability on the solution may also increase with decreasing , since the discretized difference operators will generate increasingly large estimates for spatial derivatives in the vicinity of the sonic point. In addition, the smoothing operation is always performed using nearest-neighbors, so the smoothing radius physically shrinks with resolution, diminishing the impact of the smoothing.

Figure 15: Logarithm of scaled, independent residuals of the Hamiltonian constraint (20) and slicing condition (21) for three levels of resolutions calculated from solutions in the self-similar regime. The dotted (dashed) lines are from a run which used () 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 () in order to make the convergence of the solution more apparent. Each distribution is from a solution that has been tuned to with respect to the value of for each resolution, and function values at every tenth grid point are shown. The physical velocity of the fluid for the run is shown in the bottom frame in order to facilitate comparison of features in the independent residuals to those in the solution.

In order to verify that the code is converging in the self-similar regime, we computed independent residuals (i.e. applied discretizations distinct from those used in the scheme used to compute the solution) of the Hamiltonian constraint (20) and slicing condition (21) for the three levels of resolution discussed previously (Fig. 15). The overlap of the scaled residuals seen in the figure indicates convergence. Note that the smoothing procedure has not been used to calculate the solutions shown here. We see that the scaled residuals have similar magnitudes in all regions except those that have been processed by shocks, namely . Because the self-similar solutions are converging at the expected rate, we surmise that the variations observed in for the three resolutions does not indicate a problem with convergence, but demonstrates the effect of truncation error and/or the smoothing procedure on the scaling behavior. With only three levels of resolution, it is hard to make definite claims as to whether is or is not converging to a particular value. Even so, the standard deviation of determined from the three evolutions is about of their mean, suggesting that the variation is not significant. In fact, it is comparable to the standard deviation found from the simpler ultrarelativistic perfect fluid studies of (7).

Figure 16: Scaling behavior in for several families of initial data. The “Original” date (solid line) is as before, the dotted line with squares shows the scaling behavior for runs that used a different initial velocity profile, , and the dashed line with triangles was made from runs with a different TOV solution, with . The scaling exponents, , for these runs are listed in Table 3

The final comparison entails varying the physical initial conditions of the system to investigate the universality of the critical phenomena computed with the ideal-gas EOS. The primary constituents of our model are the initial star solution and the form of the perturbation with which we drive the star to collapse. Hence, we choose to perform sets of runs to measure the scaling law using 1) a different initial star solution and 2) a different functional form of the initial velocity profile. The scaling behaviors of versus for these different configurations are compared to the results from the original configuration in Fig. 16. For the data computed using a star of central density , the only different aspect is the initial star solution. This particular value of is chosen since it is near the transition from Type II to Type I phenomena discussed in detail in (12). The second initial data set uses (27) for the initial profile of the coordinate velocity. Naturally, we see that the three data sets are vertically shifted relative to one another since each set evolved from significantly different profiles of mass-energy—the details of the initial data set the scale for for specific values of . However, only the slopes of the curves are relevant for estimating .

From the values listed in Table 3, we see that varies more significantly with the particular star solution used than with the form of the velocity profile. In fact, we are able to tune closer to the critical solution with the more compact star, a trend that can also be seen in Table 2. Nonetheless, the scaling exponents computed using the two distinct initial star configurations agree to within .

The change in the initial velocity profile only affects the computed value of by . This suggests that other methods of perturbation would also yield close to the same value. The concordance of results from these three different families of initial data imply that universality of critical solutions is maintained for perfect fluids governed by an ideal-gas EOS, at least for the case . It would be interesting to see whether these properties hold with even more realistic equations of state, as well as for other values of .

iv.2 Final Determination of

Using the calculated values of from the various methods, floor sizes, and grid resolutions, we are able to provide an estimate of the systematic error inherent in our numerical model. Further, by assuming that the universality is strictly true, we can even use the variation in computed from the different initial data families for this estimate. Taking the average and calculating the standard deviation from all of the values for the ideal-gas EOS listed in Table 3, we estimate a scaling exponent value of


This is in agreement with the value of computed from the black hole mass scaling fit shown in Fig. 7.

Roe 0 0.94 0.4687537
Roe 0 0.94 0.4687535
Roe 0 0.95 0.4687516
Roe 1 0.92 0.4682903
Roe 2 0.93 0.4682461
Roe 0 0.94 0.4299032
Roe 0 0.92 0.4482047
Marquina 0 0.94 0.4687682
Ultra-rel. - - - 0.97 -
Table 3: Scaling exponents and critical parameters computed from fits to the expected scaling behavior in . Scaling exponents were obtained from runs using initial star solutions with different central densities () and using different floor magnitudes (), levels of refinement (), and velocity profiles (). The runs labeled “Roe” use the approximate Roe solver with smoothing, the “Marquina” run used the Marquina flux formula, and the “Ultra-rel.” scaling exponent was computed from our results involving the collapse of Gaussian profiles of ultrarelativistic fluid.

In addition, we can compare our final estimate of to values previously found for the ultrarelativistic fluid. As already mentioned, was measured at three different refinement levels in (7), and a value


was quoted.

Instead of solving the full set of PDEs, can also be found by solving the eigenvalue problem that results from performing first order perturbation theory about the CSS solution. This was done in two ways in Brady et al. (2002): using the common shooting method, and solving the eigenvalue problem directly after differencing the equations to second order. The scaling exponents calculated were, respectively, and .

iv.3 Possible Presence of a Kink Instability

As mentioned in Sec. III.2, we witness an instability near the sonic point of solutions tuned close to the threshold. After thorough numerical experimentation, we are still left uncertain about its genesis. One possibility is that it has physical origin. For instance, Harada (52) reported the presence of an unstable kink mode for spherically symmetric ultrarelativistic perfect fluids when . The unstable kink mode manifests itself as an ever-steepening discontinuity in at the sonic point, , that diverges in finite proper time. A mild, seed discontinuity at is necessary for the mode to grow, but—since it diverges in finite time—minute discontinuities inevitable in discretized solutions, for instance, are expected to be sufficient to manifest the instability. Physically, the “seed” discontinuity could be the scale at which the continuum approximation of hydrodynamics breaks down, i.e. at the particle scale.

In our nonlinear PDE solutions that use the ideal-gas EOS, we do in fact find a growing discontinuity in and at (Fig. 4), and this is precisely where our instability develops.

However, as was the case in (7), our solutions of the PDEs for an explicitly ultrarelativistic fluid do not develop instabilities in the vicinity of for . One possible reason for this could be the fact that in this case the flow is never transonic, i.e. there is no sonic point since , and the flow can never attain this velocity. For the ideal-gas fluid, there is a sonic point, as (albeit arbitrarily close to ) since will always be non-zero during a numerical evolution. However, the results of (7) include ultrarelativistic near-threshold solutions for . This case does allow for the presence of the sonic point, so it may be that sonic-point kinks were not seen in this instance because never became sufficiently steep at to excite the kink mode.

Figure 17: Scaling behavior in for different values of . The “Original” (solid line with circles) was made from runs with and as before. Both the (dashed line with triangles) and (dotted line with squares) runs used with . The scaling exponents are given in the text.

In Fig. 17, we plot fits of the Type II scaling behavior of stars with three different values of . The distribution is our fiducial system, while the other two——use a different progenitor TOV solution having instead of . The change in initial state was necessary since use of for and would not necessarily lead to Type II behavior using velocity induced collapse. In addition, the values and bound by a difference of the critical value, , above which the kink mode becomes unstable (52). The scaling exponents derived from the two new sets of computations are and . The scaling exponent is the same as that calculated with the ultrarelativistic PDEs in (7), while the scaling exponent is consistent with the value obtained for , the value of closest to 1.88 for which was computed in (7).

We find that there is no consistent behavior in code stability as the value is crossed. In fact, we find the opposite of the expected behavior: we are able to tune the kink-unstable () data closer to than the kink-stable () data. We therefore find it unlikely that the kink mode is the cause of our sonic point instability.

Even if the kink mode is unstable for our EOS, tuning toward the CSS solution while in the presence of another unstable mode is not without precedent. For example, in a study of the spherically-symmetric general relativistic harmonic-map (nonlinear sigma model), Liebling (53) discovered that by judicious choice of an initial data family, he could tune to a critical solution that had been shown to have two unstable modes. However, it is unclear what relation this work might have to our current study, since the type of initial data that we are studying does not seem to have been chosen in any particularly special way (i.e. we suspect that the initial data that we have used is generic, whereas that used by Liebling to tune to the two-mode-unstable solution was, by construction, non-generic). What is clear is that this issue requires further investigation, but we will leave that to future studies.

V Conclusion

In this work, we simulated spherically-symmetric relativistic perfect fluid flow in the strong-field regime of general relatively. Specifically, a perfect fluid that admits a length scale, for example one that follows a relativistic ideal-gas law, was used to investigate the dynamics of compact, stellar objects. These stars were modeled as neutron stars by using a stiff equation of state, approximating the behavior of some realistic state equations. Our models were then used to study the dynamics of neutron stars so far out of equilibrium that they are driven to gravitational collapse.

Since these systems entail highly-relativistic fluid motions and strong, nonlinear effects from the fluid-gravitational interaction, a numerical treatment is challenging. To achieve stable evolutions in near-luminal flows, while using high resolution shock capturing techniques, the primitive variable solver required improvements. In addition, an instability was found to develop in calculations near the threshold of black hole formation, and this necessitated the use of new computational methods, that were only partially successful in stabilizing the calculations.

We find our value for the scaling exponent, , given in Eq. (37) agrees well with those found in Brady et al. (2002), and agrees with the value of computed in  (7) to within the uncertainty quoted in that work. We note that a discrepancy in the values of computed using ideal-gas and ultrarelativistic fluids was observed in (7), and this is also the case for our calculations. Our ultrarelativistic value, , agrees well with the value calculated in (7), but deviates by an estimated 3 standard deviations from the value extracted from our ideal-gas calculations. It is somewhat interesting, yet probably coincidental, that our results from the ideal-gas system of equations lead to estimates of that agree with the perturbation calculations better than those values found from the ultrarelativistic PDE calculations.

Our findings thus do not support some of the results found, and claims made by Novak (19) for the case of fluid collapse with an ideal-gas EOS and . This previous work suggested that the Type II behavior observed in such a case was not well approximated by a universal (with respect to initial data) ultrarelativistic limit. However, using different stars and velocity profiles, and by varying other aspects of the numerical model, we have found scaling behavior that is insensitive to approximations made in the numerical solution, and which does appear to be universal with respect to families of initial data. Moreover, we have found that the scaling exponent and critical solution for the collapse governed by the ideal-gas EOS agrees well with their ultrarelativistic counterparts.

Ultimately, it is our goal to expand the model a great deal, making the matter description more realistic and eliminating symmetry. As a first step, we wish to develop adaptive mesh refinement procedures for conservative systems that will be required to study critical phenomena of stellar objects in axial-symmetry (54).

It remains to be seen whether the universal scaling behavior we have observed is also seen with more realistic state equations such as the one Novak used. Since accurate measurements of have only been found for equations of state with constant adiabatic index , and since seems to only depend on for perfect fluids, it will be interesting to investigate in detail what the scaling behavior—if any—will be like for realistic state equations with variable . However, to the extent that any given realistic EOS admits a unique ultrarelativistic limit, characterized by a single value of , we can expect to see universal Type II behavior of the sort discussed in this paper, for at least some collapse scenarios.

We wish to acknowledge financial support from CIFAR, NSERC, and NSF PHY 02-05155. SCN wishes to thank I. Olabarrieta for many helpful comments and his ENO code. All numerical calculations were performed on UBC’s vn PIII and vnp4 clusters (supported by CFI and BCKDF).

Appendix A Finite Difference Equations and Approximate Riemann Solvers

We employ High-Resolution Shock Capturing (HRSC) algorithms to solve the equations of motion for the fluid (14). Such methods have become increasingly popular in the field of relativistic hydrodynamics since they are flux conservative, and ensure that discontinuities are well resolved and propagate at the correct speeds in the continuum limit. A key ingredient to these schemes is their use of solvers for the Riemann problem at every cell interface. This is crucial for the conservative nature of these schemes since the solution to the Riemann problem is always a weak solution of the hyperbolic conservation laws. The “high-resolution” aspect of the algorithms denotes that in regions where the grid functions are smooth, the integration procedure is at least accurate. Many of the HRSC methods used in this paper have been used in previous works such as (26), (19), and (25) to name only a few relevant sources. Also, excellent references on conservative methods for general systems of hyperbolic conservation equations have been written by LeVeque (45); (46).

Unlike finite difference methods, finite volume or conservative methods calculate the cell-averages of grid functions instead of the grid functions themselves. The difference equations for conservative methods are not derived from Taylor-series approximations of derivatives, but from differences of integrals. For instance, we difference the fluid EOM (14) in the following manner:


At first glance, this equation seems no different than a finite difference approximation of Eq. (14). However, the difference between the two approaches becomes apparent when we examine how specific quantities in (39) are defined. is the spatial average of over the cell centered at , is the spatio-temporal average of the source function centered at , and the numerical flux is the time average of from to . In practice, we approximate as the source of the averages, .

The techniques for determining the numerical flux are especially important since they set the spatial accuracy of the overall scheme and are primarily responsible for how well shocks are resolved by the method. Unless otherwise stated, all results presented in this paper were produced using an approximate Roe-type solver as outlined in (26) for calculating the numerical flux. The Roe solver (35) approximately solves the Riemann problem at each cell interface by casting the conservation equation (14) into quasi-linear form. We approximate the Roe matrix as




and where and will be defined shortly. After solving the linear Riemann problem, the numerical flux of is computed using the following expression (45):


Here, and are the eigenvalues and right eigenvectors, respectively, of , and are, respectively, the values of to the left and right of the cell boundary, and are the decomposed values of the jumps in the space of characteristic values:


In order to calculate all quantities associated with , such as and , we use the average of the left and right states, , defined by (41).

As far as we know, the eigenvectors for our new formulation, (14) and (17), have not previously been published, so for completeness we present them here. Since the transformation from to is linear, the eigenvalues remain the same:


Using Maple and assuming the the ideal-gas EOS (3), we have calculated the left and right eigenvectors; the Marquina flux formula requires the right eigenvectors (36). Using the typical normalization for the eigenvectors (), leads to a very complicated set of eigenvectors. Hence, we used the normalizations


which leads to significant simplification. The right eigenvectors become:


while the associated left eigenvectors are