{'fraction_non_alphanumeric': 0.03865334916532793, 'fraction_numerical': 0.015338408884542852, 'mean_word_length': 4.392542372881356, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Axisymmetric general relativistic hydrodynamics: Long-term evolution of neutron stars and stellar collapse to neutron stars and black holes

We report a new implementation for axisymmetric simulation in full general relativity. In this implementation, the Einstein equations are solved using the Nakamura-Shibata formulation with the so-called cartoon method to impose an axisymmetric boundary condition, and the general relativistic hydrodynamic equations are solved using a high-resolution shock-capturing scheme based on an approximate Riemann solver. As tests, we performed the following simulations: (i) long-term evolution of non-rotating and rapidly rotating neutron stars, (ii) long-term evolution of neutron stars of a high-amplitude damping oscillation accompanied with shock formation, (iii) collapse of unstable neutron stars to black holes, and (iv) stellar collapses to neutron stars. The tests (i)--(iii) were carried out with the $\Gamma$-law equation of state, and the test (iv) with a more realistic parametric equation of state for high-density matter. We found that this new implementation works very well: It is possible to perform the simulations for stable neutron stars for more than 10 dynamical time scales, to capture strong shocks formed at stellar core collapses, and to accurately compute the mass of black holes formed after the collapse and subsequent accretion. In conclusion, this implementation is robust enough to apply to astrophysical problems such as stellar core collapse of massive stars to a neutron star and black hole, phase transition of a neutron star to a high-density star, and accretion-induced collapse of a neutron star to a black hole. The result for the first simulation of stellar core collapse to a neutron star started from a realistic initial condition is also presented.


I. INTRODUCTION
In the 1980s, one of the most important issues in the field of numerical relativity involved performing simulations of rotating stellar collapse with the assumption of axial symmetry. Simulations of rotating stellar collapse in full general relativity were first performed by Nakamura and collaborators [1,2]. Using the (2+1)+1 formalism developed by Maeda et al. [3], they succeeded in performing simulations of a rotating collapse of massive stars to black holes. They used cylindrical coordinates (̟, z) with a nonuniform grid spacing and with at most (42,42) grid resolution for (̟, z) because of restricted computational resources.
To compute gravitational waves emitted during gravitational collapse to black holes, Stark and Piran [4] subsequently performed simulations similar to those of Nakamura et al., adopting spherical polar coordinates with a typical grid size (100, 16) for (r, θ). The distinguishing feature of their work is that they adopted the Bardeen-Piran formalism [5], which is well suited for computation of gravitational waves in the wave zone. As a result of this choice of formalism, they succeeded in computing gravitational wave forms, and clarified that the wave forms are characterized by the quasinormal mode of rotating black holes formed after gravitational collapse and that the total radiated energy of gravitational waves is at most 0.1 % of the gravitational mass of the system [4].
Since the completion of their works, no new work in this field was done for the next 15 years [6]. Although several questions that they originally wished to answer have been answered by their simulations, it was not feasible to perform sophisticated astrophysical simulations in the 1980s. This is likely due to the fact that the computational resources were severely restricted, and, in addition, that techniques in numerical relativity such as methods to provide realistic initial conditions and to perform the long-term simulations were not sufficiently developed. As a result, there still remain many unsolved issues in astrophysics and general relativity that can be studied with axisymmetric hydrodynamic simulations in full general relativity. Among them, realistic simulations of rotating core collapse of massive stars, which thereby become black holes or protoneutron stars in full general relativity, have not yet been performed. Stellar collapse is a common phenomenon in the universe, and hence, understanding the formation mechanism of black holes and neutron stars in nature is one of the most important issues in astrophysics. Actually, the study of the formation of rapidly rotating black holes with surrounding accretion disks in stellar core collapse is currently one of the hot topics in connection with a hypothetical scenario for the central engine of γ-ray bursts [7]. To date, simulations of a rotating collapse of a massive stellar core have been done in the Newtonian gravity [8][9][10][11][12][13] or in an approximate general relativistic gravity [14] using the so-called conformal flatness approximation (or the Isenberg-Wilson-Mathews approximation). In rotating stellar core collapses, general relativity plays an important role. As demonstrated in [14], general relativistic effects modify the collapse, bounce, and amplitude of gravitational waves emitted significantly, even in the formation of neutron stars. Of course, general relativity plays a crucial role in the formation of black holes. Thus, general relativistic simulation is inevitable to precisely understand the nature of stellar core collapses.
One long-standing issue for axisymmetric simulations in full general relativity has been to develop methods in which the accuracy and stability for a long-term simulation can be preserved. In axisymmetric simulations, we have in general used cylindrical and/or spherical polar coordinate systems, which have coordinate singularities at the origin and along the symmetric axis ̟ = 0. At such coordinate singularities, the finite differencing scheme has to be changed, resulting often in numerical instabilities. To stabilize computation, we have often been required to add artificial viscosities around the coordinate singularities to stabilize the numerical system [15].
Recently, the Potsdam numerical relativity group has proposed the so-called cartoon method by which a robust numerical relativity implementation for axisymmetric systems can be made [16]. The essence of their idea is that the Cartesian coordinates (x, y, z) could be used even for simulations of axisymmetric systems if the Einstein field equations are solved only for the y = 0 (or x = 0) plane, using the boundary condition at y = ±∆y (or x = ±∆x) provided by the axial symmetry. (Here, ∆x, ∆y, and ∆z denote the grid spacing.) Since the field equations are written in the Cartesian coordinate system, we neither have singular terms nor do we have to change the finite differencing scheme anywhere, except at the outer boundaries. Thus, it is possible to perform a stable and accurate long-term simulation without any prescription or artificial viscosities, but only by a minor modification of a three-dimensional implementation that has already been developed [17][18][19][20].
Other important progress has been made regarding computational resources. Current large-scale supercomputers that we can use are typically of several hundred Gbytes memory. Necessary memory in an axisymmetric simulation with double precision, with N 2 grid points, and with N v variables is where N v is ∼ 200 in our general relativistic implementation. This implies that the memory of current supercomputers is large enough to carry out an axisymmetric numerical simulation with N ∼ several thousands. Using N of order 10 3 , it is feasible to carry out a well-resolved simulation and a careful convergence test, changing the grid resolution for a wide range from N ∼ 100 to 1000. This situation is in contrast with that of three-dimensional numerical relativity, since it is still very difficult to carry out a three-dimensional simulation with N ∼ 10 3 , for which the required computational memory is of order TByte. Motivated by the status mentioned above, we recently started a project in axisymmetric numerical relativity. In [21], we reported a numerical hydrodynamic implementation for axisymmetric spacetimes that is made using the cartoon method and incorporating a hydrodynamic implementation in the cylindrical coordinates. In that paper, we presented numerical results for simulations of rotating stellar collapse adopting simple initial conditions and simple equations of state to investigate the effects of rotation on the criteria for prompt collapse to black holes in an idealized setting. We demonstrated that the axisymmetric implementation and current computational resources are well suited to systematically perform stable and well-resolved hydrodynamic simulations in axisymmetric numerical relativity. That implementation has also been applied to a study of collapse of rotating supermassive stars to supermassive black holes [22].
We have recently remade a hydrodynamic implementation using a high-resolution shock-capturing scheme based on a Godunov-type scheme [23][24][25][26][27][28]. Although the previous one [21] works well for problems in which shocks are weak, such as evolution of single rotating stars and collapse of neutron stars and supermassive stars to a black hole, it is expected that such implementation cannot produce an accurate numerical result for problems in which shocks are strong. During rotating core collapses to a neutron star or a black hole, strong shocks are likely to be accompanied. Therefore, implementing a high-resolution shock-capturing scheme, such as that adopted in [28], is a promising strategy. To check that the new implementation works well, we have performed a wide variety of test simulations. In this paper, we present the numerical results, paying particular attention to long-term numerical simulations of neutron stars as done in, e.g., [18,28], and to stellar collapse in which strong shocks are accompanied. Finally, we present the first numerical results of stellar core collapse of a neutron star for which the simulation is started from a realistic initial condition. In addition to presenting the successful numerical results, we address the advantage of axisymmetric simulations in testing a new general relativistic hydrodynamic implementation, since we can study in detail the convergence of numerical results in the test simulations changing the grid number for a wide range to a well-resolved level (e.g., N ∼ several hundreds), which is still difficult in three-dimensional simulations because of restricted computational resources.
The paper is organized as follows. In Sec. II, we describe the formulation that we adopt. In Sec. III, we define global quantities of the system and describe the calibration method for the numerical results. In Sec. IV, we present the numerical results. Section V is devoted to a summary and discussion. Throughout this paper, we use the geometrical units in which G = 1 = c, where G and c denote the gravitational constant and speed of light. We use Cartesian coordinates, x k = (x, y, z), as the spatial coordinates, with r = x 2 + y 2 + z 2 , ̟ = x 2 + y 2 , and ϕ = tan −1 (y/x). t denotes the coordinate time. Greek indices µ, ν, · · · denote x, y, z, and t, small Latin indices i, j, · · · denote x, y, and z, and capital Latin indices A, B, · · · denote x and z.

A. Formulation for the Einstein equation
The Einstein equation is solved in the (3+1) formulation, in which the line element is written in the form where α, β k , and γ ij are the lapse function, shift vector, and three-metric. The three-metric is defined by where n µ is a timelike, unit-normal vector that is orthogonal to a spacelike hypersurface, and its components are written as (1/α, −β k /α). The extrinsic curvature is defined as where L − n is the Lie derivative with respect to n µ , and D i is the covariant derivative with respect to γ jk . The Einstein field equations are solved using the same formulation, gauge conditions, and outer boundary conditions as in previous papers [17][18][19][20]29,30,21]: We adopt the so-called Nakamura-Shibata formulation [2,31] with some modification from the original version (see [20], to which the reader may refer for basic equations and gauge conditions in the latest version). In this formalism, we evolve the following geometric variables using a free evolution code: The Hamiltonian and momentum constraint equations are solved at t = 0, and used to check the accuracy of numerical solutions during computation. The slicing and spatial gauge conditions for determining α and β k are basically the same as those adopted in our previous series of papers [17][18][19][20]29,30,21], i.e., we impose an "approximate" maximal slicing condition (K ≈ 0) and an "approximate" minimum distortion (AMD) gauge condition [D i (∂ tγ ij ) ≈ 0, whereD i is the covariant derivative with respect toγ ij ]. However, in contrast with previous papers [17,29,19,20], we do not modify the spatial gauge condition even in the formation of black holes, since in the axisymmetric simulation, a sufficient number of grid points can be taken to resolve black hole formation and subsequent evolution even using the AMD gauge condition without modification.
We impose the axially symmetric condition to the geometric variables using the so-called cartoon method proposed by Alcubierre et al. [16]. First, we define the computational domain as 0 ≤ x, z ≤ L and −∆y ≤ y ≤ ∆y, where L denotes the location of the outer boundaries, and reflection symmetry with respect to the z = 0 plane is assumed. With this computational domain, we need only three points in the y direction, 0 and ±∆y. We determine here that the Einstein equation is solved only in the y = 0 plane. Then, the boundary conditions at y = ±∆y that are necessary in evaluating y derivatives are supplied from the assumption of axial symmetry as where and ϕ(x) = tan −1 [±∆y/ x 2 + (∆y) 2 ]. Q ij , Q i , and Q denote (γ ij ,Ã ij ), (F i , β i ), and (φ, K, α), respectively, and Q (0) ij , Q (0) i and Q (0) are the values of Q ij , Q i , and Q at ( x 2 + (∆y) 2 , 0, z), which are interpolated using Lagrange's formula [32] with three nearby grid points along the x direction (i.e., x ± ∆x and x). At x = L, we use only two points, x − ∆x and x, for the extrapolation.
To impose the gauge conditions, as well as to solve the constraint equations in preparing the initial conditions, we solve scalar and vector elliptic-type equations of the form [18] where ∆ flat denotes the Laplacian in the flat three-dimensional space, and S and S i denote the source terms. Using the interpolation mentioned above, ∂ yy Q and ∂ yy Q i are evaluated in the finite differencing as On the other hand, the finite differencing in the x and z directions, ∂ xx Q i and ∂ zz Q i , is written in the standard form as Thus, in the finite differencing form for each component of Eq. (13), only one component of Q i is included, implying that each component of the vector elliptic-type equation is solved independently, as in the case of the scalar elliptic equation.
Finally, we note a necessary modification for numerically handling Einstein's evolution equations in the axisymmetric case. In our formalism, the evolution equations are written in the form [31,17,20] where Q denotes one of the geometric variablesγ ij ,Ã ij , φ, K, and F k . In the three-dimensional case, we apply an upwind scheme to numerically handle the transport term β k ∂ k Q for all the components [17]. In the axisymmetric case, the same method is used for β x ∂ x Q and β z ∂ z Q, but is not for β y ∂ y Q, since it is not appropriate. (Remember that in the hydrodynamic equations in the axisymmetric case, there is no transport term for the rotational direction.) For β y ∂ y Q, we use the following schemes: For Q = φ, K,γ zz , andÃ zz , we set β y ∂ y Q = 0 because of symmetry. For other variables, we simply use the cell-centered (second-order) finite difference.

