5 Results

5.1 Models selected for polytropes

We have started our comparison project by selecting several representative polytropic models, the parameters of which are shown in Tables 1-2. We have chosen

1.

models of very low central density (nearly Newtonian) with slow and rapid rotation,

2.

models of high central density (relativistic) with slow and rapid rotation, and

3.

models at the maximum mass for each EOS.

Note that, the maximum mass model almost coincides with the maximum angular velocity model, unless there is a large phase transition at densities close to the central density of the maximum mass star (cf. Cook et al. 1994b and Stergioulas & Friedman 1995). For all EOSs in this comparison the two models almost coincide.

In order to evaluate the performance of our numerical codes for models with discontinuous density distribution, we also compare a number of homogeneous models which cover both highly relativistic and Newtonian, rapidly rotating and nonrotating cases, as shown in Tables 3 and 4 (the contents of these tables will be described in Sect. 5.4).

 

Table 1:   Polytropic models

 

Table 2:   Polytropic models (continued)

 

Table 3:   Spherical constant density models

 

Table 4:   Rotating constant density models

5.2 Models selected for realistic equations of state

As discussed in the previous section, we use six representative realistic EOSs: C, G, L, WFF3+FPS, WFF3+NV and FPS. In addition, we use the causal limit EOS CLES. For each equation of state, we compute several models as shown Tables 5 - 7. The models correspond to the maximum mass model, a fast rotating $1.4\ M_\odot$ model and a nonrotating model for each EOS.

 

Table 5:   Spherical models with realistic EOSs

 

Table 6:   Rotating models with realistic EOSs

5.3 Computed quantities

5.3.1 Grid and physical quantities

For this comparison project, KEH(OR) and KEH(SF) have used grids with ($71 \times 201)$ and ($261 \times 401$) (angular$\times$radial) grid-points. In the equatorial plane, half of the radial grid-points are inside the star. BGSM uses 21 $\times$ 41 or 33 $\times$ 65 grid points (note that the notion of "grid points'' is not very significant for a spectral method; the above numbers should better be referred to as the numbers of basis functions employed in the expansions of the physical fields).

Here we summarize the notation of computed physical quantities:

$\varepsilon_{\rm c}$ Central energy density
$r_{\rm p}/r_{\rm e}$ Ratio of polar to equatorial radii
$\Omega$ Angular velocity of the star
M₀ Baryon mass
$M\rm _p$ Proper mass
M Gravitational mass
$R_{\rm circ}$ Equatorial circumferential radius
$r_{\rm e}$ Equatorial coordinate radius
J Total angular momentum
I Moment of inertia about the rotation axis
T Rotational energy
W Gravitational energy
$v_{\rm e}$ Velocity of comoving observer at the equator
relative to the locally non-rotating observer
$Z\rm _p$ Polar redshift
$Z\rm _c$ Central redshift
$Z_{\rm eq}^{\rm b}$ Equatorial redshift in the backward direction
$Z_{\rm eq}^{\rm f}$ Equatorial redshift in the forward direction
e Intrinsic eccentricity of the star's surface
GRV2 Two dimensional virial identity
GRV3 Three dimensional virial identity.
Some of the quantities in the above list can be expressed as follows:

$$\begin{eqnarray} Z_{\rm p} & = & {\rm e}^{-2 \nu_{\rm p}} - 1, \\ Z_{\rm eq}^{\r... ...m e} {\rm e}^{(\nu_{\rm e}-\beta_{\rm e})/2} \omega_{\rm e}} - 1, \end{eqnarray}$$ (26) (27) (28)

where subscripts $\rm _p$ and $\rm _e$ denote values at the pole and the equatorial surface, respectively.

$$M_0 = 2 \pi \int \rho {{\rm e}^{2 \alpha + \beta} \over \sqrt{ 1- v^2}} r^2 \sin \theta {\rm d}r {\rm d}\theta,$$ (29)

$$M_{\rm p} = 2 \pi \int \varepsilon {{\rm e}^{2 \alpha + \beta} \over \sqrt{ 1- v^2}} r^2 \sin \theta {\rm d}r {\rm d}\theta,$$ (30)

