2.7 Numerical schemes

All available methods for solving the system of equations describing the equilibrium of rotating relativistic stars are numerical, as no analytical self-consistent solution for both the interior and exterior spacetime has been found. The first numerical solutions were obtained by Wilson [323] and by Bonazzola and Schneider [48]. Here, we will review the following methods: Hartle’s slow rotation formalism, the Newton–Raphson linearization scheme due to Butterworth and Ipser [55Jump To The Next Citation Point], a scheme using Green’s functions by Komatsu et al. [183Jump To The Next Citation Point184Jump To The Next Citation Point], a minimal surface scheme due to Neugebauer and Herold [234], and spectral-method schemes by Bonazzola et al. [47Jump To The Next Citation Point46Jump To The Next Citation Point] and by Ansorg et al. [18Jump To The Next Citation Point]. Below we give a description of each method and its various implementations (codes).

2.7.1 Hartle

To order 𝒪 (Ω2 ), the structure of a star changes only by quadrupole terms and the equilibrium equations become a set of ordinary differential equations. Hartle’s [143148] method computes rotating stars in this slow rotation approximation, and a review of slowly rotating models has been compiled by Datta [82]. Weber et al. [317Jump To The Next Citation Point319] also implement Hartle’s formalism to explore the rotational properties of four new EOSs.

Weber and Glendenning [318Jump To The Next Citation Point] improve on Hartle’s formalism in order to obtain a more accurate estimate of the angular velocity at the mass-shedding limit, but their models still show large discrepancies compared to corresponding models computed without the slow rotation approximation [260Jump To The Next Citation Point]. Thus, Hartle’s formalism is appropriate for typical pulsar (and most millisecond pulsar) rotational periods, but it is not the method of choice for computing models of rapidly rotating relativistic stars near the mass-shedding limit.

2.7.2 Butterworth and Ipser (BI)

The BI scheme [55Jump To The Next Citation Point] solves the four field equations following a Newton–Raphson-like linearization and iteration procedure. One starts with a nonrotating model and increases the angular velocity in small steps, treating a new rotating model as a linear perturbation of the previously computed rotating model. Each linearized field equation is discretized and the resulting linear system is solved. The four field equations and the hydrostationary equilibrium equation are solved separately and iteratively until convergence is achieved.

Space is truncated at a finite distance from the star and the boundary conditions there are imposed by expanding the metric potentials in powers of 1∕r. Angular derivatives are approximated by high-accuracy formulae and models with density discontinuities are treated specially at the surface. An equilibrium model is specified by fixing its rest mass and angular velocity.

The original BI code was used to construct uniform density models and polytropic models [5554]. Friedman et al. [113114Jump To The Next Citation Point] (FIP) extend the BI code to obtain a large number of rapidly rotating models based on a variety of realistic EOSs. Lattimer et al. [193Jump To The Next Citation Point] used a code that was also based on the BI scheme to construct rotating stars using “exotic” and schematic EOSs, including pion or kaon condensation and strange quark matter.

2.7.3 Komatsu, Eriguchi, and Hachisu (KEH)

In the KEH scheme [183Jump To The Next Citation Point184Jump To The Next Citation Point], the same set of field equations as in BI is used, but the three elliptic-type field equations are converted into integral equations using appropriate Green’s functions. The boundary conditions at large distance from the star are thus incorporated into the integral equations, but the region of integration is truncated at a finite distance from the star. The fourth field equation is an ordinary first order differential equation. The field equations and the equation of hydrostationary equilibrium are solved iteratively, fixing the maximum energy density and the ratio of the polar radius to the equatorial radius, until convergence is achieved. In [18318495Jump To The Next Citation Point] the original KEH code is used to construct uniformly and differentially rotating stars for both polytropic and realistic EOSs.

Cook, Shapiro, and Teukolsky (CST) improve on the KEH scheme by introducing a new radial variable that maps the semi-infinite region [0,∞ ) to the closed region [0,1]. In this way, the region of integration is not truncated and the model converges to a higher accuracy. Details of the code are presented in [69Jump To The Next Citation Point] and polytropic and realistic models are computed in [71Jump To The Next Citation Point] and [70Jump To The Next Citation Point].

Stergioulas and Friedman (SF) implement their own KEH code following the CST scheme. They improve on the accuracy of the code by a special treatment of the second order radial derivative that appears in the source term of the first order differential equation for one of the metric functions. This derivative was introducing a numerical error of 1 – 2% in the bulk properties of the most rapidly rotating stars computed in the original implementation of the KEH scheme. The SF code is presented in [294Jump To The Next Citation Point] and in [292Jump To The Next Citation Point]. It is available as a public domain code, named rns, and can be downloaded from [291Jump To The Next Citation Point].

