Recent advances in numerical relativity have enabled the longterm dynamical evolution of rotating stars. Several interesting phenomena, such as dynamical instabilities, pulsation modes, and neutron star and black hole formation in rotating collapse have now been studied in full general relativity. The current studies are limited to relativistic polytropes, but new 3D simulations with realistic equations of state should be expected in the near future.
The goal of this article is to present a summary of theoretical and numerical methods that are used to describe the equilibrium properties of rotating relativistic stars, their oscillations and their dynamical evolution. It focuses on the most recently available preprints, in order to rapidly communicate new methods and results. At the end of some sections, the reader is directed to papers that could not be presented in detail here, or to other review articles. As new developments in the field occur, updated versions of this article will appear.
2 The Equilibrium Structure of Rotating Relativistic Stars
2.1 Assumptions
A relativistic star can have a complicated structure (such as a solid crust, magnetic field, possible superfluid interior, possible quark core, etc.). Still, its bulk properties can be computed with reasonable accuracy by making several simplifying assumptions.
The matter can be modeled to be a perfect fluid because observations of pulsar glitches have shown that the departures from a perfect fluid equilibrium (due to the presence of a solid crust) are of order
${10}^{5}$
(see [
111]
). The temperature of a cold neutron star does not affect its bulk properties and can be assumed to be 0 K, because its thermal energy (
$\ll 1MeV\sim {10}^{10}K$
) is much smaller than Fermi energies of the interior (
$>60MeV$
). One can then use a zerotemperature, barotropic equation of state (EOS) to describe the matter:
$$\begin{array}{c}\varepsilon =\varepsilon \left(P\right),\end{array}$$ 
(1)

where
$\varepsilon $
is the energy density and
$P$
is the pressure. At birth, a neutron star is expected to be rotating differentially, but as the neutron star cools, several mechanisms can act to enforce uniform rotation. Kinematical shear viscosity is acting against differential rotation on a timescale that has been estimated to be [
100,
101,
77]
$$\begin{array}{c}\tau \sim 18{\left(\frac{\rho}{{10}^{15}gc{m}^{3}}\right)}^{5/4}{\left(\frac{T}{{10}^{9}K}\right)}^{2}\left(\frac{R}{{10}^{6}cm}\right)yr,\end{array}$$ 
(2)

where
$\rho $
,
$T$
and
$R$
are the central density, temperature, and radius of the star. It has also been suggested that convective and turbulent motions may enforce uniform rotation on a timescale of the order of days [
152]
. In recent work, Shapiro [
266]
suggests that magnetic braking of differential rotation by Alfvén waves could be the most effective damping mechanism, acting on short timescales of the order of minutes.
Within roughly a year after its formation, the temperature of a neutron star becomes less than
${10}^{9}K$
and its outer core is expected to become superfluid (see [
226]
and references therein).
Rotation causes superfluid neutrons to form an array of quantized vortices, with an intervortex spacing of
$$\begin{array}{c}{d}_{n}\sim 3.4\times {10}^{3}{\Omega}_{2}^{1/2}cm,\end{array}$$ 
(3)

where
${\Omega}_{2}$
is the angular velocity of the star in
${10}^{2}{s}^{1}$
. On scales much larger than the intervortex spacing, e.g. on the order of 1 cm, the fluid motions can be averaged and the rotation can be considered to be uniform [
284]
. With such an assumption, the error in computing the metric is of order
$$\begin{array}{c}{\left(\frac{1cm}{R}\right)}^{2}\sim {10}^{12},\end{array}$$ 
(4)

assuming
$R\sim 10km$
to be a typical neutron star radius.
The above arguments show that the bulk properties of an isolated rotating relativistic star can be modeled accurately by a uniformly rotating, zerotemperature perfect fluid. Effects of differential rotation and of finite temperature need only be considered during the first year (or less) after the formation of a relativistic star.
2.2 Geometry of spacetime
In general relativity, the spacetime geometry of a rotating star in equilibrium can be described by a stationary and axisymmetric metric
${g}_{ab}$
of the form
$$\begin{array}{c}d{s}^{2}={e}^{2\nu}d{t}^{2}+{e}^{2\psi}(d\phi \omega dt{)}^{2}+{e}^{2\mu}(d{r}^{2}+{r}^{2}d{\theta}^{2}),\end{array}$$ 
(5)

where
$\nu $
,
$\psi $
,
$\omega $
and
$\mu $
are four metric functions that depend on the coordinates
$r$
and
$\theta $
only (see e.g. Bardeen and Wagoner [
25]
). Unless otherwise noted, we will assume
$c=G=1$
. In the exterior vacuum, it is possible to reduce the number of metric functions to three, but as long as one is interested in describing the whole spacetime (including the sourceregion of nonzero pressure), four different metric functions are required. It is convenient to write
${e}^{\psi}$
in the the form
$$\begin{array}{c}{e}^{\psi}=rsin\theta B{e}^{\nu},\end{array}$$ 
(6)

where
$B$
is again a function of
$r$
and
$\theta $
only [
23]
.
One arrives at the above form of the metric assuming that
i) the spacetime has a timelike Killing vector field
${t}^{a}$
and a second Killing vector field
${\phi}^{a}$
corresponding to axial symmetry, ii) the spacetime is asymptotically flat, i.e.
${t}_{a}{t}^{a}=1$
,
${\phi}_{a}{\phi}^{a}=+\infty $
and
${t}_{a}{\phi}^{a}=0$
at spatial infinity.
According to a theorem by Carter [
56]
, the two Killing vectors commute and one can choose coordinates
${x}^{0}=t$
and
${x}^{3}=\phi $
(where
${x}^{a}$
,
$a=0,...,3$
are the coordinates of the spacetime), such that
${t}^{a}$
and
${\phi}^{a}$
are coordinate vector fields. If, furthermore, the source of the gravitational field satisfies the circularity condition (absence of meridional convective currents), then another theorem [
57]
shows that the 2surfaces orthogonal to
${t}^{a}$
and
${\phi}^{a}$
can be described by the remaining two coordinates
${x}^{1}$
and
${x}^{2}$
. A common choice for
${x}^{1}$
and
${x}^{2}$
are quasiisotropic coordinates, for which
${g}_{r\theta}=0$
and
${g}_{\theta \theta}={r}^{2}{g}_{rr}$
(in spherical polar coordinates), or
${g}_{\varpi z}=0$
and
${g}_{zz}={r}^{2}{g}_{\varpi \varpi}$
(in cylindrical coordinates). In the slow rotation formalism by Hartle [
142]
, a different form of the metric is used, requiring
${g}_{\theta \theta}={g}_{\phi \phi}/{sin}^{2}\theta $
.
The three metric functions
$\nu $
,
$\psi $
and
$\omega $
can be written as invariant combinations of the two Killing vectors
${t}^{a}$
and
${\phi}^{a}$
, through the relations
$$\begin{array}{ccc}{t}_{a}{t}^{a}& =& {g}_{tt},\end{array}$$ 
(7)

$$\begin{array}{ccc}{\phi}_{a}{\phi}^{a}& =& {g}_{\phi \phi},\end{array}$$ 
(8)

$$\begin{array}{ccc}{t}_{a}{\phi}^{a}& =& {g}_{t\phi},\end{array}$$ 
(9)

while the fourth metric function
$\mu $
determines the conformal factor
${e}^{2\mu}$
that characterizes the geometry of the orthogonal 2surfaces.
There are two main effects that distinguish a rotating relativistic star from its nonrotating counterpart: The shape of the star is flattened by centrifugal forces (an effect that first appears at second order in the rotation rate), and the local inertial frames are dragged by the rotation of the source of the gravitational field. While the former effect is also present in the Newtonian limit, the latter is a purely relativistic effect. The study of the dragging of inertial frames in the spacetime of a rotating star is assisted by the introduction of the local ZeroAngularMomentumObservers (ZAMO) [
22,
23]
. These are observers whose worldlines are normal to the
$t=const.$
hypersurfaces, and they are also called Eulerian observers. Then, the metric function
$\omega $
is the angular velocity of the local ZAMO with respect to an observer at rest at infinity. Also,
${e}^{\nu}$
is the time dilation factor between the proper time of the local ZAMO and coordinate time
$t$
(proper time at infinity) along a radial coordinate line. The metric function
$\psi $
has a geometrical meaning:
${e}^{\psi}$
is the proper circumferential radius of a circle around the axis of symmetry. In the nonrotating limit, the metric ( 5 ) reduces to the metric of a nonrotating relativistic star in isotropic coordinates (see [
321]
for the definition of these coordinates).
In rapidly rotating models, an
ergosphere can appear, where
${g}_{tt}>0$
. In this region, the rotational framedragging is strong enough to prohibit counterrotating timelike or null geodesics to exist, and particles can have negative energy with respect to a stationary observer at infinity.
Radiation fields (scalar, electromagnetic, or gravitational) can become unstable in the ergosphere [
107]
, but the associated growth time is comparable to the age of the universe [
67]
.
The asymptotic behaviour of the metric functions
$\nu $
and
$\omega $
is
$$\begin{array}{c}\begin{array}{ccc}\nu & \sim & \frac{M}{r}+\frac{Q}{{r}^{3}}{P}_{2}(cos\theta ),\\ \omega & \sim & \frac{2J}{{r}^{3}},\end{array}\end{array}$$ 
(10)

