A&A 381, L49-L52 (2002)
DOI: 10.1051/0004-6361:20011643

Highly accurate calculation of rotating neutron stars

M. Ansorg - A. Kleinwächter - R. Meinel

Theoretisch-Physikalisches Institut, University of Jena, Max-Wien-Platz 1, 07743 Jena, Germany

Received 6 November 2001/ Accepted 20 November 2001

Abstract
A new spectral code for constructing general-relativistic models of rapidly rotating stars with an unprecedented accuracy is presented. As a first application, we reexamine uniformly rotating homogeneous stars and compare our results with those obtained by several previous codes. Moreover, representative relativistic examples corresponding to highly flattened rotating bodies are given.

Key words: stars: rotation - stars: neutron - gravitation - relativity - methods: numerical


1 Introduction

The study of relativistic, axisymmetric and stationary, uniformly rotating perfect fluid bodies is motivated by extraordinarily compact astrophysical objects like neutron stars. Several numerical codes have been developed in order to calculate the structure and the gravitational field of these bodies (Bonazzola & Schneider 1974; Wilson 1972; Butterworth & Ipser 1975,1976; Friedman et al. 1986,1989; Lattimer et al. 1990; Neugebauer & Herold 1992; Herold & Neugebauer 1992; Komatsu et al. 1989a,1989b; Eriguchi et al. 1994; Stergioulas & Friedman 1995; Bonazzola et al. 1993; for reviews see Friedman 1998, and Stergioulas 1998). While they obtain an accuracy of up to 5 digits for sufficiently smooth equations of state, these methods yield fewer than 4 digits in the case of stiff equations of state (e.g., for constant density), which is due to particular Gibbs phenomena at the star's surface. In order to avoid these Gibbs phenomena, Bonazzola et al. (1998) used a multi-domain spectral method with which they were able to achieve an accuracy of 12 digits for the Maclaurin sequence of homogeneous Newtonian bodies.

In this letter we introduce a new numerical code, which is based on a multi-domain spectral method for representing all metric functions. We intend to use this method to investigate neutron stars with realistic equations of state. In particular, our multi-domain method lends itself to considering several layers inside the star, which are characterized by different equations of state. As we will outline below, we obtain a hitherto unobtainable accuracy which permits its application even in limiting cases such as the mass-shedding limit. Moreover, we are able to study extremely flattened, homogeneous Einsteinian bodies. Such bodies were the subject of interesting conjectures made by Bardeen (1971).

As a first application of our method, we reexamine a particular example of a uniformly rotating homogeneous star that was used by Nozawa et al. (1998) to compare three different codes. We give our results, which possess a substantially higher accuracy. Moreover, we discuss representative, relativistic examples corresponding to highly flattened rotating bodies.

In what follows, units are used in which the velocity of light as well as Newton's constant of gravitation are equal to 1.

2 Metric tensor and field equations

The line element for an axisymmetric, stationary, uniformly rotating perfect fluid body assumes in Lewis-Papapetrou coordinates $(\rho,\zeta,\varphi,t)$ the following form:

