3 Codes

We will compare three different codes, which follow the KEH scheme, developed by Komatsu et al. (1989a,b) for general relativistic polytropes with uniform and differential rotation and improved by Cook et al. (1992, 1994a,b) and the BGSM scheme due to Bonazzola et al. (1993).

Concerning the KEH scheme, for the present comparison, we will use the original KEH code, KEH(OR), and the KEH code by Stergioulas & Friedman (1995), KEH(SF), which follows the Cook et al. (1992, 1994a,b) modification of the KEH scheme.

3.1 A short description of each code

In this section the three different numerical codes are briefly described. Details of these codes can be found in Komatsu et al. (1989a,b) and in Eriguchi et al. (1994) for the KEH(OR) code, in Stergioulas & Friedman (1995) for the KEH(SF) code and in Bonazzola et al. (1993) for the BGSM code.

3.1.1 The KEH(OR) code

Komatsu et al. (1989a) have developed a new scheme for solving rapidly rotating relativistic stars. The Einstein equations for three metric potentials $\nu, \beta$ and $\omega$ are transformed into integral equations by using appropriate Green's functions for the elliptical type differential operators. In principle, one can choose Green's functions which decrease as 1/r or more rapidly at large distances. Consequently boundary conditions at infinity, i.e. asymptotically flat conditions can be easily included in the integral equations. It is noted that in this integral representation the integrand contains the metric and the matter quantities such as the energy density and the pressure. The fourth metric $\alpha$ obeys a first order partial differential equation which can be easily integrated, if the other metric potentials are known. The domain of integration is truncated at a finite distance from the star (roughly twice the equatorial radius) and the metric potentials in the integrands are assumed to vanish at that finite distance (instead of at infinity).

The KEH scheme is the extended version of the self-consistent-field (SCF) scheme which was developed for solving Newtonian rotating stars (Ostriker & Mark 1968; Hachisu 1986) and applied to relativistic rotating stars by Bonazzola & Schneider (1974) with a choice of metric functions different from that of the codes considered in this article. In the SCF method, the iteration proceeds as follows. If one assumes initial guesses for the matter quantities and the metric potentials, new (and better) values for $\nu, \beta$ and $\omega$ can be obtained using the integral representations for the metric potentials The fourth metric potential $\alpha$ can be easily solved as mentioned before. By using newly obtained metric potentials, a new density and a new pressure can be computed from the hydrostationary equilibrium Eq. (8). One needs to repeat the same procedure until the relative differences between the newly obtained quantities and the old ones become less than a certain small number, typically 10^-5.

In the original KEH code, the ratio of the central pressure ($p\rm _c$) to the central energy density ($\varepsilon \rm _c$),

$$\kappa \equiv \frac{p\rm _c} {\rm \varepsilon_c},$$ (9)

and the ratio of the polar radius ($r_{\rm p}$) to the equatorial radius ($r_{\rm e}$), $r_{\rm p}/r_{\rm e}$, are chosen as two independent model parameters which specify the model uniquely. The KEH(OR) code used in this comparison differs from the code used in Komatsu et al. (1989a) only in one aspect, that is mentioned in the next section.

3.1.2 The KEH(SF) code

The KEH(SF) code (Stergioulas & Friedman 1995) differs from the original KEH scheme in two ways. First, it follows Cook et al. (1992) in using a redefined radial variable

$$s \equiv \frac{r} {r + r\rm _e}.$$ (10)

where $r\rm _e$ is the value of the coordinate r at the equator. By this transformation, the region $[0, \infty]$ in the r coordinate is transformed to the region [0, 1] in the s coordinate. Consequently, the domain of integration of the field equations does not have to be truncated at a finite distance from the star (as in the KEH(OR) code) but covers all space. With this choice, the boundary conditions can be applied accurately at infinity.

Second, Stergioulas & Friedman found that the choice of coordinates in the original KEH scheme results in the metric potential $\alpha$oscillating in the radial direction. The oscillation is especially pronounced inside the star and introduces an error of $1-2\%$ in the mass, radius and other quantities. The problem was fixed by using a finite difference formula for the second order radial derivative that uses twice the grid-spacing. Although this formula is, in principle, of lower accuracy, the oscillations are damped completely, resulting in a more accurate stellar model.

The KEH(OR) code used in this comparison has been modified so as to use the same second order derivative formula as Stergioulas & Friedman, to smooth out the oscillations in the metric potential $\alpha$.

3.1.3 The BGSM code

Bonazzola et al. (1993) have developed a new formulation based on the 3 + 1 formalism which has been used in hydrodynamics in general relativity. Their choice of slicing and gauge in the 3 + 1 formalism results in the same form of the metric usually chosen for stationary problems, i.e. (2). Consequently the Einstein equations are reduced to the same differential equations of elliptic type as used by other schemes.

The main part of the BGSM formulation is similar to that of the KEH formulation, except for the metric coefficient $\zeta = \alpha + \nu$ for which a second-order (elliptic) equation is used instead of a first-order equation in KEH. The Einstein equations are reorganized so as to "pick" out the Laplacian operators in two and three dimensional flat spaces and regard all other remaining terms as "source terms" in the Poisson-like equations in two and three dimensional flat spaces. Concerning the matter, essentially the same equation as (8) is used for the hydrostationary equation.