B. Formulation for the hydrodynamic equations in general relativity
The hydrodynamic equations in general relativity are written as where ∇ µ is the covariant derivative with respect to the spacetime metric g µν , ρ is the baryon rest-mass density, u µ is the four-velocity, and T µν = ρhu µ u ν + P g µν .
Here, P is the pressure, h ≡ 1 + ε + P/ρ is the enthalpy, and ε is the specific internal energy. Equations (15) and (16) are the continuity equation and the equations of motion, respectively. We adopt the so-called high-resolution shock-capturing scheme in numerically handling the transport terms of hydrodynamic equations. To use such a scheme, the hydrodynamic equations should be of a conservative form as where e ≡ e 6φ ρ * T µν n µ n ν = hw − P ρw , and η is a determinant in curvilinear coordinates; in the cylindrical coordinates, η = ̟. We note that subscripts i, j, · · · here denote the components in curvilinear spatial coordinates. Equations (18), (19), and (20) are the continuity, Euler, and energy equations. The Euler and energy equations are derived from γ ν j ∇ µ T µ ν = 0 and n ν ∇ µ T µ ν = 0, respectively. We solve the hydrodynamic equations in the assumption of axial symmetry. Thus, we first write equations in the cylindrical coordinates (̟, ϕ, z). However, the Einstein equations are solved in the y = 0 plane with the Cartesian coordinates. Hence, we rewrite the hydrodynamic equations in the Cartesian coordinates using relations such as ̟ = x and u ϕ = xu y for y = 0. Then, the explicit forms of the equations can be written as where a subscript A denotes x or z, and i, j, · · · here denote x, y, and z. For numerically handling the transport terms as ∂ x (· · ·) and ∂ z (· · ·), we apply an approximate Riemann solver with third-order (piecewise parabolic) spatial interpolation. Other terms are regarded as the source terms. No artificial viscosity is added, in contrast with our previous axisymmetric implementation [21]. The time integration is done with the second-order Runge-Kutta method as explained in [18]. Detailed numerical methods with respect to the treatment of the transport terms are also described in Appendix A. We note that from Eqs. (26) and (28), conservation of baryon rest-mass and angular momentum is derived. However, we write these equations as nonconservative forms and, hence, these conserved quantities are not precisely conserved in numerical computation. To suppress the growth of violation of the conservations in an acceptable level (e.g., within 1%), we should be careful in the grid resolution (see Sec. IV).
In every time step of computation, w at each grid point is obtained by solving the following equation which is derived from the normalization relation of the four-velocity as Here, P = P (ρ, ε) = P [ρ * /(we 6φ ),ê] (see Sec. II C) and ρ = ρ * /(we 6φ ). Thus for a givenγ ij , φ,û i ,ê, and ρ * , Eq. (30) constitutes an algebraic equation for w, which can be solved by standard numerical techniques [32]. After w is obtained, ρ, P , ε, h, and v i can be updated. We note that this procedure is essentially the same as that used in our previous papers (see [18] for details).

