format: text

text:
# Relativistic simulations of rotational core collapse. I. Methods, initial models, and code tests

## 1 Introduction

The advent of gravitational wave astronomy – major interferometric gravitational wave detectors have already started taking data – aggravates the need for gravitational wave templates computed from theoretical models of potential sources of gravitational radiation (Pradier et al. 2001). Among the most promising sources are aspherical core collapse supernovae. However, reliable core collapse simulations, in particular those of rotational core collapse, require a relativistic treatment of the dynamics because of the counteracting stabilizing effect of rotation and the destabilizing effect of relativistic gravity. In addition, the simulations should incorporate a realistic equation of state (EoS) and an accurate treatment of neutrino transport.

Due to these complexities, all previous investigations have either neglected or approximated one or all of these requirements (see, e.g. Müller (1998) and references therein). Up to now there exist no successful multidimensional numerical simulations of rotational core collapse to neutron stars in general relativity, even for simple matter models. On the other hand, simulations which include better microphysics in the form of realistic (nuclear) equations of state or neutrino transport have either been confined to spherical symmetry or have used Newtonian gravity.

There is an obvious need to improve this situation if one wants to answer questions such as: What is the quantitative influence of general relativity on the dynamics of core collapse compared to Newtonian gravity? How do higher densities and velocities expected in relativistic simulations change the properties of the rotating proto-neutron star? To what extent do relativistic effects alter the gravitational wave signal? How does the rotation rate evolve in a relativistic simulation, and what consequences does that have on the development of dynamical or secular instabilities in the neutron star?

To this end we have started a project to compute reliable gravitational wave templates from rotational core collapse and nonspherical core collapse supernovae. In a first step we have improved the treatment of the dynamics beyond the still commonly used Newtonian approach and computed the general relativistic hydrodynamics of rotational core collapse consistently with the evolution of the axisymmetric dynamic spacetime. The computational tool developed by us for this purpose is described in this publication, while its application to the gravitational collapse of rotating stellar cores approximated as polytropes is devoted to an accompanying paper (Dimmelmeier et al. 2002). We note that in Dimmelmeier et al. (2001) we have already presented preliminary results of our investigation. General relativistic simulations of improved models including both a realistic equation of state and an approximate treatment of weak interactions and neutrino transport are planned in the near future.

The asumptions we adopt in the present and accompanying paper are as follows: Our matter model obeys an ideal gas EoS with the pressure $`P`$ consisting of a polytropic and a thermal part (Janka et al. 1993; Zwerger & Müller 1997). Since we are mostly interested in the gravitational radiation emission, which is controlled by the bulk motion of the fluid, we neglect any other microphysics like electron capture and neutrino transport. The initial configurations are rotating relativistic polytropes in equilibrium (Komatsu et al. 1989a, b; Stergioulas & Friedman 1995) which are marginally stable. Their collapse is initiated by decreasing the effective adiabatic index $`\gamma _\mathrm{p}`$ to some prescribed fixed value. Our initial data solver allows us to construct both uniformly and differentially rotating polytropes.

We solve the general relativistic hydrodynamic equations numerically. The metric equations are approximated by adopting the so-called conformal flatness condition for the three-metric $`\gamma _{ij}`$ (Wilson et al. 1996). This results in a constrained system of elliptic metric equations which are much simpler than the full system of Einstein equations. The dimensionality of the problem is reduced by assuming both axisymmetry with respect to the rotation axis, and equatorial plane symmetry. The numerical code uses a spherical grid in spherical polar $`(r,\theta )`$ coordinates.

This paper is organized as follows: In Section 2 we present the mathematical framework used in our approach, namely the conformally flat metric equations and the relativistic hydrodynamic equations. Section 3 is devoted to describing the initial models. Next, in Section 4 we describe all relevant issues concerning our numerical implementation. Section 5 presents a comprehensive number of tests we have used to validate our code. The paper ends with a summary in Section 6.

Unless otherwise stated we use geometrized units ($`G=c=1`$). Greek (Latin) indices run from 0 to 3 (1 to 3). Latin indices are raised and lowered using the three-metric.

## 2 Mathematical framework

### 2.1 Metric equations

We adopt the ADM 3 + 1 formalism (Arnowitt et al. 1962) to foliate the spacetime into a set of non-intersecting spacelike hypersurfaces. The line element reads

$$ds^2=\alpha ^2dt^2+\gamma _{ij}(dx^i+\beta ^idt)(dx^j+\beta ^jdt),$$
(1)
with $`\alpha `$ being the lapse function, which describes the rate of advance of time along a timelike unit vector $`n^\mu `$ normal to a hypersurface, $`\beta ^i`$ being the spacelike shift three-vector which describes the motion of coordinates within a surface, and $`\gamma _{ij}`$ being the spatial three-metric.

In the 3 + 1 formalism, the Einstein equations split into evolution equations for the three-metric $`\gamma _{ij}`$ and the extrinsic curvature $`K_{ij}`$, and constraint equations which must be fulfilled at every spacelike hypersurface:

$`_t\gamma _{ij}`$ $`=`$ $`2\alpha K_{ij}+_i\beta _j+_j\beta _i,`$ (2)
$`_tK_{ij}`$ $`=`$ $`_i_j\alpha +\alpha (R_{ij}+KK_{ij}2K_{ik}K_j^k)`$ (3)
$`+\beta ^k_kK_{ij}+K_{ik}_j\beta ^k+K_{jk}_i\beta ^k`$
$`8\pi \alpha (S_{ij}{\displaystyle \frac{\gamma _{ij}}{2}}(S\rho _\mathrm{H})),`$
$`0`$ $`=`$ $`R+K^2K_{ij}K^{ij}16\pi \rho _\mathrm{H},`$ (4)
$`0`$ $`=`$ $`_i(K^{ij}\gamma ^{ij}K)8\pi S^j.`$ (5)
In these equations $`_i`$ is the covariant derivative with respect to the three-metric $`\gamma _{ij}`$, $`K_{ij}`$ is the extrinsic curvature, $`R_{ij}`$ is the Ricci tensor, $`R`$ is the scalar curvature and $`K=K_i^i`$ is the trace of the extrinsic curvature. The matter fields appearing in the above equations, $`S_{ij}`$, $`S^j`$ and $`\rho _\mathrm{H}`$, are the spatial components of the stress-energy tensor, the momenta and the total energy, respectively.

Following Wilson et al. (1996) (see also Flanagan (1999); Mathews & Wilson (2000)) we approximate the general metric $`g_{\mu \nu }`$ by replacing its spatial three-metric $`\gamma _{ij}`$ with the conformally flat (CF) three-metric (conformal flatness condition – CFC hereafter):

$$\gamma _{ij}=\varphi ^4\widehat{\gamma }_{ij},$$
(6)
where $`\widehat{\gamma }_{ij}`$ is the flat metric ($`\widehat{\gamma }_{ij}=\delta _{ij}`$ in Cartesian coordinates). In general, the conformal factor $`\varphi `$ depends on the coordinates $`(t,r,\theta )`$. Therefore, at all times during a numerical simulation we assume that all off-diagonal components of the three-metric are zero, and the diagonal elements have the common factor $`\varphi ^4`$. Note that all metric quantities with a hat are defined with respect to the flat three-metric $`\widehat{\gamma }_{ij}`$.

Within this approximation the ADM equations reduce to a set of five coupled elliptic (Poisson-like) equations for the metric components,

$`\widehat{\mathrm{\Delta }}\varphi `$ $`=`$ $`2\pi \varphi ^5\left(\rho hW^2P+{\displaystyle \frac{K_{ij}K^{ij}}{16\pi }}\right),`$ (7)
$`\widehat{\mathrm{\Delta }}(\alpha \varphi )`$ $`=`$ $`2\pi \alpha \varphi ^5\left(\rho h(3W^22)+5P+{\displaystyle \frac{7K_{ij}K^{ij}}{16\pi }}\right),`$ (8)
$`\widehat{\mathrm{\Delta }}\beta ^i`$ $`=`$ $`16\pi \alpha \varphi ^4S^i+2\widehat{K}^{ij}\widehat{}_j\left({\displaystyle \frac{\alpha }{\varphi ^6}}\right){\displaystyle \frac{1}{3}}\widehat{}^i\widehat{}_k\beta ^k,`$ (9)
where $`\widehat{}`$ and $`\widehat{\mathrm{\Delta }}`$ are the flat space Nabla and Laplace operator, respectively. The transformation behavior between the extrinsic curvature defined on $`\gamma _{ij}`$ and $`\widehat{\gamma }_{ij}`$ is as follows:

$$K_{ij}=\varphi ^2\widehat{K}_{ij},K^{ij}=\varphi ^{10}\widehat{K}^{ij}.$$
(10)
The metric equations (79) couple to each other via their right hand sides, and in case of the equations for $`\beta ^i`$ via the operator $`\widehat{\mathrm{\Delta }}`$ acting on the vector $`\beta ^i`$. The equations are dominated by the source terms involving the hydrodynamic quantities $`\rho `$, $`P`$, and $`v^i`$, whereas the nonlinear coupling through the other, purely metric, source terms becomes only important for strong gravity. On each time slice the metric is solely determined by the instantaneous hydrodynamic state, i.e. the distribution of matter in space.

Applying the CFC to the traceless part of Eq. (2) yields the following relation between the conformal factor and the shift vector components:

$$_t\varphi =\frac{\varphi }{6}_k\beta ^k.$$
(11)
With this, Eq. (2) transforms into a definition for the extrinsic curvature components:

$$K_{ij}=\frac{1}{2\alpha }\left(_i\beta _j+_j\beta _i\frac{2}{3}\gamma _{ij}_k\beta ^k\right).$$
(12)
The above approximation is exact in spherical symmetry. Nevertheless we have used it for our simulations of rotational core collapse. The rationale behind this choice is the well-justified assumption that the energy emitted in form of gravitational waves is only about $`10^6`$ of the energy released during gravitational collapse in a Type II supernova event (Müller 1982; Mönchmeyer et al. 1991; Zwerger & Müller 1997). Furthermore, during core collapse deviations from spherical symmetry are modest even for rapidly rotating cores compared to those occuring in a neutron star merger (Wilson et al. 1996)). This fact is also reflected in the strength of the gravitational wave signal from core collapse, which is three orders of magnitude smaller than that resulting from the plunge phase of a compact binary. In Section 5.5 we will discuss the quality of the CFC in the context of rotational core collapse in more detail.

