Third-Order Corrections and Mass-Shedding Limit of Rotating Neutron Stars Computed By a Complex-Plane Strategy ()
1. Introduction
In a recent paper [1], we have applied the so-called “complex-plane strategy” (CPS), originally developed and used for computing classical polytropic models in rapid rotation (see e.g. [2,3]), to compute rapidly rotating neutron stars simulated by general-relativistic polytropic models, i.e. neutron stars obeying the well-known polytropic “equation of state” (EOS) (see e.g. [1], Section 2.1, Equations (5)-(9)). In this study, we implement Hartle’s perturbation method ([4-6]) in order to compute 1) the structure of a rotating neutron star up to terms of third order in the angular velocity, and 2) the mass-shedding limit, i.e. the angular velocity above which the gravitational attraction, compared to the centrifugal force, is not sufficient to keep matter bound to the surface ([7], Section 6.5.2; [8], Section 5.2.2; [9], Section 5). Here, we will try to avoid, as much as possible, rewriting and repeating issues from papers referred in the text; readers interested in this subject can find full details in the particular references.
The third-order corrections involve the functions and ([6], Equations (3.1) and (3.2); [10], Equations (1)-(5)). The function represents a third-order contribution to the angular momentum, moment of inertia, rotational kinetic energy, and gravitational potential energy. The function affects the massshedding velocity and, accordingly, the mass-shedding limit ([10], Section 2A; see also [11], Section 3). Both and contribute to the dragging of the inertial frames.
Unless stated otherwise, the physical quantities involved in this study are expressed in gravitational units (see e.g. [1], Section 1.2).
2. The Perturbed Metric and the Mass-Shedding Limit
The perturbation of the metric due to rotation, up to third-order terms in the angular velocity, is given in [10] (Equation (1); see also [12], Equation (25)). The functions and, involved in this perturbation, satisfy the second-order differential equations (A29) and (A41) of [10], respectively. The third-order correction to the angular momentum is given in [10] (Equation (28)); the third-order correction to the moment of inertia is given in [10] (Equation (31)); the third-order correction to the rotational kinetic energy is given in [12] (Equation (60)); and the third-order correction to the gravitational potential energy is given in [12] (Equation (61)).
Hartle’s perturbation method has been developed for treating slowly rotating neutron stars, in the sense that the angular velocity is considered to be small when compared to the “maximum angular velocity” , also called “critical angular velocity”, at which mass shedding occurs at the equator (see e.g. [5], Equation (2); see also [13], Equation (86). The quantities and are the mass and radius of the nonrotating model. corresponds to the Newtonian balance of centrifugal and gravitational forces. This Newtonian upper bound of the angular velocity is a significantly overestimated limit as far as neutron stars are concerned.
An absolute upper bound on neutron star uniform rotation is given by the “Keplerian angular velocity” or, equivalently, “mass-shedding limit”, which is the maximum allowed angular velocity of a particle in Keplerian orbit at the equator. If the angular velocity will be slightly greater than, then mass shedding would occur at the equator of a neutron star. So, is the relativistic analog of. Computing has attracted the attention of several investigators, and several numerical methods have been developed towards this task (for a discussion on this matter, see [13], Section 3.7).
3. The System of Differential Equations in the Framework of Hartle’s Perturbation Method
The main target of this investigation—which is in fact the continuation of the numerical treatment presented in [1]—is to solve numerically in the complex plane the system of the differential equations arising in the framework of Hartle’s perturbation method when taking into account up to third-order terms in the angular velocity and then to compute the mass-shedding limit.
The system under consideration consists of the following differential equations.
(01)-(02) The two first-order “Oppenheimer-Volkoff” (OV) equations of hydrostatic equilibrium and of massenergy ([1], Equations (8) and (9) with initial conditions (15) and (16)).
(03) The first-order differential equation governing the gravitational potential ([1], Equation (29) with boundary condition (30); see also the discussion following this equation).
(04)-(05) The second-order “frame dragging equation”, which can be reformulated as an equivalent system of two first-order differential equations for the angular velocity in the local inertial frame and its derivative ([1], Equations (32) and (33) with initial conditions (35a) and (35b); see also the discussion following Equation (33)).
(06)-(07) The two first-order differential equations for the mass perturbation function and the pressure perturbation function describing the spherical deformation of the star ([1], Equations (37) and (38) with initial conditions (39a) and (39b)).
(08)-(09) The two first-order differential equations for the functions and describing the quadrupole deformation of the star ([1], Equations (42) and (43) with initial conditions (44a) and (44b)).
(10)-(11) The two first-order homogeneous differential equations for the functions and ([1], Equations (47) and (48) with initial conditions (49a) and (49b)).
(12)-(15) The two second-order differential equations for the functions and ([10], Equations (A29) and (A41)) describing the third-order perturbative corrections, which can be reformulated as two equivalent systems, each consisting of two first-order differential equations ([10], Equations (A32a)-(A32b) with zero initial conditions, and (A43a)-(A43b) with zero initial conditions; details are given in [10], Appendix 3).
(16)-(19) The two second-order homogeneous differential equations for the functions and ([10], Equations (A29) and (A41) with,), which can be reformulated as two equivalent systems, each consisting of two first-order homogeneous differential equations ([10], Equations (A32a)-(A32b) with, and initial conditions (A34a)-(A34b), and (A43a)-(A43b) with and initial conditions (A45a)-(A45b), respectively).
Thus, from the point of view of numerical analysis, our task is to solve numerically in the complex plane an “initial value problem” (IVP) defined on the system of the 19 first-order “ordinary differential equations” (ODE, ODEs) referred above.
4. Solving the IVP with the ATOMFT System
To solve the IVP discussed in Section 3, we use the ATOMFT System ([14-20]; details on this software package are given in [1], Section 3.4). To avoid any singularities and/or indeterminate forms in the real axis, especially at the points and, we apply. The numerical integration is performed along a contour (i.e. a complex path), lying in the complex plane, being parallel to the real axis, and distancing a small imaginary distance from it, i.e. along the straight line
(1)
joining the points and in. For convenience, we use the abbreviations for the real part of the complex radial coordinate, and for the imaginary part of. In this study, for the parameters, , and we use the values, , and. For the initial condition of the complex rest-mass density,
(2)
we use the relation.
As ATOMFT advances the solution of the ODEs on successive points along the complex path, the numerical output of this method is used to construct the table (i.e. the table (57) of [1], now extended with all the third-order functions),
(3)
Then, using standard numerical methods, we interpolate the real variables of this table in terms of the real variable. Thus we construct by numerical interpolation an extended family of real-valued functions of a real variable, which can be used in the place of the original complex-valued functions of a complex variable for various purposes.
5. The Numerical Procedure after Having Solved the IVP
ATOMFT produces an output file containing the solution of the complex IVP. Next, a program takes control over this file and performs the following tasks.
1) Reads the output file, composes the extended table, and constructs by numerical interpolation the extended family of real-valued functions of the real variable.
2) Calculates the radius of the star as the first root of the algebraic equation ([1], Equation (22)), the gravitational mass ([1], Equation (24)), the baryonic mass ([1], Equation (108)), and the proper mass ([1], Equation (109)).
3) Calculates the angular momentum, moment of inertia, rotational kinetic energy, and gravitational potential energy ([12], Equations (21)-(23), and (24), respectively).
4) Calculates the changes in the radius, in the gravitational mass, in the baryonic mass, and in the proper mass ([12], Equation (43), (44), (48), and (49), respectively), all due to spherical deformation ([12], Section 3.3).
5) Calculates the third-order correction to the angular momentum [10], Equation (28)), the third-order correction to the moment of inertia ([12], Equation (59)), the third-order correction to the rotational kinetic energy ([12], Equation (60)), and the third-order correction to the gravitational potential energy ([12], Equation (61)).
Note that the real parts of the complex quantities referred in 2)-5) are interpreted as the corresponding familiar physical quantities.
6. Computing the Mass-Shedding Limit
In this study we use the procedure described in [10] (Section 2A) for calculating the mass-shedding limit. In particular, we consider a fluid element belonging to the star and located at the equatorial surface. This element moves with a velocity given by Equation (15) of [10]. On the other hand, the velocity of a particle moving on circular orbit in the corotating direction just outside the equator is given by Equation (17) of [10].
In the framework of, the velocities and are complex-valued functions of the complex variable, with real parts
(4)
To compute the mass-shedding limit, we need to construct sequences of models with constant baryonic mass and variable angular velocity (to be discussed at the end of this section). For a starting value
, the velocities and calculated at the equatorial radius are different,
(5)
For a second selected value, we have
(6)
We find, however, that; thus these two velocities converge to a limit as the selected gets increasing, i.e.. So, we gradually increase in order to calculate this limit. Eventually, for a value for which where is a prescribed tolerance, the fluid element on the surface will not be bound anymore and the model will start losing matter from its equator. The value is the massshedding limit within the prescribed accuracy, i.e..
In this study, we generate sequences of constant baryonic mass by using the procedure described in [10] (Sec IIB), with just a few slight changes due to the use of both the polytropic EOS and the method. This procedure is as follows.
Step 1. For each value of the polytropic index and for the central mass-energy density of the corresponding “nonrotating model of maximum mass”, we calculate the central rest-mass density of this model ([1], Equation (7)),
(7)
and we solve the OV equations ([1], Equations (8) and (9)) to find the radius of the star ([1], Equations (17) and (22)).
Step 2. We calculate the baryonic mass by solving numerically the integral (108) of [1]. Thus each model is identified by a constant baryonic mass
(8)
Step 3. For an assigned value, we solve the IVP of Sec. 3 for a sufficiently large number of values; in fact, for the corresponding values of the central rest-mass density obtained by Equation 7.
Step 4. For each rest-mass density, we calculate the correction to the baryonic mass ([12], Equation (48)) and the total baryonic mass of the corresponding model,
(9)
Step 5. Among all these rotating models, we select the one with baryonic mass,
(10)
Step 6. Having computed the model with angular velocity and total baryonic mass, we proceed with the calculation of the velocities and.
To achieve convergence, we repeat Steps 3-6 with gradually increasing angular velocities up to a value fulfilling the condition.
7. Numerical Results
We compute general-relativistic polytropic models of maximum mass simulating neutron stars in rapid uniform rotation (on models of maximum mass, see [9], Section 4; [12], Section 6; [13], Section 5.2). The procedure followed to compute such models is described in [1] (Section 4). We then compute the corresponding uniformly rotating models with angular velocities equal to their Keplerian angular velocities; in this study, we take as “reference values” of those computed by the well-known RNS package [21].
Next, we calculate the angular momentum, moment of inertia, rotational kinetic energy, gravitational potential energy, mass spherical deformation, baryonic mass spherical deformation, proper mass spherical deformation, and radius spherical deformation. Table 5 lists the percentage differences, , between the results of the present investigation, , and the corresponding ones, , in Tables 3 and 7 of [12]. Furthermore, we calculate the third-order corrections, , and. Knowing the real parts of these quantities, we then find the corresponding total quantities,
(11)
(12)
(13)
(14)
We give indicative numerical results in Tables 1-4 for, respectively. Comparison of the results in Tables 2 and 4 with corresponding ones in Tables 3 and 7 of [12] shows that they are compatible numerical results. The corresponding percentage differences are listed in Table 5. Note that the numerical results for the model, listed in Table 3, have been computed with different to that used in Table 5 of [12]; thus these two sets of results are not comparable.
Finally, Tables 6-8 give the values of the constant baryonic mass, equatorial radius, velocity, mass-shedding limit, massshedding limit calculated by the RNS package, and the percentage differences between the results of and.
Table 1. Numerical results for the n = 1.0 polytropic model of maximum mass; uniform rotation with Keplerian angular velocity, ΩK = 3.333 × 10−7 cm−1.
Table 2. Numerical results for the n = 1.5 polytropic model of maximum mass; uniform rotation with Keplerian angular velocity, ΩK = 2.658 × 10−7 cm−1.
Table 3. Numerical results for the n = 2.0 polytropic model of maximum mass; uniform rotation with Keplerian angular velocity, ΩK = 1.844 × 10−7 cm−1.
Table 4. Numerical results for the n = 2.5 polytropic model of maximum mass; uniform rotation with Keplerian angular velocity, ΩK = 1.141 × 10−7 cm−1.
Table 5. Percentage dierences between results of this investigation and the corresponding ones of Tables 3 and 7 in [12], for the n = 1.5 and n = 2.5 polytropic models.
Table 6. Mass-shedding limit for the n = 1.0 polytropic model.
Table 7. Mass-shedding limit for the n = 1.5 polytropic model.
Table 8. Mass-shedding limit for the n = 2.5 polytropic model.
8. Conclusion
The numerical results of Tables 1-5 show that gives results compatible with those presented in [12], on the full extent of Hartle’s perturbation theory, i.e. up to terms of third order in the angular velocity. It seems therefore that is an accurate and reliable numerical method for treating rapidly rotating neutron stars. On the other hand, for the mass-shedding limit(s) computed in this study, we have found that the differences with respect to corresponding results of a nonperturbative method like RNS, listed in Tables 6-8, are about 20%. As discussed in [10] (Section 3), such differences are eventually expected within the framework of Hartle’s perturbation method.
9. Acknowledgements
The authors acknowledge the use of the ATOMFT System.