where
$M$
,
$J$
and
$Q$
are the gravitational mass, angular momentum and quadrupole moment of the source of the gravitational field (see Section 2.5 for definitions). The asymptotic expansion of the dragging potential
$\omega $
shows that it decays rapidly far from the star, so that its effect will be significant mainly in the vicinity of the star.
2.3 The rotating fluid
When sources of nonisotropic stresses (such as a magnetic field or a solid state of parts of the star), viscous stresses, and heat transport are neglected in constructing an equilibrium model of a relativistic star, then the matter can be modeled as a perfect fluid, described by the stressenergy tensor
$$\begin{array}{c}{T}^{ab}=(\varepsilon +P){u}^{a}{u}^{b}+P{g}^{ab},\end{array}$$ 
(11)

where
${u}^{a}$
is the fluid's 4velocity. In terms of the two Killing vectors
${t}^{a}$
and
${\phi}^{a}$
, the 4velocity can be written as
$$\begin{array}{c}{u}^{a}=\frac{{e}^{\nu}}{\sqrt{1{v}^{2}}}({t}^{a}+\Omega {\phi}^{a}),\end{array}$$ 
(12)

where
$v$
is the 3velocity of the fluid with respect to a local ZAMO, given by
$$\begin{array}{c}v=(\Omega \omega ){e}^{\psi \nu},\end{array}$$ 
(13)

and
$\Omega \equiv {u}^{\phi}/{u}^{t}=d\phi /dt$
is the angular velocity of the fluid in the coordinate frame, which is equivalent to the angular velocity of the fluid as seen by an observer at rest at infinity. Stationary configurations can be differentially rotating, while uniform rotation (
$\Omega =$
const.) is a special case (see Section 2.5 ).
2.4 Equations of structure
Having specified an equation of state of the form
$\varepsilon =\varepsilon \left(P\right)$
, the structure of the star is determined by solving four components of Einstein's gravitational field equations
$$\begin{array}{c}{R}_{ab}=8\pi \left({T}_{ab}\frac{1}{2}{g}_{ab}T\right),\end{array}$$ 
(14)

(where
${R}_{ab}$
is the Ricci tensor and
$T={{T}_{a}}^{a}$
) and the equation of hydrostationary equilibrium.
Setting
$\zeta =\mu +\nu $
, one common choice for the gravitational field equations is [
54]
$$\begin{array}{ccc}\nabla \cdot (B\nabla \nu )& =& \frac{1}{2}{r}^{2}{sin}^{2}\theta {B}^{3}{e}^{4\nu}\nabla \omega \cdot \nabla \omega \end{array}$$ 
(15)

$$\begin{array}{ccc}& & +4\pi B{e}^{2\zeta 2\nu}\left[\frac{(\varepsilon +P)(1+{v}^{2})}{1{v}^{2}}+2P\right],\end{array}$$ 
(16)

$$\begin{array}{ccc}\nabla \cdot ({r}^{2}{sin}^{2}\theta {B}^{3}{e}^{4\nu}\nabla \omega )& =& 16\pi rsin\theta {B}^{2}{e}^{2\zeta 4\nu}\frac{(\varepsilon +P)v}{1{v}^{2}},\end{array}$$ 
(17)

$$\begin{array}{ccc}\nabla \cdot (rsin\theta \nabla B)& =& 16\pi rsin\theta B{e}^{2\zeta 2\nu}P,\end{array}$$ 
(18)

supplemented by a first order differential equation for
$\zeta $
(see [
54]
). Above,
$\nabla $
is the 3dimensional derivative operator in a flat 3space with spherical polar coordinates
$r$
,
$\theta $
,
$\phi $
.
Thus, three of the four gravitational field equations are elliptic, while the fourth equation is a first order partial differential equation, relating only metric functions. The remaining nonzero components of the gravitational field equations yield two more elliptic equations and one first order partial differential equation, which are consistent with the above set of four equations.
The equation of hydrostationary equilibrium follows from the projection of the conservation of the stressenergy tensor normal to the 4velocity
$({\delta}_{b}^{c}+{u}^{c}{u}_{b}){\nabla}_{a}{T}^{ab}=0$
, and is written as
$$\begin{array}{c}{P}_{,i}+(\varepsilon +P)\left[\nu {,}_{i}+\frac{1}{1{v}^{2}}\left(vv{,}_{i}+{v}^{2}\frac{\Omega {,}_{i}}{\Omega \omega}\right)\right]=0,\end{array}$$ 
(19)

where a comma denotes partial differentiation and
$i=1,...,3$
. When the equation of state is barotropic then the hydrostationary equilibrium equation has a first integral of motion
$$\begin{array}{c}{\int}_{0}^{P}\frac{dP}{\varepsilon +P}ln\left({u}^{a}{\nabla}_{a}t\right)+{\int}_{{\Omega}_{c}}^{\Omega}F(\Omega )d\Omega =constant.=\nu {}_{pole},\end{array}$$ 
(20)

where
$F(\Omega )={u}_{\phi}{u}^{t}$
is some specifiable function of
$\Omega $
only, and
${\Omega}_{c}$
is the angular velocity on the symmetry axis. In the Newtonian limit, the assumption of a barotropic equation of state implies that the differential rotation is necessarily constant on cylinders, and the existence of the integral of motion ( 20 ) is a direct consequence of the Poincaré–Wavre theorem (which implies that when the rotation is constant on cylinders, the effective gravity can be derived from a potential; see [
301]
).
2.5 Rotation law and equilibrium quantities
A special case of rotation law is uniform rotation (uniform angular velocity in the coordinate frame), which minimizes the total massenergy of a configuration for a given baryon number and total angular momentum [
48,
146]
. In this case, the term involving
$F(\Omega )$
in ( 20 ) vanishes.
More generally, a simple choice of a differentialrotation law is
$$\begin{array}{c}F(\Omega )={A}^{2}({\Omega}_{c}\Omega )=\frac{(\Omega \omega ){r}^{2}{sin}^{2}\theta {e}^{2(\beta \nu )}}{1(\Omega \omega {)}^{2}{r}^{2}{sin}^{2}\theta {e}^{2(\beta \nu )}},\end{array}$$ 
(21)

where
$A$
is a constant [
183,
184]
. When
$A\to \infty $
, the above rotation law reduces to the uniform rotation case. In the Newtonian limit and when
$A\to 0$
, the rotation law becomes a socalled
$j$
constant rotation law (specific angular momentum constant in space), which satisfies the Rayleigh criterion for local dynamical stability against axisymmetric disturbances (
$j$
should not decrease outwards,
$dj/d\Omega <0$
). The same criterion is also satisfied in the relativistic case [
184]
. It should be noted that differentially rotating stars may also be subject to a shear instability that tends to suppress differential rotation [
335]
.
The above rotation law is a simple choice that has proven to be computationally convenient.
More physically plausible choices must be obtained through numerical simulations of the formation of relativistic stars.
circumferential radius

$R={e}^{\psi}$

gravitational mass

$M=\int ({T}_{ab}\frac{1}{2}{g}_{ab}T){t}^{a}{\hat{n}}^{b}dV$

baryon mass

${M}_{0}=\int \rho {u}_{a}{\hat{n}}^{a}dV$

internal energy

$U=\int u{u}_{a}{\hat{n}}^{a}dV$

proper mass

${M}_{p}={M}_{0}+U$

gravitational binding energy

$W=M{M}_{p}T$

angular momentum

$J=\int {T}_{ab}{\phi}^{a}{\hat{n}}^{b}dV$

moment of inertia

$I=J/\Omega $

kinetic energy

$T=\frac{1}{2}J\Omega $



Table 1
: Equilibrium properties.
Equilibrium quantities for rotating stars, such as gravitational mass, baryon mass, or angular momentum, for example, can be obtained as integrals over the source of the gravitational field. A list of the most important equilibrium quantities that can be computed for axisymmetric models, along with the equations that define them, is displayed in Table 1 . There,
$\rho $
is the restmass density,
$u=\varepsilon \rho {c}^{2}$
is the internal energy density,
${\hat{n}}^{a}={\nabla}_{a}t/{\nabla}_{b}t{\nabla}^{b}t{}^{1/2}$
is the unit normal vector field to the
$t=const.$
spacelike hypersurfaces, and
$dV=\sqrt{{}^{3}g}{d}^{3}x$
is the proper 3volume element (with
${}^{3}g$
being the determinant of the 3metric). It should be noted that the moment of inertia cannot be computed directly as an integral quantity over the source of the gravitational field. In addition, there exists no unique generalization of the Newtonian definition of the moment of inertia in general relativity and
$I=J/\Omega $
is a common choice.
2.6 Equations of state
2.6.1 Relativistic polytropes
An analytic equation of state that is commonly used to model relativistic stars is the adiabatic, relativistic polytropic EOS of Tooper [
311]
:
$$\begin{array}{c}P=K{\rho}^{\Gamma},\end{array}$$ 
(22)

$$\begin{array}{c}\varepsilon =\rho {c}^{2}+\frac{P}{\Gamma 1},\end{array}$$ 
(23)

where
$K$
and
$\Gamma $
are the polytropic constant and polytropic exponent, respectively. Notice that the above definition is different from the form
$P=K{\varepsilon}^{\Gamma}$
(also due to Tooper [
310]
) that has also been used as a generalization of the Newtonian polytropic EOS. Instead of
$\Gamma $
, one often uses the polytropic index
$N$
, defined through
$$\begin{array}{c}\Gamma =1+\frac{1}{N}.\end{array}$$ 
(24)