$$\begin{eqnarray} M &=& 2 \pi \int \Biggl[ {\rm e}^{2 \alpha + \beta} \left\{ {(... ...p) v \over 1-v^2} \Biggr] r^2 \sin \theta {\rm d}r {\rm d}\theta,\end{eqnarray}$$ (31)

$$J = 2 \pi \int {\rm e}^{2\alpha + 2 \beta} {(\varepsilon + p) v \over 1-v^2} r^3 \sin^2 \theta {\rm d}r {\rm d}\theta,$$ (32)

$$T = {1 \over 2} \int \Omega {\rm d}J = 2 \pi \int {\rm e}^{2... ...v \over 1-v^2} \Omega r^3 \sin^2 \theta {\rm d}r {\rm d}\theta,$$ (33)

$$W = M_{\rm p} c^2 + T - M c^2,$$ (34)

and

$$I = { J \over \Omega}.$$ (35)

The eccentricity of the meridional cross section is defined by the following procedure (Friedman et al. 1986). If the surface of the star is defined by

$$r = r_{\rm s}(\theta),$$ (36)

the metric of the stellar surface can be expressed as

$${\rm d} \sigma^2 = {\rm e} ^{2 \beta} r^2 \sin^2 \theta {\rm... ...heta} \right)^2 + r_{\rm s}(\theta)^2 \right] {\rm d} \theta^2.$$ (37)

If we embed this surface in the flat three dimensional space, it is expressed as

$$R = R_{\rm s}(z),$$ (38)

in cylindrical coordinates $(R, \varphi, z)$.The 2-metric of this surface is

$${\rm d} \sigma^2 = \left[ \left({ {\rm d} R_{\rm s} \over {\... ...ight)^2 + 1 \right] {\rm d}z^2 + R_{\rm s}^2 {\rm d} \varphi^2.$$ (39)

Comparing these two equations, we have the following relations, if they express the same surface geometry:

$$R_{\rm s}(\theta) = {\rm e}^{\beta} r \sin \theta,$$ (40)

and

$$\begin{eqnarray} z_{\rm s}(\theta) = \int_\theta^{\pi/2} {\rm d} \theta & \Biggl... ...{\rm d} R_{\rm s} \over {\rm d} \theta} \right)^2 \Biggr\}^{1/2}.\end{eqnarray}$$ (41)

Using these quantities the eccentricity is defined as  

$${\rm e} \equiv \sqrt{ 1 - \left({ z_{\rm s}(\theta = 0) \over R_{\rm s}(\theta=\pi/2)} \right)^2 }.$$ (42)

For polytropes, it is convenient to express quantities in dimensionless form, by using K^N/2 as a fundamental length scale, as was done in Cook et al. (1994a). In geometrized units (c = G = 1), dimensionless quantities are define as follows:

$$\begin{eqnarray} \bar r & \equiv & K^{-N/2} r, \\ \bar R_{\rm circ} & \equiv & K... ...\bar M & \equiv & K^{-N/2} M, \\ \bar M_0 & \equiv & K^{-N/2} M_0.\end{eqnarray}$$ (43) (44) (45) (46) (47) (48) (49) (50) (51) (52)

 

Table 7:   Rotating models with realistic EOSs (continued)

5.3.2 Virial theorem

Equilibrium configurations in Newtonian gravity satisfy the following relation:

$$2 T + 3 (\gamma -1) U + W = 0 ,$$ (53)

where U is the internal energy. This relation is called the virial relation and has been used to check the accuracy of numerically obtained solutions (see e.g. Ostriker & Mark 1968; Tassoul 1978).

In general relativity, similar relations were first found by Bonazzola (1973). Recently, two virial identities in general relativity have been discovered by Gourgoulhon & Bonazzola (1994) and Bonazzola & Gourgoulhon (1994). Those identities are valid for a general asymptotically flat spacetime. We can use these identities to estimate the numerical error. Let us define two quantities $\lambda_2$ and $\lambda_3$ as follows:

$$\lambda_2 \equiv \frac{ \displaystyle 8 \pi \int_0^{+\infty}... ...\partial \omega \right] \, r \, {\rm d}r \, {\rm d}\theta },$$ (54)

and

$$\begin{eqnarray} \lambda_3 &\equiv& 4 \pi \int_0^{+\infty} \int_0^{\pi} \left[... ...g] {\rm e}^{\psi} \, r \, {\rm d}r\, {\rm d}\theta \Bigg\} ^{-1},\end{eqnarray}$$ (55)