2.7.4 Bonazzola et al. (BGSM)

In the BGSM scheme [47Jump To The Next Citation Point], the field equations are derived in the 3+1 formulation. All four chosen equations that describe the gravitational field are of elliptic type. This avoids the problem with the second order radial derivative in the source term of the ODE used in BI and KEH. The equations are solved using a spectral method, i.e., all functions are expanded in terms of trigonometric functions in both the angular and radial directions and a Fast Fourier Transform (FFT) is used to obtain coefficients. Outside the star a redefined radial variable is used, which maps infinity to a finite distance.

In [260Jump To The Next Citation Point261Jump To The Next Citation Point] the code is used to construct a large number of models based on recent EOSs. The accuracy of the computed models is estimated using two general relativistic virial identities, valid for general asymptotically flat spacetimes [132Jump To The Next Citation Point43Jump To The Next Citation Point] (see Section 2.7.7).

While the field equations used in the BI and KEH schemes assume a perfect fluid, isotropic stress-energy tensor, the BGSM formulation makes no assumption about the isotropy of Tab. Thus, the BGSM code can compute stars with a magnetic field, a solid crust, or a solid interior, and it can also be used to construct rotating boson stars.

2.7.5 Lorene/rotstar

Bonazzola et al. [46] have improved the BGSM spectral method by allowing for several domains of integration. One of the domain boundaries is chosen to coincide with the surface of the star and a regularization procedure is introduced for the divergent derivatives at the surface (that appear in the density field when stiff equations of state are used). This allows models to be computed that are nearly free of Gibbs phenomena at the surface. The same method is also suitable for constructing quasi-stationary models of binary neutron stars. The new method has been used in [133Jump To The Next Citation Point] for computing models of rapidly rotating strange stars and it has also been used in 3D computations of the onset of the viscosity-driven instability to bar-mode formation [129Jump To The Next Citation Point].

2.7.6 Ansorg et al. (AKM)

A new multi-domain spectral method has been introduced in [18Jump To The Next Citation Point19]. The method can use several domains inside the star, one for each possible phase transition. Surface-adapted coordinates are used and approximated by a two-dimensional Chebyshev expansion. Requiring transition conditions to be satisfied at the boundary of each domain, the field and fluid equations are solved as a free boundary value problem by a Newton–Raphson method, starting from an initial guess. The field equations are simplified by using a corotating reference frame. Applying this new method to the computation of rapidly rotating homogeneous relativistic stars, Ansorg et al. achieve near machine accuracy, except for configurations at the mass-shedding limit (see Section 2.7.8)! The code has been used in a systematic study of uniformly rotating homogeneous stars in general relativity [264].

2.7.7 The virial identities

Equilibrium configurations in Newtonian gravity satisfy the well-known virial relation

2T + 3(Γ − 1)U + W = 0. (28 )
This can be used to check the accuracy of computed numerical solutions. In general relativity, a different identity, valid for a stationary and axisymmetric spacetime, was found in [39Jump To The Next Citation Point]. More recently, two relativistic virial identities, valid for general asymptotically flat spacetimes, have been derived by Bonazzola and Gourgoulhon [132Jump To The Next Citation Point43Jump To The Next Citation Point]. The 3-dimensional virial identity (GRV3) [132] is the extension of the Newtonian virial identity (28View Equation) to general relativity. The 2-dimensional (GRV2) [43Jump To The Next Citation Point] virial identity is the generalization of the identity found in [39] (for axisymmetric spacetimes) to general asymptotically flat spacetimes. In [43], the Newtonian limit of GRV2, in axisymmetry, is also derived. Previously, such a Newtonian identity had only been known for spherical configurations [59].

The two virial identities are an important tool for checking the accuracy of numerical models and have been repeatedly used by several authors [47Jump To The Next Citation Point260261235Jump To The Next Citation Point18Jump To The Next Citation Point].

2.7.8 Direct comparison of numerical codes

The accuracy of the above numerical codes can be estimated, if one constructs exactly the same models with different codes and compares them directly. The first such comparison of rapidly rotating models constructed with the FIP and SF codes is presented by Stergioulas and Friedman in [294Jump To The Next Citation Point]. Rapidly rotating models constructed with several EOSs agree to 0.1 – 1.2% in the masses and radii and to better than 2% in any other quantity that was compared (angular velocity and momentum, central values of metric functions, etc.). This is a very satisfactory agreement, considering that the BI code was using relatively few grid points, due to limitations of computing power at the time of its implementation.

