Quasi-radial modes of pulsating neutron stars: numerical results for general-relativistic rigidly rotating polytropic models

In this paper we compute general-relativistic polytropic models simulating rigidly rotating, pulsating neutron stars. These relativistic compact objects, with a radius of $\sim 10 \, \mathrm{km}$ and mass between $\sim 1.4$ and $3.2$ solar masses, are closely related to pulsars. We emphasize on computing the change in the pulsation eigenfrequencies owing to a rigid rotation, which, in turn, is a decisive issue for studying stability of such objects. In our computations, we keep rotational perturbation terms of up to second order in the angular velocity.


Introduction
According to "Hartle's perturbation method" (also called "Hartle-Thorne perturbation method"; [1], [2]), we treat rotating relativistic neutron stars as perturbative solutions of the Einstein's field equations, which describe a static spherically symmetric relativistic model. Our aim here is to compute, based on the work of Hartle et al. [3] and also of Hartle and Friedmann [4], the zerothand second-order eigenfrequencies of the lowest three modes of radial pulsation for general-relativistic polytropic models. The zeroth-order eigenfrequencies, [σ 2 ] (0) , are the eigenfrequencies of the nonrotating model, while the secondorder ones, [σ 2 ] (2) , are the rotation-induced changes in the eigenfrequencies.
In this study, we use "gravitational units" (abbreviated "gu"; see e.g. [5] and [6]), also called "geometrized units". In these units the speed of light, c, and the gravitational constant, G, are equal to unity, c gu = G gu = 1, and the length is the unique "base unit", measured in centimeters (abbreviated "cm"). On the other hand, the well-known "centimeter-gram-second units" (abbreviated "cgs") have three base units: the length measured in cm, the mass measured in grams (abbreviated "g"), and the time measured in seconds (abbreviated "s").
To treat numerically the system (4)-(5), we need an equation relating P to E, P = P (E), that is, an "equation of state" (EOS). In this study, we consider the polytropic EOS (see e.g. [7], Equation (3); see also [8], Section II) The parameters K and n are the so-called "polyropic constant" and "polytropic index", respectively. The so-called "adiabatic index", Γ, is defined as (cf. [3], and, in general, is associated with perturbations about equilibrium under constant entropy. On the other hand, the so-called "adiabatic index associated with the equation of state", γ, is defined as (cf. [3], Equation (3.12)) In the polytropic EOS (7), ̺(r) is the rest-mass density, related to the massenergy density via the equation ( [6], Equation (9)) To solve the system (4)-(5), we write the two equations as ( [6], Equations (8) and (9), respectively) The initial conditions are now To complete the study of the nonrotating model, we need to solve the differential equation for the gravitational potential Φ ( [1], Equation (29b)) obeying the boundary condition at the surface of the star (cf. [6], Equation (30))

Rigid Rotation
When a star is rotating, its shape deviates from sphericity. Furthermore, when the star is in equilibrium state, there is a balance between the pressure forces, the gravitational forces, and the centrifugal forces. Assuming that the star is rotating rigidly and slowly (see however the starting remark in Section 7), the perturbed metric is given by ( [2], Equation (4)) where P l = P l (cos θ) are Legendre polynomials; the perturbation functions m 0 , h 0 , m 2 , h 2 , and υ 2 are radial functions proportional to Ω 2 ; the centrifugal forces depend on the angular velocity Ω relative to a distant observer, as well as on the angular velocity ̟ relative to the local inertial frame, these two angular velocities connected via the relation ( [2], Equation (6)) where ω is the angular velocity of the local inertial frame. To calculate ̟, we solve the so-called "frame-dragging equation" ( [2], Eq. (9)) where ( [2], Equation (10)) Outside the star, ̟ takes the form ( [2], Equation (13b)) where J is the total angular momentum of the star, given by ( [2], Equation (13a)) where R is the radius of the star. To solve Equation (18), we integrate from the center of the star outwards, imposing the initial conditions where the constant ̟ c is chosen arbitrarily. Since, due to Equation (20), the angular velocity Ω arb corresponding to the arbitrarily chosen initial value ̟ c is equal to and since, in general, the prescribed angular velocity Ω is different to Ω arb , we must rescale the solution ̟(r) of Equation (18) in order for the prescribed value Ω to be applied to our model (for details on this matter, see [5], Section 5.2),

Spherical Deformation
The spherical deformations due to rigid rotation are the mass perturbation function, m 0 , and the pressure perturbation function, p 0 . They are calculated by integrating the "l = 0 equations of hydrostatic equilibrium" ([1], Section VII), which become after a long algebra ( [2], Equations (15a) and (15b), respectively) with initial conditions m 0 = 0 and p 0 = 0 at r = 0. The perturbation function h 0 (to be used in Section 5), involved in the metric (16), is defined outside the star as ( [2], Equation (17a)) and inside the star as ( [2], Equation (17b)) The constant h 0c is determined so that h 0 be continuous across the surface (see e.g. [5], Equation (47)) In these equations, M is the mass-energy of the star (also called gravitational mass); and δM is the increase in the gravitational mass due to spherical deformation, given by ( [5], Equation (43))

Quasi-Radial Pulsation
When a star is rotating rigidly with a small angular velocity Ω (see however the starting remark in Section 7), the squared eigenfrequencies σ 2 of its pulsation modes can be expanded in powers of Ω ( [4], Equation (2.1)), where the superscript "(0)" denotes terms of zeroth order in Ω, and the superscript "(2)" terms of second order in Ω. The aim of the present study is to compute the zeroth-and the second-order eigenfrequencies of the lowest three pulsation modes.
To compute the eigenvalue(s) [σ 2 ] (0) , we work as follows. We start the numerical integration for a trial value σ 2 and initial conditions as above. We integrate towards the surface and then check if the resulting solution U satisfies the boundary condition Γ P U (R) ′ = 0. From the point of view of numerical analysis, this boundary condition can be treated as an algebraic equation of the form f σ 2 = 0; thus our task, to compute the root(s) [σ 2 ] (0) of this equation, can be achieved by a standard numerical method (e.g. the bisection method).

The Computations
In this study, to compute nonrotating models (Section 2), rigid rotations (Section 3), and spherical deformations (Section 4), we use the corresponding numerical methods described in very detail in [5] (Sections 5.1 and 5.2). We then combine these methods with the numerical framework described in Section 5 for computing the zeroth-and second-order eigenfrequencies of pulsation.
To implement all the required methods, and thus to compute the results presented here, we have written and used a Mathematica program.

Numerical Results and Discussion
To begin with, it is worth remarking here that, as it has been verified by several authors (see e.g. [10], Sections 4 and 7; see also [11], Section 5.3; for differentially rotating neutron stars, see [5], Section 6), Hartle's perturbation method gives remarkably accurate results even when applied to rapidly rotating neutron stars, although this method has been developed as a slow-rotation perturbation method.
As discussed in [4] (Section I), pulsars are identified as rotating neutron stars and, therefore, there is a strong interest in studying the influence of a rigid rotation on the properties of such relativistic objects. In particular, it is of great interest to compute the pulsation frequencies of the quasi-radial modes (i.e. modes which would be radial in the absence of rotation) for several models and thus to have a measure of the effect of general relativity on these frequencies. To that purpose, we have computed, and present in this section, relevant numerical results.
Regarding our computations, we resolve four nonrotating general-relativistic polytropic models with polytropic indices n = 1.0, 1.5, 2.0, and 2.5. Each model is resolved for five central mass-energy densities: E c = 10 13 , 3.16 × 10 13 , 10 14 , 3.16 × 10 14 , and 10 15 cgs. These values have been chosen to be below and relatively close to the "maximum-mass densities" of the corresponding models, being in fact the more interesting ones when considering neutron stars. It is worth mentioning here that the total mass M of a model, treated as a function of the central density E c , M = M (E c ), obtains a maximum value M max for a specific value E max c ; such a model is called "maximum-mass model", and the central density of this model is called "maximum-mass density". The maximummass densities of our models, computed by a method described in [6] (Section 4), are E max c = 3.793 × 10 15 , 4.890 × 10 15 , 4.656 × 10 15 , and 3.489 × 10 15 cgs, respectively ( [6], Tables 2, 3, 4, and 5, respectively). All models, studied here, have E max c < 5 × 10 15 cgs. Accordingly, in our computations the sequence of central mass-energy densities is terminated at E c = 10 15 cgs.
Next, each density case is resolved for the three lowest modes of pulsation: Mode 0, Mode 1, and Mode 2. For each mode, we compute a rigidly rotating configuration with angular velocity equal to the corresponding Keplerian angular velocity, Ω K . Hartle's perturbation method uses proper expansions in the rotation parameter ǫ = Ω/Ω max , where Ω max = G M/R 3 is the angular velocity for which mass shedding starts occuring at the equator of a star. Thus Ω max describes the Newtonian balance between centrifugal and gravitational forces. However, this Newtonian upper bound appears to be a rather overestimated limit for neutron stars. For such relativistic objects, the appropriate upper bound is Ω K . Hence, if the angular velocity is slightly greater than Ω K , then mass shedding occurs at the equator of a neutron star. Ω K can be computed by several methods (for a discussion on this matter, see [11]; for a detailed description of such a method, see [12]; see also [13], and references therein, for results concerning general-relativistic polytropic models). In this study, the Keplerian angular velocities given in Tables 1, 2, 3, and 4 have been computed by using the well-known "Rotating Neutron Stars Package" (RNS) [14].
Furthermore, we assume that the adiabatic indices Γ (Equation (8)) and γ (Equation (9)) do coincide for the polytropic models under consideration, In addition, a second step towards simplification is to assume, as several authors do (see e.g. [4]; for a different view, on the other hand, see e.g. [9]), that Γ is associated with the polytropic index n via the polytropic relation (cf. [6], Equation (5)) In all tables of the present paper, parenthesized (positive or negative) integers following numerical values denote powers of ten. For example, the entry 3.16(13) is equal to 3.16 × 10 13 and the entry 3.51(−15) is equal to 3.51 × 10 −15 .
For the five density cases with polytropic index n = 1.5, we can compare our results in Table 5 with respective results in Table 3 of [4]. To the purpose of such comparisons, we have written Table 5 in exactly the same format with that of Table 3 of [4]. We find excellent agreement between respective results, except for two eigenfrequencies [σ 2 ] (2) belonging to Mode 2. The first one arises when E c = 10 14 cgs and leads to a difference ∼ 4% (our result is "-15.95", while the result in [4] is "-16.6"); and the second one occurs when E c = 10 15 cgs and leads to a difference ∼ 3% (our result is "-14.15", while that in [4] is "-13.7"). Since all other results do almost coincide, it seems that these two differences are of rather minor significance.
Our main remarks on the numerical results presented in Tables 5, 6, 7, and 8 have as follows. First, in all cases examined, the eigenfrequencies [σ 2 ] (0) and [σ 2 ] (2) increase in absolute value with the central density. Equivalently, since (for the cases examined) increasing central density implies increasing gravitational mass, all the eigenfrequencies increase in absolute value with the gravitational mass. In addition, since increasing central density does also imply increasing Keplerian angular velocity for a rotating configuration, all the eigenfrequencies [σ 2 ] (2) increase in absolute value with the Keplerian angular velocity.
Second, all eigenfrequencies [σ 2 ] (0) are positive, thus representing stable nonrotating pulsating configurations. Likewise, the eigenfrequencies [σ 2 ] (2) are positive for Mode 0 and for the soft polytropic EOSs n = 2.5 and n = 2.0. Consequently, the effect of rotation is to stabilize such rotationally perturbed configurations. On the other hand, the eigenfrequencies [σ 2 ] (2) are negative for Modes 1 and 2 of the soft EOSs n = 2.5 and n = 2.0, as well as for the three modes of the stiff EOSs n = 1.5 and n = 1.0, thus turning to destabilize the corresponding rotating configurations. It is worth mentioning here that, among the members of a collection of EOSs, the EOS deriving the larger value P for a given E is the "stiffest" EOS in the collection; while the EOS leading to the smaller value P for the same E is the "softest" EOS in this collection. Note that, for increasing polytropic index n, the polytropic EOSs are getting softer; equivalently, for decreasing n, the polytropic EOSs become stiffer.
Finally, in all cases examined, the zeroth-order eigenfrequencies [σ 2 ] (0) are ∼one order of magnitude greater than the respective squared Keplerian angular velocities Ω 2 K ; this inequality is in fact a necessary condition for the perturbation theory developed in [3] and [4] to be valid, as discussed in detail in [4] (Section II). Consequently, all the cases examined lie within the domain of applicability of the theory, and, especially, of Equation (31) for computing the second-order eigenfrequencies [σ 2 ] (2) which represent the rotationally induced changes in the pulsation eigenfrequencies.     (14) 2.68(4) 1.29(6) 6.38(-8) Table 5: Eigenfrequencies of the lowest three modes of a nonrotating and a rotating general-relativistic polytropic model with polytropic index n = 1.5, and polytropic constant and Keplerian angular velocities as in Table 2.
Mode 0  Table 6: Eigenfrequencies of the lowest three modes of a nonrotating and a rotating general-relativistic polytropic model with polytropic index n = 1.0, and polytropic constant and Keplerian angular velocities as in Table 1.  (13) Table 7: Eigenfrequencies of the lowest three modes of a nonrotating and a rotating general-relativistic polytropic model with polytropic index n = 2.0, and polytropic constant and Keplerian angular velocities as in Table 3.  (13) Table 8: Eigenfrequencies of the lowest three modes of a nonrotating and a rotating general-relativistic polytropic model with polytropic index n = 2.5, and polytropic constant and Keplerian angular velocities as in Table 4.  (13)