\begin{displaymath}{\rm d}s^2={\rm e}^{-2U}[{\rm e}^{2k}({\rm d}\rho^2+{\rm d}\z... ...2 {\rm d}\varphi^2] -{\rm e}^{2U}({\rm d}t+a{\rm d}\varphi)^2.\end{displaymath}

To define the coordinates $(\rho,\zeta)$ uniquely, we require that the metric coefficients and their first derivatives be continous at the surface of the body.

In the vacuum region, there emerge three field equations of second order to determine the potentials U, a and W. The function k follows from the other potentials by a line integral[*].

In the interior of the body we use the metric functions valid in the comoving frame of reference. Here, the only new coordinate is $\varphi'=\varphi-\Omega t$, where $\Omega $ is the angular velocity of the body. The corresponding line element also assumes the above form with potentials U', a', k' and W', which are given by

\begin{displaymath}{\rm e}^{2U'}={\rm e}^{2U}[(1+\Omega a)^2-\Omega^2 W^2 {\rm e}^{-4U}],\end{displaymath}


\begin{displaymath}(1-\Omega a'){\rm e}^{2U'} = (1+\Omega a){\rm e}^{2U},\end{displaymath}


\begin{displaymath}k'-U'=k-U \quad \mbox{and} \quad W'=W.\end{displaymath}

Since in the comoving frame the energy-momentum tensor reads

\begin{displaymath}T^{ik}=(\mu+p)u^iu^k+pg^{ik},\quad u^k={\rm e}^{-U'}\delta_4^k,\end{displaymath}

where $\mu$ is the mass-energy density and p the pressure, the field equations assume a particularly simple form, see, e.g., Kramer et al. (1980), Eqs. (19.35a-c). For a given equation of state, $p=p(\mu) $ or $\mu=\mu(p)$, the relativistic Euler equations $T^{ik}_{\quad;k}=0$ yield

\begin{displaymath}{\rm e}^{U'}\exp\left[\int_0^p\frac{{\rm d}p'}{\mu(p')+p'}\right] ={\rm e}^{V_0}=\mbox{const.}\end{displaymath}

Hence, pressure and density can be expressed in terms of U'. At the surface B of the star, the pressure vanishes which leads to a constant surface potential U'=V0. In particular, for homogeneous stars ( $\mu=\mbox{const.}$) we obtain

\begin{displaymath}p=\mu\left({\rm e}^{V_0-U'}-1\right).\end{displaymath}

Taking formula (19.35d)[*] and the condition (19.37) of Kramer et al. (1980) we may express k' by the potentials a', U' and W via a line integral. Thus again a system of three field equations emerges.

All potentials satisfy regularity conditions at infinity and along the axis of rotation ($\rho=0$) and possess moreover reflectional symmetry with respect to the plane $\zeta=0$ (see Meinel & Neugebauer 1995).

3 The numerical scheme

The numerical scheme to solve the field equations with respect to boundary and transition conditions is based on a multi-domain spectral method. After imposing reflectional symmetry, the set of all relevant $(\rho,\zeta)$-values, $\{(\rho,\zeta)\mbox{:}\quad\rho\geq 0,\,\zeta\geq0\}$, is divided into several subregions for physical reasons. In the simplest case, we only take two subregions, the interior and the exterior of the body. However, if we consider several layers inside the star, which are characterized by different equations of state, we will be forced to allow for more than two regions. In this letter, we will restrict ourselves to only two subregions.

Each of these subregions is mapped onto the square I2=[0,1]2. In order to do this we introduce a function $y_{\rm B}$ defined on the interval I=[0,1] with

\begin{displaymath}y_{\rm B}(0)=1\,,\quad y_{\rm B}(1)=0, \end{displaymath}

which describes the surface of the body by

\begin{displaymath}B=\left\{(\rho,\zeta)\mbox{:}\quad\rho^2=r_{\rm e}^2\hspace*{... ...p}^2\hspace*{0.5mm}y_{\rm B}(t)\,, \quad 0\leq t\leq 1\right\},\end{displaymath}

where $r_{\rm e}$ and $r_{\rm p}$ are the equatorial and polar coordinate radii of the body respectively. As a particular example for the mapping in question, we used for the interior the transformation

\begin{displaymath}\rho^2=r_{\rm e}^2\hspace*{0.1mm}s\hspace*{0.1mm}t\,,\quad \... ...hspace*{0.1mm}s\hspace*{0.1mm}y_{\rm B}(t)\,,\quad (s,t)\in I^2\end{displaymath}

and for the exterior

\begin{displaymath}\rho^2=\frac{r_{\rm e}^2\hspace*{0.1mm}t}{s^2}\,,\quad \zeta... ...\rm p}^2\hspace*{0.5mm}y_{\rm B}(t)}{s^2}\,,\quad (s,t)\in I^2.\end{displaymath}

In this manner, the axis $\rho=0$ and the plane $\zeta=0$ are mapped to the coordinate boundaries t=0 and t=1 respectively. Furthermore, the surface B of the body is characterized by s=1. For the interior and exterior transformation, the point s=0 corresponds to the origin and to infinity respectively.

We assume all potentials to be smooth functions on I2 such that we may approximate them well by two-dimensional Chebyshev-expansions with respect to the coordinates s and t. In the same manner, we represent the unknown boundary function $y_{\rm B}$ as well as the boundary values $a'_{\rm B}$ and $W_{\rm B}$ of the potentials a' and W in terms of (one-dimensional) Chebyshev-polynomials with respect to the coordinate t. If these three one-dimensional functions were given, we would have to solve a particular interior and exterior boundary value problem[*] of the respective field equations. However, we have to deal with a free boundary value problem, where these three functions are not known, but have to be determined such that the normal derivatives of the potentials U', a' and W behave continuously at the surface B [*]. Taking only a finite number of Chebyshev-coefficients into account for the interior potentials U', a' and W, the exterior potentials U, a and W and the surface quantities $y_{\rm B}$, $a'_{\rm B}$ and $W_{\rm B}$, our numerical scheme consists in determining all unknown Chebyshev-coefficients by satisfying the interior and exterior field equations at a number of grid points in I2 and, moreover, requiring the above transition conditions at the surface. The total number of unknown coefficients equals the total number of equations. The system is solved by a Newton-Raphson method, where the initial guess for the entire solution in the case of constant density can be taken from the analytical Newtonian Maclaurin solution.

4 Results

4.1 Comparison with previous codes

Nozawa et al. (1998) compared three different codes for various choices of the equation of state. We take one particular example in which for constant density they prescribed the normalized central pressure $\bar{p}_{\rm c}=p_{\rm c}/\mu=1$ and the ratio $r_{\rm p}/r_{\rm e}=7\times 10^{-1}$. The results of this comparison are given in Table 11 of Nozawa et al. (1998). Here we calculate the same quantities with an accuracy of 12 digits and list them in Col. 2 of Table 1[*]. Columns 3-5 refer to the codes by Komatsu et al. (1989a, 1989b) and Eriguchi et al. (1994) [abbreviated by KEH(OR)], Stergioulas & Friedman (1995) [abbreviated by KEH(SF)] and Bonazzola et al. (1993) [abbreviated by BGSM] and give the relative error of the quantity in question, i.e. for example, $\vert\bar{M}[\mbox{KEH(OR)}]-\bar{M}\vert/\bar{M}\approx0.023$. The fact that $r_{\rm p}/r_{\rm e}$ was not exactly 0.7 in the BGSM calculation does not affect the comparison substantially. In Nozawa et al. (1998), the general-relativistic virial identities GRV2 and GRV3 (derived by Bonazzola & Gourgoulhon 1994 and Gourgoulhon & Bonazzola 1994) were calculated to check the accuracy of the numerical solution. For our result, this check yields $1.6\times 10^{-13}$ for GRV2 and $4.4\times 10^{-13}$ for GRV3. Note that we used 23 Chebyshev-polynomials for each dimension in this calculation. As an additional test of accuracy, we calculated the angular momentum and the gravitational mass in two different ways: (i) from the asymptotic behaviour of the metric and (ii) by means of integrals over the matter distribution (cf. Bardeen & Wagoner 1971, Eqs. (II.24), (II.26) and (II.23), (II.25) respectively). We get a relative deviation of $1.7\times 10^{-14}$ for the mass and $6.2\times 10^{-14}$ for the angular momentum.

 

 
Table 1: Results for the constant density model with $\bar{p}_{\rm c}~=~1$, $r_{\rm p}/r_{\rm e}=7\times 10^{-1}$, and relative errors of the codes (1): KEH(OR), (2): KEH(SF) and (3): BGSM. $\bar{\Omega}=\Omega/\mu^{1/2}$, $\bar{M}=M \mu^{1/2}$, $\bar{M}_0={M}_0 \mu^{1/2}$, $\bar{R}_{\rm circ}={R}_{\rm circ}\,\mu^{1/2}$ and $\bar{J}=J\mu$ are normalized values of the angular velocity $\Omega $, gravitational mass M, rest mass M0, equatorial circumferential radius ${R}_{\rm circ}$ and angular momentum J. $Z_{\rm p}$ is the polar redshift, $Z_{\rm eq}^{\rm f}$ and $Z_{\rm eq}^{\rm b}$ are equatorial redshifts for photons emitted in the forward and backward direction.
  (1) (2) (3)
$\bar{p}_{\rm c}$ 1      
$r_{\rm p}/r_{\rm e}$ 0.7     0.11%
$\bar{\Omega}$ 1.41170848318 1.1% 0.32% 0.97%
$\bar{M}$ 0.135798178809 2.3% 0.19% 0.86%
$\bar{M}_0$ 0.186338658186 0.17% 0.32% 1.4%
$\bar{R}_{\rm circ}$ 0.345476187602 0.098% 0.053% 0.27%
$\bar{J}$ 0.0140585992949 1.6% 0.045% 2.3%
$Z_{\rm p}$ 1.70735395213 6.1% 0.013% 2.1%
$Z_{\rm eq}^{\rm f}$ -0.162534082217 1.7% 1.9% 4.4%
$Z_{\rm eq}^{\rm b}$ 11.3539142587 17% 0.10% 8.1%


4.2 Highly flattened rotating bodies

Chandrasekhar (1967) has shown that the post-Newtonian Maclaurin spheroids become singular at the eccentricity $\epsilon =\epsilon_1 = 0.98522$..., the point of the first axisymmetric secular instability. Bardeen (1971) confirmed this result and discussed the possibility of two Newtonian axisymmetric sequences bifurcating from the Maclaurin spheroid at $\epsilon=\epsilon_1$, a first one that should finally evolve towards the Dyson-Wong rings (Dyson 1892,1893; Wong 1974; Kowalewsky 1895; Poincaré 1885a,1885b,1885c; see also Lichtenstein 1933) and a second one, which he called the "central bulge configuration'', that should end in a mass-shedding limit. Eriguchi & Sugimoto (1981) indeed found the first sequence, and called it the "one-ring sequence''. The properties of the bifurcation have been analyzed by Christodoulou et al. (1995).

By means of a Newtonian version of our code we were able to find the second sequence as well. It appears as a continuation of the one-ring sequence across the Maclaurin sequence and is characterized by a (bi-convex) "lens shape'' of the solutions.

The point $\epsilon=\epsilon_1$ of the Maclaurin sequence corresponds to the value R=R1=0.27320... of the rotation parameter $R=\bar{J}^2/\bar{M}_0^{10/3}$. Bardeen (1971) speculated that there should be a gap in the R-values of relativistic solutions around R1 and that the relativistic solutions should show some properties of the Newtonian "central bulge'' or Dyson ring solutions when approaching this gap from the sphere end or from the disc end of the Maclaurin sequence.


 

 
Table 2: Three highly flattened, constant density models. The corresponding surface shapes can be found in Fig. 1. In these calculations we used 23 Chebyshev-polynomials for each dimension. The relative error in the calculated parameters given above is less than 10-6.
  (a) (b) (c)
$\bar{p}_{\rm c}$ 0.002 0.004 0.003
$r_{\rm p}/r_{\rm e}$ 0.2 0.1 0.1
$\bar{\Omega}$ 1.089468e+00 1.004798e+00 9.305304e-01
$\bar{M}$ 8.314919e-04 1.099261e-02 5.603298e-03
$\bar{M}_0$ 8.353271e-04 1.124879e-02 5.677891e-03
$\bar{R}_{\rm circ}$ 1.028186e-01 2.474602e-01 2.249448e-01
$\bar{J}$ 3.663135e-06 3.093789e-04 1.072577e-04
$Z_{\rm p}$ 1.585397e-02 9.132868e-02 5.381392e-02
$Z_{\rm eq}^{\rm f}$ -9.890231e-02 -1.912324e-01 -1.720749e-01
$Z_{\rm eq}^{\rm b}$ 1.308359e-01 3.827610e-01 2.825350e-01
       
R 0.2444465 0.3001205 0.3522896
       
GRV2 1.1e-09 3.1e-11 1.9e-08
GRV3 1.6e-09 3.2e-11 2.3e-08



  \begin{figure} \par\includegraphics[width=6.2cm,clip]{Dk067.eps} \end{figure} Figure 1: Meridional cross-sections of the solutions specified in Table 2.
Open with DEXTER

In order to check these conjectures, we have calculated three relativistic models in the vicinity of R=R1. These models are characterized by parameters as shown in Table 2 and their (coordinate) shapes are depicted in Fig. 1. Note that the high accuracy of our code is crucial when investigating the subtle behaviour of the relativistic solutions in this region - an earlier attempt by Butterworth (1979) could not finally clarify these questions. Solution (a) shows indeed a "lens shape'', whereas solution (b) has a "one-ring tendency'' as does solution (c), albeit far less pronounced. The solution (a) is rather close to the Newtonian "lens shape'' sequence. The situation becomes more and more complex when the bifurcation points R2=0.36633..., R3=0.43527..., etc. (corresponding to $\epsilon_2=0.99375$..., $\epsilon_3=0.99657$..., etc.) of the Newtonian two-ring, three-ring, etc. sequences are approached (for the two-ring sequence see Eriguchi & Hachisu 1982). Our solution (c) already has a value for R close to R2, see Table 2.

A detailed analysis of highly flattened, rotating Newtonian as well as Einsteinian bodies including a discussion of stability aspects and of the route to infinitesimally thin, relativistic discs (Bardeen & Wagoner 1969,1971; Neugebauer & Meinel 1993,1995) will be the subject of a future publication.

Acknowledgements

This work was supported by the German Deutsche Forschungsgemeinschaft (DFG-project ME 1820/1).

References

 

Copyright ESO 2002