with the abridged notation

$$\partial \mu \partial \psi \equiv {\partial \mu \over \parti... ...er \partial \theta} {\partial \psi \over \partial \theta} \ .$$ (56)

Then, we define:

$$\begin{eqnarray} GRV2 & \equiv & \vert 1 - \lambda_2\vert, \\ GRV3 & \equiv & \vert 1 - \lambda_3\vert.\end{eqnarray}$$ (57) (58)

If the Einstein equations are satisfied, these quantities satisfy the following virial identities:

$$\begin{eqnarray} GRV2 & = & 0, \\ GRV3 & = & 0.\end{eqnarray}$$ (59) (60)

Since exact solutions for the stationary problems satisfy the above relations, we can choose GRV2 and GRV3 as the error indicators for numerically obtained solutions. Note that due to its three dimensional character, GRV3 gives a larger weight to the external layers of the star. GRV3 is a relativistic generalization of the Newtonian virial theorem Eq. (43).

In practice, however, it should be noted that the virial identities in the above form are not always close to the accuracy of numerical results. In particular, for GRV2 the integration is done by integrating in the $\theta$ coordinate as seen from the definition of GRV2. If one does not use the $\theta$ coordinate for solving equilibrium structures, one needs to change variables and in that procedure accuracy may be lost and the resultant values may become worse than the "real" accuracy before the variable change.

Concerning the quantity GRV3, the metric potentials in the vacuum region contribute to the integral considerably. It implies that if only the finite regions are treated, as in the KEH(OR) code, a large portion of the integrand cannot be taken into account, However, the expressions for the GRV2 and GRV3 are not unique because we can do the integral by part and replace the second derivatives with the matter terms, using Einstein's equations. In this way, the contribution far away from the star becomes less important. In the KEH(OR) code, GRV2 and GRV3 are evaluated through

$$\begin{eqnarray} \lambda_2' &\equiv& 8 \pi \int_0^{r_{\rm max}}\int_0^{\pi} \le... ...a-2\nu} r^3\sin^2\theta\, {\rm d}r\, {\rm d}\theta \Bigg\} ^{-1},\end{eqnarray}$$ (61)

$$\begin{eqnarray} \lambda_3' &\equiv& 4 \pi \int_0^{\pi} \int_0^{r_{\rm max}} \l... ...r^2{\rm e}^{\beta}\sin\theta {\rm d}r {\rm d}\theta \Bigg\} ^{-1},\end{eqnarray}$$ (62)

with another abridged notation