C. Equations of state
We adopt two equations of state. One is the so-called Γ-law equation of state of the form where Γ is an adiabatic constant. In using Eq. (31), we always give an initial condition using the polytropic equation of state P = K P ρ Γ of the identical Γ, where K P is a polytropic constant. In this paper, we set the adiabatic constant as Γ = 2 [i.e., the polytropic index n is given by n = 1/(Γ − 1) = 1] as a qualitative approximation of moderately stiff equations of state for neutron stars. We note that if we prepare a polytropic star as an initial condition and the system evolves in an adiabatic manner (with no shock, cooling, and heating), the equation of state is preserved in the polytropic form even using Eq. (31); i.e., the value P/ρ Γ ≡ K P (x µ ) for any fluid element remains a constant (= K P ). The second one is a parametric equation of state that has been used by Yamada and Sato [11], and by Müller and his collaborators [12,14] for the simulation of a rotating stellar core collapse. In this equation of state, we assume that the pressure consists of the sum of polytropic and thermal parts as The polytropic part is in general given as P P = K P (ρ)ρ Γ(ρ) , where K P and Γ are not constants but functions of density ρ. In this paper, we follow [14] for the choice of K P (ρ) and Γ(ρ): For density smaller than the nuclear density ρ nuc ≡ 2 × 10 14 g/cm 3 , Γ = Γ 1 (=const) is set to be < ∼ 4/3, and for ρ ≥ ρ nuc , Γ = Γ 2 (= const) ≥ 2. Thus, where K 1 and K 2 are constants. Since P P should be continuous, we demand that the relation, K 2 = K 1 ρ Γ1−Γ2 , should be satisfied. Following [12,14], we set K 1 = 5 × 10 14 cgs, because we can well approximate the polytropic part of the equation of state for ρ < ρ nuc in which the degenerate pressure of electrons is dominant. Taking into account that the specific internal energy should also be continuous at ρ = ρ nuc , the polytropic specific internal energy ε P is written as With these settings, we mimic a realistic equation of state for high-density, cold nuclear matter. The thermal part of pressure plays a role in the case that shocks are generated. Here, we write it as where ε th ≡ ε − ε P . Following [12,14], we set Γ th = 1.5 in this paper. We performed simulations of rotating stellar collapses using this parametric equation of state. In choosing this, we always give equilibrium stars as initial conditions using the polytropic equation of state where K 0 is a constant. Following [12,14], we set K 0 = 5 × 10 14 cm 3 /s 2 /gr 1/3 , with which a soft equation of state governed by the electron degenerate pressure is well approximated [33]. Here, K 0 and K 1 are related by , where ρ 0 = 1 g/cm 3 .

D. Adding atmosphere
In using high-resolution shock-capturing schemes, we have to add an atmosphere of small density outside stars, since ρ and P have to be nonzero. At t = 0, we put an atmosphere of uniform density and specific internal energy in the computational domain of ρ = 0, according to the following methods. For the Γ-law equation of state with Γ = 2, the uniform density of the atmosphere is set as ρ a = 10 −6 ρ max , where ρ max denotes the maximum density of a star. The specific internal energy is given using the polytropic constant as K P /4. For parametric equations of state (32), the uniform density of the atmosphere is set as ρ a ≈ 1 g/cm 3 . In this case, the specific internal energy is given using the polytropic constant as K 0 .
For the Γ-law equation of state with Γ = 2, the density decreases steeply around the surface of a neutron star. In such a case, numerical instability could often turn on around the stellar surface, if the density of the atmosphere is too low. This is the reason that we attach the atmosphere of relatively high density. On the other hand, a small value of ρ a is acceptable in parametric equations of state.

III. GLOBAL QUANTITIES AND METHOD FOR CALIBRATION
We monitor the conservation of the total baryon rest-mass M * , ADM mass M , and angular momentum J, which are computed in the y = 0 plane as where E = ρhw 2 −P andR is the Ricci scalar with respect toγ ij . M * should be conserved in any system. Because of the axial symmetry, J should also be conserved. On the other hand, M is not conserved in general because of gravitational radiation. However, the total radiated energy of gravitational waves is likely to be quite small in the axisymmetric spacetime, so that we can consider M as an approximately conserved quantity. In our axisymmetric hydrodynamic implementation, M * and J are not guaranteed to be conserved precisely. Thus, monitoring the conservation of them is a good check of numerical accuracy. In addition to the mass and angular momentum, we also check the conservation of the specific angular momentum spectrum [37], where j is the specific angular momentum computed as xû y (= hu ϕ ) and j 0 denotes a particular value for j. The numerical accuracy is also checked monitoring the violation of the Hamiltonian constraint, which is written as where ψ ≡ e φ , and∆ denotes the Laplacian with respect toγ ij . In this paper, we define the averaged violation according to where Namely, we use ρ * as the weight factor for the average. The reason that we introduce this average is as follows. In using high-resolution shock-capturing schemes, we add an atmosphere of small density outside neutron stars and/or collapsing stars. In the atmosphere, a small error in the metric results in a large violation of the Hamiltonian constraint because E is a very small value. Furthermore, the volume fraction occupied by the atmosphere in the whole computational domain is larger than that for main bodies. Thus, if we simply compute the volume integral of |V |, it is close to unity irrespective of the grid resolution. However, the numerical accuracy in the atmosphere is not very important for evolution of the main bodies and for global evolution of the system in which we are interested. Therefore, to monitor whether the main bodies (neutron stars and collapsing stars) are accurately computed or not, this type of weight factor is necessary.

IV. NUMERICAL RESULTS
In the numerical simulations reported in Secs. IV A-C below, we adopted a fixed uniform grid, in which the grid spacing ∆x = ∆y = ∆z is constant, with grid size (N + 1, N + 1) for (x, z) to cover a computational domain as 0 ≤ x, z ≤ L, where L = N ∆x. In the simulation reported in Sec. IV D, we varied the grid spacing during the computation, but still used the uniform grid in which ∆x = ∆y = ∆z. To check the convergence of the numerical results for ∆x → 0, numerical computations were carried out with three levels of the grid resolution while fixing L. All the computations were done on the FACOM VPP5000 machine in the data processing center of the National Astronomical Observatory of Japan. The memory and CPU time in one run with a grid size N = 600 and with four processors are about 1 GByte and 12 CPU hours for ∼ 30000 time steps.
We note that for a simulation with N = 600 in three spatial dimensions, we would need ∼ 300 Gbytes memory and it takes ∼ 1000 CPU hours for 30000 time steps using 32 processors [20]. Such simulation is pragmatically impossible, because the computational costs are too high (note that under normal circumstances, we can use at most ∼ 1000 CPU hours per year). However, taking N = 600 in an axisymmetric simulation is an easy task with the current computational resources.