As a consequence of setting the off-diagonal elements of $`\gamma _{ij}`$ to zero, the degrees of freedom representing gravitational waves are removed from the spacetime. Therefore, gravitational wave emission is calculated in a post-processing step using the quadrupole formula. A more detailed description of the wave extraction methods is given in the accompanying paper (Dimmelmeier et al. 2002).

### 2.2 General relativistic hydrodynamic equations

The hydrodynamic evolution of a relativistic perfect fluid with rest-mass current $`J^\mu =\rho u^\mu `$ and energy-momentum tensor $`T^{\mu \nu }=\rho hu^\mu u^\nu +Pg^{\mu \nu }`$ in a (dynamic) spacetime $`g^{\mu \nu }`$ is determined by a system of local conservation equations, which read:

$$_\mu J^\mu =0,_\mu T^{\mu \nu }=0.$$
(13)
Here $`_\mu `$ is the covariant derivative and $`u^\mu `$ is the four-velocity of the fluid. Its three-velocity, as measured by an Eulerian observer at rest in a spacelike hypersurface $`\mathrm{\Sigma }_t`$ is

$$v^i=\frac{u^i}{\alpha u^t}+\frac{\beta ^i}{\alpha }.$$
(14)
Following the work of Banyuls et al. (1997) we now introduce the following set of conserved variables in terms of the primitive (physical) hydrodynamic variables $`(\rho ,v_i,ϵ)`$, with $`\rho `$ and $`ϵ`$ being the rest-mass density and specific internal energy density, respectively:

$`D`$ $`=`$ $`J^\mu n_\mu =\rho W,`$ (15)
$`S^i`$ $`=`$ $`_\nu ^iT^{\mu \nu }n_\mu =\rho hW^2v^i,`$ (16)
$`\tau `$ $`=`$ $`T^{\mu \nu }n_\mu n_\nu J^\mu n_\mu =\rho hW^2PD.`$ (17)
In the above expressions, $`n_\mu `$ is the unitary vector normal to the slice and $`_\nu ^i`$ denotes the projection operator. Furthermore $`W`$ is the Lorentz factor defined as $`W=\alpha u^t`$ and satisfying the relation $`W=1/\sqrt{1v_iv^i}`$. In addition, $`h=1+ϵ+P/\rho `$ is the specific enthalpy and $`P`$ is the pressure.

The local conservation laws, Eq. (13), can be written as a first-order, flux-conservative hyperbolic system of equations,

$$\frac{1}{\sqrt{g}}\left[\frac{\sqrt{\gamma }𝑼}{x^0}+\frac{\sqrt{g}𝑭^i}{x^i}\right]=𝑸.$$
(18)
with the state vector, flux vector and source vector given by

$`𝑼`$ $`=`$ $`[D,S_j,\tau ]`$ (19)
$`𝑭^i`$ $`=`$ $`[D(v^i{\displaystyle \frac{\beta ^i}{\alpha }}),S_j(v^i{\displaystyle \frac{\beta ^i}{\alpha }})+\delta _j^iP,`$ (20)
$`\tau (v^i{\displaystyle \frac{\beta ^i}{\alpha }})+Pv^i],`$
$`𝑸`$ $`=`$ $`[0,T^{\mu \nu }({\displaystyle \frac{g_{\nu j}}{x^\mu }}\Gamma _{\mu \nu }^\lambda g_{\lambda j}),`$ (21)
$`\alpha (T^{\mu 0}{\displaystyle \frac{\mathrm{ln}\alpha }{x^\mu }}T^{\mu \nu }\Gamma _{\mu \nu }^0)].`$
Here $`\sqrt{g}=\alpha \sqrt{\gamma }`$, with $`g=det(g_{\mu \nu })`$ and $`\gamma =det(\gamma _{ij})`$ being the determinant of the four-metric and three-metric, respectively. In addition, $`\Gamma _{\mu \nu }^\lambda `$ are the Christoffel symbols.

## 3 Initial models for rotational core collapse

### 3.1 Equation of state

For a rotating iron core before collapse the polytropic relation between the pressure $`P`$ and the rest mass density $`\rho `$,

$$P=K\rho ^\gamma ,$$
(22)
with $`\gamma =\gamma _{\mathrm{ini}}=4/3`$ and $`K4.897\times 10^{14}`$ (in cgs units) is a fair approximation of the density and pressure stratification (Zwerger & Müller 1997).

To model the physical processes which lead to the onset of collapse, we reduce the effective adiabatic index from $`\gamma _{\mathrm{ini}}`$ to $`\gamma _1`$ on the initial time slice. During the infall phase of core collapse, we also assume the matter of the core to obey a polytropic EoS, Eq. (22), which is consistent with the ideal gas EoS for a compressible inviscid fluid,

$$P=(\gamma 1)\rho ϵ,$$
(23)
where $`ϵ`$ is the specific internal energy.

During and after core bounce, when the shock starts to propagate out, the matter in the post-shock region is heated, i.e. kinetic energy is dissipated into internal energy. Therefore, in this stage we demand that the total pressure consists of a polytropic and a thermal part, $`P=P_\mathrm{p}+P_{\mathrm{th}}`$. We further assume that the polytropic index jumps from $`\gamma _1`$ to $`\gamma _2`$ for densities larger than nuclear matter density $`\rho _{\mathrm{nuc}}`$. This approximates the stiffening of the EoS. Requiring that the pressure and internal energy are continuous at the transition density $`\rho _{\mathrm{nuc}}`$, one finds for the polytropic internal energy