$$\partial^2 \psi \equiv {\partial^2 \psi \over {\partial r}... ...{1\over r^2\tan\theta}{\partial \psi \over \partial \theta} \ ,$$ (63)

and $r_{\rm max}$ is the distance of the truncated point beyond which actual numerical computations are not carried out in the KEH(OR) code.

This rewrite does not break the mathematical identity. In a sense, it may make "identity" more trivial, and then what information the identities provide us becomes unclear. However, as far as the same expression of the identity in the same code is used, they can play a role as indicators of accuracy among models solved by each code.

5.4 Tables of models and comparison

  The physical parameters of 45 models, computed in this comparison project, are displayed in Tables 1-7. All quantities are displayed to three significant figures. The lower and the upper bounds on each quantity, as obtained by comparing results from the three codes, are shown in the upper and lower rows for each corresponding quantity, respectively. It follows that quantities expressed in a single row can be regarded as "exact", to three significant figures.

Tables 1 and 2 display results for polytropes with index $N = 0.5, \ 0.75, \ 1.0$and 1.5. For each value of the polytropic index N we compute the following models:

1.

a spherical Newtonian model (denoted by the symbol sn),

2.

a rapidly rotating Newtonian model (rn),

3.

a nearly spherical, relativistic model (sr),

4.

the maximum mass model (mr) and

5.

a rapidly rotating relativistic model (rr).

For the constant density case (N = 0), the spherical Newtonian and spherical relativistic models are displayed in Table 3 and the rapidly rotating Newtonian and rapidly rotating relativistic models in Table 4. While all other models are specified by the central energy density $\rm \varepsilon_c$ and the ratio of the polar radius to the equatorial radius $r_{\rm p}/r_{\rm e}$, constant density models are specified by their central pressure and $r_{\rm p}/r_{\rm e}$.

For realistic equations of state, spherical models are shown in Table 5 and rotating models with $1.4\ M_\odot$ (br: binary pulsar mass and relativistic) and maximum mass models are shown in Tables 6 and 7. For the equation of state L, model L(L)mr uses four-point Lagrange interpolation, while L(H)mr is the same model but computed using cubic Hermite interpolation.

From the tables displaying polytropic models one can see that the three codes have a good agreement on most quantities especially for soft polytropes. For stiff polytropes (N<1.0) the agreement is somewhat smaller. For constant density models, the relative differences between the three codes become several percent. More sensitive quantities are the three redshifts and the eccentricity. It should be noted that redshift factors are local quantities which reflect the metric potentials at each point. This implies that local values of the metric potentials do not have the same agreement between different numerical codes, as integrated global quantities. For the eccentricity, one needs to compute the length along the surface of the star (see the definition of the eccentricity (42)). Since in the KEH codes $\mu=\cos \theta$ is used as the angular variable there arise numerical errors near the pole region, i.e. $\theta \simeq 0$. Thus, the differences in the values of the eccentricity also reflect this numerical error due to the choice of coordinates. This causes differences of up to a few percent in the eccentricity for rapidly rotating models. On the other hand, global quantities such as angular velocity, mass, radius and angular momentum agree quite well among results of different codes.

From Tables 5-7, similar tendencies can be observed for realistic equations of state. Models for the most EOSs, except EOS CLES, have a good agreement between the three codes, although the agreement is not as good as for polytropic models. By comparing models constructed with EOSs WFF(FPS) and WFF(NV), it is evident that the choice of the low density EOSs affects very little the structure of the star.

The main reason for the large differences in the constant density case is that the discontinuous density distribution is creating Gibbs phenomena near the surface and this affects all three codes. The reason for the smaller agreement for realistic EOSs, compared to polytropic EOSs, is that the necessary interpolation between tabulated data affects the accuracy with which the equation of hydrostationary equilibrium is satisfied. For EOS L, the choice of the interpolation scheme also affects the accuracy of the computed models, with the cubic Hermite scheme being a better choice compared to a four-point Lagrange interpolation (see the discussion in a later section). For the other realistic EOSs the choice of the interpolation scheme had a negligible effect on the accuracy of computed models.

 

Table 8:   Detailed comparison of polytropic models

 

Table 9:   Detailed comparison of polytropic models (continued)

 

Table 10:   Detailed comparison of constant density models

 

Table 11:   Detailed comparison of constant density models (continued)

 

Table 12:   Detailed comparison of models with realistic EOSs

 

Table 13:   Detailed comparison of models with realistic EOSs (continued)

 

Table 14:   Detailed comparison of models with realistic EOSs (continued)

5.5 Detailed comparison

In order to investigate further the differences among numerical results obtained by the three codes, we show more detailed results for models: N15sn, N15mr and N05sn in Table 8, N05mr in Table 9, N00sn, N00rn and N00sr in Table 10, N00rr in Table 11, Gmr, Lsr and L(L)mr in Table 12, L(H)mr, WFF(FPS)sr and WFF(FPS)br in Table 13 and CLESmr in Table 14. In these tables, values to eight figures for each physical quantity are shown, as well as the relative differences among results of the three codes.

The three relative differences ${\rm diff1}$, ${\rm diff2}$ and ${\rm diff3}$are defined as

$$\begin{eqnarray} {\rm diff1} & \equiv & {{\rm BGSM} - {\rm KEH(OR)} \over {\rm K... ...iff3} & \equiv & {{\rm KEH(SF)} - {\rm BGSM} \over {\rm BGSM}} \ .\end{eqnarray}$$ (64) (65) (66)

From Tables 8 and 9, we see that the agreement between the KEH(SF) and BGSM codes for the rapidly rotating, relativistic polytropic models models N15mr and N05mr is between 10^-4 and $3 \ 10^{-4}$ in all computed quantities (global and local), except for the eccentricity (due to reasons we explained in the previous section). This very good agreement shows that both the BGSM and KEH schemes (the latter when applying boundary conditions exactly at infinity) are suitable for the construction of highly accurate initial-data configurations of rapidly rotating relativistic stars, modeled as polytropes with typical index N>0.5.

From the same tables, we see that the agreement between the KEH(OR) and BGSM codes is similar to the agreement between KEH(SF) and BGSM in the global quantities of model N15mr but to within 10^-3 for the local quantities of this model. For the stiffer polytrope N05mr the agreement is 10^-3 and 10^-2 for global and local quantities, respectively. This difference in accuracy between KEH(SF) and KEH(OR) is expected, since in KEH(OR) boundary conditions are applied only approximately at a finite distance from the star. Considering how close to the star the domain of integration is truncated, the KEH(OR) code performs very well. This is explained as follows:

Since the integration is performed over only a finite region, the truncated part of the integral, $I_{\rm tr}$, can be expressed as

$$I_{\rm tr}({\bf r}) \equiv \int_{\rm out} S({\bf r^{'}}) G({\bf r}, {\bf r^{'}}) {\rm d}^3 {\bf r^{'}},$$ (67)

where S and G are the source term and the Green's function, respectively, and the subscript out denotes that the integration covers the region with $r \ge R_{\rm max}$. Here $R_{\rm max}$ is the maximum radius of the region of integration. Although the source terms of the integrals for the metric potentials in the KEH(OR) scheme contain both matter and the metric terms, only the metric terms contribute to the integral the outside of the star. Since we impose asymptotically flat conditions for the metric potentials, the radial dependence of the source term can be considered to be

$$S \sim {1 \over r^k} ,$$ (68)

where $k \ge 4$, because of the behavior of metric functions at large radii. Consequently, $I_{\rm tr}$ can be roughly expressed as

$$\begin{eqnarray} I_{\rm tr}({\bf r}) & \sim & {\rm constant} , \qquad (r < r_{\r... ...over R_{\rm max}} \right) . \qquad (r_{\rm e} < r < R_{\rm max}).\end{eqnarray}$$ (69)

From the above equation, we can see the following: First, for the stellar interior, i.e. $r < r_{\rm e}$,$I_{\rm tr} \sim$ constant, so that the relative error $I_{\rm tr}/I_0$can be very small where I₀ is the exact value of the integral. Second, for the outer region of the stars, i.e. $r_{\rm e} \le r \le R_{\rm max}$, the value of $I_{\rm tr}$ becomes larger for larger values of r because of the logarithmic term. Third, for the same region, i.e. $r_{\rm e} \le r \le R_{\rm max}$, if the source term depends weakly on the radial coordinate, i.e., for smaller values of k, the contribution of the truncated part becomes larger. This occurs for the potentials $\nu$ and $\omega$. However, for the metric functions B and $\zeta$, since values of k are large, the relative differences are rather small.

In the extreme case of constant density models (Tables 10 and 11), the three codes agree on the computed physical quantities typically only within a few percent and this is caused by the sharp density discontinuity at the surface of the star. The numerical schemes in this comparison assume that the density distribution is a smooth function of coordinates, thus, in the case of density discontinuities, this assumption is violated and Gibbs phenomena appear, resulting in low accuracy of the computed models.

From Tables 12-14 it follows that the agreement of the KEH(SF) code to the BGSM code is between 10^-3 and 10^-4 for realistic EOSs, except for EOS CLES, where the agreement is an order of magnitude smaller. KEH(OR) and BGSM agree on the realistic models within 10^-2 and 10^-3, i.e. similar to the agreement for the N=0.5 polytrope. The somewhat lesser agreement for realistic EOSs is due to the use of interpolation between the tabulated equation of state data (see the discussion next).

In Tables 8 to 14, we also display the virial quantities GRV2 and GRV3. In the ideal case, these should exactly vanish, so the smaller the values for GRV2 and GRV3 are, the better is the accuracy of the computed model. The opposite is not always true, i.e. in some models the computed values for GRV2 or GRV3 do not reflect an overall better agreement in physical quantities among the different codes. This indicates that the computation of GRV2 and GRV3 may itself be prone to numerical error. This seems to be the case for GRV2 in rapidly rotating models computed with the KEH(SF) code, where one first has to interpolate data between $\cos \theta$ and $\theta$ grids to be able to compute GRV2. Note that the displayed values of the two virial quantities for the KEH(OR) code correspond to the modified virial identities (61) and (62) and not to the original identities. The computation of GRV3 for the KEH(OR) code is affected significantly by the truncation of the domain of integration.