Quasi-equilibrium models for triaxially deformed rotating compact stars
Quasi-equilibrium models of rapidly rotating triaxially deformed stars are computed in general relativistic gravity, assuming a conformally flat spatial geometry (Isenberg-Wilson-Mathews formulation) and a polytropic equation of state. Highly deformed solutions are calculated on the initial slice covered by spherical coordinate grids, centered at the source, in all angular directions up to a large truncation radius. Constant rest mass sequences are calculated from nearly axisymmetric to maximally deformed triaxial configurations. Selected parameters are to model (proto-) neutron stars; the compactness is for polytropic index and for . We confirmed that the triaxial solutions exist for these parameters as in the case of Newtonian polytropes. However, it is also found that the triaxial sequences become shorter for higher compactness, and those may disappear at a certain large compactness for the case. In the scenario of the contraction of proto-neutron stars being subject to strong viscosity and rapid cooling, it is plausible that, once the viscosity driven secular instability sets in during the contraction, the proto-neutron stars are always maximally deformed triaxial configurations, as long as the compactness and the equation of state parameters allow such triaxial sequences. Detection of gravitational waves from such sources may be used as another probe for the nuclear equation of state.
Rapidly rotating compact objects are expected to be formed as new born neutron stars after stellar core collapses, or as differentially rotating hypermassive neutron stars after binary neutron star mergers. Accretion onto neutron stars in X-ray binaries can also lead to rapid rotation. All of these have been extensively studied as strong sources of gravitational waves for the ground based laser interferometers LIGO/GEO600/VIRGO/TAMA (See e.g. reviews () and references therein).
Classical models of rotating stars are a class of ellipsoidal figures of equilibrium; self-gravitating rotating incompressible fluids in Newtonian gravity. Such solutions include sequences of axisymmetric Maclaurin ellipsoids, or non-axisymmetic Jacobi, Dedekind, or Riemann S-type ellipsoids Ch69 (). These models are used to study the secular evolutions of rapidly rotating stars due to the viscosity and the radiation back-reaction of gravitational waves ellipsoidal (). Lai and Shapiro LS95 () have developed an ellipsoidal approximation to the rotating polytropes, and applied the model to clarify the secular evolution of rapidly rotating neutron stars in detail, and more recently focused on the viscosity driven secular instability STU04 ().
As discussed in reviews (); ellipsoidal (); LS95 (), and shown by a number of numerical simulations of rapidly rotating compact stars, core collapse, and binary neutron star mergers SimDI (); SimCC (); SimBNS (), a transient triaxially deformed compact object may survive within a secular time scale. In this paper, we consider uniformly rotating models of such triaxially deformed compact objects, an extension of the Jacobi ellipsoid in general relativity. In Newtonian gravity, such solutions exist for rotating polytropes with polytropic index James64 (). Here, the polytropic equation of state (EOS) relates the pressure with the baryon rest mass density .
In general relativity, such configurations are not in equilibrium due to the back-reaction of gravitational radiation. However as in the case of quasi-equilibrium initial data of binary neutron stars, triaxially deformed uniformly rotating stars are in quasi-equilibrium, as long as the gravitational luminosity is small enough that the energy radiated away within a rotational period is small compared to the binding energy of the star, which is always the case, and if the viscosity is strong enough for the flow field to become uniformly rotating during the evolution. Therefore, as an important application, a sequence of uniformly rotating quasi-equilibrium solutions may model a secular evolution from the proto-neutron star to the neutron star in the strong viscosity limit, and each solution may serve as the initial data for the general relativistic hydrodynamic simulations of such objects.
Models of rapidly rotating neutron stars have been extensively studied as stationary, axisymmetric, perfect-fluid spacetimes Nick03 (), while less attention has been paid to uniformly rotating triaxial solutions, not only because they are not exact equilibria due to gravitational radiation reaction, but also because, in early calculations FIP86 (), models of the EOS above nuclear density did not allow large enough values of where a triaxial sequence is expected to bifurcate from an axisymmetric sequence. Here, is the ratio of kinetic energy to gravitational potential energy . However, as seen, for example, in Haensel (), a value (an effective adiabatic index ) may be possible for recent models of EOS for high density nuclear matter, above , and, as we will see below, triaxial quasi-equilibrium solutions do exist even in strong gravity for relatively small polytropic indexes, such as or , as in the Newtonian case.
About a decade ago, Nozawa succeeded in computing uniformly rotating triaxial polytropes in general relativistic gravity in his thesis Nozawa (), although his calculations were limited by the computational resources to low resolutions. A few studies approximating the fluid as an ellipsoidal configuration in general relativistic gravity have been made GRellpisoid (), and perturbative analyses have located the bifurcation point suggesting the existence of solutions having triaxial bar-mode deformations BFGall (); Dorota (); Saijo06 (); SL96 (); YE97 (); SF98 (). Our computations of triaxially deformed stars can be used to locate the instability points on the axisymmetric sequence as well as to estimate the gravitational wave amplitude and luminosity from such objects.
In this paper, we present our first results on triaxial configurations of rapidly rotating general relativistic stars as models of neutron stars in extreme rotation. We assume a conformally flat spatial slice, and solve the constraints and spatial trace of the Einstein equation (Isenberg-Wilson-Mathews (IWM) formulation) ISEN78 (); WM89 () 111The validity of the IWM formulation for axisymmetric configurations is discussed in CST96 (); WM_TEST (). This is different from the formulation used by Nozawa Nozawa () in which the line element is chosen to be the same form as that of stationary axisymmetric spacetime, but an azimuthal dependence is allowed (see usui99 () for the same formulation). The formulations and the code are described in the next section. For testing the code, selected axisymmetric solutions are compared with the results in the literature, and the bifurcation points of axisymmetric and triaxial sequences in weak gravity are examined. Then, we present the results of deformation sequences of constant rest mass systematically in the range of two parameters, the compactness and the polytropic index , appropriate for realistic neutron stars. Applications of such triaxial solutions in the contraction of a newly born proto-neutron star are briefly discussed in the final section. Throughout the paper, we use units such that . For our tensor notation, we adopt the use of Greek letters for spacetime indices, and Latin letters for spatial indices.
Ii Formulation and numerical method
ii.1 IWM formulation
The IWM formulation for computing spatially conformally flat initial data we use for computing non-axisymmetric quasi-equilibrium of rotating compact star is briefly described. The spacetime is foliated by the family of spacelike hypersurfaces . The future-pointing timelike normal to is related to the timelike vector , which is tangent to a curve , by , where is the lapse, and the shift which satisfies . Restricting the projection tensor , a spatial metric is defined on . In the IWM formulation, the spatial metric is chosen to be conformally flat, , where is a flat metric on each slice, and is a conformal factor. Then the metric takes the form
in a chart . The extrinsic curvature of the foliation is defined by
In the IWM formulation, the Einstein equation is decomposed with respect to the foliation, and the following 5 components Eqs.(3)-(5) are chosen to be solved for the five metric coefficients on the initial slice :
where the first and second equations are the constraints. These equations are written in the form of elliptic equations with the non-linear source terms, respectively,
where is the trace of , its tracefree part, the conformally weighted shift defined by and , is the flat Laplacian and is the covariant derivative with respect to the flat three-metric . The source terms of matter are defined by , , and . Expressions of the sources in terms of the metric potentials and fluid variables are given in Appendix A.
We choose a maximally embedded slice . Because the spatial metric is conformally flat, does not involve time derivatives of the spatial metric,
The elliptic equations and are separated, and the former is substituted to the latter:
The potentials are solved for simultaneously, and the shift is reconstructed from Eq.(15).
ii.2 Formulation for a relativistic fluid in equilibrium
A perfect fluid is described by the stress-energy tensor
where is the 4-velocity of the fluid, its pressure, and the energy density. As a consequence of the Bianchi identity, the stress-energy tensor is covariantly conserved:
When the fluid is close to equilibrium, one can obtain a simpler set of equations. Introducing the specific enthalpy defined by , where is the baryon rest mass density, the left hand side of Eq. (20) can be written
where is the specific entropy. In the derivation, the local first law of thermodynamics was used. In local thermodynamic equilibrium, one can also assume the conservation of baryon mass,
Consequently, the conservation of specific entropy along the fluid world line,
and, the relativistic Euler equations,
are obtained. Assuming the flow field to be isentropic everywhere inside the neutron star matter, , we have a one-parameter equations of state (EOS) .
We assume a stationary state in the rotating frame for the fluid source. Imposing symmetry along the helical vector where is a constant angular velocity of a rotating frame, we have
with interpreted as the scalar . For a corotational flow, , the rest mass conservation becomes trivial, and the relativistic Euler equation has the first integral
where is the injection energy. From the normalization of the four velocity , one obtains
As a first step in the calculation of a highly deformed triaxial compact star, we assume a simple polytropic EOS,
where is a constant, and is the polytropic index. Then is related to by
We also refer to the polytropic exponent defined by .
ii.3 Numerical computation
The Poisson solver and the iteration scheme used to solve the system of elliptic equations with non-linear source terms, are similar to the ones used in a previously developed initial data code for binary black holes and neutron stars uryu (); tsokaros (). However, the code itself has been completely rewritten, so that further extensions can be incorporated easily. One of the revisions of the code is that no symmetry is a priori assumed on the spatial slice ; that is, the spherical coordinate grids centered at the source cover all angular directions, up to a certain large truncation radius. Hence, for example, asymmetric magnetic fields may be later included without major modifications to the code. Computation of binary solutions using the same coordinate setup is also possible. The other major change is a simpler, more robust choice of finite differencing. In this section, we briefly describe the necessary steps for constructing the code, which are 1. Spherical coordinates and the length scale, 2. Summary of variables and equations for coding, 3. Poisson solver, 4. Grid spacing, 5. Finite differencing and iteration, 6. Computation of a sequence of solutions.
ii.3.1 Spherical coordinates and the length scale
The slice is covered by a spherical coordinate patch . For a single star calculation, the radial coordinate extends from the center of the star to the asymptotic radius , where is the radius of the neutron star along the semi-major axis, defined by the and lines. We also refer to Cartesian coordinates whose positive , and directions are along , and , respectively.
The quantity is introduced as an additional parameter in the formulation used in our code, normalizing the radial coordinate as
For a polytropic EOS, one can rescale the length using the polytropic constant as , or simply setting (see e.g. uryu ()). As a result, we have three parameters in our formulation.
Furthermore, we introduce surface fitted coordinates on which the fluid variables are defined. Assuming that the surface of a neutron star can be described by a function of the angular coordinate , the surface fitted coordinates are defined by
where is defined in a region .
ii.3.2 Summary of variables and equations for coding
As mentioned in Sec. II.1, the field equations (11), (13), (17) and (18) are solved for the metric potentials and, as in Sec. II.2 and II.3.1, a comoving fluid in equilibrium is characterized by one fluid variable, which is chosen to be the relativistic enthalpy , and three parameters .
The field equations are normalized to have the following form; representing each of the metric potentials by ,
where the flat Laplacian corresponds now to the normalized coordinate . The source term includes the metric potentials and their derivatives, while also includes the matter variables and the parameters , while the dependence on the length scale is explicitly separated in Eq.(32).
The fluid variable is determined by Eq.(26) coupled to the EOS (28), and the relations (27) and (29). The three parameters are determined by the following three quantities: the surface radii along two of the three semi-major axes, and the value of the central density. These quantities are used to impose three conditions on Eq.(26), which are solved with respect to the three parameters in each iteration cycle.
ii.3.3 Poisson solver
where and are positions, . We choose the Green function without boundary,
and perform a multipole expansion in associated Legendre functions,
where the radial Green function is defined by
and the coefficients are equal to for , and for .
ii.3.4 Grid spacing
The field equations in the integral form (33) are discretized on the spherical grids, and iterated until convergence is achieved. Our code allows us to use any non-equidistant grid spacing in all the spatial coordinates, , , , and . The radial grid points are equidistant in the region and non-equidistant in as follows: writing , we have
where the constant is determined from the relation
We choose equidistant grid spacing for and , that is, , and . Our notations for the grid points are summarized in Table 1.
|:||Radial coordinate where the grid starts.|
|:||Radial coordinate where the grid ends.|
|:||Radial coordinate between and where the|
|grid changes from equidistant to non-equidistant.|
|:||Total number of intervals between and .|
|:||Number of intervals between and .|
|:||Total number of intervals for .|
|:||Total number of intervals for .|
|:||Total number of intervals for .|
ii.3.5 Finite differencing and iteration
For the numerical integration of Eq.(33) we select the mid-point rule. Accordingly, source terms are evaluated at the middle of successive grid points. The linear interpolation formula and the second order Lagrange formula are applied for computing the source term fields and their derivatives respectively, at the mid-points of the , and grids.
The reason for selecting a rather low (second) order finite difference scheme is the following: When the field quantities vary rapidly, such as at a density discontinuity in a neutron star, higher order interpolating formulas as well as finite difference formulas tend to overshoot near the region, and may cause a non-convergent iteration. To overcome this behavior, one may either (i) separate the computing regions at the discontinuity, or (ii) use lower order polynomial approximations. With the first approach, pseudo-spectral methods have been successfully implemented by Meudon () and achieved an evanescent error. We select the second idea to keep the code as simple and flexible as possible, and improve the accuracy by simply increasing the number of grid points.
where the softening parameter is chosen to be for accelerated convergence. Then we check the relative difference of successive cycles
as a criteria for the convergence. We typically stop the iteration when this quantity becomes less than .
ii.3.6 Constructing a sequence of solutions
We compute constant rest mass sequences of isentropic equilibrium solutions in the IWM formulation. Constant entropy is modeled by setting the parameter in the EOS to a constant. In a quasi-equilibrium evolution of a rotating star, the angular velocity remains constant when the viscosity of matter is dominant. Then the solution sequence approximately models an evolution, with an error that includes neglecting the increase of entropy due to viscosity.
Each solution is computed by setting the central value of to , the ratio of the stellar radii along the and axes for the axisymmetric configuration, and the ratio of the and axes for the triaxial configuration. To compute a constant rest mass sequence, one iterates changing until converges to the specified value. We use a discrete Newton-Raphson iteration for the rest mass.
A sequence of rotating star solutions with a certain EOS is labeled by the compactness of a non-rotating spherical star having the same rest mass . We denote the gravitational mass of this spherical star by , and the compactness by . This labeling for each sequence is possible as long as it is a normal sequence that has a stable spherical star in the limit that the angular velocity goes to zero, which is not the case for a supermassive sequence. In the following, we focus on the normal sequences whose compactness is close to its value for a neutron star, around .
Formulas used for computing the rest mass , ADM mass , Komar mass , as well as the total angular momentum are presented in Appendix A. For polytropic EOS, one can normalize these quantities by a certain power of the polytropic constant , as shown in the same Appendix. Hence we choose units to present solution sequences.
A sequence of solutions with constant rest mass is considered as an evolutionary track of adiabatic changes in quasi-equilibrium. Under this assumption, the solutions in each sequence are parameterized by the angular velocity , and the first-law relation
is satisfied, as proved in FUS ().
Iii Code test
iii.1 Axisymmetric solutions
Axisymmetric solutions calculated by our new code are compared with the results in the literature CST96 (); NSGE98 (). We show the results of comparisons for models presented in Table I of Cook, Shapiro and Teukolsky CST96 () (hereafter CST), which correspond to a solution sequence with constant rest mass for the case with the polytropic index . This value of is close to the maximum rest mass of a non-rotating spherical solution; the gravitational mass and the compactness of the same non-rotating solution are and . In Table 3, selected solutions calculated with the highest resolution I-5 in Table 2 are compared with the results shown in Table I of CST96 (). Fractional errors in any quantities are less than .
As discussed in the section II.3.5, our choice of finite difference approximations is second order. The rate of convergence of our code is checked using different resolutions, whose setups of coordinate grids are shown in Table 2. The grid spacing of each coordinate is proportionally scaled as , , , , from type I-1 to I-5. Here, we show the results of the convergence test with respect to the resolutions, fixing the maximum number of multipoles as shown in Table 2. 222For the convergence tests with respect to the order of the Legendre expansion, see tsokaros ().
When a sufficient number of multipoles is kept, the differences between numerically computed quantities with different resolutions and their exact value are written
where () is a quantity computed using one of the resolution types in Table 2, is its exact value, represents the grid spacing associated with the type setup, and is a constant. Then, keeping the leading term, differences between different resolutions become
To see the order in a log-log plot, we select the combinations of different resolutions that give the same ratio , and in our choice, these are , , and . In Fig.1, these combinations normalized by of selected quantities are plotted against the grid spacing, where represents the grid spacing in arbitrary units. It is clearly seen that the local quantities, here and , converge to , and that integral (global) quantities also approach second order convergence as the resolution increases.
We also checked the convergence with the different sets of resolutions type II-1, 2 and 3. These setups have fewer grid points in the angular coordinates and . We found that the highest resolution type II-3 agrees well with the results of the higher resolutions of type I-4 or 5 for the axisymmetric solutions.
iii.2 Triaxial solutions with
We calculate triaxial sequences for small compactness , which is in the Newtonian regime, to check the value of at the bifurcation point of the triaxial sequence from the axisymmetric sequence. The triaxial and axisymmetric solution sequences for and are plotted in Fig. 2. Extrapolating the triaxial sequences to corresponding axisymmetric sequences, values at the branch points are determined approximately as and for and , respectively. This value may be compared with the Newtonian results such as the ellipsoidal approximation for STU04 ().
Iv Triaxial solutions
iv.1 Accuracy of the sequences of solutions
Triaxially deformed solutions are calculated for selected values of the polytropic index, and . Models with are calculated for and with for .
We noticed that it is necessary to increase the numbers of grid points as much as in type I-5 in Table 2 to have a smoothly changing sequence of triaxial solutions. For lower resolutions, the sequences appear to be less smooth especially for the plot of and for the part of the sequences closer to axisymmetric solutions. One of the reasons for this may be that, when one compares neighboring solutions of deformed sequences for mass, binding energy, or angular momentum, the change in these quantities for triaxial sequences is much smaller than that for axially symmetric sequences of about the same amount of deformation.
iv.2 Properties of the triaxial sequences
As a sample of calculated solutions, the density contours in the , , and planes and the surface plot of the model with parameters and and the largest deformation are presented in Fig.3. The solution corresponds to the last row of data shown in Table 5 in Appendix B.
In Figs. 4 and 5, and are plotted for and respectively, against the eccentricity for the constant rest mass sequences shown in the same Table 5. For uniformly rotating Newtonian polytropes, at the bifurcation point weakly depends on the difference of the EOS parameter, the polytropic index, whose value is about . For highly differential rotations, may vary largely Saijo06 (); Yoshida02 (). The definitions of and in general relativity are given in Appendix A.
Our results for the solution sequences of uniformly rotating relativistic polytropes with and suggest that the value of at the bifurcation point strongly depends on compactness . In Table 4, approximate values of quantities at the bifurcation point of each model are shown, which are evaluated by linearly extrapolating the triaxial sequence to the corresponding axisymmetric sequence. The value of at the bifurcation point becomes for the compact model , and it will certainly increase for a more compact sequence.
As seen in the plot of Fig. 5, the triaxial solution sequence for becomes shorter as increases. In fact, we were not able to find a triaxial solution sequence for ; the triaxial sequence may disappear at a certain value of between 0.14-0.2. We discuss an interesting consequence of the disappearance of triaxial sequences for high compactness in the last section.
V Discussion: proto-neutron star contraction
As a result of massive stellar core collapses, proto-neutron stars are formed and contract to more compact neutron stars within the time scale of cooling of a few tens of seconds protoNS (). Even for the small rotation rate of the collapsing stellar core, the ratio of the proto-neutron star becomes much higher than the value where the axisymmetric solution becomes secularly unstable against the viscosity driven bar mode instability LS95 (); STU04 (). Therefore, uniformly rotating triaxial solutions discussed in this paper may describe a quasi-stationary model of proto-neutron star contraction in the range of , assuming the following: (1) a certain mechanism of strong viscosity operates during the contraction, (2) the time scale rapid cooling is shorter than that of gravitational radiation reaction, (3) the effective polytropic (adiabatic) index of the EOS for the realistic neutron star matter is small (large) enough to allow uniformly rotating triaxial solutions, and (4) those triaxial solutions are dynamically stable.
Such an evolutionary track of a proto-neutron star contraction has been considered using a compressible ellipsoidal model333In their work, changes in the entropy during the evolution is modeled by the changes in the adiabatic constant of the polytropic EOS. STU04 (). Our results add two further important features to this. First, the sequences of triaxial solutions terminate at the maximally deformed models, at the mass-shedding limits, and the changes in ADM mass or total angular momentum are small along the triaxial sequences from the bifurcation points to the termination points, even for the stiffer EOS such as , as seen in Table 5; the triaxial sequences are not very long at all. Secondly, the triaxial sequences may become shorter and disappear as the compactness becomes larger for a relatively less stiff EOS such as .
The angular velocity near the braking limit is estimated as . Therefore once the secular bar mode instability sets in, conserving and , the proto-neutron star evolves towards a maximally deformed triaxial configuration as it contracts, say, from to . And then, it is likely that it always evolves along the sequence of maximally deformed configurations during the contraction as long as such triaxial solutions exist and are dynamically stable in the parameter region of the effective and . Excess angular momentum arising from contraction may be transported outward by mass ejected from the Lagrange point at the cusp of the longest semi-major axis. (Note that the time scale of the mass ejection may be that of cooling, which is much longer than the dynamical time scale.) Furthermore, if the effective satisfies , the triaxial solution may disappear when the solution reaches a certain value of the compactness and higher. Dynamically stability of such uniformly rotating solutions are not known, but it is unlikely that the dynamical instability appears within such short triaxial sequences, along which the ratio is nearly constant.
It is estimated that the amplitude of gravitational wave (GW) signals from such objects may be detectable using the ground based laser interferometric detectors, if the source is within a few tens of Mpc LS95 (); Saijo06 (). Detection of the persistent GW signals even after the proto-neutron star contraction phase suggests a large effective , while the shutdown of the signal during the contraction implies the relatively smaller . Detection of such GW signal may set another constraint on the EOS parameter of high density matter. Source modeling for constructing the wave templates may be straightforward because one can concentrate on calculating the maximally deformed configurations. Our next plan is to include more realistic nuclear EOS in the code, then to estimate the gravitational wave amplitude for those EOSs that allow the triaxial solutions.
Acknowledgements.We would like to thank John Friedman for discussions and warm encouragement. KU thanks Yoshiharu Eriguchi for discussions and for providing a reprint of PhD thesis by Tetsuo Nozawa, and Shin Yoshida for discussions. This work was supported by NSF grants No. PHY0071044, PHY0503366, NASA Grant No. NNG05GB99G, the Greek State Scholarships Foundation, and JSPS Grant-in-Aid for Scientific Research(C) 20540275.
Appendix A Formulas for mass and angular momentum
Definitions of the quantities shown in tables and figures that characterize each solution of a rotating relativistic star, and their expressions in terms of the metric potentials in the IWM formulation, are summarized in this Appendix.
The rest mass of the star is written as
where and .
The ADM mass becomes
where and , and coincides with at spatial infinity.
The Komar mass associated with a timelike Killing field is written
and, in the IWM formulation, we have
where was used. The above derivation holds if the global timelike Killing field exists. For the spacetime of a triaxially deformed rotating star, no such timelike Killing field exists. Instead, an asymptotic Komar mass can be written
In SUF04 (), we have derived sufficient conditions of the fall off of the 3-metric and extrinsic curvature and their time derivative for the relation to be satisfied. In the IWM formulation the fall off of each field is sufficiently fast to have the equality. And also in this case, the above two definitions for agree.
The total angular momentum calculated in the asymptotics is written
The relativistic analog of the kinetic energy is defined by
therefore, for uniform rotation we have . Also the relativistic analog of the gravitational potential energy is defined by
where is the proper mass defined by