$$ϵ_\mathrm{p}=\{\begin{array}{cc}\frac{K}{\gamma _11}\rho ^{\gamma _11}\hfill & \\ \frac{K}{\gamma _21}\rho ^{\gamma _21}\rho _{\mathrm{nuc}}^{\gamma _1\gamma _2}+\frac{(\gamma _2\gamma _1)K}{(\gamma _21)(\gamma _11)}\rho _{\mathrm{nuc}}^{\gamma _11}\hfill & \end{array}$$
(24)
for $`\rho \rho _{\mathrm{nuc}}`$ and $`\rho >\rho _{\mathrm{nuc}}`$, respectively (see Zwerger (1995)). The thermal part of the pressure, which is due to shock heating, is given by $`P_{\mathrm{th}}=(\gamma _{\mathrm{th}}1)\rho ϵ_{\mathrm{th}}`$, where $`\gamma _{\mathrm{th}}=1.5`$ and $`ϵ_{\mathrm{th}}=ϵϵ_\mathrm{p}`$. This describes a mixture of relativistic ($`\gamma =4/3`$) and nonrelativistic ($`\gamma =5/3`$) components of an ideal fluid.

From these considerations one can construct an EoS for which both the total pressure $`P`$ and the individual contributions $`P_\mathrm{p}`$ and $`P_{\mathrm{th}}`$ are continuous at $`\rho _{\mathrm{nuc}}`$:

$`P`$ $`=`$ $`{\displaystyle \frac{\gamma \gamma _{\mathrm{th}}}{\gamma 1}}K\rho _{\mathrm{nuc}}^{\gamma _1\gamma }\rho ^\gamma {\displaystyle \frac{(\gamma _{\mathrm{th}}1)(\gamma \gamma _1)}{(\gamma _11)(\gamma _21)}}K\rho _{\mathrm{nuc}}^{\gamma _11}\rho `$ (25)
$`+(\gamma _{\mathrm{th}}1)\rho ϵ.`$
### 3.2 Rotating relativistic stars in equilibrium

Our initial models are obtained using Hachisu’s self-consistent field method (Komatsu et al. 1989a) which is a general relativistic method for solving the hydrostatic equilibrium equations for rotating polytropic matter distributions. The general metric to describe rotating axisymmetric relativistic matter configurations in equilibrium is

$`ds^2`$ $`=`$ $`e^{2\widehat{\nu }}dt^2+e^{2\widehat{\alpha }}(dr^2+r^2d\theta ^2)`$ (26)
$`+e^{2\widehat{\beta }}r^2\mathrm{sin}^2\theta (d\phi \omega dt)^2,`$
with the metric potentials $`\widehat{\nu }`$, $`\widehat{\alpha }`$, $`\widehat{\beta }`$ and $`\omega `$.

The only nonvanishing three-velocity component is $`v^3`$, the velocity of the fluid as measured by a ZAMO (zero angular momentum observer). Thus, the Lorentz factor is given by $`W=1/\sqrt{1v_3v^3}`$. With the definitions $`u^0=We^{\widehat{\nu }}`$, and $`u^3=\mathrm{\Omega }We^{\widehat{\nu }}`$, where $`\mathrm{\Omega }`$ is the angular velocity of the fluid as measured from infinity, the four-velocity reads:

$$u^\mu =We^{\widehat{\nu }}(1,0,0,\mathrm{\Omega }).$$
(27)
Thus, one gets for the $`\phi `$-component of the three-velocity:

$`v^3`$ $`=`$ $`{\displaystyle \frac{u^3}{u^0e^{\widehat{\nu }}}}{\displaystyle \frac{\omega }{e^{\widehat{\nu }}}}=e^{\widehat{\nu }}(\mathrm{\Omega }\omega ),`$ (28)
$`v_3`$ $`=`$ $`e^{2\widehat{\beta }\widehat{\nu }}r^2\mathrm{sin}^2\theta (\mathrm{\Omega }\omega ),`$ (29)
$`v_\phi `$ $``$ $`\sqrt{v_3v^3}=e^{\widehat{\beta }\widehat{\nu }}r\mathrm{sin}\theta (\mathrm{\Omega }\omega ).`$ (30)
Here $`v_\phi `$ is the proper rotation velocity with respect to the ZAMO. The specific angular momentum of the fluid $`ju^0u_3`$ is given by

$$j=W^2e^{2(\widehat{\beta }\widehat{\nu })}r^2\mathrm{sin}^2\theta (\mathrm{\Omega }\omega )=\frac{v_3v^3}{(1v_3v^3)(\mathrm{\Omega }\omega )}.$$
(31)
Using the Einstein equations one can derive the equation of hydrostatic equilibrium for an axisymmetric mass distribution rotating with the angular velocity $`\mathrm{\Omega }=\mathrm{\Omega }(r)`$:

$$\frac{P}{\rho h}+\widehat{\nu }\frac{v_\phi }{1v_\phi ^2}v_\phi +j\mathrm{\Omega }=0.$$
(32)
For barotropes where $`p`$ depends only on $`\rho `$ the integrability of Eq. (32) requires that $`j`$ is a function of $`\mathrm{\Omega }`$, only. Hence, the simplest choice for a rotation law is

$$j=A^2(\mathrm{\Omega }_\mathrm{c}\mathrm{\Omega }),$$
(33)
where $`\mathrm{\Omega }_\mathrm{c}`$ is the value of $`\mathrm{\Omega }`$ at the coordinate center and $`A`$ is a positive constant (Komatsu et al. 1989a). In the Newtonian limit the rotation law (33) reduces to

$$\frac{\mathrm{\Omega }}{\mathrm{\Omega }_\mathrm{c}}=\frac{A^2}{A^2+d^2}=\{\begin{array}{ccc}1\hfill & \text{for}A\mathrm{}\hfill & \\ \frac{A^2}{d^2}\hfill & \text{for}A0\hfill & \end{array}$$
(34)
where $`d=r\mathrm{sin}\theta `$ is the distance from the rotation axis. The upper case in Eq. (34) corresponds to a rigid rotator and the lower one to a configuration with constant specific angular momentum.

Using the rotation law specified in Eq. (33) one can integrate the equation of hydrostatic equilibrium (32) to obtain (Komatsu et al. 1989a)

$$\mathrm{ln}h+\widehat{\nu }+\frac{1}{2}\mathrm{ln}(1v_\phi ^2)\frac{1}{2}A^2(\mathrm{\Omega }\mathrm{\Omega }_\mathrm{c})^2=\mathrm{const},$$
(35)
The four-metric components $`\widehat{\alpha }`$, $`\gamma =\widehat{\beta }+\widehat{\nu }`$, $`\delta =\widehat{\nu }\widehat{\beta }`$ and $`\omega `$ are given by the four coupled partial differential equations (Komatsu et al. 1989a)

$`\mathrm{\Delta }\left(\delta e^{\gamma /2}\right)`$ $`=`$ $`S_\delta ,`$ (36)
$`\left(\mathrm{\Delta }+{\displaystyle \frac{1}{r}}_r{\displaystyle \frac{\mu }{r^2}}_\mu \right)\gamma e^{\gamma /2}`$ $`=`$ $`S_\gamma ,`$ (37)
$`\left(\mathrm{\Delta }+{\displaystyle \frac{2}{r}}_r{\displaystyle \frac{2\mu }{r^2}}_\mu \right)\omega e^{(\gamma 2\delta )/2}`$ $`=`$ $`S_\omega ,`$ (38)
$`_\mu \widehat{\alpha }`$ $`=`$ $`S_{\widehat{\alpha }},`$ (39)
where $`\mathrm{\Delta }`$ is the Laplacian and $`\mu \mathrm{cos}\theta `$. The source terms $`S_\delta `$, $`S_\gamma `$, $`S_\omega `$, and $`S_{\widehat{\alpha }}`$ depend in general on all four metric functions. In order to solve this system of equations the first three PDEs are converted into integral equations using Green functions which in turn are expanded into Legendre polynomials. In discretized form, these equations are iterated until convergence is obtained. During each iteration, the discretized equation (39) for $`\widehat{\alpha }`$ is integrated.

### 3.3 Initial models and collapse parameters

For a given adiabatic index $`\gamma _{\mathrm{rot}\mathrm{ini}}`$, the initial models are determined by three parameters: The central density $`\rho _{\mathrm{c}\mathrm{ini}}`$, the degree of differential rotation $`A`$ (see Eq. (33)), and the rotation rate $`\beta _{\mathrm{rot}\mathrm{ini}}`$, which is given by the ratio of rotational energy and the absolute value of the gravitational binding energy. The central density of an iron core is about $`10^{10}\mathrm{g}\mathrm{cm}^3`$. We adopt this value for $`\rho _{\mathrm{c}\mathrm{ini}}`$ in all initial models.

Due to the lack of consistent models of rotating iron cores from stellar evolution calculations, we follow Zwerger & Müller (1997) and treat $`A`$ and $`\beta _{\mathrm{rot}\mathrm{ini}}`$ as free parameters of the initial models. The only other free parameters are $`\gamma _1`$, the adiabatic index at subnuclear densities ($`\rho <\rho _{\mathrm{nuc}}`$), and $`\gamma _2`$, the adiabatic index at supranuclear densities ($`\rho \rho _{\mathrm{nuc}}`$). As in Zwerger & Müller (1997), we set $`\gamma _2=2.5`$ and $`\gamma _1`$ equal to 1.325, 1.320, 1.310, 1.300 or 1.280, respectively. In total, we have evolved 26 models (see Table 1). One of the models (A3B2G4<sub>soft</sub>) has also been run with $`\gamma _2=2.0`$ in order to test the influence of a softer supranuclear EoS on the collapse dynamics. The nuclear density is set to $`\rho _{\mathrm{nuc}}=2.0\times 10^{14}\mathrm{g}\mathrm{cm}^3`$ in all simulations.

### 3.4 Atmosphere treatment

One common problem of many hydrodynamic codes is their inability to handle regions where the density is zero (or very small compared to the typical densities of interest). Such regions are encountered, for example, when simulating pulsating stars on an Eulerian grid, or when simulating rotating stars on a spherical grid. In the first case a very low density “atmosphere” is required outside the stellar surface into which the star can expand during its pulsations. In the second case the stellar surface is flattened due to centrifugal forces, i.e. part of the computational grid extends outside the flattened surface and must be filled with some low density atmosphere. Thus, when computing initial models of rotating cores we introduce a low density atmosphere outside the stellar matter distribution (see Table 2). An atmosphere is assumed in all zones where the density $`\rho `$ is below a prescribed threshold value which is practically defined as a fraction $`f_{\mathrm{atm}\mathrm{thr}}`$ of the central density of the initial model, i.e. the density of the atmosphere $`\rho _{\mathrm{atm}}=f_{\mathrm{atm}}\rho _{\mathrm{c}\mathrm{ini}}`$. During a simulation the atmosphere is reset after each timestep. This ensures that it adapts to the time-dependent “shape” of the stellar surface.

For the values specified in Table 2 and an initial central density $`\rho _{\mathrm{c}\mathrm{ini}}=10^{10}\mathrm{g}\mathrm{cm}^3`$ the homogeneous atmosphere has a density $`\rho _{\mathrm{atm}}=1\mathrm{g}\mathrm{cm}^3`$. For all practical purposes this can be considered as a low density atmosphere, which has no influence on the dynamics of core collapse, the only tradeoff being a slight loss of angular momentum to the atmosphere from regions of the star which are close to the boundary (see Section 5.4). The EoS in the atmosphere is assumed to be polytropic, too.

## 4 Numerical implementation

### 4.1 Computational grid

We use spherical polar coordinates $`(t,r,\theta ,\phi )`$, and assume axial symmetry with respect to the rotation axis, i.e. none of the variables depends on the azimuthal coordinate $`\phi `$. However, motions in the $`\phi `$-direction are allowed, and the associated velocity $`v_3`$ and momentum $`S_3`$ are nonzero.

We further assume symmetry with respect to the equator at $`\theta =\pi /2`$. Therefore, $`\theta `$ spans the interval between 0 and $`\pi /2`$. The angular grid consisting of $`n_\theta `$ zones is equally spaced, with $`\theta _j`$ denoting the coordinates of the cell centers in angular direction ($`1jn_\theta `$). The $`n_r`$ radial zones (zone centers located at $`r_i`$ ) are logarithmically spaced. Each computational cell $`(r_i,\theta _j)`$ is bounded by interfaces at $`r_{i\frac{1}{2}}`$, $`r_{i+\frac{1}{2}}`$, $`\theta _{j\frac{1}{2}}`$, and $`\theta _{j+\frac{1}{2}}`$, respectively. The angular size of a zone is $`\mathrm{\Delta }\theta =\frac{1}{2}\pi /n_\theta `$, while the radial size is given by $`\mathrm{\Delta }r_i=r_{i+\frac{1}{2}}r_{i\frac{1}{2}}`$. The radial size of the three innermost cells is enlarged in order to ease the timestep restriction; for $`i>3`$, $`\mathrm{\Delta }r_{i+1}/\mathrm{\Delta }r_i=\mathrm{const}`$. Note that $`r_{\frac{1}{2}}=0`$ is the origin, and $`r_{n_r+\frac{1}{2}}=R`$ defines the outer boundary of the grid. The rotation axis coincides with $`\theta _{\frac{1}{2}}=0`$, and $`\theta _{n_\theta +\frac{1}{2}}=\frac{1}{2}\pi `$ denotes the equatorial plane. We use $`n_r=200`$ radial zones and $`n_\theta =30`$ angular zones, which corresponds to an angular resolution of $`3^{}`$. The radial size $`\mathrm{\Delta }r_i`$ of the innermost cells is about 500 m. Convergence tests show that this grid resolution is sufficient to resolve the important flow features of core collapse dynamics.

### 4.2 Boundary conditions

We impose symmetry conditions for both hydrodynamic and metric quantities at the center ($`r=0`$) and the rotation axis ($`\theta =0`$), and at the equatorial plane ($`\theta =\pi /2`$). Quantities are assumed to be either symmetric or antisymmetric across a boundary. We combine these two options in the following notation for a quantity $`q_{i,j}`$:

$`q_{\frac{1}{2},j}`$ $`=`$ $`\pm q_{\frac{1}{2},j},`$ (40)
$`q_{i,n_\theta +\frac{1}{2}}`$ $`=`$ $`\pm q_{i,n_\theta \frac{1}{2}},`$ (41)
$`q_{i,\frac{1}{2}}`$ $`=`$ $`\pm q_{i,\frac{1}{2}},`$ (42)
where

$$\pm =\{\begin{array}{c}+\text{symmetric},\hfill \\ \text{antisymmetric}.\hfill \end{array}$$
(43)
The symmetry conditions for the scalar quantities $`\rho `$, $`ϵ`$, $`\varphi `$ and $`\alpha `$ are straightforward. Like their Newtonian counterparts $`\rho `$, $`ϵ`$ and $`\mathrm{\Phi }`$ they have to be continuous across all boundaries. Thus, the scalar quantities are assumed symmetric (+) with respect to all boundaries. Contrary to Newtonian hydrodynamics, the symmetry conditions for the vector quantities are nontrivial, as they occur both in covariant and contravariant form.

By defining “physical” velocities (the same can be done for the shift vector components $`\beta ^i`$),

$`v_r`$ $``$ $`\sqrt{v_1v^1}=\varphi ^2v_1=\varphi ^2v^1,`$ (44)
$`v_\theta `$ $``$ $`\sqrt{v_2v^2}=\varphi ^2r^1v_2=\varphi ^2rv^2,`$ (45)
$`v_\phi `$ $``$ $`\sqrt{v_3v^3}=\varphi ^2r^1\mathrm{sin}^1\theta v_3=\varphi ^2r\mathrm{sin}\theta v^3,`$ (46)
which are bounded by the speed of light ($`0|v_r|,|v_\theta |,|v_\phi |1`$), we can specify the symmetry conditions from Newtonian hydrodynamics. Keeping in mind the symmetry behaviour of the geometrical terms in the three-metric components $`\gamma _{ij}`$ and $`\gamma ^{ij}`$, the symmetry conditions for the “physical velocity vector” and the shift vector read:

|  | $`v_r`$ | $`v_\theta `$ | $`v_\phi `$ | $`\beta ^1`$ | $`\beta ^2`$ | $`\beta ^3`$ |
| --- | --- | --- | --- | --- | --- | --- |
| center | $``$ | $``$ | $``$ | $``$ | $`+`$ | $`+`$ |
| pole | $`+`$ | $``$ | $``$ | $`+`$ | $``$ | $`+`$ |
| equator | $`+`$ | $``$ | $`+`$ | $`+`$ | $``$ | $`+`$ |

Here we demand that (i) all velocities have to vanish at the center (no mass flow across the center), that (ii) the meridional velocity $`v_\theta `$ is zero along the rotation axis and in the equatorial plane (no mass flow across these boundaries), and that (iii) the azimuthal velocity $`v_\phi `$ is zero along the rotation axis. Hence, these entries are set to ($``$), while all others, where the velocities are continuous, are set to (+). The shift vector components $`\beta ^i`$ are treated according to the contravariant three-velocity vector $`v^i`$, as the shift vector corresponds to a “coordinate velocity”.

At the outer radial boundary, the hydrodynamic quantities are extrapolated from the interior to the boundary zones. However, this procedure cannot be used for the metric quantities, as they are determined by elliptic equations which define a boundary value problem. We thus determine their boundary conditions by matching the interior metric to an exterior spacetime metric. We assume that the exterior spacetime is given by the Schwarzschild solution for a spherically symmetric vacuum spacetime. The quality of this approximation is sufficiently good for all practical purposes provided the spacetime at the surface of the star does not deviate too much from spherical symmetry. This is true for iron core initial data, and also holds for rotational core collapse, as the outer boundary of the radial grid is fixed (at the equatorial radius of the initial model) at about $`1500\mathrm{km}`$ and as the monopole term always dominates the gravitational field.

The Schwarzschild solution in isotropic coordinates is

$`ds^2`$ $`=`$ $`{\displaystyle \frac{\left(1\frac{M_{\mathrm{grav}}}{2r}\right)^2}{\left(1+\frac{M_{\mathrm{grav}}}{2r}\right)^2}}dt^2`$ (47)
$`+\left(1+{\displaystyle \frac{M_{\mathrm{grav}}}{2r}}\right)^4(dr^2+r^2d\mathrm{\Omega }^2).`$
For a vanishing shift vector this metric can be exactly matched to the conformally flat interior metric:

$$\varphi =1+\frac{M_{\mathrm{grav}}}{2r},\alpha =\frac{\left(1\frac{M_{\mathrm{grav}}}{2r}\right)}{\left(1+\frac{M_{\mathrm{grav}}}{2r}\right)},\beta ^i=0.$$
(48)
As the gravitational mass $`M_{\mathrm{grav}}`$ is conserved and as $`\beta ^i`$ is always close to zero at the outer boundary, this matching to the isotropic Schwarzschild metric provides an acceptable outer boundary condition for the metric coefficients $`\varphi `$, $`\alpha `$ and $`\beta ^i`$.

We stress that these boundary conditions are only an approximation to the exact vacuum spacetime solution, since no rotating spacetime can be described by a conformally flat metric. By the above matching we neglect effects like frame dragging ($`\beta ^30`$), which would in general be present in the vacuum spacetime region outside the computational domain.

### 4.3 The hydrodynamic solver

The hydrodynamic solver of our code performs the numerical integration of system (18) using a so-called Godunov-type (high-resolution shock-capturing – HRSC) scheme. Such schemes are specifically designed to solve nonlinear hyperbolic systems of conservation laws (see, e.g., Toro (1997) for definitions). In a Godunov-type method the knowledge of the characteristic structure of the equations is essential to design a solution procedure based upon either exact or approximate Riemann solvers. These solvers compute, at every cell-interface of the numerical grid, the solution of local Riemann problems. Therefore, they automatically guarantee the proper capturing of discontinuities which may arise naturally in the solution space of a nonlinear hyperbolic system.

The time update of system (18) from $`t^n`$ to $`t^{n+1}`$ is performed according to the following conservative algorithm:

$`𝑼_{i,j}^{n+1}=𝑼_{i,j}^n`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }r_i}}(\widehat{𝑭}_{i+1/2,j}^r\widehat{𝑭}_{i1/2,j}^r)`$ (49)
$``$ $`{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }\theta }}(\widehat{𝑭}_{i,j+1/2}^\theta \widehat{𝑭}_{i,j1/2}^\theta )`$
$`+`$ $`\mathrm{\Delta }t𝑸_{i,j}.`$
The numerical fluxes $`\widehat{𝑭}^r`$ and $`\widehat{𝑭}^\theta `$ are computed by means of Marquina’s approximate flux-formula (Donat & Marquina 1996). The fluxes are obtained independently for each direction. The time update of the state-vector $`𝑼`$ is done simultaneously in both spatial directions using a method of lines in combination with a second-order (in time) conservative Runge–Kutta scheme. Moreover, in order to set up a family of local Riemann problems at every cell-interface we use the piecewise parabolic method (PPM) of Colella & Woodward (1984) for the reconstruction of cell interface values, which provides third-order accuracy in space.

Time derivatives of all metric quantities ($`\varphi `$, $`\alpha `$, and $`\beta ^i`$) are required to compute the Christoffel symbols, and hence the source term $`𝑸`$. As the equations for the metric in the CFC approximation do not provide explicit expressions for the time derivatives of these quantities, we approximate them numerically by discretized derivatives using values from the two time slices $`\mathrm{\Sigma }_{t^n}`$ and $`\mathrm{\Sigma }_{t^{n1}}`$, e.g.

$$\left(\frac{\varphi }{t}\right)_{i,j}^n=\frac{\varphi _{i,j}^n\varphi _{i,j}^{n1}}{\mathrm{\Delta }t^{n1}}.$$
(50)
Finally we note that after the hydrodynamics update we need to recover the primitive variables $`(\rho ,v_i,ϵ)`$ from the conserved ones $`(D,S_i,\tau )`$. Since the relation between these two sets is only given implicitely, we have to resort to a Newton-Raphson iterative method. This procedure is discussed in Appendix A.

### 4.4 The metric solver

One of the standard methods to solve an elliptic linear system like Poisson’s equation on a numerical grid is to discretize the equation. Let the vector of unknowns be $`𝒖=(u_{i,j})`$, and the vector of sources be $`𝒇=(f_{i,j})`$, where $`(i,j)`$ labels the position $`(r_i,\theta _j)`$ on the grid. This discretization leads to a linear matrix equation with the right hand side being the source vector of the original equation:

$$𝑨𝒖=𝒇.$$
(51)
In general the matrix $`𝑨`$ is sparse, i.e. the number of nonzero elements is much smaller than the total number of elements. For solving sparse linear systems there exist a variety of efficient numerical methods and software packages for various computer architectures.

Although the CFC metric equations (79) do not form a linear system like Poisson’s equation, they can still be reduced to a linear system. Our strategy is as follows: We first define a vector of unknowns, whose components are the five metric quantities:

$$\widehat{u}=(u^k)=(\varphi ,\alpha \varphi ,\beta ^1,\beta ^2,\beta ^3).$$
(52)
Then the metric equations can be written as

$$\widehat{f}(\widehat{u})=0,$$
(53)
with $`\widehat{f}`$ denoting the vector of the metric equations for $`\widehat{u}`$.

The explicit form of the metric equations (53) is given in Appendix B. We discretize these equations on the $`(r,\theta )`$-grid, denoting the metric components $`u^k`$ ($`k=1,\mathrm{},5`$) at the cell-center $`(r_i,\theta _j)`$ by $`u_{i,j}^k`$. Their values at cell-interfaces (needed for the hydrodynamic solver) are obtained by interpolation. The radial and angular derivatives are approximated by central finite differences. A 9-point stencil is required to formulate the discretized equations. Hence, the system of metric equations, when discretized, gives rise to a nonlinear system of dimension $`5\times n_r\times n_\theta `$,

$$\widehat{𝒇}(\widehat{𝒖})=0,$$
(54)
with the vector of discretized equations $`\widehat{𝒇}=\widehat{f}_{i,j}=f_{i,j}^k`$ for the discretized unknowns $`\widehat{𝒖}=\widehat{u}_{l,m}=u_{l,m}^n`$. For this system we have to find the roots $`u_{i,j}^l`$.

We use a multi-dimensional Newton–Raphson iteration method to solve the nonlinear system, i.e. the nonlinear system of dimension $`5\times n_r\times n_\theta `$ is reduced to a linear one of the same dimension for the Jacobi matrix of $`\widehat{𝒇}`$. This linear matrix equation has to be solved once in each iteration step. For our set of equations the Jacobi matrix of the nonlinear system is in general both sparse and diagonally dominated, which reduces the complexity of the linear problem significantly.

The performance of our axisymmetric general relativistic code is mostly determined by the computation time spent in the metric solver, and thus by the computational efficiency of the linear solver. Therefore, the linear problem has to be solved as efficiently as possible. After extensive investigation we have decided to use a block tridiagonal sweeping method (Potter 1973). As the matrix which defines the linear problem has nine bands of blocks rather than three (due to the 9-point stencil; see above), we combine the left three, the middle three, and the right three bands, collecting them in three bands of submatrices – left, diagonal and right. Although these submatrices are sparse, the block tridiagonal sweeping algorithm leads to dense submatrices which have to be inverted. This is achieved using a standard $`LU`$ decomposition scheme for dense matrices as implemented in the LAPACK numerical library.

As far as computational performance is concerned, the block tridiagonal sweeping method exhibits an $`n_r\times n_\theta ^2`$ behavior, and is thus superior to both the Gauss–Jordan elimination and the bi-conjugate gradient stabilized methods which we have also tried. In particular, as it is a direct solver it does not suffer from convergence problems.

The tolerance measure we use to control convergence of the Newton–Raphson method is the maximum increment of the solution vector on the whole computational grid:

$$\mathrm{\Delta }\widehat{u}_{\mathrm{max}}^s=\mathrm{max}(\mathrm{\Delta }\widehat{𝒖}^s)=\mathrm{max}(\mathrm{\Delta }u_{i,j}^{ks}).$$
(55)
We are able to meet this convergence criterion for a tolerance level of $`\mathrm{\Delta }\widehat{u}_{\mathrm{max}}^s10^{15}`$ in all our computations.

In Fig. 1 we show how $`\mathrm{\Delta }\widehat{u}_{\mathrm{max}}^s`$ decreases with the number of iterations $`s`$ in the metric computation of a typical core collapse model around the time of maximum density. The convergence measures exhibits an inverse exponential behavior until saturation is reached. This occurs when the increment of $`\widehat{𝒖}^s`$ approaches machine precision.

In addition, we monitor the residual of the function vector $`r_{\widehat{𝒇}}=\widehat{𝒇}(\widehat{𝒖}^s))`$ on the entire computational grid. If $`\widehat{𝒖}^s`$ is the exact solution of the discretized equations, then $`r_{\widehat{𝒇}}=0`$. In practice we obtain values for $`r_{\widehat{𝒇}}`$ which are as small as roughly $`10^9`$ of the individual terms of $`\widehat{𝒇}(\widehat{𝒖}^s))`$. This measure ensures convergence to the correct solution of the metric equations.

### 4.5 Metric extrapolation

We have performed comparisons between the execution time for one metric solution (for different linear solver methods) and one hydrodynamic time step. The comparison shows that, particularly for grid sizes appropriate for our collapse simulations, the computation is completely dominated by the metric solver: One metric computation can be as time-consuming as about 100 hydrodynamic time steps!

On the other hand, as long as the metric quantities do not change too rapidly with time, it is a fair approximation to solve for the metric not at every time step. We introduce the metric resolution parameter $`\mathrm{\Delta }n_\mathrm{m}`$, such that $`\mathrm{mod}(n,\mathrm{\Delta }n_\mathrm{m})=0`$. Then $`\mathrm{\Delta }n_\mathrm{m}`$ expresses the number of “hydrodynamic” time steps between two “metric” time slices where the metric is calculated.

During a simulation, the metric resolution parameter should vary with time in a way which adapts to the specific phase of the collapse scenario. We split the core collapse evolution into three different phases: The pre-bounce phase lasts from the start of the evolution until the maximum density on the grid $`\rho _{\mathrm{max}}`$ reaches a threshold value, which is a specified fraction of the nuclear density: $`\rho _{\mathrm{thr}}f_{\mathrm{ref}}\rho _{\mathrm{nuc}}`$. During that phase the metric is calculated every $`n_\mathrm{m}^{(1)}`$ “hydrodynamic” time step. The subsequent bounce phase lasts for a time $`\mathrm{\Delta }t_{\mathrm{ref}}`$, and the metric is solved on every $`n_\mathrm{m}^{(2)}`$th time slice. In the subsequent post-bounce phase, the metric resolution parameter is set to $`n_\mathrm{m}^{(3)}`$.

In models showing multiple bounces, i.e. where the maximum density has multiple distinct peaks in time often separated by several 10 ms, the evolution is split into a series of consecutive bounce phases interspersed with phases where the density is below the threshold marked by $`\rho _{\mathrm{thr}}`$. During the bounce phases the metric resolution is given by $`n_\mathrm{m}^{(2)}`$, and before the first bounce and in between bounces by $`n_\mathrm{m}^{(1)}`$.

The collapse phases introduced above are sketched in Fig. 2, while the actual values for $`f_{\mathrm{ref}}`$ and $`\mathrm{\Delta }t_{\mathrm{ref}}`$ used in our simulations are summarized for the two different collapse types in Table 3.

In order to approximate the actual evolution of the metric on those “hydrodynamic” time slices where it is not calculated, the state of the metric at old “metric” time slices can be used to extrapolate the metric quantities forward in time. Fig. 3 shows the approximation of the actual metric calculated at every time slice by a metric which is calculated at every 10th time step, and is kept constant or is linearly or parabolically extrapolated in between. Tests with collapse models demonstrate that our choice of values for the metric resolution parameters is appropriate to obtain the desired accuracy. Due to the superior accuracy of the 3rd order parabolic metric extrapolation scheme we use this method in all our simulations.

## 5 Code tests

We present now the tests we have used to calibrate our code. It is worth pointing out that the present code also includes the possibility of performing simulations of rotational core collapse employing Newtonian gravity. Many routines and numerical algorithms are common to both the relativistic and the Newtonian computer code, and by comparing with previous Newtonian simulations we have also checked the correct implementation of those routines. For this purpose we have run some of the models of the comprehensive sample of Zwerger & Müller (1997) finding excellent agreement for both the dynamics of the collapse and the gravitational waveforms. Therefore, whenever we present Newtonian results below, it is assumed that those have been computed with our current code. We focus hence on checking our code in the relativistic regime.

### 5.1 Relativistic shock tube tests

Shock tubes can efficiently test the ability of a code to handle shock fronts. In a shock tube test, initial data are given for a fluid system which has two states of constant pressure, density and zero velocity, initially separated by a boundary. The time-evolution of these initial data yields a combination of constant states separated by shocks, contact discontinuities, or rarefaction waves. These tests allow to compare the numerical results against the analytic solution (Martí & Müller 1994) in a straightforward way.

We have simulated the two relativistic shock tube problems analyzed by Martí & Müller (1996). The initially constant states are given in Table 4. Both cases are computed in a flat one-dimensional Minkowski spacetime using Cartesian coordinates in the interval $`x[0,1]`$ spanned by 500 equidistant zones. The EoS is assumed to be that of an ideal fluid with adiabatic index $`\gamma =5/3`$, and the initial discontinuity is placed at $`x=0.5`$. The results for the most stringent problem 2 (a relativistic blast wave) are displayed in Fig. 4, after an evolution time of $`t_{\mathrm{final}}=0.25`$. The profiles of the density $`\rho `$, the pressure $`P`$, and the velocity $`v_x`$ accurately reproduce the analytic solution.

### 5.2 Rotating neutron stars

An important test that any axisymmetric hydrodynamics code should pass is the computation of the evolution of an equilibrium model of a rotating stellar configuration. If the initial model is computed accurately, and if the evolution scheme is implemented properly, the matter and metric profiles of the evolved model should at most oscillate slightly around their initial equilibrium states, i.e. the model should remain in equilibrium. In particular, for a rotating model, the rotation profile has to be preserved for as many rotation periods as possible.

In order to enhance the relativistic effects we have used a rapidly and uniformly rotating neutron star model with a polytropic matter distribution (see Table 5). This test allows for an independent check of the metric and the hydrodynamics parts of the code. If desired, the code can solve only the fluid equations while keeping the metric fixed (Cowling approximation), or evolve both the hydrodynamics and the spacetime in a fully coupled evolution. In both cases we observe that the neutron star remains in equilibrium to high accuracy for several 10 milliseconds, which corresponds to 10 rotation periods. In this test run, the metric equations are solved every 200th hydrodynamic time step, which corresponds to a time interval between two metric time slices of $`0.06\mathrm{ms}`$. This interval is small compared to both the neutron star’s rotation period and the frequencies of its strongest radial oscillations.

In Fig. 5 we show the time evolution of the density $`\rho `$, the lapse $`\alpha `$ and the radial velocity $`v_r=\sqrt{\gamma ^{11}}v_1`$ at a radius of 8.77 km in the equatorial plane for a fully coupled evolution. The small amplitude oscillations are triggered by numerical discretization errors. The fact that these oscillations are hardly damped reflects the low numerical viscosity of the HRSC scheme. By performing a Fourier transform of these evolved data one can obtain the frequencies of the different modes of pulsation to high accuracy (Font et al. 2000, 2001a). Comparisons performed with the results in the Cowling approximation of Font et al. (2001a) yield a very good agreement for the mode frequencies.

As demonstrated in Fig. 5, the code can keep the equilibrium initial data to a high degree. The small secular drift observed in $`\rho `$ and $`\alpha `$ is an artifact of the numerical scheme and has also been observed by Font et al. (2000). This drift can largely be minimized using a polytropic zero-temperature EoS, which is a fair assumption for studying small amplitude pulsations of polytropic stars. The effect of different EoS on the drift is discussed in Font et al. (2001b).

In Fig. 6 we plot the radial profiles of the density $`\rho `$ and the rotation velocity $`v_\phi =\sqrt{\gamma ^{33}}v_3`$ along the equator for the same neutron star model. Even after 5 rotation periods, the rotation profile is very close to its initial distribution. Only at the stellar surface, the angular velocity slightly deviates from its initial shape due to the interaction with the atmosphere. As shown in Font et al. (2000) the use of high-order schemes such as PPM is essential in preserving the rotational profile.

Since the neutron star model specified in Table 5 is rapidly rotating and thus deviates strongly from spherical symmetry, its evolution provides a convincing way to test the ability of the CFC to yield a good approximation of the exact spacetime (see also Section 5.5). The initial model has been calculated using the exact metric, Eq. (26), without assuming conformal flatness. The initial data are then matched to the CFC spacetime. The fact that the evolution code maintains stability for many rotation periods already proves the validity of the CFC for rotating neutron star spacetimes. This observation is supported by Fig. 7 which shows how well the conformal factor $`\varphi ^{\mathrm{CFC}}`$ of the CF metric approximates the corresponding metric component $`\varphi ^{\mathrm{ex}}`$ of the exact metric. The excellent approximation of the azimuthal component of the exact shift vector $`\beta ^{3\mathrm{ex}}`$ by the CFC shift vector component $`\beta ^{3\mathrm{CFC}}`$ is shown in the lower panels. Our results are in good agreement with those of Cook et al. (1996).

### 5.3 Spherical core collapse

The problem of relativistic core collapse for an ideal fluid in spherical symmetry was first considered by May & White (1966). Their code used a Lagrangian formulation, where the coordinates label mass shells rather than spatial positions and thus comove with the fluid. Artificial viscosity terms were included to damp out spurious numerical oscillations behind shock fronts. In this Section we compare our Eulerian code against a Lagrangian finite difference code (Dimmelmeier 1998) for a spherically symmetric core collapse.

For this test we have set up a nonrotating equilibrium star as initial data on a two-dimensional $`(r,\theta )`$-grid, with the same EoS and central density as in the (rotating) initial models listed in Table 1. During the evolution, we do not use the hybrid EoS (25), but a simple ideal fluid EoS given by Eq. (23).

The effective stiffening of the EoS at supranuclear matter densities is modeled by a variable adiabatic index given by the following relation (van Riper 1979; Romero et al. 1996):

$$\gamma =\{\begin{array}{cc}\gamma _1\hfill & \text{for }\rho <\rho _{\mathrm{nuc}},\hfill \\ \gamma _1+\frac{1}{2}\mathrm{log}\frac{\rho }{\rho _{\mathrm{nuc}}}\hfill & \text{for }\rho \rho _{\mathrm{nuc}},\hfill \end{array}$$
(56)
with $`\gamma _1=1.320`$. Besides comparing our results with those of an independent spherically symmetric code, this test also allows us to assess the capability of the code to maintain spherical symmetry throughout the collapse.

The differences in the ability of both codes to handle shocks are demonstrated in Fig. 8. In the upper panel we plot the evolution of the central density around the time of bounce. We note that due to the relatively soft supranuclear EoS used for this run the core dives deeply into the potential well reaching densities as high as $`5\rho _{\mathrm{nuc}}`$. Even during this epoch, where the highest densities are reached, the evolution computed with our code follows very accurately the corresponding one obtained with the May–White code. However, due to its superior effective radial resolution, the May–White code resolves the drop in the central density after the maximum peak better, and also produces a smoother density evolution. Up to the formation of the shock front, both codes yield equally good results for the radial infall velocity profiles. Such profiles are plotted in the lower panel of Fig. 8 at $`t=60\mathrm{ms}`$ (before bounce) and $`t=75\mathrm{ms}`$ (after bounce), respectively. With the HRSC scheme of our code the shock front is steeper and confined to fewer grid zones than with the artificial viscosity scheme implemented in the May–White code. Additionally, the velocity behind the shock is higher and the shock propagates faster, which is a result of much less numerical dissipation. We note that if the artificial viscosity is reduced in the May–White code to cure these negative effects, unphysical spurious oscillations grow behind the shock front (see inset in the lower panel of Fig. 8).

### 5.4 Integral quantities – Conservation of rest-mass and angular momentum

In the 3 + 1 formalism of general relativity, the definitions for the (baryon) rest mass $`M_{\mathrm{rest}}`$, the proper mass $`M_{\mathrm{proper}}`$, the gravitational mass $`M_{\mathrm{grav}}`$, the angular momentum $`J_{\mathrm{rot}}`$, and the rotational mass (energy) $`M_{\mathrm{rot}}`$ are (Friedman et al. 1986):

$`M_{\mathrm{rest}}`$ $`=`$ $`{\displaystyle \rho u_\mu n^\mu dV},`$ (57)
$`M_{\mathrm{proper}}`$ $`=`$ $`{\displaystyle \rho (1+ϵ)u_\mu n^\mu dV},`$ (58)
$`M_{\mathrm{grav}}`$ $`=`$ $`{\displaystyle \left(2T_{\mu \nu }g_{\mu \nu }T\right)t^\mu n^\nu dV},`$ (59)
$`J_{\mathrm{rot}}`$ $`=`$ $`{\displaystyle T_{\mu \nu }s^\mu n^\nu dV},`$ (60)
$`M_{\mathrm{rot}}`$ $`=`$ $`{\displaystyle T_{\mu \nu }s^\mu n^\nu \frac{s_\lambda u^\lambda }{2n_\kappa u^\kappa }dV},`$ (61)
where $`t^\mu `$ and $`s^\mu `$ are a time-like and a space-like Killing vector, respectively. The rotation axis is perpendicular to $`s^\mu `$. For the CF metric in axisymmetry, where $`s^\mu =\widehat{e}^3`$ (i.e. the unit vector in the $`\phi `$-direction), the above integral quantities have the following form:

$`M_{\mathrm{rest}}`$ $`=`$ $`2\pi {\displaystyle r^2\mathrm{sin}\theta \varphi ^6\rho W\mathrm{d}r\mathrm{d}\theta },`$ (62)
$`M_{\mathrm{proper}}`$ $`=`$ $`2\pi {\displaystyle r^2\mathrm{sin}\theta \varphi ^6\rho (1+ϵ)Wdrd\theta },`$ (63)
$`M_{\mathrm{grav}}`$ $`=`$ $`2\pi {\displaystyle }r^2\mathrm{sin}\theta \varphi ^6(\alpha (2\rho hW^2+2P\rho h)`$ (64)
$`2\rho hW^2v_i\beta ^i)\mathrm{d}r\mathrm{d}\theta ,`$
$`J_{\mathrm{rot}}`$ $`=`$ $`2\pi {\displaystyle r^2\mathrm{sin}\theta \varphi ^6\rho hW^2v_3\mathrm{d}r\mathrm{d}\theta },`$ (65)
$`M_{\mathrm{rot}}`$ $`=`$ $`2\pi {\displaystyle r^2\mathrm{sin}\theta \varphi ^6\rho hW^2\frac{v_3(\alpha v^3\beta ^3)}{2}\mathrm{d}r\mathrm{d}\theta }.`$ (66)
Note that for $`c=1`$, the energies have the same units as the corresponding masses. The potential (mass) energy is given by $`M=M_{\mathrm{grav}}M_{\mathrm{proper}}M_{\mathrm{rot}}`$.

This formulation of the above integral quantities is an extension of the masses and the angular momentum defined in Komatsu et al. (1989a) for nonzero radial and angular velocity and shift vector components, applied to a conformally flat metric. It is obvious that the integrand of the rest mass and the angular momentum can be written as $`\sqrt{\gamma }D`$ and $`\sqrt{\gamma }S_3`$, respectively. These are the first and fourth component of Eq. (18) and correspond to the continuity equation and the angular momentum conservation equation, respectively. We have used these equations to check the conservation properties of our numerical code.

First, we have performed tests where we have evolved the hydrodynamic state vector in a fixed spacetime metric. Additionally, the fluxes at the outer boundary are explicitely set to zero and no artifical atmosphere is used. In this situation, we obtain conservation of the rest mass and the angular momentum of

$$dM_{\mathrm{rest}}\left|\frac{M_{\mathrm{rest}}}{M_{\mathrm{rest}\mathrm{\hspace{0.17em}0}}}1\right|10^{15},$$
(67)
$$dJ_{\mathrm{rot}}\left|\frac{J_{\mathrm{rot}}}{J_{\mathrm{rot}\mathrm{\hspace{0.17em}0}}}1\right|10^{15},$$
(68)
i.e. up to machine precision. Here $`M_{\mathrm{rest}\mathrm{\hspace{0.17em}0}}`$ and $`J_{\mathrm{rot}\mathrm{\hspace{0.17em}0}}`$ are the initial total rest mass and total angular momentum, respectively. In Newtonian gravity, the same level of conservation is reached.

In a dynamic core collapse simulation the above restrictions have to be abandoned. The nonstationarity of the metric and (in the case of rotation) the artificial atmosphere prevent a highly accurate conservation of the integral quantities. In Fig. 9 we plot the behavior of the mass and the angular momentum in the rotational core collapse model A3B2G4 (see Table 1). The dashed lines correspond to a Newtonian simulation and the solid ones are the relativistic counterparts. In Newtonian gravity the only source for mass and angular momentum “accretion” is the artificial atmosphere. This becomes visible as a small monotonous increase in $`M_{\mathrm{rest}}`$ (middle panel) and $`J_{\mathrm{rot}}`$ (bottom panel). Until bounce, indicated by the vertical dotted line, the mass increase in relativistic gravity exhibits the same behavior. However, due to the nonzero source term in the relativistic angular momentum equation, $`dJ_{\mathrm{rot}}`$ is larger than its Newtonian counterpart already during the infall phase.

Around the time of bounce another effect which prohibits exact numerical conservation of total rest mass and angular momentum becomes important. The conservation equations (18) are formulated for $`\sqrt{\gamma }𝑼`$. But in our numerical scheme we evolve only the state vector $`𝑼`$ and, for numerical reasons, we treat the geometric volume term $`\sqrt{\gamma }`$ as a source term. Therefore, for any strong temporal variation of the metric, e.g. during bounce, $`dM_{\mathrm{rest}}`$ and $`dJ_{\mathrm{rot}}`$ vary considerably, too. This is clearly visible in Fig. 9. When the central density has settled down to an almost constant value, and the core has reached a new equilibrium configuration, the rates of increase of $`M_{\mathrm{rest}}`$ and $`J_{\mathrm{rot}}`$ are again similar in relativistic and in Newtonian gravity.

The effects of the artificial atmosphere on the conservation of total mass and angular momentum can be infered from Fig. 10. In the highly nonspherical rotational core collapse model A4B2G2, the shock propagates outward at different speeds in the polar and the equatorial direction. When the shock reaches the surface of the star at the pole ($`t85\mathrm{ms}`$), a considerable amount of mass, but not much angular momentum is lost to the atmosphere. The equator is reached by the shock at a later time, since the shock propagates slower in this direction and the equatorial radius is larger. Around that time at $`t92.5\mathrm{ms}`$ the mass loss curve steepens and significantly more angular momentum is lost, as the relative contributions of fluid elements to the total angular momentum is higher near the equatorial plane than at the pole. Fig. 10 also shows the variation of $`dM_{\mathrm{rest}}`$ and $`dJ_{\mathrm{rot}}`$ due to rapid changes of the volume element $`\sqrt{\gamma }`$ around peak maximum density.

### 5.5 Quality of the conformally flat condition during rotational core collapse

As mentioned in Section 5.2, the accuracy of the CFC approximation has been considered by Cook et al. (1996), who find remarkably good results for rapidly rotating relativistic stars in equilibrium, typical errors for different variables being smaller than 5%. However, the quality of the CFC approximation degrades in the case of extremely relativistic nonspherical configurations, e.g. for rigidly rotating infinitesimally thin disks of dust, which have been investigated by Kley & Schäfer (1999). These authors also point out that general relativity and the truncated CFC approach are identical up to the first post-Newtonian order. Although the above comparisons have been made for stationary configurations, we want to emphasize that the CFC approximation does not make any explicit assumptions about the spacetime being static or dynamic. Due to the elliptic character of the metric equations (79), the spacetime metric is solely determined by the instantaneous hydrodynamic state. Therefore the deviation of the CF metric from the exact spacetime metric at a particular time will only be due to the current deviation of the matter distribution from spherical symmetry.

In this Section we analyze the validity of the CFC in simulations of rotational core collapse. There are several ways to quantify the violations of the exact metric equations. Firstly, one could use the CF metric obtained by solving Eqs. (79) to compute the Einstein tensor $`G_{\mu \nu }`$ and to check the Einstein equations $`G_{\mu \nu }=8\pi T_{\mu \nu }`$. However, this is a difficult approach due to the huge number of terms involved.

On the other hand, the ADM metric equations (25) are equivalent to the Einstein equations. As we only use a subset of the ADM equations to calculate the CF metric, we can use the remaining ones to test the quality of the CFC approximation. Therefore, we re-insert the CF metric components into the individual ADM metric equations and compare the left and right hand sides of the respective equations. We note that such a procedure has also been used by Gourgoulhon et al. (2002) in validating initial data for binary black holes. Comparisons with Damour’s analytic post-Newtonian expansion show satisfactory agreement in this case up to third post-Newtonian order (Gourgoulhon 2001).

If the inserted metric is an exact solution of the Einstein equations, the numerical inequality between the left and right hand sides of the ADM equations has to converge to zero with increasing resolution. For an approximate metric like the CF metric there will remain a nonzero residual

$$r_q\left|\frac{\mathrm{lhs}q}{\mathrm{rhs}q}1\right|>0$$
(69)
in the ADM equations for at least some metric quantities $`q`$ irrespective of numerical resolution. Here $`\mathrm{lhs}q`$ and $`\mathrm{rhs}q`$ stand for the left and right hand side of the corresponding ADM metric equation, respectively.

We want to point out that the metric equations (79) are analytically equivalent to the ADM constraint equations. Any solution of the CF metric equations will also be a solution of the ADM constraint equations, at least up to the precision of the numerical scheme.

The evolution equations for the three-metric $`\gamma _{ij}`$, Eq. (2), are also identically fulfilled for the off-diagonal elements of $`\gamma _{ij}`$, because for a CF metric they transform into the definition of the extrinsic curvature, Eq. (12). As the three-metric $`\gamma _{ij}`$ is conformally flat, of the six evolution equations (2) there remains only the evolution equation for the conformal factor $`\varphi `$, Eq. (11) which must be monitored. We have found that the residual of this equation,

$$r_\varphi =\left|\frac{\mathrm{lhs}\varphi }{\mathrm{rhs}\varphi }1\right|=\left|\frac{_t\varphi }{\frac{1}{6}\varphi _k\beta ^k}1\right|,$$
(70)
is very close to zero in various test situations including rotational core collapse. For model A3B4G3 the radial profile of $`r_\varphi `$ is plotted in Fig. 11 shortly before maximum density is reached. Except where a pole develops due to a zero denominator in Eq. (70), $`r_\varphi `$ is everywhere less than 0.2%, particularly in the dense interior. Hence, we infer that the ADM evolution equations for the diagonal elements of $`\gamma _{ij}`$ are solved to high degree of accuracy by the CF metric as well.

The ADM evolution equations for the extrinsic curvature provide another way of checking the accuracy of the CFC approximation. As a consequence of the maximal slicing condition, $`K=0`$, the sum of the left and right hand sides of the equation for the diagonal elements of the extrinsic curvature are identical up to numerical discretization errors at all grid points in any evolution. Furthermore, in spherical symmetry, the left and right hand sides of all six evolution equations for $`K_{ij}`$ agree, because the CF metric is an exact solution of the spacetime in this case.

A violation of the evolution equations for the diagonal elements of $`K_{ij}`$, Eq. (3), occurs in all highly nonspherical and relativistic spacetimes, e.g. in the core collapse model A3B4G3 (Fig. 12). The agreement between both sides of the equations is good at low densities at the beginning of the infall phase. However, when high densities are encountered in the center during and after bounce the misalignment becomes as prominent as for the rotating equilibrium neutron star model of Section 5.2 (bottom panel of Fig. 12). On the other hand, the evolution equations for off-diagonal components are solved to high accuracy. This is demonstrated by he very good agreement of the left and right hand sides of the equation for $`K_{13}`$ even in the center of the core (upper panel of Fig. 12).

These observations raise the question why the left and right hand sides of some of the evolution equations for the extrinsic curvature disagree considerably for nonspherical matter distributions. The reason becomes obvious when analyzing the structure of those equations. In comprehensive tests computing many dynamic collapse models we have found that the terms in Eq. (3 containing extrinsic containing extrinsic curvature components are negligibly small compared to the other three terms $`_i_j\alpha `$, $`\alpha R_{ij}`$ and $`8\pi \alpha (S_{ij}\frac{1}{2}\gamma _{ij}(S\rho _\mathrm{H}))`$ everywhere on the grid and during all epochs. The relevant terms depend only on the metric quantities $`\alpha `$ (which is fixed by a gauge condition) and $`\varphi `$, but not on the shift vector components $`\beta ^i`$. This points towards the inability of the CF three-metric $`\gamma _{ij}`$, which has only one nonzero component $`\varphi `$, to exactly represent the physical spacetime. On the other hand, we note that the quality of the shift vector components, which are calculated directly according to the ADM momentum constraint equations (5) (which are equivalent to the metric equation (9)), is not measured by Eq. (3).

Therefore, the misalignment of the evolution equation for the extrinsic curvature is not a result of strong deviations of the shift vector components, but is instead due to the restriction imposed onto the three-metric $`\gamma _{ij}`$ by the CFC approximation. This conjecture is additionally supported by the fact that, as demonstrated in Fig. 7 in Section 5.2, even for high density, rapidly rotating neutron stars in equilibrium the value of the shift vector component $`\beta ^3`$ calculated by the CFC equation (9) remains close to the exact value.

This also explains the perfect matching of the left and right hand sides of Eq. (3) for $`K_{13}`$ even for rapidly rotating and collapsing models. For this component the problematic terms $`_i_j\alpha `$ and $`\alpha R_{ij}`$ vanish. However, for $`K_{12}`$ and $`K_{23}`$ the deviations are more pronounced, since the polar velocity $`v_2`$ is very close to zero in our collapse models. As a result the left and right hand sides of these evolution equations are dominated by truncation errors.

## 6 Conclusions

We have presented a detailed description of a new axisymmetric general relativistic code for rotational core collapse which solves the coupled system of Einstein equations and fluid equations. Those equations are integrated adopting the ADM 3 + 1 formalism to foliate the spacetime into a set of non-intersecting spacelike hypersurfaces. In order to simplify the metric equations we approximate the general metric using the conformal flatness condition (CFC) introduced by Wilson et al. (1996), whereby the spatial components of the metric are set equal to the flat three-metric times a conformal factor which depends on the coordinates. The CFC implies a considerable simplification of the ADM 3 + 1 metric equations which reduce to a set of five coupled non-linear elliptic equations for the metric components. We have carefully analyzed the applicability and quality of the CFC, showing that this approximation is appropriate for simulations of rotational core collapse.

We have also demonstrated that our general relativistic hydrodynamics solver accurately passes a set of difficult test problems including special relativistic shock tubes, the preservation of the rotation profile and of the equilibrium of rapidly and differentially rotating relativistic polytropes, and the spherical collapse of a relativistic polytrope. The solver, which is based on a high-resolution shock capturing scheme with an approximate Riemann solver, also conserves accurately rest-mass and angular momentum in dynamic spacetimes.

The code described in this work was developed with the aim to simulate the general relativistic axisymmetric collapse of rotating stellar cores. To this end we have computed and presented here a large set of relativistic rotating precollapse equilibrium models with different amounts and distributions of angular momentum. Using these initial models we have performed a series of collapse simulations and calculated the resulting gravitational wave signal. The results of this first relativistic study of rotational core collapse and of its gravitational wave signature are discussed in an accompanying paper (Dimmelmeier et al. 2002).

###### Acknowledgements.

We are grateful to Y. Eriguchi for providing us with a code for constructing the intial models. It is a pleasure to acknowledge J. M$`^{\underset{¯}{\text{a}}}`$ Ibáñez and G. Schäfer for very useful comments and suggestions. We also want to thank M. Rampp for insightful discussions about the sweeping method for sparse linear problems. J.A.F. acknowledges support from a Marie Curie fellowship by the European Union (MCFI-2001-00032) and from the Spanish Ministerio de Ciencia y Tecnología (AYA2001-3490-C02-01). All computations were performed on the NEC SX-5/3C supercomputer at the Rechenzentrum Garching.

## Appendix A Recovery of the primitive variables

For an ideal gas EoS a suitable method to recover the primitive quantities from the state vector was introduced by Martí et al. (1991). It has been described in detail in Martí & Müller (1996). We have extended this method to be adequate for our hybrid EoS (25).

We calculate a function $`f(P^{})P(\rho ^{},ϵ^{})P^{}`$, where $`\rho ^{}`$ and $`ϵ^{}`$ depend on the conserved quantities and $`P^{}`$ only (in all relations involving the pressure we replace $`P`$ by $`P^{}`$):

$`f(P^{})`$ $`=`$ $`K\left(D\sqrt{X}\right)^\gamma \left(1{\displaystyle \frac{\gamma _{\mathrm{th}}1}{\gamma 1}}\right)`$ (71)
$`+(\gamma _{\mathrm{th}}1)X\left[\tau +D\left(1{\displaystyle \frac{1}{\sqrt{X}}}\right)+P^{}\left(1{\displaystyle \frac{1}{X}}\right)\right]`$
$`{\displaystyle \frac{(\gamma _{\mathrm{th}}1)(\gamma \gamma _1)}{(\gamma _11)(\gamma _21)}}K\rho _{\mathrm{nuc}}^{\gamma _11}D\sqrt{X}P^{},`$
where we use the auxiliary quantity $`X1S^2/(\tau +D+P^{})^2`$, and $`S=S_iS^i`$. The new pressure is iteratively computed using the Newton–Raphson method:

$$P^{\mathrm{new}}=P^{}f(P^{})\left(\frac{\mathrm{d}f(P^{})}{\mathrm{d}P^{}}\right)^1,$$
(72)
where the derivative of $`f(P^{})`$ is given by

$`{\displaystyle \frac{\mathrm{d}f(P^{})}{\mathrm{d}P^{}}}`$ $`=`$ $`{\displaystyle \frac{K\left(D\sqrt{X}\right)^\gamma S^2\left(1\frac{\gamma _{\mathrm{th}}1}{\gamma 1}\right)}{X(\tau +D+P^{})^3}}`$ (73)
$`+{\displaystyle \frac{2(\gamma _{\mathrm{th}}1)S^2\left[\tau +D\left(1\frac{1}{\sqrt{X}}\right)+P^{}\left(1\frac{1}{X}\right)\right]}{(\tau +D+P^{})^3}}`$
$`+(\gamma _{\mathrm{th}}1)X[{\displaystyle \frac{DS^2}{X^{3/2}(\tau +D+P^{})^3}}+1`$
$`{\displaystyle \frac{1}{X}}+{\displaystyle \frac{2P^{}S^2}{X^2(\tau +D+P^{})^3}}]`$
$`{\displaystyle \frac{(\gamma _{\mathrm{th}}1)(\gamma \gamma _1)K\rho _{\mathrm{nuc}}^{\gamma _11}DS}{(\gamma _11)(\gamma _21)\sqrt{X}(\tau +D+P^{})^3}}1.`$
During each iteration, the primitive quantities and the Lorentz factor, which are needed in the EoS, are updated using the following relations:

$`v_i^{}`$ $`=`$ $`{\displaystyle \frac{S_i}{\tau +D+P^{}}},`$ (74)
$`W^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1v_i^{}v^i}}},`$ (75)
$`\rho ^{}`$ $`=`$ $`{\displaystyle \frac{D}{W^{}}},`$ (76)
$`ϵ^{}`$ $`=`$ $`{\displaystyle \frac{\tau +D(1W^{})+P^{}(1W^2)}{DW^{}}}.`$ (77)
The iteration is continued until convergence is reached.

## Appendix B Explicit form of the metric equations

Using the notation of Eqs. (52) and (53), the metric equations read (with $`u_{,f}=_fu`$, $`u_{,fg}=_f_gu`$):

$`f^1(u^k)`$ $`=`$ $`u_{,rr}^1+{\displaystyle \frac{2u_{,r}^1}{r}}+{\displaystyle \frac{u_{,\theta \theta }^1}{r^2}}+{\displaystyle \frac{\mathrm{cot}\theta }{r^2}}u_{,\theta }^1`$ (78)
$`+2\pi (u^1)^5(\rho hW^2P)`$
$`+{\displaystyle \frac{(u^1)^7}{48(u^2)^2}}[3r^2\mathrm{sin}^2\theta (u_{,r}^5)^2+{\displaystyle \frac{3(u_{,\theta }^3)^2}{r^2}}+{\displaystyle \frac{4(u^3)^2}{r^2}}.`$
$`+4(u_{,r}^3)^2+3r^2(u_{,r}^4)^2+3\mathrm{sin}^2\theta (u_{,\theta }^5)^2{\displaystyle \frac{8u^3u_{,r}^3}{r}}`$
$`+6u_{,\theta }^3u_{,r}^4+{\displaystyle \frac{4u^3u_{,\theta }^4}{r}}4u_{,r}^3u_{,\theta }^4+4(u_{,\theta }^4)^2`$
$`+{\displaystyle \frac{4\mathrm{cot}\theta u^3u^4}{r}}4\mathrm{cot}\theta u^4u_{,r}^3`$
$`.4\mathrm{cot}\theta u^4u_{,\theta }^4+4\mathrm{cot}^2\theta (u^4)^2]=0,`$
$`f^2(u^k)`$ $`=`$ $`u_{,rr}^2+{\displaystyle \frac{2u_{2,r}^2}{r}}+{\displaystyle \frac{u_{,\theta \theta }^2}{r^2}}+{\displaystyle \frac{\mathrm{cot}\theta }{r^2}}u_{,\theta }^2`$ (79)
$`2\pi (u^1)^4u^2(\rho h(3W^22)+5P)`$
$`{\displaystyle \frac{7(u^1)^7}{48(u^2)^2}}[3r^2\mathrm{sin}^2\theta (u_{,r}^5)^2+{\displaystyle \frac{3(u_{,\theta }^3)^2}{r^2}}+{\displaystyle \frac{4(u^3)^2}{r^2}}.`$
$`+4(u_{,r}^3)^2+3r^2(u_{,r}^4)^2+3\mathrm{sin}^2\theta (u_{,\theta }^5)^2{\displaystyle \frac{8u^3u_{,r}^3}{r}}`$
$`+6u_{,\theta }^3u_{,r}^4+{\displaystyle \frac{4u^3u_{,\theta }^4}{r}}4u_{,r}^3u_{,\theta }^4+4(u_{,\theta }^4)^2`$
$`+{\displaystyle \frac{4\mathrm{cot}\theta u^3u^4}{r}}.4\mathrm{cot}\theta u^4u_{,r}^3`$
$`4\mathrm{cot}\theta u^4u_{,\theta }^4+4\mathrm{cot}^2\theta (u^4)^2]=0,`$
$`f^3(u^k)`$ $`=`$ $`u_{,rr}^3+{\displaystyle \frac{2u_{,r}^3}{r}}+{\displaystyle \frac{u_{,\theta \theta }^3}{r^2}}+{\displaystyle \frac{\mathrm{cot}\theta }{r^2}}u_{,\theta }^3`$ (80)
$`{\displaystyle \frac{2u^3}{r^2}}{\displaystyle \frac{2\mathrm{cot}\theta u^4}{r}}{\displaystyle \frac{2u_{,\theta }^4}{r}}{\displaystyle \frac{16\pi u^2S_1}{u^1}}`$
$`{\displaystyle \frac{1}{u^2}}[(u_{,\theta }^2{\displaystyle \frac{7u^2u_{,\theta }^1}{u^1}})(u_{,r}^4+{\displaystyle \frac{u_{,\theta }^3}{r^2}}).`$
$`+{\displaystyle \frac{2}{3}}(u_{,r}^2{\displaystyle \frac{7u^2u_{,r}^1}{u^1}})(2u_{,r}^3u_{,\theta }^4{\displaystyle \frac{2u^3}{r}}.`$
$`..\mathrm{cot}\theta u^4)]+{\displaystyle \frac{1}{3}}(u_{,rr}^3+u_{,r\theta }^4+{\displaystyle \frac{2u_{,r}^3}{r}}.`$
$`.{\displaystyle \frac{2u^3}{r^2}}+\mathrm{cot}\theta u_{,r}^4)=0,`$
$`f^4(u^k)`$ $`=`$ $`u_{,rr}^4+{\displaystyle \frac{4u_{,r}^4}{r}}+{\displaystyle \frac{u_{,\theta \theta }^4}{r^2}}+{\displaystyle \frac{\mathrm{cot}\theta }{r^2}}u_{,\theta }^4+{\displaystyle \frac{2u_{,\theta }^3}{r^3}}`$ (81)
$`+{\displaystyle \frac{(1\mathrm{cot}\theta )u^4}{r^2}}{\displaystyle \frac{16\pi u^2S_2}{r^2u^1}}`$
$`{\displaystyle \frac{1}{u^2}}[(u_{,r}^2{\displaystyle \frac{7u^2u_{,r}^1}{u^1}})(u_{,r}^4+{\displaystyle \frac{u_{,\theta }^1}{r^2}}).`$
$`+{\displaystyle \frac{2}{3r^2}}(u_{,\theta }^2{\displaystyle \frac{7u^2u_{,\theta }^1}{u^1}})(u_{,r}^3+2u_{,\theta }^4.`$
$`..+{\displaystyle \frac{u^3}{r}}\mathrm{cot}\theta u^4)]+{\displaystyle \frac{1}{3r^2}}(u_{,r\theta }^3+u_{,\theta \theta }^4.`$
$`.+{\displaystyle \frac{2u_{,\theta }^3}{r}}{\displaystyle \frac{u^4}{\mathrm{sin}^2\theta }})=0,`$
$`f^5(u^k)`$ $`=`$ $`u_{,rr}^5+{\displaystyle \frac{4u_{,r}^5}{r}}+{\displaystyle \frac{u_{,\theta \theta }^5}{r^2}}+{\displaystyle \frac{3\mathrm{cot}\theta }{r^2}}u_{,\theta }^5`$ (82)
$`{\displaystyle \frac{16\pi u^2S_3}{r^2\mathrm{sin}^2\theta u^1}}{\displaystyle \frac{1}{u^2}}[(u_{,r}^2{\displaystyle \frac{7u^2u_{,r}^1}{u^1}})u_{,r}^5`$
$`+{\displaystyle \frac{1}{r^2}}(u_{,\theta }^2{\displaystyle \frac{7u^2u_{,\theta }^1}{u^1}})u_{,\theta }^5]=0.`$