For the above equation of state, the quantity
${c}^{(\Gamma 2)/(\Gamma 1)}\sqrt{{K}^{1/(\Gamma 1)}/G}$
has units of length. In gravitational units (
$c=G=1$
), one can thus use
${K}^{N/2}$
as a fundamental length scale to define dimensionless quantities. Equilibrium models are then characterized by the polytropic index
$N$
and their dimensionless central energy density. Equilibrium properties can be scaled to different dimensional values, using appropriate values for
$K$
. For
$N<1.0$
(
$N>1.0$
) one obtains stiff (soft) models, while for
$N\sim 0.5\text{\u2013}1.0$
, one obtains models with bulk properties that are comparable to those of observed neutron star radii and masses.
Notice that for the above polytropic EOS, the polytropic index
$\Gamma $
coincides with the adiabatic index of a relativistic isentropic fluid
$$\begin{array}{c}\Gamma ={\Gamma}_{ad}\equiv \frac{\varepsilon +P}{P}\frac{dP}{d\varepsilon}.\end{array}$$ 
(25)

This is not the case for the polytropic equation of state
$P=K{\varepsilon}^{\Gamma}$
, which satisfies ( 25 ) only in the Newtonian limit.
2.6.2 Hadronic equations of state
The true equation of state that describes the interior of compact stars is, still, largely unknown.
This comes as a consequence of our inability to verify experimentally the different theories that describe the strong interactions between baryons and the manybody theories of dense matter, at densities larger than about twice the nuclear density (
i.e. at densities larger than about
$5\times {10}^{14}gc{m}^{3}$
).
Many different socalled realistic EOSs have been proposed to date that all produce neutron star models that satisfy the currently available observational constraints. The two most accurate constraints are that the EOS must admit nonrotating neutron stars with gravitational mass of at least
$1.44{M}_{\odot}$
and allow rotational periods at least as small as 1.56 ms (see [
242,
186]
). Recently, the first direct determination of the gravitational redshift of spectral lines produced in the neutron star photosphere has been obtained [
73]
. This determination (in the case of the lowmass Xray binary EXO 0748676) yielded a redshift of
$z=0.35$
at the surface of the neutron star, corresponding to a mass to radius ratio of
$M/R=0.23$
(in gravitational units), which is compatible with most normal nuclear matter EOSs and incompatible with some exotic matter EOSs. The theoretically proposed EOSs are qualitatively and quantitatively very different from each other. Some are based on relativistic manybody theories while others use nonrelativistic theories with baryonbaryon interaction potentials. A classic collection of early proposed EOSs was compiled by Arnett and Bowers [
19]
, while recent EOSs are used in Salgado et al. [
260]
and in [
83]
. A review of many modern EOSs can be found in a recent article by Haensel [
137]
. Detailed descriptions and tables of several modern EOSs, especially EOSs with phase transitions, can be found in Glendenning's book [
121]
.
High density equations of state with pion condensation have been proposed by Migdal [
227]
and Sawyer and Scalapino [
263]
. The possibility of kaon condensation is discussed by Brown and Bethe [
50]
(but see also Pandharipande et al. [
240]
). Properties of rotating Skyrmion stars have been computed in [
236]
.
The realistic EOSs are supplied in the form of an energy density vs. pressure table and intermediate values are interpolated. This results in some loss of accuracy because the usual interpolation methods do not preserve thermodynamical consistency. Swesty [
300]
devised a cubic Hermite interpolation scheme that does preserve thermodynamical consistency and the scheme has been shown to indeed produce higher accuracy neutron star models in Nozawa et al. [
235]
.
Usually, the interior of compact stars is modeled as a onecomponent ideal fluid. When neutron stars cool below the superfluid transition temperature, the part of the star that becomes superfluid can be described by a twofluid model and new effects arise. Andersson and Comer [
8]
have recently used such a description in a detailed study of slowly rotating superfluid neutron stars in general relativity, while the first rapidly rotating models are presented in [
247]
.
2.6.3 Strange quark equations of state
Strange quark stars are likely to exist, if the ground state of matter at large atomic number is in the form of a quark fluid, which would then be composed of about equal numbers of up, down, and strange quarks together with electrons, which give overall charge neutrality [
37,
97]
. The strangeness per unit baryon number is
$\simeq 1$
. The first relativistic models of stars composed of quark matter were computed by Ipser, Kislinger, and Morley [
157]
and by Brecher and Caporaso [
49]
, while the first extensive study of strange quark star properties is due to Witten [
325]
.
The strange quark matter equation of state can be represented by the following linear relation between pressure and energy density:
$$\begin{array}{c}P=a(\varepsilon {\varepsilon}_{0}),\end{array}$$ 
(26)

where
${\varepsilon}_{0}$
is the energy density at the surface of a bare strange star (neglecting a possible thin crust of normal matter). The MIT bag model of strange quark matter involves three parameters, the bag constant,
$\mathcal{\mathcal{B}}={\varepsilon}_{0}/4$
, the mass of the strange quark,
${m}_{s}$
, and the QCD coupling constant,
${\alpha}_{c}$
.
The constant
$a$
in ( 26 ) is equal to
$1/3$
if one neglects the mass of the strange quark, while it takes the value of
$a=0.289$
for
${m}_{s}=250MeV$
. When measured in units of
${\mathcal{\mathcal{B}}}_{60}=\mathcal{\mathcal{B}}/(60MeVf{m}^{3})$
, the constant
$B$
is restricted to be in the range
$$\begin{array}{c}0.9821<{\mathcal{\mathcal{B}}}_{60}<1.525,\end{array}$$ 
(27)