The characteristic features of the BGSM code can be found in the numerical solving method, i.e. the pseudo-spectral method (Gottlieb & Orszag 1977; Bonazzola et al. 1996, 1997, 1998b). In the spectral method all functions are expanded in terms of certain base functions and algebraic equations for coefficients which appear in the expansion are solved. Therefore there are two distinct procedures in this method: one is obtaining coefficients from the functions and the other is constructing functions by using the coefficients. Since these two steps need to be, in general, performed many times, it is highly desirable to use a fast algorithm. In the spectral method of the BGSM code, Bonazzola et al. (1993) have adopted trigonometric functions for the angle variable and the Chebyshev polynomials for the radial variable. Consequently for the angle part of any function, the fast fourier transform (FFT) can be employed. Concerning the radial variable, a new variable which is related to the radial variable by a simple equation is introduced so that the Chebyshev polynomials are expressed by the trigonometric functions. After this transformation, one can use the FFT also for the radial variable.

The BGSM code can handle the region extended to infinity as is done by the KEH(SF) code. This can be performed by introducing a new radial variable u as follows:

$$u \equiv \frac{R_0} {r}, \qquad {\ \rm for \ outside \ of \ the \ matter}$$ (11)

where R₀ is the maximum value of the stellar radius. The boundary conditions at infinity, are easily applied at u = 0.

It may be fair to note that in the Newtonian rotating star problems a similar expansion was used by Ostriker & Mark (1968) in the SCF method, although Ostriker and Mark used the integral form of the Newtonian potential instead of solving the Poisson equation directly and they did not use the FFT.

3.2 Relations among the three different codes

Here we summarize similarities and difference between the three codes:

A) Common features through all three codes:

1.

Quasi-isotropic coordinates are used. It implies that the metric components are essentially the same, although background views deriving the metric are not the same.

2.

Integral form of the hydrostationary equation for the matter is used.

3.

Poisson-like operators are "picked" out from the Einstein equations and the other terms are treated as "source terms".

B) The KEH(OR) code differs from the other two codes in that the boundary conditions are not applied at infinity, but approximate boundary conditions are applied at a finite distance from the star.

C) A difference between the BGSM code and the other codes is the use of a second-order (elliptic)equation for $\zeta$ in BGSM versus a first-order equation for $\alpha=\zeta-\nu$ in KEH.

We reorganized our codes so as to make the differences as small as possible. The codes agree exactly on:

1.

the physical model parameters by which we can specify a model uniquely,

2.

the values of physical constants, and

3.

the equation of state of matter.

Since the three codes use different grids and/or numerical methods for solving the field equations, there will always be a residual difference in the results, even after this reorganization. This residual difference is what we want to determine.

3.3 Starting point of computations

3.3.1 Constants

Values used in this paper for the velocity of light c, the gravitational constant G, the mass of the sun $M_{\odot}$ and the baryon mass $m_{\rm B}$ are as follows:

$$\begin{eqnarray} c & = & 2.9979 \ 10^{10}\ {\rm cm \ s}^{-1}, \nonumber \\ G & ... ... \nonumber \\ m_{\rm B} & = & 1.66 \ 10^{-24}\ {\rm g}. \nonumber\end{eqnarray}$$

Note that some of the above constants differ slightly from the ones used in the papers where the three codes were first presented.

3.3.2 Interpolation of tabulated equation-of-state data

For realistic equations of state, the energy density, pressure and other thermodynamical quantities are given in tables. Intermediate values need to be obtained by a method of interpolation. We will use two different interpolation schemes, the four-point Lagrange interpolation (hereafter LI) and the cubic Hermite interpolation (HI) (Swesty 1996):

A) Lagrange interpolation.

Let us assume that there is a table which relates the variable x to the variable y at n points, i.e. a set of values (x_i, y_i) for $i = 1, \dots, n$ are tabulated. For the LI scheme the interpolated formula $y_{\rm LI}$ can be expressed as

$$y_{\rm LI}(x) \equiv \sum_{i = 1}^n y_i \frac{p_i(x)} {p^{'}(x)\vert _{x=x_i}},$$ (12)

where

$$\begin{eqnarray} p(x) & \equiv & (x-x_1) (x-x_2) \dots (x-x_n), \\ p_i(x) & \equiv & \frac{p(x)} {(x-x_i)}.\end{eqnarray}$$ (13) (14)

The prime denotes the differentiation with respect to the argument. This scheme does not guarantee the values of derivatives at the points in the table. In this paper we use the four point Lagrange interpolation, i.e. n = 4.

B) Hermite interpolation.

In the Hermite interpolation, the interpolated formula $y_{\rm HI}$for $x_i \le x \le x_{i+1}$ is expressed as

$$\begin{eqnarray} y_{\rm HI}(x) & \equiv & y_i h_0(w) + y_{i+1} h_0(1-w) \nonumbe... ...\frac{{\rm d}y}{{\rm d}x} \right)_{i+1} (x_{i+1} - x_i) h_1(1-w),\end{eqnarray}$$ (15)

where

$$w \equiv \frac {x - x_i} {x_{i+1} - x_i}.$$ (16)

Here h₀ and h₁ are the cubic Hermite functions defined by

$$\begin{eqnarray} h_0(w) & \equiv & 2 w^3 - 3 w^2 + 1, \\ h_1(w) & \equiv & w^3 - 2 w^2 + w .\end{eqnarray}$$ (17) (18)

It is noted that these cubic Hermite functions are unique polynomials of degree three satisfying the following relations:

$$\begin{eqnarray} & & h_0(0) = 1, \\ & & h_0(1) = h'_0(0) = h'_0(1) = 0, \\ & & h'_1(0) = 1, \\ & & h_1(0) = h_1(1) = h'_1(1) = 0. \end{eqnarray}$$ (19) (20) (21) (22)

Contrary to the LI, in the Hermite interpolation the values as well as their first derivatives at mesh points are exactly reproduced by the interpolation formula. The main advantage of the Hermite interpolation is that it preserves the thermodynamical consistency (Swesty 1996).