A. Spherical neutron stars
In this subsection, we focus on the long-term evolution of spherical neutron stars. The initial conditions are given using the polytropic equation of state with Γ = 2, and during the time evolution, the Γ-law equation of state (31) is used. Although we did the same test simulations for the previous hydrodynamic implementation and obtained successful outputs from it [18], we repeated the tests again in the present new implementation to demonstrate that it also works well. A difference in the present tests from the previous ones is that we have performed much longer-term simulations than those in the previous tests, since it is computationally inexpensive and pragmatically possible to do in the axisymmetric case.
In the polytropic equations of state, the polynomial relation c 3−n K n/2 P G −3/2 has dimension of mass. With this property, all the quantities can be scaled to be nondimensional, if we multiply an appropriate combination of c, G, and K P . Thus, we will only show the nondimensional quantities in using this equation of state. In other words, we adopt the units of c = G = K P = 1. In these units, the maximum ADM mass and baryon rest-mass of spherical neutron stars are ≈ 0.164 and 0.180, respectively, with the central density ρ c ≈ 0.318. Recovering K P , the mass and density in dimensional units can be written as In Table I, we list several quantities for five models of spherical neutron stars that we pick up in this subsection. Below, we refer to these models as models (S1) -(S5). Models (S1) -(S4) are stable against gravitational collapse, while (S5) is marginally stable. In Figs    by 0.2%. We note that whenever we superimpose a perturbation to an equilibrium configuration, we reinforce the Hamiltonian and momentum constraints at t = 0. Numerical results are shown for N = 90, 120, and 180 with a fixed value of L, to demonstrate that the convergence is achieved. The simulations continued for ∼ 30 dynamical time scales (see also Fig. 2) until the crash of the run, irrespective of the grid resolution, although the accuracy deteriorates gradually with time. Here, we refer to the period of the fundamental radial (and quasiradial) oscillation as the dynamical time scale.
L is chosen as ∼ 3R s , where R s denotes the coordinate radius of the neutron star. For simulating spherical systems, we imposed the outer boundary conditions as (For nonspherical problems, we impose an outgoing boundary condition forγ ij andÃ ij .) These boundary conditions are adequate but not physically perfect. As a result, for the choice of a too small value of L as ∼ R s , the numerical solution is affected by the spurious effects of the outer boundaries, resulting in an earlier crash of the run. However, for L > 2R s , the results are not significantly modified by the spurious effects.
As we mentioned in Sec. III, the baryon rest-mass is not numerically conserved strictly. However, the violation does not seriously affect the numerical results. Indeed, the averaged values of the central density and lapse remain constants, as they should. The numerical results indicate that if we want to suppress the violation of the mass conservation within 1% (2%) after 10 dynamical time scales, the radius of the neutron star should be covered by more than 40 (30) grid points. The error of mass conservation converges to zero with improving the grid resolution at approximately secondorder. The averaged violation of the Hamiltonian constraint also indicates approximate second-order convergence. Therefore, we can conclude that the numerical solution converges to the exact solution in the limit ∆x → 0. We note, however, that the convergence is only approximately at second-order, because near stellar surfaces, the gradients of hydrodynamic variables are so steep that transport terms are often computed with first-order accuracy in space. The convergence may also become first-order if shocks are generated during numerical computations (see Sec. IV C), since near the shocks, the hydrodynamic computations are done with first-order accuracy. A similar tendency is reported by Miller et al. in the simulations of a head-on collision of two neutron stars [34].
In the last figure of Fig. 1(a), it is shown that K c is not zero exactly but relaxes to a finite value. This indicates that even in solving the equation for the maximal slicing condition, K deviates from zero as long as finite differencing methods are used. Indeed, K c converges to zero at second-order with improving the grid resolution. Thus the maximal slicing condition K = 0 cannot be precisely imposed in numerical computation even in a well-resolved simulation, if we adopt finite differencing schemes. This tells us that we should follow the evolution of K and should not a priori set K = 0 in numerical computation, even in choosing the maximal slicing condition.
Long-term simulations for more compact stable stars (S3) and (S4) were also carried out. In Fig. 2, we display the results for models (S2), (S3), and (S4) together. For (S3), the simulation continued for ∼ 20 dynamical time scales until the crash of the run. However, for (S4), the star starts collapsing to a black hole after about five oscillation periods. The reason for this consequence is clear. The ADM mass of model (S4) is ≈ 99% of the maximum allowed value. Thus, with a slight increase of the mass as a result of the accumulation of numerical error, the mass exceeded the maximum allowed value for stable stars, resulting in eventual gravitational collapse. It is interesting to note that in this case, the computation was able to be continued until the formation of a black hole of the apparent horizon mass [see Eq. (49) for definition] ∼ M , where M is the initial ADM mass. Thus, the increase of the central density and the decrease of α c [see Fig. 2 (a)] do not imply that the computation crashed. To avoid the collapse to a black hole and to make the oscillating time longer, we need to take more grid numbers to improve the grid resolution. Actually, we have checked that we can increase the oscillation cycles in numerical computation with improvement of the grid resolution.
The duration of the simulation to crash for model (S3) was shorter than that of (S2). This fact indicates that with increasing the compactness of neutron stars, the computations crash earlier. We note here that the duration of the simulation for these models depends only weakly on the grid resolution. Thus, it does not depend on the numerical accuracy, but seems to depend on certain factors associated with formulation or gauge conditions or outer boundaries. The same tendency is found in the simulations of rotating neutron stars. We will argue this point in Sec. IV B again.
Although the simulations can be continued only for a finite time scale, the duration of > 10 dynamical time scales seems to be sufficiently long. Indeed, we can accurately extract the oscillation frequencies of the fundamental radial mode from these simulations. In Fig. 2(b), we display the Fourier spectra of ρ c (t), which is defined as where t f is chosen as ∼ 70ρ The Fourier spectra indicate that the oscillation period of the f mode of the radial oscillation is ≈ 5.0ρ . These values agree well with those derived from Chandrasekhar's semianalytic formula [35,18]. Furthermore, the result for (S2) is in good agreement with that for a spherical star of ρ c = 0.128 reported in [28]. Thus, we conclude that the computation can be continued for a sufficiently long time to accurately obtain the oscillation frequencies of even extremely relativistic neutron stars. The simulation would be able to be carried out to study nonspherical oscillations of neutron stars, as was done in [18,36]. In Fig. 3, we display the time evolution of several quantities for collapse of a marginally stable spherical neutron star (S5). To induce the collapse, we initially reduced the pressure by 0.5%. With this setting, the neutron star collapses to a black hole in ∼ 10ρ −1/2 c,init . We have checked that even with a 0.2% decrease of the pressure, the star collapses to a black hole in ∼ 20ρ −1/2 c,init , which is longer than that for the case of 0.5%. Numerical results are presented for N = 90, 120, and 180 and demonstrate that the convergence is achieved. Note that with lower grid resolution, it takes a slightly longer time to form a black hole. This is because the dissipation that prevents the increase of density is larger with lower grid resolution.
In the final phase of the collapse, the grid resolution around the black hole forming region became so bad that the computation crashed. It is interesting to note that the magnitude of the averaged violation of the Hamiltonian constraint, ERROR, relaxes to ∼ 0.2 irrespective of the grid resolution at t ∼ 10ρ −1/2 c,init . This verifies that the computation crashes when ERROR is ∼ 0.2. We can also observe that as the accuracy deteriorates, (i) the central density stops increasing and instead starts decreasing, and (ii) exponential decrease of α c that is a feature in the maximal slicing condition is modified. The reason for (i) is as follows. In our implementation, ρ * is a fundamental quantity to evolve, and ρ is computed from ρ * /(e 6φ w). In numerical computations, ρ * and φ monotonically increase, but since φ around the origin is too large [of O(1)] in the late phase of the collapse, a small error in φ leads to a large error in ρ. As a result, ρ decreases in the late phase. The reason for (ii) is simply that the computation crashed. Indeed, the time at which the behavior of α c starts changing agrees with that at which the magnitude of the averaged violation of the Hamiltonian constraint saturates to ∼ 0.2.
Although the accuracy deteriorates in the final phase of the collapse, the simulation can be carried out at least until the formation of an almost static black hole. To confirm the black hole formation, apparent horizons were located during the simulations. In Fig. 3(b), we display the time evolution of the mass of the apparent horizon in units of the ADM mass of the system. Here, the mass of the apparent horizon is defined as In Fig. 4(a), we show the time evolution of several quantities for an oscillating spherical neutron star of a high amplitude. In this simulation, we picked up a low-mass spherical star (S1), and to induce an oscillation of a high amplitude, we initially reduced the pressure by 40%. The numerical results are shown for N = 90, 120, and 180, and demonstrate that convergence is achieved. In each case, the neutron star is initially covered by 39, 52, and 78 grid points, respectively.
Because of a significant decrease of the pressure, the radius of the neutron star decreases by a factor of ∼ 2 soon after the simulation starts. However, the magnitude of the pressure decrease is not large enough for the star to collapse to a black hole. Instead, the star bounces when the central density becomes about six times of initial value, and repeats oscillations subsequently. As the density approaches the maximum, shocks are formed around the stellar surface, and as a result, outer envelopes explode. To illustrate that shock heating indeed occurs, we display ε/(K P ρ Γ−1 ) at t ≈ 0 and after about one oscillation period in Fig. 4(b). As mentioned in Sec. II C, in the absence of shocks, this quantity does not change from the initial value, but in the presence of shock heating it increases. Figure 4(b) clearly shows that in the outer envelope, the shock heating is significant. On the other hand, a negligible effect of shocks can be seen around the central region. Since the shocks are generated only in the atmosphere, the averaged violation of the Hamiltonian constraint still converges approximately at second-order with improving the grid resolution. Figure 4(a) indicates that the amplitude of the oscillation gradually decreases, and after several oscillation periods, it settles down approximately to a constant. This illustrates that the kinetic energy of the oscillation is dissipated by the shocks gradually. Similar results are reported in [28] for a simulation of a migrating neutron star. The lower figure of Fig. 4(b) shows the density profiles at t = 0 and t ∼ one oscillation period. This indicates that the density profile is modified to a more centrally condensed state as a result of the shock dissipation.
We emphasize that this simulation is highly dynamical and general relativistic. Even in such a case, the simulation was able to be continued for more than 20 oscillation periods. This illustrates the robustness of our implementation for dynamical problems in general relativity.  c,init , J/M 2 , and |T /W | of rotating neutron stars at mass shedding limits with Γ = 2 that we pick up in this paper. Here T and W are the rotational kinetic energy and gravitational potential energy, and R denotes the circumference radius at the equator. All the quantities are shown in units of c = G = KP = 1. The star denoted with † is unstable. In the last column, numerical results of the quasiradial oscillation period for the f mode, Posc, are presented in units of ρ −1/2 c,init .