assuming
${m}_{s}=0$
. The lower limit is set by the requirement of stability of neutrons with respect to a spontaneous fusion into strangelets, while the upper limit is determined by the energy per baryon of
${}^{56}$
Fe at zero pressure (930.4 MeV). For other values of
${m}_{s}$
the above limits are modified somewhat.
A more recent attempt to describe deconfined strange quark matter is the Dey
et al. EOS [
86]
, which has asymptotic freedom built in. It describes deconfined quarks at high densities and confinement at zero pressure. The Dey et al. EOS can be approximated by a linear relation of the same form as the MIT bag model strange star EOS ( 26 ). In such a linear approximation, typical values of the constant
$a$
are 0.45–0.46 [
127]
.
Going further A review of strange quark star properties can be found in [
320]
. Hybrid stars that have a mixedphase region of quark and hadronic matter, have also been proposed (see e.g. [
121]
).
A study of the relaxation effect in dissipative relativistic fluid theories is presented in [
199]
.
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 selfconsistent 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 [
47]
. Here, we will review the following methods: Hartle's slow rotation formalism, the Newton–Raphson linearization scheme due to Butterworth and Ipser [
54]
, a scheme using Green's functions by Komatsu et al. [
183,
184]
, a minimal surface scheme due to Neugebauer and Herold [
234]
, and spectralmethod schemes by Bonazzola et al. [
46,
45]
and by Ansorg et al. [
17]
.
Below we give a description of each method and its various implementations (codes).
2.7.1 Hartle
To order
$\mathcal{O}\left({\Omega}^{2}\right)$
, the structure of a star changes only by quadrupole terms and the equilibrium equations become a set of ordinary differential equations. Hartle's [
142,
147]
method computes rotating stars in this slow rotation approximation, and a review of slowly rotating models has been compiled by Datta [
81]
. Weber et al. [
317,
319]
also implement Hartle's formalism to explore the rotational properties of four new EOSs. Weber and Glendenning [
318]
improve on Hartle's formalism in order to obtain a more accurate estimate of the angular velocity at the massshedding limit, but their models still show large discrepancies compared to corresponding models computed without the slow rotation approximation [
260]
.
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 massshedding limit.
2.7.2 Butterworth and Ipser (BI)
The BI scheme [
54]
solves the four field equations following a Newton–Raphsonlike 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 highaccuracy 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 [
54,
53]
.
Friedman
et al. [
112,
113]
(FIP) extend the BI code to obtain a large number of rapidly rotating models based on a variety of realistic EOSs. Lattimer et al. [
194]
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 [
183,
184]
, the same set of field equations as in BI is used, but the three elliptictype 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 [
183,
184,
94]
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 semiinfinite region
$[0,\infty )$
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 [
68]
and polytropic and realistic models are computed in [
70]
and [
69]
.
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 [
294]
and in [
292]
. It is available as a public domain code, named rns, and can be downloaded from [
291]
.
2.7.4 Bonazzola et al. (BGSM)
In the BGSM scheme [
46]
, 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 [
260,
261]
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 [
131,
42]
(see Section 2.7.7 ).
While the field equations used in the BI and KEH schemes assume a perfect fluid, isotropic stressenergy tensor, the BGSM formulation makes no assumption about the isotropy of
${T}_{ab}$
. 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. [
45]
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 quasistationary models of binary neutron stars. The new method has been used in [
132]
for computing models of rapidly rotating strange stars and it has also been used in 3D computations of the onset of the viscositydriven instability to barmode formation [
128]
.
2.7.6 Ansorg et al. (AKM)
A new multidomain spectral method has been introduced in [
17,
18]
. The method can use several domains inside the star, one for each possible phase transition. Surfaceadapted coordinates are used and approximated by a twodimensional 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 massshedding 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 wellknown virial relation
$$\begin{array}{c}2T+3(\Gamma 1)U+W=0.\end{array}$$ 
(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 [
38]
. More recently, two relativistic virial identities, valid for general asymptotically flat spacetimes, have been derived by Bonazzola and Gourgoulhon [
131,
42]
. The 3dimensional virial identity (GRV3) [
131]
is the extension of the Newtonian virial identity ( 28 ) to general relativity. The 2dimensional (GRV2) [
42]
virial identity is the generalization of the identity found in [
38]
(for axisymmetric spacetimes) to general asymptotically flat spacetimes. In [
42]
, the Newtonian limit of GRV2, in axisymmetry, is also derived. Previously, such a Newtonian identity had only been known for spherical configurations [
58]
.
The two virial identities are an important tool for checking the accuracy of numerical models and have been repeatedly used by several authors [
46,
260,
261,
235,
17]
.
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 [
294]
. 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 [
294]
, 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. [
94]
, resulted from the fact that Eriguchi et al. and FIP used different versions of a tabulated EOS. Nozawa et al. [
235]
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 threepoint finitedifference 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. [
17]
. Their multidomain 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 [
17]
, which includes results obtained by the Lorene/rotstar code in [
128]
and by the SF code with higher resolution than the resolution used in [
235]
, 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 doubleprecision accuracy, while the Lorene/rotstar code has a typical relative accuracy of
$2\times {10}^{4}$
to
$7\times {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 [
235]
.
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 massshedding 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 [
41]
. A formulation for nonaxisymmetric, uniformly rotating equilibrium configurations in the second postNewtonian approximation is presented in [
21]
.

AKM

Lorene/

SF

SF

BGSM

KEH





rotstar

(
$260\times 400$
)

(
$70\times 200$
)



${\overline{p}}_{c}$

1

.

0

${r}_{p}/{r}_{e}$

0

.

7




$1\times {10}^{3}$


$\overline{\Omega}$

1

.

41170848318

$9\times {10}^{6}$

$3\times {10}^{4}$

$3\times {10}^{3}$

$1\times {10}^{2}$

$1\times {10}^{2}$

$\overline{M}$

0

.

135798178809

$2\times {10}^{4}$

$2\times {10}^{5}$

$2\times {10}^{3}$

$9\times {10}^{3}$

$2\times {10}^{2}$

${\overline{M}}_{0}$

0

.

186338658186

$2\times {10}^{4}$

$2\times {10}^{4}$

$3\times {10}^{3}$

$1\times {10}^{2}$

$2\times {10}^{3}$

${\overline{R}}_{circ}$

0

.

345476187602

$5\times {10}^{5}$

$3\times {10}^{5}$

$5\times {10}^{4}$

$3\times {10}^{3}$

$1\times {10}^{3}$

$\overline{J}$

0

.

0140585992949

$2\times {10}^{5}$

$4\times {10}^{4}$

$5\times {10}^{4}$

$2\times {10}^{2}$

$2\times {10}^{2}$

${Z}_{p}$

1

.

70735395213

$1\times {10}^{5}$

$4\times {10}^{5}$

$1\times {10}^{4}$

$2\times {10}^{2}$

$6\times {10}^{2}$

${Z}_{eq}^{f}$

$0$

.

162534082217

$2\times {10}^{4}$

$2\times {10}^{3}$

$2\times {10}^{2}$

$4\times {10}^{2}$

$2\times {10}^{2}$

${Z}_{eq}^{b}$

11

.

3539142587

$7\times {10}^{6}$

$7\times {10}^{5}$

$1\times {10}^{3}$

$8\times {10}^{2}$

$2\times {10}^{1}$

$GRV3$

$4\times {10}^{13}$

$3\times {10}^{6}$

$3\times {10}^{5}$

$1\times {10}^{3}$

$4\times {10}^{3}$

$1\times {10}^{1}$



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
$\times $
radial) grid points. See [235] for the definition of the various equilibrium quantities.
2.8 Analytic approximations to the exterior spacetime
The exterior metric of a rapidly rotating neutron star differs considerably from the Kerr metric.
The two metrics agree only to lowest order in the rotational velocity [
148]
. At higher order, the multipole moments of the gravitational field created by a rapidly rotating compact star are different from the multipole moments of the Kerr field. There have been many attempts in the past to find analytic solutions to the Einstein equations in the stationary, axisymmetric case, that could describe a rapidly rotating neutron star. An interesting solution has been found recently by Manko et al. [
219,
220]
. For nonmagnetized sources of zero net charge, the solution reduces to a 3parameter solution, involving the mass, specific angular momentum, and a parameter that depends on the quadrupole moment of the source. Although this solution depends explicitly only on the quadrupole moment, it approximates the gravitational field of a rapidly rotating star with higher nonzero multipole moments. It would be interesting to determine whether this analytic quadrupole solution approximates the exterior field of a rapidly rotating star more accurately than the quadrupole,
$\mathcal{O}\left({\Omega}^{2}\right)$
, slow rotation approximation.
The above analytic solution and an earlier one that was not represented in terms of rational functions [
218]
have been used in studies of energy release during disk accretion onto a rapidly rotating neutron star [
278,
279]
. In [
276]
, a different approximation to the exterior spacetime, in the form of a multipole expansion far from the star, has been used to derive approximate analytic expressions for the location of the innermost stable circular orbit (ISCO). Even though the analytic solutions in [
276]
converge slowly to an exact numerical solution at the surface of the star, the analytic expressions for the location and angular velocity at the ISCO are in good agreement with numerical results.
2.9 Properties of equilibrium models
2.9.1 Bulk properties of equilibrium models
Neutron star models constructed with various realistic EOSs have considerably different bulk properties, due to the large uncertainties in the equation of state at high densities. Very compressible (soft) EOSs produce models with small maximum mass, small radius, and large rotation rate. On the other hand, less compressible (stiff ) EOSs produce models with a large maximum mass, large radius, and low rotation rate.
The gravitational mass, equatorial radius, and rotational period of the maximum mass model constructed with one of the softest EOSs (EOS B) (
$1.63{M}_{\odot}$
, 9.3 km, 0.4 ms) are a factor of two smaller than the mass, radius, and period of the corresponding model constructed by one of the stiffest EOSs (EOS L) (
$3.27{M}_{\odot}$
, 18.3 km, 0.8 ms). The two models differ by a factor of 5 in central energy density and by a factor of 8 in the moment of inertia!
Not all properties of the maximum mass models between proposed EOSs differ considerably, at least not within groups of similar EOSs. For example, most realistic hadronic EOSs predict a maximum mass model with a ratio of rotational to gravitational energy
$T/W$
of
$0.11\pm 0.02$
, a dimensionless angular momentum
$cJ/G{M}^{2}$
of
$0.64\pm 0.06$
, and an eccentricity of
$0.66\pm 0.04$
[
111]
.
Hence, within the set of realistic hadronic EOSs, some properties are directly related to the stiffness of the EOS while other properties are rather insensitive to stiffness. On the other hand, if one considers strange quark EOSs, then for the maximum mass model
$T/W$
can become a factor of about two larger than for hadronic EOSs.
Figure 1
: 2D surface of equilibrium models for EOS L. The surface is bounded by the nonrotating (
$J=0$
) and massshedding (
$\Omega ={\Omega}_{K}$
) limits and formed by constant
$J$
and constant
${M}_{0}$
sequences (solid lines). The projection of these sequences in the
$J$
–
$M$
plane are shown as longdashed lines. Also shown are the axisymmetric instability sequence (shortdashed line). The projection of the 2D surface in the
$J$
–
$M$
plane shows an overlapping (see dotted lines). (Figure 7 of Stergioulas and Friedman [
294]
.)
Compared to nonrotating stars, the effect of rotation is to increase the equatorial radius of the star and also to increase the mass that can be sustained at a given central energy density. As a result, the mass of the maximum mass rotating model is roughly 15–20% higher than the mass of the maximum mass nonrotating model, for typical realistic hadronic EOSs. The corresponding increase in radius is 30–40%. The effect of rotation in increasing the mass and radius becomes more pronounced in the case of strange quark EOSs (see Section 2.9.8 ).
The deformed shape of a rapidly rotating star creates a distortion, away from spherical symmetry, in its gravitational field. Far from the star, the dominant multipole moment of the rotational distortion is measured by the quadrupolemoment tensor
${Q}_{ab}$
. For uniformly rotating, axisymmetric, and equatorially symmetric configurations, one can define a scalar quadrupole moment
$Q$
, which can be extracted from the asymptotic expansion of the metric function
$\nu $
at large
$r$
, as in Equation ( 10 ).
Laarakkers and Poisson [
187]
numerically compute the scalar quadrupole moment
$Q$
for several equations of state, using the rotating neutron star code rns [
291]
. They find that for fixed gravitational mass
$M$
, the quadrupole moment is given as a simple quadratic fit,
$$\begin{array}{c}Q=a\frac{{J}^{2}}{M{c}^{2}},\end{array}$$ 
(29)

where
$J$
is the angular momentum of the star and
$a$
is a dimensionless quantity that depends on the equation of state. The above quadratic fit reproduces
$Q$
with remarkable accuracy. The quantity
$a$
varies between
$a\sim 2$
for very soft EOSs and
$a\sim 8$
for very stiff EOSs, for
$M=1.4{M}_{\odot}$
neutron stars. This is considerably different from a Kerr black hole, for which
$a=1$
[
304]
.
For a given zerotemperature EOS, the uniformly rotating equilibrium models form a 2dimensional surface in the 3dimensional space of central energy density, gravitational mass, and angular momentum [
294]
, as shown in Figure 1 for EOS L. The surface is limited by the nonrotating models (
$J=0$
) and by the models rotating at the massshedding (Kepler) limit, i.e. at the maximum allowed angular velocity (above which the star sheds mass at the equator). Cook et al. [
68,
70,
69]
have shown that the model with maximum angular velocity does not coincide with the maximum mass model, but is generally very close to it in central density and mass. Stergioulas and Friedman [
294]
show that the maximum angular velocity and maximum baryon mass equilibrium models are also distinct. The distinction becomes significant in the case where the EOS has a large phase transition near the central density of the maximum mass model; otherwise the models of maximum mass, baryon mass, angular velocity, and angular momentum can be considered to coincide for most purposes.
Going further Although rotating relativistic stars are nearly perfectly axisymmetric, a small degree of asymmetry (e.g. frozen into the solid crust during its formation) can become a source of gravitational waves. A recent review of this can be found in [
165]
.
2.9.2 Massshedding limit and the empirical formula
Massshedding occurs when the angular velocity of the star reaches the angular velocity of a particle in a circular Keplerian orbit at the equator, i.e. when
$$\begin{array}{c}\Omega ={\Omega}_{K},\end{array}$$ 
(30)

where
$$\begin{array}{c}{\Omega}_{K}=\frac{{\omega}^{\prime}}{2{\psi}^{\prime}}+{e}^{\nu \psi}{\left[{c}^{2}\frac{{\nu}^{\prime}}{{\psi}^{\prime}}+{\left(\frac{{\omega}^{\prime}}{2{\psi}^{\prime}}{e}^{\psi \nu}\right)}^{2}\right]}^{1/2}+\omega .\end{array}$$ 
(31)

In differentially rotating stars, even a small amount of differential rotation can significantly increase the angular velocity required for massshedding. Thus, a newlyborn, hot, differentially rotating neutron star or a massive, compact object created in a binary neutron star merger could be sustained (temporarily) in equilibrium by differential rotation, even if a uniformly rotating configuration with the same rest mass does not exist.
In the Newtonian limit the maximum angular velocity of uniformly rotating polytropic stars is approximately
${\Omega}_{max}\simeq (2/3{)}^{3/2}(GM/{R}^{3}{)}^{1/2}$
(this is derived using the Roche model, see [
267]
). For relativistic stars, the empirical formula [
141,
113,
108]
$$\begin{array}{c}{\Omega}_{max}=0.67\sqrt{\frac{G{M}_{max}}{{R}_{max}^{3}}}\end{array}$$ 
(32)

gives the maximum angular velocity in terms of the mass and radius of the maximum mass nonrotating model with an accuracy of 5–7%, without actually having to construct rotating models.
A revised empirical formula, using a large set of EOSs, has been computed in [
140]
.
The empirical formula results from universal proportionality relations that exist between the mass and radius of the maximum mass rotating model and those of the maximum mass nonrotating model for the same EOS. Lasota
et al. [
192]
find that, for most EOSs, the coefficient in the empirical formula is an almost linear function of the parameter
$$\begin{array}{c}{\chi}_{s}=\frac{2G{M}_{max}}{{R}_{max}{c}^{2}}.\end{array}$$ 
(33)

The Lasota et al. empirical formula
$$\begin{array}{c}{\Omega}_{max}=\mathcal{C}\left({\chi}_{s}\right)\sqrt{\frac{G{M}_{max}}{{R}_{max}^{3}}},\end{array}$$ 
(34)

with
$\mathcal{C}\left({\chi}_{s}\right)=0.468+0.378{\chi}_{s}$
, reproduces the exact values with a relative error of only 1.5%.
Weber and Glendenning [
317,
318]
derive analytically a similar empirical formula in the slow rotation approximation. However, the formula they obtain involves the mass and radius of the maximum mass rotating configuration, which is different from what is involved in ( 32 ).
2.9.3 Upper limits on mass and rotation: Theory vs. observation
The maximum mass and minimum period of rotating relativistic stars computed with realistic hadronic EOSs from the Arnett and Bowers collection [
19]
are about
$3.3{M}_{\odot}$
(EOS L) and 0.4 ms (EOS B), while
$1.4{M}_{\odot}$
neutron stars, rotating at the Kepler limit, have rotational periods between 0.53 ms (EOS B) and 1.7 ms (EOS M) [
69]
. The maximum, accurately measured, neutron star mass is currently still
$1.44{M}_{\odot}$
(see e.g. [
314]
), but there are also indications for
$2.0{M}_{\odot}$
neutron stars [
167]
.
Core collapse simulations have yielded a bimodal mass distribution of the remnant, with peaks at about
$1.3{M}_{\odot}$
and
$1.7{M}_{\odot}$
[
309]
(the second peak depends on the assumption for the highdensity EOS – if a soft EOS is assumed, then black hole formation of this mass is implied). Compact stars of much higher mass, created in a neutron star binary merger, could be temporarily supported against collapse by strong differential rotation [
30]
.
When magnetic field effects are ignored, conservation of angular momentum can yield very rapidly rotating neutron stars at birth. Recent simulations of the rotational core collapse of evolved rotating progenitors [
150,
118]
have demonstrated that rotational core collapse can easily result in the creation of neutron stars with rotational periods of the order of 1 ms (and similar initial rotation periods have been estimated for neutron stars created in the accretioninduced collapse of a white dwarf [
211]
). The existence of a magnetic field may complicate this picture. Spruit and Phinney [
287]
have presented a model in which a strong internal magnetic field couples the angular velocity between core and surface during most evolutionary phases. The core rotation decouples from the rotation of the surface only after central carbon depletion takes place. Neutron stars born in this way would have very small initial rotation rates, even smaller than the ones that have been observed in pulsars associated with supernova remnants. In this model, an additional mechanism is required to spin up the neutron star to observed periods. On the other hand, Livio and Pringle [
212]
argue for a much weaker rotational coupling between core and surface by a magnetic field, allowing for the production of more rapidly rotating neutron stars than in [
287]
. A new investigation by Heger et al., yielding intermediate initial rotation rates, is presented in [
151]
. Clearly, more detailed computations are needed to resolve this important question.
The minimum observed pulsar period is still 1.56 ms [
186]
, which is close to the experimental sensitivity of most pulsar searches. New pulsar surveys, in principle sensitive down to a few tenths of a millisecond, have not been able to detect a submillisecond pulsar [
51,
80,
74,
93]
. This is not too surprising, as there are several explanations for the absence of submillisecond pulsars. In one model, the minimum rotational period of pulsars could be set by the occurrence of the
$r$
mode instability in accreting neutron stars in Low Mass Xray Binaries (LMXBs) [
11]
. Other models are based on the standard magnetospheric model for accretioninduced spinup [
322]
or on the idea that gravitational radiation (produced by accretioninduced quadrupole deformations of the deep crust) balances the spinup torque [
34,
312]
. It has also been suggested [
52]
that the absence of submillisecond pulsars in all surveys conducted so far could be a selection effect: Submillisecond pulsars could be found more likely only in close systems (of orbital period
${P}_{orb}\sim 1hr$
), however the current pulsar surveys are still lacking the required sensitivity to easily detect such systems.
The absence of submillisecond pulsars in wide systems is suggested to be due to the turningon of the accreting neutron stars as pulsars, in which case the pulsar wind is shown to halt further spinup.
Going further A review by J.L. Friedman concerning the upper limit on the rotation of relativistic stars can be found in [
109]
.
2.9.4 The upper limit on mass and rotation set by causality
If one is interested in obtaining upper limits on the mass and rotation rate, independently of the proposed EOSs, one has to rely on fundamental physical principles. Instead of using realistic EOSs, one constructs a set of schematic EOSs that satisfy only a minimal set of physical constraints, which represent what we know about the equation of state of matter with high confidence. One then searches among all these EOSs to obtain the one that gives the maximum mass or minimum period. The minimal set of constraints that have been used in such searches are that

1.
the high density EOS matches to the known low density EOS at some matching energy density
${\varepsilon}_{m}$
, and

2.
the matter at high densities satisfies the causality constraint (the speed of sound is less than the speed of light).
In relativistic perfect fluids, the speed of sound is the characteristic velocity of the evolution equations for the fluid, and the causality constraint translates into the requirement
$$\begin{array}{c}dP/d\varepsilon \le 1.\end{array}$$ 
(35)

(see [
119]
). It is assumed that the fluid will still behave as a perfect fluid when it is perturbed from equilibrium.
For nonrotating stars, Rhoades and Ruffini showed that the EOS that satisfies the above two constraints and yields the maximum mass consists of a high density region as stiff as possible (
i.e. at the causal limit,
$dP/d\varepsilon =1$
), that matches directly to the known low density EOS. For a chosen matching density
${\varepsilon}_{m}$
, they computed a maximum mass of
$3.2{M}_{\odot}$
. However, this is not the theoretically maximum mass of nonrotating neutron stars, as is often quoted in the literature.
Hartle and Sabbadini [
145]
point out that
${M}_{max}$
is sensitive to the matching energy density and Hartle [
143]
computes
${M}_{max}$
as a function of
${\varepsilon}_{m}$
:
$$\begin{array}{c}{M}_{max}=4.8{\left({\frac{2\times {10}^{14}gc{m}^{3}}{\varepsilon}}_{m}\right)}^{1/2}{M}_{\odot}.\end{array}$$ 
(36)

In the case of rotating stars, Friedman and Ipser [
110]
assume that the absolute maximum mass is obtained by the same EOS as in the nonrotating case and compute
${M}_{max}$
as a function of matching density, assuming the BPS EOS holds at low densities. A more recent computation [
185]
uses the FPS EOS at low densities, arriving at a similar result as in [
110]
:
$$\begin{array}{c}{M}_{max}^{rot}=6.1{\left({\frac{2\times {10}^{14}gc{m}^{3}}{\varepsilon}}_{m}\right)}^{1/2}{M}_{\odot},\end{array}$$ 
(37)

where
$2\times {10}^{14}gc{m}^{3}$
is roughly nuclear saturation density for the FPS EOS. A first estimate of the absolute minimum period of uniformly rotating, gravitationally bound stars was computed by Glendenning [
124]
by constructing nonrotating models and using the empirical formula ( 32 ) to estimate the minimum period. Koranda, Stergioulas, and Friedman [
185]
improve on Glendenning's results by constructing accurate, rapidly rotating models; they show that Glendenning's results are accurate to within the accuracy of the empirical formula.
Furthermore, they show that the EOS satisfying the minimal set of constraints and yielding the minimum period star consists of a high density region at the causal limit (CL EOS),
$P=(\varepsilon {\varepsilon}_{c})$
, (where
${\varepsilon}_{c}$
is the lowest energy density of this region), which is matched to the known low density EOS through an intermediate constant pressure region (that would correspond to a first order phase transition). Thus, the EOS yielding absolute minimum period models is as stiff as possible at the central density of the star (to sustain a large enough mass) and as soft as possible in the crust, in order to have the smallest possible radius (and rotational period).
The absolute minimum period of uniformly rotating stars is an (almost linear) function of the maximum observed mass of nonrotating neutron stars,
$$\begin{array}{c}{P}_{min}=0.28ms+0.2({M}_{max}^{nonrot}/{M}_{\odot}1.44)ms,\end{array}$$ 
(38)

and is rather insensitive to the matching density
${\varepsilon}_{m}$
(the above result was computed for a matching number density of
$0.1f{m}^{3}$
). In [
185]
, it is also shown that an absolute limit on the minimum period exists even without requiring that the EOS matches to a known low density EOS, i.e. if the CL EOS,
$P=(\varepsilon {\varepsilon}_{c})$
, terminates at a surface energy density of
${\varepsilon}_{c}$
. This is not so for the causal limit on the maximum mass. Thus, without matching to a lowdensity EOS, the causality limit on
${P}_{min}$
is lowered by only 3%, which shows that the currently known part of the nuclear EOS plays a negligible role in determining the absolute upper limit on the rotation of uniformly rotating, gravitationally bound stars.
The above results have been confirmed in [
138]
, where it is shown that the CL EOS has
${\chi}_{s}=0.7081$
, independent of
${\varepsilon}_{c}$
, and the empirical formula ( 34 ) reproduces the numerical result ( 38 ) to within 2%.
2.9.5 Supramassive stars and spinup prior to collapse
Since rotation increases the mass that a neutron star of given central density can support, there exist sequences of neutron stars with constant baryon mass that have no nonrotating member.
Such sequences are called supramassive, as opposed to normal sequences that do have a nonrotating member. A nonrotating star can become supramassive by accreting matter and spinningup to large rotation rates; in another scenario, neutron stars could be born supramassive after a core collapse.
A supramassive star evolves along a sequence of constant baryon mass, slowly losing angular momentum. Eventually, the star reaches a point where it becomes unstable to axisymmetric perturbations and collapses to a black hole.
In a neutron star binary merger, prompt collapse to a black hole can be avoided if the equation of state is sufficiently stiff and/or the equilibrium is supported by strong differential rotation. The maximum mass of differentially rotating supramassive neutron stars can be significantly larger than in the case of uniform rotation. A detailed study of this massincrease has recently appeared in [
214]
.
Cook
et al. [
68,
70,
69]
have discovered that a supramassive relativistic star approaching the axisymmetric instability will actually spin up before collapse, even though it loses angular momentum. This potentially observable effect is independent of the equation of state and it is more pronounced for rapidly rotating massive stars. Similarly, stars can spin up by loss of angular momentum near the massshedding limit, if the equation of state is extremely stiff or extremely soft.
If the equation of state features a phase transition,
e.g. to quark matter, then the spinup region is very large, and most millisecond pulsars (if supramassive) would need to be spinning up [
288]
; the absence of spinup in known millisecond pulsars indicates that either large phase transitions do not occur, or that the equation of state is sufficiently stiff so that millisecond pulsars are not supramassive.
2.9.6 Rotating magnetized neutron stars
The presence of a magnetic field has been ignored in the models of rapidly rotating relativistic stars that were considered in the previous sections. The reason is that the observed surface dipole magnetic field strength of pulsars ranges between
${10}^{8}G$
and
$2\times {10}^{13}G$
. These values of the magnetic field strength imply a magnetic field energy density that is too small, compared to the energy density of the fluid, to significantly affect the structure of a neutron star. However, one cannot exclude the existence of neutron stars with higher magnetic field strengths or the possibility that neutron stars are born with much stronger magnetic fields, which then decay to the observed values (of course, there are also many arguments against magnetic field decay in neutron stars [
242]
).
In addition, even though moderate magnetic field strengths do not alter the bulk properties of neutron stars, they may have an effect on the damping or growth rate of various perturbations of an equilibrium star, affecting its stability. For these reasons, a fully relativistic description of magnetized neutron stars is desirable and, in fact, Bocquet
et al. [
36]
achieved the first numerical computation of such configurations. Following we give a brief summary of their work.
A magnetized relativistic star in equilibrium can be described by the coupled Einstein–Maxwell field equations for stationary, axisymmetric rotating objects with internal electric currents. The stressenergy tensor includes the electromagnetic energy density and is nonisotropic (in contrast to the isotropic perfect fluid stressenergy tensor). The equilibrium of the matter is given not only by the balance between the gravitational force, centrifugal force, and the pressure gradient; the Lorentz force due to the electric currents also enters the balance. For simplicity, Bocquet
et al.
consider only poloidal magnetic fields that preserve the circularity of the spacetime. Also, they only consider stationary configurations, which excludes magnetic dipole moments nonaligned with the rotation axis, since in that case the star emits electromagnetic and gravitational waves. The assumption of stationarity implies that the fluid is necessarily rigidly rotating (if the matter has infinite conductivity) [
46]
. Under these assumptions, the electromagnetic field tensor
${F}^{ab}$
is derived from a potential fourvector
${A}_{a}$
with only two nonvanishing components,
${A}_{t}$
and
${A}_{\phi}$
, which are solutions of a scalar Poisson and a vector Poisson equation respectively. Thus, the two equations describing the electromagnetic field are of similar type as the four field equations that describe the gravitational field.
For magnetic field strengths larger than about
${10}^{14}G$
, one observes significant effects, such as a flattening of the equilibrium configuration. There exists a maximum value of the magnetic field strength of the order of
${10}^{18}G$
, for which the magnetic field pressure at the center of the star equals the fluid pressure. Above this value no stationary configuration can exist.
A strong magnetic field allows a maximum mass configuration with larger
${M}_{max}$
than for the same EOS with no magnetic field and this is analogous to the increase of
${M}_{max}$
induced by rotation.
For nonrotating stars, the increase in
${M}_{max}$
due to a strong magnetic field is 13–29%, depending on the EOS. Correspondingly, the maximum allowed angular velocity, for a given EOS, also increases in the presence of a strong magnetic field.
Another application of general relativistic E/M theory in neutron stars is the study of the evolution of the magnetic field during pulsar spindown. A detailed analysis of the evolution equations of the E/M field in a slowly rotating magnetized neutron star has revealed that effects due to the spacetime curvature and due to the rotational framedragging are present in the induction equations, when one assumes finite electrical conductivity (see [
251]
and references therein). Numerical solutions of the evolution equations of the E/M have shown, however, that for realistic values of the electrical conductivity, the above relativistic effects are small, even in the case of rapid rotation [
336]
.
Going further An
$\mathcal{O}(\Omega )$
slow rotation approach for the construction of rotating magnetized relativistic stars is presented in [
136]
.
2.9.7 Rapidly rotating protoneutron stars
Following the gravitational collapse of a massive stellar core, a protoneutron star (PNS) is born.
Initially it has a large radius of about 100 km and a temperature of 50–100 MeV. The PNS may be born with a large rotational kinetic energy and initially it will be differentially rotating. Due to the violent nature of the gravitational collapse, the PNS pulsates heavily, emitting significant amounts of gravitational radiation. After a few hundred pulsational periods, bulk viscosity will damp the pulsations significantly. Rapid cooling due to deleptonization transforms the PNS, shortly after its formation, into a hot neutron star of
$T\sim 10MeV$
. In addition, viscosity or other mechanisms (see Section 2.1 ) enforce uniform rotation and the neutron star becomes quasistationary. Since the details of the PNS evolution determine the properties of the resulting cold NSs, protoneutron stars need to be modeled realistically in order to understand the structure of cold neutron stars.
Figure 2
: Isoenergy density lines of a differentially rotating protoneutron star at the massshedding limit, of rest mass
${M}_{0}=1.5{M}_{\odot}$
. (Figure 5a of Goussard, Haensel and Zdunik [
134]
; used with permission.)
Hashimoto et al. [
149]
and Goussard et al. [
133]
construct fully relativistic models of rapidly rotating, hot protoneutron stars. The authors use finitetemperature EOSs [
238,
195]
to model the interior of PNSs. Important (but largely unknown) parameters that determine the local state of matter are the lepton fraction
${Y}_{l}$
and the temperature profile. Hashimoto et al. consider only the limiting case of zero lepton fraction,
${Y}_{l}=0$
, and classical isothermality, while Goussard et al.
consider several nonzero values for
${Y}_{l}$
and two different limiting temperature profiles – a constant entropy profile and a relativistic isothermal profile. In both [
149]
and [
238]
, differential rotation is neglected to a first approximation.
The construction of numerical models with the above assumptions shows that, due to the high temperature and the presence of trapped neutrinos, PNSs have a significantly larger radius than cold NSs. These two effects give the PNS an extended envelope which, however, contains only roughly 0.1% of the total mass of the star. This outer layer cools more rapidly than the interior and becomes transparent to neutrinos, while the core of the star remains hot and neutrino opaque for a longer time. The two regions are separated by the “neutrino sphere”.
Compared to the
$T=0$
case, an isothermal EOS with temperature of 25 MeV has a maximum mass model of only slightly larger mass. In contrast, an isentropic EOS with a nonzero trapped lepton number features a maximum mass model that has a considerably lower mass than the corresponding model in the
$T=0$
case and, therefore, a stable PNS transforms to a stable neutron star. If, however, one considers the hypothetical case of a large amplitude phase transition that softens the cold EOS (such as a kaon condensate), then
${M}_{max}$
of cold neutron stars is lower than the
${M}_{max}$
of PNSs, and a stable PNS with maximum mass will collapse to a black hole after the initial cooling period. This scenario of delayed collapse of nascent neutron stars has been proposed by Brown and Bethe [
50]
and investigated by Baumgarte et al. [
31]
.
An analysis of radial stability of PNSs [
126]
shows that, for hot PNSs, the maximum angular velocity model almost coincides with the maximum mass model, as is also the case for cold EOSs. Because of their increased radius, PNSs have a different massshedding limit than cold NSs. For an isothermal profile, the massshedding limit proves to be sensitive to the exact location of the neutrino sphere. For the EOSs considered in [
149]
and [
133]
, PNSs have a maximum angular velocity that is considerably less than the maximum angular velocity allowed by the cold EOSs. Stars that have nonrotating counterparts (i.e. that belong to a normal sequence) contract and speed up while they cool down. The final star with maximum rotation is thus closer to the massshedding limit of cold stars than was the hot PNS with maximum rotation. Surprisingly, stars belonging to a supramassive sequence exhibit the opposite behavior. If one assumes that a PNS evolves without losing angular momentum or accreting mass, then a cold neutron star produced by the cooling of a hot PNS has a smaller angular velocity than its progenitor. This purely relativistic effect was pointed out in [
149]
and confirmed in [
133]
.
It should be noted here that a small amount of differential rotation significantly affects the massshedding limit, allowing more massive stars to exist than uniform rotation allows. Taking differential rotation into account, Goussard
et al. [
134]
suggest that protoneutron stars created in a gravitational collapse cannot spin faster than 1.7 ms. A similar result has been obtained by Strobel et al. [
297]
. The structure of a differentially rotating protoneutron star at the massshedding limit is shown in Figure 2 . The outer layers of the star form an extended disklike structure.
The above stringent limits on the initial period of neutron stars are obtained assuming that the PNS evolves in a quasistationary manner along a sequence of equilibrium models. It is not clear whether these limits will remain valid if one studies the early evolution of PNS without the above assumption. It is conceivable that the thin hot envelope surrounding the PNS does not affect the dynamics of the bulk of the star. If the bulk of the star rotates faster than the (stationary) massshedding limit of a PNS model, then the hot envelope will simply be shed away from the star in the equatorial region (if it cannot remain bounded to the star even when differentially rotating).
Such a fully dynamical study is needed to obtain an accurate upper limit on the rotation of neutron stars.
Going further The thermal history and evolutionary tracks of rotating PNSs (in the second order slow rotation approximation) have been studied recently in [
299]
.
2.9.8 Rotating strange quark stars
Most rotational properties of strange quark stars differ considerably from the properties of rotating stars constructed with hadronic EOSs. First models of rapidly rotating strange quark stars were computed by Friedman ([
106]
, quoted in [
123,
122]
) and by Lattimer et al. [
194]
. Colpi and Miller [
65]
use the
$\mathcal{O}\left({\Omega}^{2}\right)$
approximation and find that the spin of strange stars (newlyborn, or spunup by accretion) may be limited by the CFS instability to the
$l=m=2$
$f$
mode, since rapidly rotating strange stars tend to have
$T/W>0.14$
. Rapidly rotating strange stars at the massshedding limit have been computed first by Gourgoulhon et al. [
132]
, and the structure of a representative model is displayed in Figure 3 .
Figure 3
: Meridional plane cross section of a rapidly rotating strange star at the massshedding limit, obtained with a multidomain spectral code. The various lines are isocontours of the logenthalpy
$H$
, as defined in [
132]
. Solid lines indicate a positive value of
$H$
and dashed lines a negative value (vacuum). The thick solid line denotes the stellar surface. The thick dotdashed line denotes the boundary between the two computational domains. (Figure 4 of Gourgoulhon, Haensel, Livine, Paluch, Bonazzola, and Marck [
132]
; used with permission.)
Nonrotating strange stars obey scaling relations with the constant
$\mathcal{\mathcal{B}}$
in the MIT bag model of the strange quark matter EOS (Section 2.6.3 ); Gourgoulhon et al. [
132]
also obtain scaling relations for the model with maximum rotation rate. The maximum angular velocity scales as
$$\begin{array}{c}{\Omega}_{max}=9.92\times {10}^{3}\sqrt{{\mathcal{\mathcal{B}}}_{60}}{s}^{1},\end{array}$$ 
(39)

while the allowed range of
$\mathcal{\mathcal{B}}$
implies an allowed range of
$0.513ms<{P}_{min}<0.640ms$
. The empirical formula ( 32 ) also holds for rotating strange stars with an accuracy of better than 2%.
A derivation of the empirical formula in the case of strange stars, starting from first principles, has been presented by Cheng and Harko [
61]
, who found that some properties of rapidly rotating strange stars can be reproduced by approximating the exterior spacetime by the Kerr metric.
Since both the maximum mass nonrotating and maximum mass rotating models obey similar scalings with
$B$
, the ratios
$$\begin{array}{c}\frac{{M}_{max}^{rot}}{{M}_{max}^{stat}}=1.44,\frac{{R}_{max}^{rot}}{{R}_{max}^{stat}}=1.54,\end{array}$$ 
(40)

are independent of
$\mathcal{\mathcal{B}}$
(where
${R}_{max}$
is the radius of the maximum mass model). The maximum mass increases by 44% and the radius of the maximum mass model by 54%, while the corresponding increase for hadronic stars is, at best,
$\sim 20\%$
and
$\sim 40\%$
, correspondingly. The rotational properties of strange star models that are based on the Dey et al. EOS [
86]
are similar to those of the MIT bag model EOS [
37,
325,
97]
, but some quantitative differences exist [
127]
.
Accreting strange stars in LMXBs will follow different evolutionary paths than do accreting hadronic stars in a mass vs. central energy density diagram [
342]
. When (and if ) strange stars reach the massshedding limit, the ISCO still exists [
296]
(while it disappears for most hadronic EOSs). Stergioulas, Kluźniak, and Bulik [
296]
show that the radius and location of the ISCO for the sequence of massshedding models also scales as
${\mathcal{\mathcal{B}}}^{1/2}$
, while the angular velocity of particles in circular orbit at the ISCO scales as
${\mathcal{\mathcal{B}}}^{1/2}$
. Additional scalings with the constant
$a$
in the strange quark EOS (that were proposed in [
194]
) are found to hold within an accuracy of better than
$\sim 9\%$
for the massshedding sequence
$$\begin{array}{c}M\propto {a}^{1/2},R\propto {a}^{1/4},\Omega \propto {a}^{1/8}.\end{array}$$ 
(41)

In addition, it is found that models at massshedding can have
$T/W$
as large as
$0.28$
for
$M=1.34{M}_{\odot}$
.
As strange quark stars are very compact, the angular velocity at the ISCO can become very large. If the 1066 Hz upper QPO frequency in 4U 182030 (see [
167]
and references therein) is the frequency at the ISCO, then it rules out most models of slowly rotating strange stars in LMXBs. However, in [
296]
it is shown that rapidly rotating bare strange stars are still compatible with this observation, as they can have ISCO frequencies
$<1kHz$
even for
$1.4{M}_{\odot}$
models. On the other hand, if strange stars have a thin solid crust, the ISCO frequency at the massshedding limit increases by about 10% (compared to a bare strange star of the same mass), and the above observational requirement is only satisfied for slowly rotating models near the maximum nonrotating mass, assuming some specific values of the parameters in the strange star EOS [
340,
341]
.
Moderately rotating strange stars, with spin frequencies around 300 Hz can also be accommodated for some values of the coupling constant
${\alpha}_{c}$
[
338]
(see also [
130]
for a detailed study of the ISCO frequency for rotating strange stars). The 1066 Hz requirement for the ISCO frequency depends, of course, on the adopted model of kHz QPOs in LMXBs, and other models exist (see next section).
If strange stars can have a solid crust, then the density at the bottom of the crust is the neutron drip density
${\varepsilon}_{nd}\simeq 4.1\times {10}^{11}gc{m}^{3}$
, as neutrons are absorbed by strange quark matter. A strong electric field separates the nuclei of the crust from the quark plasma. In general, the mass of the crust that a strange star can support is very small, of the order of
${10}^{5}{M}_{\odot}$
. Rapid rotation increases by a few times the mass of the crust and the thickness at the equator becomes much larger than the thickness at the poles [
341]
. Zdunik, Haensel, and Gourgoulhon [
341]
also find that the mass
${M}_{crust}$
and thickness
${t}_{crust}$
of the crust can be expanded in powers of the spin frequency
${\nu}_{3}=\nu /({10}^{3}Hz)$
as
$$\begin{array}{ccc}{M}_{crust}& =& {M}_{crust,0}(1+0.24{\nu}_{3}^{2}+0.16{\nu}_{3}^{8}),\end{array}$$ 
(42)

$$\begin{array}{ccc}{t}_{crust}& =& {t}_{crust,0}(1+0.4{\nu}_{3}^{2}+0.3{\nu}_{3}^{6}),\end{array}$$ 
(43)

where a subscript “0” denotes nonrotating values. For
$\nu \le 500Hz$
, the above expansion agrees well with the quadratic expansion derived previously by Glendenning and Weber [
125]
. In a spinning down magnetized strange quark star with crust, parts of the crust will gradually dissolve into strange quark matter, in a strongly exothermic process. In [
341]
, it is estimated that the heating due to deconfinement may exceed the neutrino luminosity from the core of a strange star older than
$\sim 1000yr$
and may therefore influence the cooling of this compact object (see also [
334]
).
2.10 Rotating relativistic stars in LMXBs
2.10.1 Particle orbits and kHz quasiperiodic oscillations
In the last few years, Xray observations of accreting sources in LMXBs have revealed a rich phenomenology that is waiting to be interpreted correctly and could lead to significant advances in our understanding of compact objects (see [
191,
313,
248]
). The most important feature of these sources is the observation of (in most cases) twin kHz quasiperiodic oscillations (QPOs). The high frequency of these variabilities and their quasiperiodic nature are evidence that they are produced in highvelocity flows near the surface of the compact star. To date, there exist a large number of different theoretical models that attempt to explain the origin of these oscillations. No consensus has been reached, yet, but once a credible explanation is found, it will lead to important constraints on the properties of the compact object that is the source of the gravitational field in which the kHz oscillations take place. The compact stars in LMXBs are spun up by accretion, so that many of them may be rotating rapidly; therefore, the correct inclusion of rotational effects in the theoretical models for kHz QPOs is important. Under simplifying assumptions for the angular momentum and mass evolution during accretion, one can use accurate rapidly rotating relativistic models to follow the possible evolutionary tracks of compact stars in LMXBs [
71,
342]
.
In most theoretical models, one or both kHz QPO frequencies are associated with the orbital motion of inhomogeneities or blobs in a thin accretion disk. In the actual calculations, the frequencies are computed in the approximation of an orbiting test particle, neglecting pressure terms. For most equations of state, stars that are massive enough possess an ISCO, and the orbital frequency at the ISCO has been proposed to be one of the two observed frequencies. To first order in the rotation rate, the orbital frequency at the prograde ISCO is given by (see Kluźniak, Michelson, and Wagoner [
170]
)
$$\begin{array}{c}{f}_{ISCO}\simeq 2210(1+0.75j)\left(\frac{1{M}_{\odot}}{M}\right)Hz,\end{array}$$ 
(44)

where
$j=J/{M}^{2}$
. At larger rotation rates, higher order contributions of
$j$
as well as contributions from the quadrupole moment
$Q$
become important and an approximate expression has been derived by Shibata and Sasaki [
276]
, which, when written as above and truncated to the lowest order contribution of
$Q$
and to
$\mathcal{O}\left({j}^{2}\right)$
, becomes
$$\begin{array}{c}{f}_{ISCO}\simeq 2210(1+0.75j+0.78{j}^{2}0.23{Q}_{2})\left(\frac{1{M}_{\odot}}{M}\right)Hz,\end{array}$$ 
(45)

where
${Q}_{2}=Q/{M}^{3}$
.
Notice that, while rotation increases the orbital frequency at the ISCO, the quadrupole moment has the opposite effect, which can become important for rapidly rotating models. Numerical evaluations of
${f}_{ISCO}$
for rapidly rotating stars have been used in [
228]
to arrive at constraints on the properties of the accreting compact object.
In other models, orbits of particles that are eccentric and slightly tilted with respect to the equatorial plane are involved. For eccentric orbits, the periastron advances with a frequency
${\nu}_{pa}$
that is the difference between the Keplerian frequency of azimuthal motion
${\nu}_{K}$
and the radial epicyclic frequency
${\nu}_{r}$
. On the other hand, particles in slightly tilted orbits fail to return to the initial displacement
$\psi $
from the equatorial plane, after a full revolution around the star. This introduces a nodal precession frequency
${\nu}_{pa}$
, which is the difference between
${\nu}_{K}$
and the frequency of the motion out of the orbital plane (meridional frequency)
${\nu}_{\psi}$
. Explicit expressions for the above frequencies, in the gravitational field of a rapidly rotating neutron star, have been derived recently by Marković [
221]
, while in [
222]
highly eccentric orbits are considered. Morsink and Stella [
229]
compute the nodal precession frequency for a wide range of neutron star masses and equations of state and (in a postNewtonian analysis) separate the precession caused by the Lense–Thirring (framedragging) effect from the precession caused by the quadrupole moment of the star. The nodal and periastron precession of inclined orbits have also been studied using an approximate analytic solution for the exterior gravitational field of rapidly rotating stars [
277]
. These precession frequencies are relativistic effects and have been used in several models to explain the kHz QPO frequencies [
290,
249,
1,
168,
4]
.
It is worth mentioning that it has recently been found that an ISCO also exists in Newtonian gravity, for models of rapidly rotating lowmass strange stars. The instability in the circular orbits is produced by the large oblateness of the star [
169,
339,
4]
.
2.10.2 Angular momentum conservation during burst oscillations
Some sources in LMXBs show signatures of type I Xray bursts, which are thermonuclear flashes on the surface of the compact star [
197]
. Such bursts show nearlycoherent oscillations in the range 270–620 Hz (see [
313,
298]
for recent reviews). One interpretation of the burst oscillations is that they are the result of rotational modulation of surface asymmetries during the burst. In such a case, the oscillation frequency should be nearly equal to the spin frequency of the star. This model currently has difficulties in explaining some observed properties, such as the oscillations seen in the tail of the burst, the frequency increase during the burst, and the need for two antipodal hot spots in some sources that ignite at the same time. Alternative models also exist (see e.g. [
248]
).
In the spinfrequency interpretation, the increase in the oscillation frequency by a few Hz during the burst is explained as follows: The burning shell is supposed to first decouple from the neutron star and then gradually settle down onto the surface. By angular momentum conservation, the shell spins up, giving rise to the observed frequency increase. Cumming
et al. [
75]
compute the expected spinup in full general relativity and taking into account rapid rotation. Assuming that the angular momentum per unit mass is conserved, the change in angular velocity with radius is given by
$$\begin{array}{c}\frac{dln\Omega}{dlnr}=2\left[\left(1\frac{{v}^{2}}{2}\frac{R}{2}\frac{\partial \nu}{\partial r}\right)\left(1\frac{\omega}{\Omega}\right)\frac{R}{2\Omega}\frac{\partial \omega}{\partial r}\right],\end{array}$$ 
(46)

where
$R$
is the equatorial radius of the star and all quantities are evaluated at the equator.
The slow rotation limit of the above result was derived previously by Abramowicz
et al. [
2]
. The fractional change in angular velocity during spinup can then be estimated as
$$\begin{array}{c}\frac{\Delta \Omega}{\Omega}=\frac{dln\Omega}{dlnr}\left(\frac{\Delta r}{R}\right),\end{array}$$ 
(47)

where
$\Delta r$
is the coordinate expansion of the burning shell, a quantity that depends on the shell's composition. Cumming et al. find that the spin down expected if the atmosphere rotates rigidly is a factor of two to three times smaller than observed values. More detailed modeling is needed to fully explain the origin and properties of burst oscillations.
Going further A very interesting topic is the modeling of the expected Xray spectrum of an accretion disk in the gravitational field of a rapidly rotating neutron star as it could lead to observational constraints on the source of the gravitational field. See e.g. [
302,
278,
279,
33,
32]
, where work initiated by Kluzniak and Wilson [
171]
in the slow rotation limit is extended to rapidly rotating relativistic stars.
3 Oscillations and Stability
The study of oscillations of relativistic stars is motivated by the prospect of detecting such oscillations in electromagnetic or gravitational wave signals. In the same way that helioseismology is providing us with information about the interior of the Sun, the observational identification of oscillation frequencies of relativistic stars could constrain the highdensity equation of state [
12]
. The oscillations could be excited after a core collapse or during the final stages of a neutron star binary merger. Rapidly rotating relativistic stars can become unstable to the emission of gravitational waves.
When the oscillations of an equilibrium star are of small magnitude compared to its radius, it will suffice to approximate them as linear perturbations. Such perturbations can be described in two equivalent ways. In the Lagrangian approach, one studies the changes in a given fluid element as it oscillates about its equilibrium position. In the Eulerian approach, one studies the change in fluid variables at a fixed point in space. Both approaches have their strengths and weaknesses.
In the Newtonian limit, the Lagrangian approach has been used to develop variational principles [