In [294Jump To The Next Citation Point], it is also shown that a large discrepancy between certain rapidly rotating models (constructed with the FIP and KEH codes) that was reported by Eriguchi et al. [95], resulted from the fact that Eriguchi et al. and FIP used different versions of a tabulated EOS.

Nozawa et al. [235Jump To The Next Citation Point] have completed an extensive direct comparison of the BGSM, SF, and the original KEH codes, using a large number of models and equations of state. More than twenty different quantities for each model are compared and the relative differences range from 10–3 to 10–4 or better, for smooth equations of state. The agreement is also excellent for soft polytropes. These checks show that all three codes are correct and successfully compute the desired models to an accuracy that depends on the number of grid points used to represent the spacetime.

If one makes the extreme assumption of uniform density, the agreement is at the level of 10–2. In the BGSM code this is due to the fact that the spectral expansion in terms of trigonometric functions cannot accurately represent functions with discontinuous first order derivatives at the surface of the star. In the KEH and SF codes, the three-point finite-difference formulae cannot accurately represent derivatives across the discontinuous surface of the star.

The accuracy of the three codes is also estimated by the use of the two virial identities. Overall, the BGSM and SF codes show a better and more consistent agreement than the KEH code with BGSM or SF. This is largely due to the fact that the KEH code does not integrate over the whole spacetime but within a finite region around the star, which introduces some error in the computed models.

A new direct comparison of different codes is presented by Ansorg et al. [18Jump To The Next Citation Point]. Their multi-domain spectral code is compared to the BGSM, KEH, and SF codes for a particular uniform density model of a rapidly rotating relativistic star. An extension of the detailed comparison in [18], which includes results obtained by the Lorene/rotstar code in [129Jump To The Next Citation Point] and by the SF code with higher resolution than the resolution used in [235Jump To The Next Citation Point], is shown in Table 2. The comparison confirms that the virial identity GRV3 is a good indicator for the accuracy of each code. For the particular model in Table 2, the AKM code achieves nearly double-precision accuracy, while the Lorene/rotstar code has a typical relative accuracy of 2 × 10–4 to 7 × 10–6 in various quantities. The SF code at high resolution comes close to the accuracy of the Lorene/rotstar code for this model. Lower accuracies are obtained with the SF, BGSM, and KEH codes at the resolutions used in [235Jump To The Next Citation Point].

The AKM code converges to machine accuracy when a large number of about 24 expansion coefficients are used at a high computational cost. With significantly fewer expansion coefficients (and comparable computational cost to the SF code at high resolution) the achieved accuracy is comparable to the accuracy of the Lorene/rotstar and SF codes. Moreover, at the mass-shedding limit, the accuracy of the AKM code reduces to about 5 digits (which is still highly accurate, of course), even with 24 expansion coefficients, due to the nonanalytic behaviour of the solution at the surface. Nevertheless, the AKM method represents a great achievement, as it is the first method to converge to machine accuracy when computing rapidly rotating stars in general relativity.

Going further   A review of spectral methods in general relativity can be found in [42]. A formulation for nonaxisymmetric, uniformly rotating equilibrium configurations in the second post-Newtonian approximation is presented in [22].


Table 2: Detailed comparison of the accuracy of different numerical codes in computing a rapidly rotating, uniform density model. The absolute value of the relative error in each quantity, compared to the AKM code, is shown for the numerical codes Lorene/rotstar, SF (at two resolutions), BGSM, and KEH (see text). The resolutions for the SF code are (angular × radial) grid points. See [235] for the definition of the various equilibrium quantities.
  AKM Lorene/ SF SF BGSM KEH
    rotstar (260 × 400) (70 × 200)    
¯pc 1.0
rp∕re 0.7 × 10–3
¯Ω 1.41170848318 × 10–6 × 10–4 × 10–3 × 10–2 × 10–2
¯ M 0.135798178809 × 10–4 × 10–5 × 10–3 × 10–3 × 10–2
¯ M0 0.186338658186 × 10–4 × 10–4 × 10–3 × 10–2 × 10–3
¯ Rcirc 0.345476187602 × 10–5 × 10–5 × 10–4 × 10–3 × 10–3
¯ J 0.0140585992949 × 10–5 × 10–4 × 10–4 × 10–2 × 10–2
Zp 1.70735395213 × 10–5 × 10–5 × 10–4 × 10–2 × 10–2
Zf eq –0.162534082217 × 10–4 × 10–3 × 10–2 × 10–2 × 10–2
b Zeq 11.3539142587 × 10–6 × 10–5 × 10–3 × 10–2 × 10–1
|GRV3 | × 10–13 × 10–6 × 10–5 × 10–3 × 10–3 × 10–1


  Go to previous page Go up Go to next page