B. Rapidly rotating neutron stars
We focus here on the long-term evolution of rigidly and rapidly rotating neutron stars at mass shedding limits for which the angular velocity at the equator is equal to the Kepler angular velocity. Following Sec. IV A, we adopt the polytropic equation of state with Γ = 2 for setting initial conditions and evolve neutron stars using the Γ-law equation of state (31). The axial ratio of the polar radius to the equatorial one is ≈ 0.58 for rotating neutron stars at mass shedding limits with Γ = 2. As in Sec. IV A, we only present the scaled dimensionless quantities with the units of c = G = K P = 1 throughout this subsection. In these units, the maximum ADM mass and baryon rest-mass of rigidly rotating neutron stars are about 0.188 and 0.207, respectively, with the central density ρ c ≈ 0.27 [38]. The central density of the marginally stable star against gravitational collapse has slightly larger density (≈ 0.295) than this value [38,29]. In Table II, we list several quantities of six rotating neutron stars that we pick up here. In the following, we refer to these neutron stars as models (R0) -(R5). Models (R0) -(R4) are stable against gravitational collapse, and (R5) is unstable and very close to the marginally stable point. The simulations continued for ∼ 10 dynamical time scales, and eventually crashed. The duration of the simulation to crash depends only weakly on the grid resolution as in the spherical cases. This implies that the crash is not triggered by accumulation of the numerical error. The duration also does not vary much even if we change the outer boundary conditions and the location of the outer boundary as L ≈ 3R e -4R e , where R e denotes the coordinate radius at the equator. Furthermore, the duration of the simulation for a rotating neutron star is shorter than that for a spherical star of identical compactness. Therefore, we deduce that the crash of computations might be associated with our choice of the spatial gauge condition or the formulation, although we do not understand the reason fully at present. There may still be room to improve the spatial gauge condition and/or the formulation, if one wants to perform an extremely long-term simulation of the duration ≫ 10 dynamical time scales. However, 10 dynamical time scales are long enough to produce scientific results for most problems, so that we do not address this problem any longer in this paper.
As in the spherical case, convergence is achieved with improvement of the grid resolution. As argued in Sec. III, the angular momentum as well as the baryon rest-mass are not conserved strictly, although they are conserved quantities. However, the violation converges to zero nearly at second-order with improving the grid resolution. The results of this convergence test indicate that the polar axis should be covered by at least ∼ 30 grid points if one wants to demand that the violation of the conservation of angular momentum and baryon rest-mass is less than a few % after ∼ 10 dynamical time scales. If the polar axis is covered by fewer than 20 grid points, the magnitude of the violation becomes larger than 10% after 10 dynamical time scales.
The long-term simulations were also performed for models (R0), (R2), (R3), and (R4). In all these simulations, we initially reduced the pressure by 0.5% and take N = 180. With this grid number, the polar (equatorial) radius is covered by 35 (60) grid points. The time evolution of the central density and central value of the lapse function (α c ) are shown together in Fig. 6(a). As in the spherical case, the duration of the simulation is shorter for more compact stars. (The moment of the crash of a run is identified with the time at which α c sharply drops.) This might be evidence that with our spatial gauge, coordinate distortion is accumulated too much for the long-term simulation, because it is likely to be accumulated more rapidly for more compact stars. As in the simulation for (S4), a high-mass rotating star (R4), for which the ADM mass is ≈ 99.5% of the maximum, collapses to a black hole in a few dynamical time scales instead of crashing, because of a slight increase of the baryon rest-mass due to numerical error. It is necessary to take N > 180 to continue computations of such high-mass stars for more than two oscillation periods. However, as long as the neutron star is not very close to the marginally stable point, the simulation can be continued for more than five dynamical time scales with N ∼ 200, and this duration is long enough to produce scientific results for most problems. For example, from these simulations, we can extract the frequency of fundamental quasi-radial oscillation modes. In Fig. 6(b), we show the Fourier spectra of the central density as in the case of Sec. IV A. The figure indicates clear peaks that denote the fundamental frequency of the quasiradial oscillation modes.
In Fig. 7, we summarize the frequencies of the radial and quasiradial oscillation modes for spherical and rapidly rotating neutron stars with Γ = 2. The filled and open circles denote the numerical results for spherical and rotating neutron stars, respectively. The filled circle along the vertical axis is plotted according to a Newtonian analysis for the spherical polytrope [39]. On the other hand, the filled and open circles along the horizontal axis are plotted by the fact that the frequency of the fundamental radial and quasiradial oscillations of the marginally stable stars is zero. We note that the frequency in dimensional units is computed from where P osc = 1/f is the oscillation period of the fundamental quasiradial mode. As mentioned in Sec. IV A, the frequencies for spherical neutron stars agree well with semianalytical results [35]. For rotating neutron stars, the frequency is slightly smaller than that for spherical neutron stars for identical ρ c . A similar tendency is found in the results of [28]. It is interesting to note that for the identical value of the central lapse function, the frequencies approximately coincide for α c < ∼ 0.6. In Fig. 8(a), we display the central density, central value of the lapse function, and averaged violation of the Hamiltonian constraint for collapse of an unstable rotating neutron star (R5). To induce the collapse, we initially reduced the pressure by 0.5%. We also carried out the simulations with the reduced factor of 0.2%, and have found that the star (R5) collapses also in this case, although it takes longer to be a black hole. Numerical results are shown for N = 240, 180, and 120 and demonstrate that the convergence is achieved well. With these grid numbers, the polar (equatorial) radius is initially covered by about 69 (120), 52 (90), and 35 (60), respectively.
During the collapse, the density (lapse function) monotonically increases (decreases) with time, and finally a black hole is formed (i.e., the apparent horizon is located) in the late time when α c < ∼ 0.03. As in the spherical case, the reason that ρ c decreases in the late stage of the collapse is that a small error in φ that is of O(1) at the formation of the black hole leads to a large error in ρ that is computed from ρ * /(we 6φ ). It is also found that at the time when the magnitude of the averaged violation of the Hamiltonian constraint becomes ∼ 0.2, the computation crashed due to the grid stretching around the horizon of a black hole, and that the magnitude of ERROR at the crash is almost independent of the grid resolution.
In Fig. 8(b), we show the time evolution of the mass of the apparent horizon defined by Eq. (49). Here, the long-dashed horizontal line denotes the expected final value S EH /(16πM 2 ), where S EH is derived from the formula for the area of the event horizon of a Kerr black hole as [33] S EH = 8πM 2  where M and J are the ADM mass and angular momentum of the collapsing neutron star. Since almost all the matter eventually falls into a black hole in this simulation, the area of the apparent horizon should settle down to S EH . The figure indicates that the area of the apparent horizon asymptotically approaches the expected value. This demonstrates that the spacetime in the final phase of our simulation almost relaxes to a stationary, Kerr black hole spacetime.

