1 Introduction

The physical state of the neutron star matter has not been fully understood yet because it is very difficult to investigate particle interactions beyond nuclear matter density ($\varepsilon_{\rm N}/c^2 \sim 2 \ 10^{14}$ g cm^-3) either from nuclear experiments or from nuclear theories, (here $\varepsilon_{\rm N}$ is the energy density of the nuclear matter and c is the velocity of light). Given this situation, one promising approach to explore the behavior of very high density matter is to make use of the macroscopic quantities of neutron stars. In particular, the mass and the rotational period of neutron stars depend crucially on the softness of the equation of state (EOS) at very high densities (see e.g. Friedman et al. 1984, 1986, 1989), thus, observational constrains, matched with theoretical models, may help in reconstructing the equation of state of very high density matter.

Given a particular equation of state, the mass of neutron stars varies with central energy density and always reaches a maximum. This implies that if the maximum mass of neutron star models constructed with a certain equation of state is smaller than the mass of observed neutron stars, that equation of state must be discarded. Currently, the largest accurately measured mass of slowly rotating neutron stars is $M_{\rm BP} = 1.44\ M_{\odot}$, where $M_{\rm BP}$ is the mass of one of the components of the binary pulsar PSR1913+16 (Taylor & Weisberg 1989) and $M_{\odot}$ is one solar mass. Individual masses of neutron stars have also been estimated in six other binary pulsars (Thorsett et al. 1993; Wolszczan 1997), as well as in six X-ray binaries (van Kerkwijk et al. 1995) but the accuracy is not as good as in PSR1913+16. Thus, equations of state which give larger masses than $M_{\rm BP}$ for slowly rotating stars, can be valid as candidates for the real equation of state at very high densities. Since the maximum mass of neutron stars is smaller for more compressible (soft) equations of state than for less compressible (stiff) equations of state, the true equation of state at high densities cannot be extremely soft.

On the other hand, stiff equations of state can be limited by considering the neutron star with the shortest rotational period, i.e. the most rapidly rotating pulsar. There exists a lower limit on the rotational period for each equation of state, because if the centrifugal force exceeds the self-gravity at the equatorial surface, no equilibrium states are allowed. The lower limit of the rotational period depends on the softness of the equation of state - the radius of neutron stars with softer equations of state is smaller, which allows for higher rotation rates. Thus, if very rapidly rotating neutron stars should be found, we could exclude most stiff equations of state. At the moment, the shortest period of observed pulsars is 1.56 ms, of PSR1937+21. Consequently, equations of state for which the shortest rotational period is larger than this value, must be excluded as candidates for the real equation of state for neutron star matter.

The discussions above require us to make use of highly accurate schemes for constructing rotating neutron star models, in order to compute precise theoretical values of masses and rotational periods. Highly accurate relativistic equilibrium models are also needed as initial data for relativistic time-evolution codes (modeling of nonlinear pulsations, collapse and generation of gravitational waves). Recently, a number of groups have succeeded in constructing models of rapidly rotating neutron stars (Friedman et al. 1984, 1986, 1988, 1989; Eriguchi et al. 1994; Salgado et al. 1994a,b; Cook et al. 1994b; Stergioulas & Friedman 1995 - for a recent review see Stergioulas 1998). However, the obtained models by those authors do not always agree with each other (see e.g. Friedman et al. 1989; Eriguchi et al. 1994; Salgado et al. 1994a; Cook et al. 1994b; Stergioulas & Friedman 1995). Although Stergioulas & Friedman (1995) have determined the cause of the discrepancy between models in Friedman et al. (1989) and Eriguchi et al. (1994), (which was due to the use of a slightly different equation of state table), the reasons of smaller differences which remain, even after using exactly the same equation of state, have not been clarified yet. This is because numerical techniques used in the different codes, such as the choice of parameters defining the model, the interpolation method, the method of integrating the field equations, a.s.o. are not the same.

In this paper, three groups using their own codes (Komatsu et al. 1989a; Eriguchi et al. 1994; Salgado et al. 1994a; Stergioulas & Friedman 1995) will decrease the differences between their results to a minimum possible, by tuning each code and using the same parameters, the same schemes of interpolation, the same equations of state, and so on. Since the basic schemes used by the three groups are different, it will be impossible to have exactly the same results and the relative differences between results are a measure of the accuracy of the codes. Models obtained with small relative differences between the three codes can be considered as "standard" models for each equation of state. Furthermore, this direct comparison allows us to investigate the effect that the choice of interpolation method, equation of state and domain of integration has on the accuracy of the codes.