C. Collapse of rotating stars with parametric equations of state
The purpose of this subsection is to demonstrate that with our implementation, it is feasible to carry out stable and accurate simulations for the collapse of rotating stars with parametric equations of state (32) that are more realistic for high-density matter than the Γ-law equation of state used in the previous two subsections. Since the equation of state for high-density matter is still not precisely known, the parametric equation of state (32) is used for several choices of Γ 1 and Γ 2 . Initial conditions are set up adopting the polytropic equation of state (36) with Γ = 4/3.
In the realistic core collapse of massive stars, the central density just before the collapse is of order 10 10 g/cm 3 [14,40]. Since the collapse leads to the formation of a neutron star of density of order 10 15 g/cm 3 or a black hole, the characteristic length scale changes by a factor of ∼ 100. This implies that we need to take N of O(10 3 ) for a well-resolved simulation in the fixed uniform grid. Although it is possible to take a large value of N as several thousands, performing such large-scale simulation is not computationally inexpensive even in the axisymmetric case. Since the main purpose of this subsection is not to present scientific results, but both to demonstrate that realistic equations of state can be adopted in our implementation and to grasp characteristic behaviors associated with new implementation with such equations of state, here we pick up more compact stars of central density ≈ 6 × 10 12 g/cm 3 as initial conditions to save the computational costs as a first step. In the next subsection, we will show a numerical result with a more realistic initial data set as illustration.
Velocity profiles of equilibrium rotating stars used as initial conditions are given according to a popular relation [41,42], where Ω a denotes the angular velocity along the z axis, and ̟ d is a constant. In the Newtonian limit, the rotational profile is written as Thus, ̟ d indicates the steepness of differential rotation. In this paper, we pick up the rigidly rotating cases in which ̟ d → ∞ and a differentially rotating case in which ̟ d /R e = 1/2, where R e is the equatorial coordinate radius. In both cases, we choose the axial ratio of polar radius to equatorial radius as 2/3. In Table III, we list several quantities for the models we adopt in the numerical computation. We refer to these models as (C1) and (C2).
The simulations were carried out for (Γ 1 , Γ 2 ) = (1.325, 2) (a), (1.3,2) (b), (1.325,2.5) (c), and (1.3, 2.5) (d). In the following, we specify our choice of Γ 1 and Γ 2 in terms of (a), (b), (c), and (d). For example, if we pick up model (C1) as the initial condition and choose Γ 1 = 1.325 and Γ 2 = 2, we refer to this model as (C1a). In Table IV, we list masses and central density for spherical neutron stars of maximum mass in the parametric, cold equation of states (33) with four choices of Γ 1 and Γ 2 , which are obtained by numerically solving the TOV equation [33]. For Γ 1 = 1.3 and Γ 2 = 2, the maximum mass is too small to be a realistic value, but with this model, we can study qualitative properties of stellar core collapses in an extremely soft equation of state.
According to a numerical result [43], rigidly rotating stars in equilibrium with Γ = 4/3 are unconditionally unstable against gravitational collapse, if the compactness M/R e is larger than ∼ 1/700. Thus, the equilibrium state of model (C1) is unstable. The criterion of the instability has not been established yet for differentially rotating stars. Thus, the stability is not clear for (C2). However, since the compactness is much larger than 1/700, the initial equilibrium of (C2) is also likely to be unstable. Hence, by decreasing Γ from 4/3 to Γ 1 < 4/3 at t = 0, the collapse is accelerated. To investigate convergence, the simulations were carried out with N = 600, 400, and 300 for all the models. In these simulations, the equatorial radius was initially covered by N grid points. (C2) has a slightly toroidal shape due to the effect of differential rotation. In Figs. 11 and 12, we also show the time evolution of several quantities for models (C1a) and (C2a). Here M * core denotes the total baryon rest-mass of high-density matter with ρ > ρ nuc . As indicated in these figures, the collapses proceed in the following manner. In the early phase in which the central density is smaller than ρ nuc , the central region of the star contracts approximately in a homologous manner because the adiabatic index is close to 4/3 [44]. After the central density exceeds ρ nuc , a core is formed at the central region, and the mass gradually increases as a result of subsequent accretion. At the time when the central density becomes ∼ 4.5ρ nuc for (C1a) and ∼ 3.5ρ nuc for (C2a), the increase of the core mass stops, and strong shocks rapidly propagate outward due to the restoring force of the core. After this moment, the core settles down toward an approximately stationary state, and finally the central density relaxes to ∼ 3ρ nuc for (C1a) and ∼ 2.5ρ nuc for (C2a). On the other hand, the shocks propagate outward, sweeping the infalling matter. In both cases (C1a) and (C2a), the baryon rest-mass of ρ ≥ ρ nuc is ∼ 0.65M ⊙ after shock propagates far from the core. Since the angular momentum at t = 0 for (C2a) is larger than that for (C1a), the baryon rest-mass of the core and central density of the final state of the core are slightly smaller for model (C2a). These facts indicate that the products after collapse depend basically on the equation of state, but are modified by the magnitude and distribution of angular momentum that are initially retained by a precollapse star. If the initial angular momentum is much larger than those for (C1a) and (C2a), the collapse would be halted before the central density reaches ρ nuc , and as a result, a neutron star would not be formed [14].
Another noticeable feature is that for the differentially rotating initial condition, a small modulation is found in the time evolution of the central density and α c after a quasistatic core is formed. This is likely because in the collapse with a differentially rotating initial condition, the collapsing star deviates highly from spherical symmetry and shocks are formed in a nonspherical manner. As a result, the nonspherical oscillation is excited in the formed core. Similar results are observed in [14].
In Figs. 11(b) and 12(b), we display the violation of the baryon rest-mass conservation, angular momentum conservation, and averaged violation of the Hamiltonian constraint. First, we note that at t ∼ 6.5 msec, the shocks reach the outer boundaries, and the matter escapes from the computational domain. This causes the rapid angular momentum decrease for t > ∼ 6.5 msec. Besides this, the errors converge to zero with improving the grid resolution, as we saw in Secs. IV A and IV B. However, there are several noticeable properties that have not been seen in the previous subsections. One is that the magnitude of the errors quickly increases when shocks are generated. The second is that the averaged violation of the Hamiltonian constraint converges to zero at about first-order (not at second-order). We cannot explain the reason for them correctly. However, we deduce it as follows: At the shocks, the hydrodynamic quantities are discontinuous, and as a result, the second partial derivatives of the metric are discontinuous and the complete analyticity for the metric is violated. This may change the global convergence property from second-order to first-order.
The magnitude of the errors for (C2a) is smaller than that for (C1a) with identical grid spacing. This is a consequence of the fact that the central density for (C2a) after shock formation is smaller than that for (C1a), and as a result, the grid resolution for (C2a) is better than that for (C1a).
To see the effects of the equations of state on the dynamics of collapse and on accumulation of the numerical errors, we display numerical results for (C1a), (C1b), (C1c), and (C1d) in Fig. 13. All the simulations were done with N = 600, by which the equatorial radius is initially covered by N grid points. Reflecting the difference of Γ 1 and Γ 2 , the products after the collapse and time evolution of the numerical errors are different. The following is a summary of the differences: (a) for Γ 1 = 1.3, the magnitude of the errors is larger than that for Γ 1 = 1.325, (b) for the identical value of Γ 1 , the central density and M * core in the final state of the core are smaller for the larger value of Γ 2 , (c) M * core is larger for the larger value of Γ 1 with the identical value of Γ 2 , (d) the magnitude of the error of the baryon rest-mass conservation always increases with time (never decreases in the absence of mass ejection from the computational domain) irrespective of Γ 1 and Γ 2 , and (e) the angular momentum increases due to the numerical error for Γ 1 = 1.325 but decreases for Γ 1 = 1.3. The reason for (a) is clear because for smaller values of Γ 1 , the central density increases by a large factor, and as a result, the grid resolution becomes worse even in the identical grid spacing. The reason for (b) is also clear because with stiffer equations of state, the central density of neutron stars is smaller. The result (c) comes from the fact that with smaller values of Γ 1 , shocks are stronger, and as a result, the amount of mass that is ejected outside the core is larger. However, the reason for (d) and (e) is not clear at all. These may be consequences of our choice of the finite differencing scheme for the hydrodynamic equations. Indeed, a similar tendency is found in the results of an approximate general relativistic simulation [14].
It is interesting to note that even for (C1b), in which the maximum allowed baryon rest-mass of neutron stars is very small (∼ 1M ⊙ ), a low-mass neutron star of M * core ∼ 0.55M ⊙ is formed after the collapse. This is because the shocks explode the infalling matter sufficiently. Thus, even in the case that the stellar core mass before collapse is much larger than the maximum allowed mass of the neutron star for a given equation of state, a neutron star instead of a black hole could be formed at least temporarily, although such a neutron star could easily collapse to a black hole by subsequent accretion of matter or by a fall-back.
In Fig. 14, the specific angular momentum spectra at selected time steps are shown for (C1d) and (C2a) as examples. In both cases, we take N = 600. The figures demonstrate that the spectral shape is well conserved throughout the simulations irrespective of the initial angular velocity profile. Therefore, we conclude that angular momentum transfer due to numerical dissipation is sufficiently small in numerical computations.
D. An example of stellar core collapse In this subsection, we present numerical results for a simulation of rotating stellar core collapse that is started from a realistic initial condition with the central density ∼ 10 10 g/cm 3 and with Γ = 4/3. We give a rigidly rotating equilibrium star that is close to the mass shedding limit as the initial condition. Several quantities of this rotating star are listed in Table V. Since it is rapidly rotating and its compactness is sufficiently small, this equilibrium star is dynamically stable against gravitational collapse even for the polytropic equation of state with Γ = 4/3 [43]. Thus, the collapse is triggered by the slight decrease of the adiabatic constant from 4/3 to Γ 1 . In the present simulation, we set Γ 1 = 1.325 and Γ 2 = 2. Since the characteristic length scale changes by a factor of ∼ 100 during the collapse, we performed the simulation changing the grid size and grid number as done in [22]. The grid size and computational domain were changed monitoring the value of the lapse function at the center (α c ), which approximately indicates the compactness of the collapsing star. Whenever we carried out regridding, we made the grid spacing half and used cubic interpolation [32] for assigning the values of variables on the finer grids. The simulation was started with N = 500, by which the equatorial radius is covered by 480 grid points initially. At t = 0, α c ≈ 0.993. We carried out the first regridding when α c was 0.975, at which the mean radius of the collapsing star became ∼ 1/4 of the initial one. In this regridding, we chose N = 900 and made the grid spacing half. The next regridding was carried out when α c = 0.95 and 0.90, and we chose N = 1500 and 2100, respectively. After α c reached 0.90, we fixed N and grid spacing. L and ∆x in the final stage are about 1050 km and 0.5 km, respectively. For this simulation, the computational time was about 100 CPU hours for ≈ 60000 time steps using eight processors of the FACOM VPP 5000 machine.
Since the computational region was reduced whenever we carried out the regridding, a small amount of mass that is outside the new computational domain was discarded. However, the magnitude of the violation of mass and angular momentum conservation is less than ∼ 0.5% and 2%, respectively [see Fig. 16(b)]. This implies that the total amount of the discarded mass is comparable to that of the numerical error associated with the finite differencing, which does not much affect the evolution of the system, as indicated in Sec. IV C.
In Fig. 15, we display the density contour curves and velocity fields at selected time steps around which shocks are formed. The time of the shock formation is ∼ 71.7 msec, which is in good agreement with that for model A1B3G1 in [14] with the correction factor which is associated with the dynamical time scale as (ρ c,init /10 10 g/cm 3 ) −1/2 . (Note that in [14], the central density of the initial condition is 10 10 g/cm 3 , while here ρ c,init ≈ 1.65 × 10 10 g/cm 3 .) This coincidence suggests that the approximate general relativistic approach adopted in [14] is indeed suited for study of axisymmetric stellar core collapse to neutron stars in general relativity.
After the shock formation, the shock fronts of prolate shape spread outward. The prolateness is produced by the fact that the shocks are stronger for the z direction due to the absence of centrifugal force [11,12,14]. In Fig. 16(a), we display the time evolution of the central density and lapse function at the center. Global features are qualitatively the same as those for the simulation of model (C1a) presented in Sec. IV C. As in that case, the central density (lapse function) monotonically increases (decreases) until it exceeds ρ nuc . When it becomes ∼ 3.5ρ nuc , the collapse is halted  and the shocks start propagating outward, while the core gradually settles down toward a quasistationary state of ρ c ∼ 2ρ nuc . Although these qualitative features are the same as those for (C1a), there are a few quantitative differences between the results of two models. First, the maximum density and central density of the formed neutron star found here are slightly smaller than those for (C1a). This is likely due to the fact that the effect of the centrifugal force plays a more important role than for (C1a). Second, the core does not relax to a static state soon, but oscillates approximately in a quasiperiodic manner for several periods. This oscillation is not conspicuous for model (C1a). These facts imply that to obtain quantitative outputs of stellar core collapse, we should start the simulations from an initial condition of a realistic density profile, although qualitative global features of the collapse can be found even using a more compact initial condition. We note that the quasiperiodic oscillation found here is also observed in [14]. As reported in [12,14], quasiperiodic gravitational waves are likely emitted associated with this oscillation. However, we have not tried to compute gravitational waves in the present work. As expected from the results in [14], the amplitude of gravitational waves is not very large, so that it would not be technically easy to extract them from the metric in which gauge modes and numerical noises are included. Developing a method for the wave extraction of a weak signal will be one of the challenging problems in the future.
In the last figure of Fig. 16(b), the time evolution of baryon rest-mass of the core (of density larger than ρ nuc ) is shown. As in the time evolution of ρ c , it reaches the maximum at which ρ c reaches the maximum value, and then settles down toward a constant ∼ 0.55M ⊙ . Thus, the temporary product after the collapse is a low mass neutron star in this simulation.

V. SUMMARY
We have presented numerical results obtained by an axisymmetric general relativistic implementation, and demonstrated that with this new implementation, it is feasible to carry out long-term simulations for spherical and rapidly rotating neutron stars, and rotating stellar collapses to a neutron star and a Kerr black hole. It is shown that the simulations for stable neutron stars can be continued for > ∼ 10 dynamical time scales until the crash of computation with N ∼ 200, even if the neutron stars are of ∼ 95% of the maximum allowed mass. The duration of the computation is long enough to obtain oscillation modes of neutron stars. The simulations are also feasible for collapse of rotating neutron stars to Kerr black holes. We have illustrated that with our implementation, the mass of a Kerr black hole formed after the collapse can be computed accurately. We have also demonstrated that it is feasible to perform the simulations of rotating stellar core collapse to a neutron star, adopting parametric equations of state that mimic realistic equations of state. This illustrates that the new implementation works well also for realistic equations of state that have not been adopted so far in fully general relativistic simulations. In conclusion, the axisymmetric numerical implementation presented here will be used for a wide variety of astrophysical simulations such as rotating core collapse of a massive star to a neutron star or a black hole, and accretion-induced collapse of a neutron star to a black hole. As the next step, we plan to perform simulations for rotating core collapse to neutron stars, and compare the results with those in an approximate relativistic approach [14]. We also consider that collapse of a massive stellar core to a black hole is one of the most interesting topics.
Because of the assumption of axial symmetry, we were able to carry out a wide variety of tests and calibrations for our new hydrodynamic implementation with low computational costs. In previous works, e.g., [18,28], tests of their hydrodynamic implementations picking up single stars have been done in three spatial dimensions. In those cases, N could be at most 100, since the computational costs were very high for N > 100. It is pragmatically very difficult to investigate the accuracy and convergence in a well-resolved simulation with N ∼ several hundreds under normal circumstances in which we can use at most ∼ 1000 CPU hours per year. In the axisymmetric case, the simulation with N ∼ several hundreds is not very expensive, and thus for detailed tests of new hydrodynamic implementations in general relativity, axisymmetric simulation has great advantages. We expect that for the testing of new gauge conditions and wave extraction techniques, it would also play an important role.
Finally, we note the following point. Although we focus only on the axisymmetric simulation in this paper, the present hydrodynamic implementation can be used for three-dimensional simulations with a slight modification, since the transport terms in the hydrodynamic equations are of the same form. Actually, we have already implemented it and checked that it works. We expect that with the same grid resolution that we adopted in this paper, the same results will be obtained (although it takes a much longer time to carry out the long-term simulations). So far, we have performed simulations for a merger of binary neutron stars only using the Γ-law equation of state [19,20]. However, with the new hydrodynamic implementation reported in this paper, we will be able to adopt a variety of equations of state. We plan to perform the simulation for a merger of binary neutron stars adopting more realistic equations of state in the future.
Thus, M A ab has information on the characteristic speed of the fluid. Following Font et al. [23,24], we calculate the Jacobian matrix from where and q c = (ρ, v x , v y , v z , ε). Explicit forms for C ab and B x ab in our notation are and where B z ab is obtained by appropriate exchanges of subscripts among x, y, z of B x ab . The eigenvalues of the matrix M A ab , λ A , correspond to the characteristic speeds of the fluid in the A-th direction, and are derived from the equation The solutions are [24,25] where and Using the eigenvalues, the spectrum decomposition for M A ab can be done in a straightforward manner as where Λ A bc is the diagonal matrix composed of λ A in the following order: For convenience of the calculation of R A ab , we define a matrix T A ab , which satisfies the relation as Since R A ab is calculated from the right eigenvectors of M A ab , T A ab is composed of vectors t that satisfy the equation Then, T A ab = [(t A a ) (1) , (t A a ) (2) , (t A a ) (3) , (t A a ) (4) , (t A a ) (5) ], and hence we obtain where After the above reconstruction of the fluid equation, the numerical fluxes in the upwind scheme are computed from where we omit the subscripts for R ab and Λ ab . Q l c and Q r c denote Q c at the left and right sides of the corresponding interfaces, and are evaluated using the third-order spatial interpolation. At the interface between the i-th and (i+1)-th cells, we define them according to where ∆ i = Q i+1 − Q i . To suppress the oscillation near shock discontinuities, we modify the interpolation using the following min-mod limiter as [46] Q l c = Q i + where For the simulations presented in this paper, we choose b = 2, since for b ≈ 1, the dissipation is so large that the envelope of neutron stars spreads outward too quickly, while for b = 4, the oscillation around shock discontinuities is too serious. The values for components of matrices R ab and Λ ab at grid interfaces are computed using the Roe-type average such as [45,46] q i+1/2 = (ρ * ) i+1 q i+1 + (ρ * ) i q i where we carry out the average for variablesû i , κ, χ, and h (i.e., q i is one of these variables). Other variables are computed from them. It should be noted that in the relativistic case, the average is not uniquely specified in contrast with the Newtonian case [45]. However, numerical results of test computations (see below) seem to indicate that this averaging is appropriate.
To confirm that our hydrodynamic implementation can capture shocks accurately, we carried out the simulations for Riemann shock tube problems and wall shock problems in the 1+1 special relativistic spacetime with (t, x) as the coordinates. In this test, we adopt the Γ-law equation of state. In both tests, we take N = 400 with ∆x = 1/N .
In Figs. 17 and 18, we compare numerical results of the Riemann shock tube problem and of the wall shock problem with analytical solution (solid curves) for the choice b = 1 and 2. As indicated in these figures, numerical results agree well with analytic solutions, in particular for b = 2.