Rotating neutron stars with exotic cores: masses, radii, stability

A set of theoretical mass-radius relations for rigidly rotating neutron stars with exotic cores, obtained in various theories of dense matter, is reviewed. Two basic observational constraints are used: the largest measured rotation frequency (716Hz) and the maximum measured mass ( 2M⊙ . The present status of measuring the radii of neutron stars is described. The theory of rigidly rotating stars in general relativity is reviewed and limitations of the slow rotation approximation are pointed out. Mass-radius relations for rotating neutron stars with hyperon and quark cores are illustrated using several models. Problems related to the non-uniqueness of the crust-core matching are mentioned. Limits on rigid rotation resulting from the mass-shedding instability and the instability with respect to the axisymmetric perturbations are summarized. The problem of instabilities and of the back-bending phenomenon are discussed in detail. Metastability and instability of a neutron star core in the case of a first-order phase transition, both between pure phases, and into a mixed-phase state, are reviewed. The case of two disjoint families (branches) of rotating neutron stars is discussed and generic features of neutron-star families and of core-quakes triggered by the instabilities are considered.


Introduction
The determination of radii R obs i of neutron stars (NS) of known masses M obs i , (i = 1, 2, . . .) would allow us to unveil the equation of state (EOS) of neutron-star cores of density significantly higher than normal nuclear density ρ 0 = 2.7 × 10 14 g cm −3 (corresponding to baryon number density n 0 = 0.16 fm −3 ). To be useful, however, the uncertainties in the values M obs i , R obs i should be sufficiently small (at the level of a few percent) and the maximum M obs i should be close to an (unknown) maximum allowable mass of NS. For the time being, a set of unevenly spaced M obs i was determined 1 [1], with a maximum value M obs max = 2.01 ± 0.04M [2] being a very strong constraint on the EOS. The precise measurement of R is still a challenge for observers (see sect. 2).
For ρ < ρ 0 constituents of matter are well established: nucleons and electrons, with a small admixture of muons at the upper subnuclear density segment where the Fermi energy of electrons exceeds the muon rest energy. However, we expect that the central density of neutron stars with Contribution to the Topical Issue on "Exotic matter in neutron stars" edited by David Blaschke, Jürgen Schaffner-Bielich, Hans-Josef Schulze. a e-mail: fortin@camk.edu.pl 1 see http://stellarcollapse.org/nsmasses. M > 1M is larger than 1.5ρ 0 -2ρ 0 . For M = 2M the star's central density may be as high as ∼ 7ρ 0 . For ρ 2ρ 0 even the actual hadronic constituents of the NS core are uncertain: are they just nucleons (zero strangeness), or more generally baryons, i.e. nucleons and hyperons (nonzero strangeness)? Maybe the density realized there is sufficient for a phase transition to quark matter? Finally, maybe in addition to baryons, real kaons or pions forming a boson condensate are present there? The uncertainty grows with increasing ρ, which is expected to be as high as 8ρ 0 -10ρ 0 at the center of the most massive NS. This uncertainty results from the lack of precise knowledge of strong interactions and the approximations (often uncontrollable) of the many-body theories of superdense hadronic matter. The uncertainty in the structure and composition of super-dense matter implies an even larger uncertainty in the EOS of NS cores. All neutron stars rotate and there are many millisecond pulsars (MSP) with rotation frequency f > 500 Hz (10 accreting X-ray pulsars and 14 radio/gamma-ray pulsars), and a radio MSP with 716 Hz is observed (sect. 2). In the present paper we restrict ourselves to rigid rotation and axisymmetric approximation for rotating NS models, nevertheless for the completeness sake we briefly recall the general picture. The approximation of rigid, axisymmetric rotation holds extremely well already a few minutes after the NS is born in the supernova core collapse.
Prior to that one expects a differentially rotating, hot and lepton-rich proto-NS with high entropy, which cannot be properly described by a cold, catalyzed matter EOS. During this period dynamo mechanism and convection may operate, increasing the interior magnetic field and leading to magneto-hydrodynamical instabilities [3]. Shortly after the birth, rigid rotation sets in due to the presence of viscosity. However, for sufficiently high rotation rates, parametrized by the kinetic energy T to potential energy W ratio, β = T /|W |, a dynamical triaxial bar-mode instability may arise in rigidly rotating stars; relativistic calculations indicate critical β 0.24. Substantial differential rotation facilitates the onset of these dynamical instabilities -they may occur at a low β ≈ 0.01. Another, secular bar-mode instability, driven by the dissipation due to viscosity or the emission of gravitational waves sets in at β 0.14 (for a review see [4]). We also expect that accreting NSs may be prone to the Rossby-type instabilities (r-modes), driven by the Coriolis force [5].
The dependence of the radius of a NS with an exotic core on its mass and the imprint of rotation on the massequatorial radius relation and the stability of rotating NS configurations is the main topic of the present paper. In our review we try to present some generic features of rotating NS with exotic (E) cores, with E being: hyperon matter, quark matter, or a baryon phase with a boson (pion or kaon) condensate. NS with E-core will be compared with standard nucleon NS models, hoping that the differences between theoretical models, confronted with observations, will help to unveil the true EOS of NS.
In order to study axisymmetric hydrostatic equilibria of rotating NS in general relativity (GR), approximate, as well as exact numerical methods of solving Einstein's equation were developed; we briefly discuss them in sect. 3. We discuss the precision which should be reached in the 2D calculations to study the mass shedding limit, the spin evolution and stability with respect to axisymmetric perturbations of rotating NS models.
We consider a relativistic star in the perfect fluid approximation. In general, the metric of space-time around a rotating NS is essentially different from that around a static star (f = 0). For f = 0 the metric depends only on the NS mass M and it is the Schwarzschild metric. For f > 0 the metric depends explicitly on the matter and pressure distribution inside a rotating star, which makes it dependent on the EOS; it also depends on f through the effect of the dragging of the inertial frames. In particular, this refers to the innermost stable orbit around an accreting MSP (sect. 3). Some properties of the circular orbits around NS and their relation to the Keplerian limit are reviewed in sect. 6.
The EOS of dense matter with a transition from a normal (N) phase to an exotic (E) one has some particularities. Generic features of such an EOS are reviewed in sect. 4. We start with the simplest case of N-E transition in full thermodynamic equilibrium, reviewed in sect. 4.1. Then we consider a more complicated case including the possibility of a metastable state and nucleation of the Ephase in the N-one in sect. 4.2.
The theoretical M -R relation (here R is the equatorial radius) depends on the rotation frequency. The region in the M -R plane allowed for rotating configurations is affected by the presence of an exotic core in massive NS. The maximum allowable NS mass is a functional of the EOS. Its value for non-rotating stars, M stat max , has to satisfy M stat max [EOS] > 2M , the largest observed mass. Rotation increases M max only by a few percent even at 716 Hz, but for the minimum mass M f min (which for f = 0 is M stat min ≈ 0.1M , see [6] and references therein), the effect of rotation is dramatic and depends indirectly also on the phase transition to an exotic high density phase (sect. 5). This property of M f min can be used to derive the EOS dependent lower bound on the mass of the fastest 716 Hz pulsar (sect. 6.3). We present in detail examples of families of rotating NS models with nucleon, hyperon, and quark cores, showing the differences between these families.
Limits on the frequency of rotation of NS are reviewed in sect. 6. There is an upper bound for the frequency of rotation for each given baryon (rest) mass of a NS M b , corresponding to a non-rotating (static) mass M s . It results from the mass-shedding instability at the equator and is called Keplerian frequency f K (M s ). This bound is quite sensitive to the EOS, because it depends on the radius of non-rotating configuration (sect. 6.1). There is also a theoretical maximum frequency for all stably rotating configurations of NS, f max = 1500 Hz-2000 Hz which depends on the EOS (sect. 6.2). The highest measured frequency of a pulsar which as for today is 716 Hz, results in an EOS dependent constraint on the mass of this fastest pulsar (sect. 6.3).
In sect. 7 we review various aspects of the spin evolution, dynamics, and stability of rotating NS with exotic cores. The softening of the EOS associated with a transition to an exotic phase can lead to a phenomenon of back bending (spin-up induced by an angular momentum loss), reviewed in sect. 7.1. In particular, we point out the possibility of the existence of unstable segments of configuration sequences, splitting the stable back bending fragment of the spin evolution track into two separate branches.
The instability induced by the softening of the EOS due to a first order phase transition into an exotic phase is discussed in sect. 7.2. A sufficiently strong softening of the EOS can lead to splitting a single one-parameter family of hydrostatic configurations of NS into two separate (disjoint) families. This feature is valid not only for static NS, but also for rotating NS models with constant f . The static criterion for the split into two branches is valid also for rigidly rotating configurations (sect. 7.2). We then review, in sect. 8, the possibility of a minicollapse of a NS due to the nucleation of the E-phase at its center during the NS evolution, and the role of metastability of the Nphase core in this process. The astrophysical signatures of a mini-collapse are described in sect. 9. In sect. 10 we review the effect of the crust formation scenario on the M -R relation.
Our conclusions are summarized in the final sect. 11.

Observational constraints from spin frequency and radius measurements
The recent discovery of two 2M pulsars [7,2] provides an important constraint on the poorly known equation of state at supra-nuclear density. In this section we summarize the current status of measurements of radius and spin frequency of NS.

Radius
The radius of a NS can in principle be extracted from the analysis of X-ray spectra emitted by the NS atmosphere (see [8] for a review). However even in the case of a nonrotating NS, due to the space-time curvature, only the apparent radius, is constrained by the modelling. It actually depends on both the radius and the mass. Measurements are complicated since they depend on the distance to the NS, its magnetic field, the composition of its atmosphere and the interstellar absorption (see e.g. [9]). On the one hand, the magnetic field of isolated NS is likely to be large (B > 10 9 G) and thus will affect their spectra, and the chemical composition of their atmosphere is unknown and difficult to determine. On the other hand, NS that undergo periods of accretion of matter from their binary companion are believed to have a low magnetic field (due to accretion-induced decay), and an atmosphere likely to be composed of light elements (H, possibly He [10,9]). Among such objects one can distinguish quiescent X-ray transients (QXT), NS in a binary system observed when the accretion has stopped or is strongly reduced, and bursting NS (BNS) i.e. NS from which recurring and very powerful bursts, so-called photospheric radius expansion (PRE) bursts, are observed. These sources are even more promising when they are located in globular clusters whose distance is likely to be accurately measured. R ∞ can also be constrained by the modelling of the shape of the X-ray pulses observed from rotation-powered radio millisecond pulsars (RP-MSP) in particular if their mass is known from radio observations. Figure 1 shows the most recent constraints on the radius R 1.4 of a 1.4M NS obtained for various types of sources (see details in [11]). The constraints QXT-1 and RP-MSP being mutually exclusive, so far no consensus on R 1.4 can be reached. However, the determination of the radius of a NS is subject to many assumptions, uncertainties and systematics effects, (see e.g. table 1 in [8]). Obtaining constraints from the PRE bursts of BNS is still subject to uncertainties and debates in particular concerning the modelling of the phenomenon itself, the selection of bursts to be used (hard state X-ray bursts vs. soft state ones) and the composition of the atmosphere (see [12][13][14][15][16][17]). As far as QXT in globular clusters are concerned, the  [22], BNS-1 [12], BNS-2 [15], QXT-1 [23], BNS+QXT [14]. Constraints QXT-2 and QXT-2 are included for discussion only, see text for details. The constraints correspond to 2-σ error bars.
composition of the atmosphere and the amount of interstellar absorption, quantified by the so-called "equivalent hydrogen column density" N H , are unknown and significantly affect results [18,19,9]. For example among the five QXT studied in [18] (constraint QXT-2 in fig. 1), one of them: NGC 6397, has a substantially smaller R ∞ compared to the value obtained for the four other sources. As a consequence the constraint QXT-2 that is derived from these five sources suggests a small NS radius. However, while for one of the five sources observations suggest a hydrogen composition for the atmosphere: ω Cen [20], the composition of the atmosphere of the four other sources, including NGC 6397, is still unknown. Using a helium atmosphere instead a hydrogen one, a larger R ∞ is obtained for NGC 6397 [9]. The constraint QXT-2 corresponds to the QXT-2 one when NGC 6397 is not included: it then favours larger radii. The quantity N H can in principle be constrained thanks to observations in various wavelengths or derived when fitting the X-ray spectra. Large discrepancies between the values derived for N H using these two approaches are however observed and as a consequence the constraint on R ∞ can vary by as much as a factor 2 (for ω Cen in [18]). Finally, the uncertainty on the distance to globular clusters can be as large as 25% [19] and further affects the constraint on the radii (see, e.g., [18,19,9]). Last but not least, taking into account NS rotation strongly complicates the analysis of the collected X-ray spectra. Both QXT and BNS are likely to rotate at a frequency of few hundreds of Hz which is expected to affect the radius determination by ∼ 10% according to [13,21].
Due to uncertainties in both the observations and the modelling of QXT, BNS and RP-MSP, no stringent constraint on the radius of NS can currently be derived. However the next generation of X-ray telescopes such as NICER [24], Athena [25] and, possibly, a LOFT-like [26] missions promise measurements of radii with an accuracy of few percents. Together with a constraint on the maximum observed NS mass, a simultaneous measurement of NS mass and radius with such a precision could enable to constrain the NS EOS.

Spin frequency
Since the discovery of the first pulsar in 1967 [27], later identified as being a NS [27,28], ∼ 3000 NS have been observed in all wavelengths, most of them as radio pulsars (see [29] for a review). Among those one can distinguish two populations: the so-called "normal pulsars" with periods of the order of few seconds and the "millisecond pulsars" (MSP) which as their name indicate have a period of the order of few milliseconds [30]. These are believed to be old NS that have been "recycled" i.e. spun-up to millisecond periods by the accretion of matter from a binary companion [31,32].
During the recycling process, a binary system can be observed as an X-ray source and its pulsar as an X-ray millisecond pulsar (XMSP). The spin frequency can be determined or estimated for three different types of XMSP [33]: -Accreting X-ray millisecond pulsars (AXMSP): X-ray pulsations due the presence of hotspots at the surface of the rotating neutron star have been observed from these sources. The spin frequency of 15 AXMSP has been measured with a great accuracy (see [34] for a review). -Nuclear X-ray millisecond pulsars (NXMSP): they exhibit oscillations during thermonuclear X-ray bursts. The frequency of the oscillations is thought to be at or close to the pulsar spin frequency, though there are still some uncertainties on the physical process that triggers the oscillations [35]. Therefore for these sources, the measurement of the spin frequency is indirect and has an uncertainty of few hertz. The spin frequency of 10 NXMSP has been determined so far; -Twin kilohertz quasi-periodic oscillations have been observed in several systems. However, their interpretation and the precise link with the rotation of the neutron star is still unclear (see, e.g., [35]). Figure 2 shows the frequencies of currently observed radio and gamma-ray pulsars (data from the ATNF Pulsar Catalogue 2 [36]) and XMSP rotating at a frequency larger than 100 Hz. Out of ∼ 2500 radio and gamma-ray pulsars with measured period, 11% of them have a spin frequency larger than 100 Hz. So far, the fastest rotating XMSP is 4U1608−522 with f = 620 Hz [37] and the fastest rotating MSP is PSR J1748−2446a in the globular cluster Terzan 5 with f = 716 Hz [38]. Oscillations at a frequency of 1122 Hz in one type I X-ray burst of XTE J1739−285 were reported [39] but not observed later [37].
Although with current observational techniques submillisecond pulsars could in principle be detected, so far all attempts were unsuccessful (see, e.g., [40][41][42]). Thus this might indicate the existence of a mechanism that 2 http://www.atnf.csiro.au/people/pulsar/psrcat. prevents accreting NS from reaching submillisecond periods. The interaction between the NS magnetic field and the accretion disk (see, e.g., [43]) or the loss of angular momentum due to the emission of gravitational radiation (see, e.g., [44,33,45]) may inhibit the recycling process and thus the formation of submillisecond pulsars. The consequences of the existence of the fastest rotating pulsar with f = 716 Hz are discussed in sect. 6.

Rotating stars in general relativity
After a few decades since the discovery of the spherically symmetric, static solution for matter distribution by Tolman,Oppenheimer and Volkoff ([46,47]) -the TOV equations-the interest of researchers on rotating, relativistic stars was revived in the 1960s Golden Era of GR. The breakthrough came just in time for the discovery of the first pulsar with the work of Hartle [48], who devised a slow-rotation approximation to an exact solution by treating rigid rotation as a small perturbation of the spherically symmetric TOV background solution. In quasi-Schwarzschild coordinates (t, r, θ, φ), such a stationary, axisymmetric spacetime is described by a generic metric 3 where the ω function is the angular velocity of the freefalling observers (which corresponds to the frame dragging of inertial frames). In the first-order expansion in terms of the star's angular spin frequency Ω, outside the star the metric function ω = 2J/r 3 , i.e., it is proportional to the total stellar angular momentum J. With respect to the metric used in the TOV solution, the O(Ω) metric differs only by the ωdt term, In addition to the usual TOV ordinary differential equations for m(r), ν(r) and pressure p(r) the off-diagonal tφ component of Einstein's equations provides a differential equation for ω(r), or equivalently, J (which then can be used to define the moment of inertia of the star I = J/Ω). Subsequently, Hartle and Thorne [49], and Sedrakyan and Chubaryan [50] obtained a second-order, O(Ω 2 ) solution of the slow-rotation approximation. Within this approximation, the star's angular momentum, the fluid velocity and the frame-dragging term are exactly the same as in the O(Ω) order (they are functions of odd powers of Ω). What is affected are however the diagonal metric terms and pressure and energy density distributions. The metric reads where h 0 (r) and m 0 (r) describe the monopole deformations, and h 2 (r), m 2 (r) and v 2 (r) describe the quadrupole deformations (the dipole term is identically zero). The function P 2 (cos θ) is the Legendre polynomial. The pressure inside the star is modified as follows: p(r, θ) = p(r) + (ρ + p) (p 0 + p 2 P 2 (cos θ)) , with ρ denoting the energy density. Similarly to the firstorder expansion, the solution consists of the TOV background solution, supplemented with an additional set of first-order ordinary differential and algebraic equations for the monopole and quadrupole terms. The gravitational mass M of a star with angular momentum J is where M and R are the mass and radius of a non-rotating configuration with the same ρ c . Rotational corrections deform the star to a spheroid shape. The θ dependence of the radius is with ξ 0 and ξ 2 being functions of p 0 (r), p 2 (r), v 2 (r), h 2 (r), as well as the equation of state and structure. The equatorial circumferential radius equals for which one can evaluate the exterior solution of g φφ at the surface; R e is the coordinate equatorial radius obtained by integrating the equations of structure. A second-order slow-rotation approximation allows for defining the star's quadrupole moment Q by comparison of the metric terms with their Newtonian analogues. Consequently it allows to characterize the exterior metric of a rotating object using its multipole moments: gravitational mass M , angular momentum J and quadrupole moment Q. This, together with the fact that slow-rotation approximation offers a direct and intuitive link to a spherically symmetric Schwarzschild configurations (by using the same system of coordinates), as well as a relative simplicity of the set of ordinary differential equations to solve is its biggest advantage.
There are drawbacks of this approach, however. Slowrotation approximation cannot be applied to all rotation rate: by definition Ω/Ω K 1, where Ω K = 2πf K is the Keplerian (mass-shedding) angular frequency. Moreover, the very definition of a spheroid shape in eq. (7) prevents the star to accurately reproduce the mass-shedding limit. Thirdly, O(Ω 2 ) definitions of the gravitational mass M and angular momentum J are in general not accurate enough to robustly indicate an instability (by means, for example, the arguments based on the turning-point theorem of [51][52][53], which is a sufficient condition for instability).
To obtain accurate results for any rotation rate, one needs to change the way the problem is posed. On these grounds highly accurate numerical schemes like [54][55][56] were developed. Among other things, it is useful to abandon the Schwarzschild coordinates in favor of e.g., quasiisotropic coordinates. The metric is then expressed as where the ν potential is as previously related to the gravitational potential of the source, and e μ and e ψ are called conformal factors. Note the different relation between the r and θ coordinates with respect to, e.g., eq. (3). The Einstein equations to be solved can be derived in a number of ways. A particularly interesting, widely accepted and successful approach in numerical relativity is the 3+1 decomposition of spacetime i.e., a specific "slicing" of the fourdimensional spacetime into spacelike three-dimensional hypersurfaces in order to deal with the three-dimensional tensor fields to obtain solutions (see, e.g., [57,58] as well as [59] and [60] for the specialized case of rotating relativistic stars). In the quasi-isotropic gauge with a choice of slicing (maximal slicing), the Einstein equations for a stationary, axisymmetric star are expressed as a system of four coupled non-linear elliptic partial differential (Poisson-like) equations: where the right-hand sides of each equation, σ i , are the source terms describing non-linear metric terms and matter via the energy-momentum tensor, usually assumed to describe the perfect fluid. The only boundary condition for this system exists in spatial infinity and is provided by the asymptotically flat metric. Global quantities, such as gravitational mass M are formulated as surface or volume integrals using the asymptotic behavior of appropriate metric functions far from the source (the so-called ADM mass [61]), or by exploiting the symmetries of the problem (the so-called Komar mass for stationary spacetimes by taking into account the existence of a time Killing vector, [62]). Similar reasoning applies to the angular momentum J. Redefining e ψ as Br sin θ, one can get the relation between the coordinate and the circumferential radius. Analogous to eq. (8), the circumferential equatorial radius of the star, i.e., the length of the equator divided by 2π, is simply with R e being the coordinate equatorial radius. In what follows we will denote the circumferential equatorial radius of the rotating star by R eq .

Accuracy of solutions
An additional important aspect of solving numerically the Einstein equations is the choice of numerical methods. In some applications the spectral methods prove to be superior over traditionally more widely used finite-differences methods. With the spectral decomposition, functions are expanded in terms of adequately chosen basis functions, and resulting algebraic equations for the expansion coefficients are solved. When properly implemented, the difference between the series expansion and the real solution vanishes like e −N , where N is the number of expansion coefficients (the evanescent error). As an example of a real implementation, the formulation of [54] is using spectral methods in a numerical library LORENE 4 in a nrotstar code; other highly accurate implementations can easily reach machine precision [56]. Note also that numerical relativity with spectral methods is particularly suitable for precision studies of instabilities, due to very low numerical viscosity of spectral methods (see e.g., [65][66][67][68] and sect. 7.1 of this article). The accuracy of a numerical solution may be checked in a number of ways. For example, for stationary asymptotically flat spacetime, the Komar mass is in principle equal to the ADM mass -their difference expresses thus the imperfection of numerical solution and is proportional to the accuracy achieved. A very sensible and widely used accuracy indicator is a relativistic generalization of the classical virial theorem by [69,70], applicable in case of asymptotically flat four-dimensional spacetimes such as in the case of rotating relativistic stars.
As an example, figs. 3 and 4 show the effect of rigid rotation on the mass-radius M (R) sequences for two recent EOS -a "standard" nucleonic DH EOS of [63] and stiff TM1C EOS that includes hyperons. The latter is a non-unified EOS (see discussion in sect. 5.1) where the 4 http://lorene.obspm.fr. 1.4M , leading to large radii of NS models for this mass range, is explained in sect. 5.1.

Fig. 4.
Non-spherical shapes of NS deformed due to rapid rotation. Isocontours denote constant fluid proper energy density; vertical direction is aligned with star's angular momentum. Thicker line correspond to the surface. Upper slice corresponds to the TM1C EOS NS with M = 1.4M , rotating at 716 Hz which for this configuration is the mass-shedding limit (Keplerian rotation; note the cusp at the equator). Lower slice corresponds to the DH EOS with the same mass and spin frequency. The configurations correspond to stars marked on fig. 3. The circumferential radius of TM1C NS is 20.9 km, whereas for DH NS it is 12.9 km. The results were obtained using the LORENE/nrotstar code.
DH EOS is used for the crust and the model by [64] for the core. Configurations were obtained with the use of LORENE/nrotstar. Specifically, fig. 4 shows how strongly the shape of a rapidly rotating relativistic star depends on the EOS. Note that the configuration at the verge of a breakup (mass-shedding limit, which is related to a cusp on the stellar surface at the equator) cannot be accurately simulated with the slow-rotation approximation.
Apart from the shape, rigid rotation changes the values of global parameters of a star; it increases its equatorial radius R eq and also the maximum gravitational mass M max with which the star can still be stable for a given central EOS parameters. The maximum increase of M max is, with respect to the non-rotating configurations, about 15-20% for dense matter described by realistic hadronic EOSs. Note that the mass increase caused by rotation cannot be therefore proposed as a general solution to the problem of maximum mass decrease caused by a substantial phase transition/softening in some exotic EOS, in order to reconcile their inconsistency with recent observational data. The maximal increase of R eq is about 30-40% for the mass-shedding configurations (see sect. 6 for more details).
Exact numerical solutions were compared with the slow-rotation approximation in a number of articles. Notably, differences between the results obtained using the early implementation of LORENE/nrotstar [54] and the results obtained by [71] are described in [72]; the differences in mass for rotating maximum mass models is reported to be of about 5%, whereas the differences in radii about 15%. The slow-rotation approximation is particularly sensitive to quantities that depend on the derivatives of the metric functions, e.g. the radius or the innermost stable circular orbit around a rotating compact star. In a comparison performed in [73], where the radii of orbits around strange quark stars were studied, the difference between the slow-rotation approximation and exact results for moderately rotating stars (f 500 Hz) at the canonical mass of 1.4M is of the order of 1 km. The discrepancy grows with stellar mass and spin frequency. 4 Transition to an exotic core: equilibrium, metastability, and instability

Thermodynamic equilibrium considerations
Some theories predict that with increasing density the NS core undergoes a transition from a normal (N) state to a new exotic (E) phase. Some of these predicted phase transitions are second-order ones (e.g. kaon condensation, pion condensation), so that density and composition of the matter are continuous at the transition point P t , n t , while the speed of sound drops discontinuously c E < c N . However, for many models the softening in the E-phase just after threshold is so strong that it results in dP/dn < 0, and therefore induces a density jump between the N-and E-phases, coexisting at some P NE (see fig. 5). In this way, we have effectively a first order phase transition at constant P = P NE between the N-phase and the E-phase. This occurs for instance for a sufficiently strong pion or kaon condensation in nucleon matter. Another example of a genuine first order phase transition is quark deconfinement in dense hadronic matter. In general, in the first order phase transition at constant pressure, N-and E-phases Possible EOS in the n-P plane with various types of transition between the N-and E-phases. If the surface tension at the N-E interface is large: a mixed state is not present, and the N-phase is stable up to P NE (solid N-line). For P > PNE and n > nN the N-phase is metastable with respect to the nucleation of the E-phase (dotted line). The rate of nucleation grows rapidly with overcompression P − P NE, and at P = P crit N , n = n crit N the E-phase nucleates. After equilibration we get the E-phase coexisting with the N-phase at pressure P NE, with the density jump nN −→ nE (thin dotted line horizontal segment, constant pressure). After further compression the EOS continues in a pure E-phase (thick solid line). If the surface tension at the N-E interface is small : the equilibration produces a mixed state of the E-and N-phases, starting at P are separated by a surface with a surface tension σ > 0. We define the baryon chemical potential in a given phase as where E is the energy density (including rest mass energy of particles). Thermodynamic equilibrium at a given P is realized by a state (phase) of the dense matter with a minimum value of μ b . However, it has been shown by Glendenning [74], that if the surface contribution to the energy is sufficiently small, the state of thermodynamic equilibrium (i.e. with minimum μ b (P )) is actually a mixture of coexisting phases E and N (mixed m state) in the pressure interval P

Central compression, metastability, and nucleation of the exotic phase
Consider an element of matter at the center of a NS, consisting of the N-phase at density n c close to, but smaller than n (m) N . The value of n c can change due to the NS evolution induced by: 1) angular momentum loss due to dipole radiation from the radio pulsar; 2) mass gain in the process of matter accretion from a companion in a close binary system. In most cases it is an increase of n c (compression) that proceeds on a certain timescale τ comp = n c /ṅ c . The increase of n c during the spinning-down of an isolated pulsar depends on the mass of the star and reaches 5-30% for evolution from the Keplerian to the non-rotating configuration [75]. This increase however is proportional to the square of f initial /f Kepler , and for an initial period ∼ 10 ms it is less than 1%. The timescales involved are longer than 1 Gyr. In accreting NS the crucial parameter is the total angular momentum transferred to the star and the magnetic torque due to the interaction with an accretion disk (for details see [75]). The central compression is small (few percent) in the absence of a magnetic torque and for accretion from the marginally stable orbit, but could be as large as 10% for B ∼ 5 × 10 8 G after accretion of 0.1M . It would take 10 Myr for a NS accreting at the rate 10 −8 M /yr in a low mass X-ray binary. The compression timescale becomes rather short (years) for a very specific class of young magnetars (see [76] and references therein).
In scenarios (1) and (2) the temperature at the NS center is < 10 9 K and therefore thermal contributions to thermodynamic quantities are small: matter is strongly degenerate, and thermal fluctuations are negligibly small compared to the energy barriers separating the N and m states. In any case, the transition to the m state has to be initiated by a droplet of the E-phase. Even neglecting the surface tension contribution to the energy of an E-drop (σ = 0), it is energetically possible only for n b > n N (at lower n b the drop decays back into the N-phase). However, to nucleate the E-phase in the N-medium, an energy barrier, created actually by the surface contribution, has to be quantum-mechanically penetrated with the energy supply coming from quantum fluctuations. After an E-drop nucleates, it grows into a bulk E-phase which coexists stably with the N-phase at P NE , with a density jump n N −→ n E at the phase interface.
Sometimes it may be convenient to visualize the evolution of the NS center in the P -μ b plane, fig. 6. Assume that the central core is being compressed, so that P c grows in one of the astrophysical processes described in the first paragraph of the present section. Compression corresponds to a trajectory in the P -μ b plane. Even after passing P c = P (m) N , the core with P c > P > P (m) N consists of a pure N-phase and grows in time, because the mixed state m cannot be reached due to the impossibility of nucleation of the E-phase because of the energy barriers (resulting from surface tension and Coulomb interaction). After reaching P c = P NE , the star's center enters the metastable (overcompressed) segment P > P NE of the N-phase EOS. The lifetime with respect to nucleation of the E-phase τ nucl decreases rapidly with growing overpressure ΔP = P c − P NE . Nucleation of quark matter in dense baryon cores was studied in [77][78][79], while the nucleation of the pion-condensed state was discussed in [80][81][82]. As soon as τ nucl ∼ τ comp (which takes place at P c = P crit ), droplets of the E-phase appear spontaneously and grow into regions of the E-phase coexisting with the N-phase. This kinetic (non-equilibrium) process implies a local pressure deficit and a collapse of the central core. If the heat release is sufficiently large, one can contemplate a redistribution and coagulation of the E-droplets in the N-phase, creating a m-state core, larger than a uniform E-core could have been.
Such a situation is schematically depicted in fig. 7. The central core is being compressed while still in the Nphase, until at P c = P crit the E-phase nucleates in the Nmedium. Assuming that thermodynamic fluctuations are sufficiently strong, one gets a mixed m state extending down to the (P (m) N , ρ (m) N ) point. This corresponds to a full thermodynamic equilibrium in the core. However, before this final state has been reached, nucleation of the E-phase at P crit implied a local pressure deficit, and a collapse of a NS into a new more compact configuration took place. After the E-and N-phases mix, a large m core is formed in configuration C * , with a mean central density ρ * c > ρ crit and P * c < P crit . The dynamics of this minicollapse process induced by a core-quake is discussed in sect. 8.

Hyperonic cores
Hyperons (baryons containing at least one strange quark) were discovered in laboratory in the early 1950s and are studied experimentally since then. In the late 1950s it was suggested that hyperons could also be present in NS cores Fig. 7. Trajectory of the neutron-star center in the ρc-Pc plane during spin-down or accretion, leading to a core-quake after the nucleation of the E phase (configuration C crit ) which implies a collapse into a configuration C * with a mixed-phase core. The baryon number A is conserved in the collapse process. C 0 is the last strictly stable configuration with a N-phase core. For further explanation see the text.
(see [6,83] for a historical perspective). Indeed, although hyperons are unstable under terrestrial conditions, at densities typical for the NS centers, the Pauli exclusion principle prevents them from decaying into nucleons.
Relatively little is known about the properties of the interactions of hyperons with other baryons from hyperon scattering (see discussion in [84,85]). On the one hand, thanks to the study of Λ-hypernuclei and Ξ-hypernuclear states in laboratory, the potential for the Λ and Ξ hyperons in symmetric nuclear matter at saturation density is found to be attractive. But on the other hand, contradicting results were found for the Σ potential. Moreover, only few double-Λ hypernuclei were studied indicating an attractive Λ-Λ potential and no other pairs of hyperons as Λ-Ξ or Ξ-Ξ were observed.
Allowing for a possible transition to hyperonic matter results in a softening of the EOS at densities n b > 2n 0 and, as a consequence in a decrease of the maximum mass. This is illustrated in fig. 8, where M -R relations at various spin frequencies are plotted for the DH EOS and for the EOS obtained within the RMF (relativistic mean field) approach with the TM1 model [86], as an example. For the latter four different EOS for the core are shown. The first model noY corresponds to a purely nucleonic composition. Three models allowing for a transition to hyperonic matter at high density are plotted. In the Y model, the vector-isoscalar hidden-strangeness φ mesons is included and vector-mesons couplings to baryons are given by the SU (6) symmetry following [87]. The φ meson which is coupled to hyperons only, yields an additional repulsion between hyperons and thus leads to a stiffening of the EOS and an increase of the maximum mass. The Yss includes in addition to the φ meson, the scalar-isoscalar hidden- strangeness σ * meson in SU (6) symmetry following [87]. This meson which also couples to hyperons only, enables reproducing a weakly attractive Λ-Λ potential [88]. The consequence of its inclusion is a mild softening of the EOS and thus a slight decrease of the maximum mass as can be seen from fig. 8 for non-rotating stars. Finally, the Yssz model corresponds to the Yss one except that SU (6) symmetry is broken following [64]. The effect of this breaking is studied in detail in, e.g., [89]. The specific choice of parameters in the Yssz model, results in a stiffening of the EOS, which becomes even stiffer than the Y model although σ * mesons are included. As a consequence the maximum mass is the highest of all the hyperonic models but nevertheless lower than the one for a purely nucleonic star.
Observations of massive neutron stars with a mass 2M [7,2] are therefore challenging for hyperonic EOS. Reconciling the possibility of a transition to hyperonic matter at high density with observations of massive NS requires solving the so-called "hyperon puzzle". In [11], a systematic study of all EOS for hyperonic matter, consistent with a 2M maximum mass, available at the time of publication, was conducted (14 EOS, all but one being RMF models). It was shown that all of them give pressures in pure neutron matter at densities close to n 0 which are too large compared with recent precise many-body calculations for pure neutron matter. These calculations were performed using two different approaches: quantum Monte Carlo method [90] and chiral effective field theory [91], and are in a remarkable mutual agreement. A large pressure for n b n 0 2n 0 in fact is needed to balance the hyperon softening at higher density, and is correlated with large radii: R > 13 km for neutron stars with masses M = 1.0-1.6M . Hyperonic EOS consistent with a 2M maximum mass and with a pressure close to n 0 consistent with [90,91] are obtained in e.g. [85,92]. In [93] a RMF model, so-called KVOR, where hadron masses and coupling constants are scaled by functions depending on the scalar field is formulated. Hyperonic NS with M ≥ 2M are obtained in two versions of this model. The MKVORH one [94] assumes a relatively low value of the nucleon effective mass at saturation and yields M max = 2.2-2.3M and R 1.4 = 12.2 km while the KVORcut3 version [95] has higher effective mass at n 0 and M max = 2.0-2.3M for R 1.4 = 13.0 km. However, even at maximum mass the strangeness per baryon is very small: 3 × 10 −2 . A solution to the "hyperon puzzle" is therefore reached in [94,95] when there are nearly no hyperons. Figure 8 also presents M -R relations for NS rotating at 700 and 900 Hz together with M -R relation for NS rotating at the mass-shedding limit f K . Let us now compare results for the DH and TM1 EOS. On the one hand, the increase of the maximum mass due to rotation is larger for the TM1 EOS, with larger radii for non-rotating configurations, than for the DH EOS. On the other hand, the increase of the minimum mass, which is located at the intersection between the f K curve and the M -R one, is larger for the TM1 EOS, for a given rotation rate. In other words the DH EOS, which has smaller radii for nonrotating models than the TM1 one, has for a given rotation rate a broader range of masses than TM1. Under the effect of rotation the M -R relation for TM1 become flatter than the DH EOS. The property that M -R relations become flatter for EOS with larger radii for non-rotating configurations than for those with smaller radii is even more dramatic for very high spin frequency, close to the Keplerian frequency [96]. For the TM1 model, the softer the high density part (e.g. when comparing the Yss EOS to the Y one), the smaller the increase of M max and the larger the increase of R(M max ) with rotation. As a consequence for a given rotation rate, the M -R curve is flatter for the softer TM1 EOS and the narrower is the stellar mass range.
There exists very few unified EOS, in the sense that the same nuclear interaction model is used to describe both the clusterized matter in the crust and the homogeneous one in the core [63,97]. Two non-unified EOS for the crust and the core are usually "glued" together by ensuring that the pressure P (n b ) and energy density ρ(n b ) are increasing functions of the baryon number density n b . However there is no unique prescription for the transition between two different EOS [92]. For example in fig. 9 three possible choices of "gluing" for the same crust and core EOS: the DH [63] and NL3 RMF models [98] are shown. The NL3l EOS corresponds to gluing DH and NL3 at the density at which the P (n b ) relations for the two EOS cross: n b = 0.046 fm −3 . For the NL3h EOS, the core EOS is glued to the crust at n b = 0.16 fm −3 . Finally, the EOS NL3u corresponds to using the same nuclear model for both the crust and the core for the NL3 parametrization. The crust model is taken from [99], where the Thomas-Fermi approach is used to describe the clusterized matter in the crust. Similarly, unified EOS for the TM1 model are shown in fig. 8. Figure 9 shows the M -R relations obtained for the three NL3 EOS: for non-rotating stars the difference in the equatorial radii between the EOS de- creases when the mass increases but can be as large as 4% of the radius for a 1.4M NS and 3% for a 1.6M NS. For a given mass the difference in radii between various matching prescriptions increases with spin frequency. For example for a NS rotating at 700 Hz the difference in radii increases to 13% and 6% of the radii for a 1.4 and 1.6M NS, respectively. Therefore calculations of unified EOS is of great importance in order to properly describe the macrophysical properties of NS [92].

Hybrid stars
Some theories of dense hadronic matter predict a deconfinement of quarks at densities achievable in the cores of massive NS. The phase transition from the baryon phase of matter (N) to the deconfined quark matter (Q) is usually assumed to be of first order. It softens the EOS due to the density jump at constant P NQ (transition between pure N and Q phases with a density jump at the interface λ = ρ Q /ρ N ), or via a mixed-phase region. In spite of this softening, the existence of 2M pulsars does not exclude quark cores in NS, but imposes rather tight constraints on the EOS with N-Q transition [100,101]. Two subsequent phase transitions through the intermediate quark phase energetically preferred in a rather narrow range of densities were also considered [102][103][104][105].
In all these cases the density jump (or two jumps) is the main source of softening of matter at the transition pressure P NE determining M (P c ), the M (R) dependence in the vicinity of the configuration C 0 with central pressure equal to the phase transition pressure P c = P NE .
It should be mentioned that the baryon-quark phase transition between two phases can proceed through a mixed phase [106,74] in which the condition of local charge neutrality (for the two phases separately) is relaxed. Although the properties of stars (mass-radius relation) are not the same for an EOS with phase transition involving a mixed or pure phases, significant differences are observed for a rather small region of central pressures (close to the transition pressure). Global properties such as the existence of instability regions, "twins" and maximum mass of hybrid star are very similar [107].
The conditions for the quark EOS resulting from the maximum mass constraint (M max > 2M ) are the following: the softening effect of the first order phase transition has to be compensated by a stiff quark EOS (larger λ corresponds to a larger sound velocity c Q of the quark phase, [104]). The EOS for the baryon phase should also be relatively stiff -the configuration C 0 cannot be too compact and close to the maximum mass for the N-phase stars. If the latter requirement is not fulfilled, the phase transition at the center leads almost immediately to the dynamical instability and gravitational collapse of a hybrid star. For phase transitions to quark matter at densities 2.5-3.5n 0 , the sound velocity in the quark phase 5 should be larger than 0.6c.
Benic et al. [109] presented the possibility of the existence of high-mass twins -two families of dense objects (baryon and hybrid stars) with masses about 2M . In their model two conditions discussed before are fulfilled: the matter in the NJL8 quark phase is stiff (c s ∼ 0.6-0.9c), and the compactness of the C 0 configuration for 2M is not large, but comparable with the compactness of a typical nucleon NS with M = 1.4M and R 11 km.
In fig. 10 we present the M -R eq relations for the model with a first order phase transition to quark matter (Q, the black curve). The effect of rotation for hybrid stars is similar to that of hyperon stars ( fig. 8), although in this case the relatively large stiffness of the quark matter EOS leads to the maximum mass larger than for a baryon star. The minimum mass of the hybrid star rotating with f = 700 Hz is larger than 1M , much larger than for the DH EOS, and comparable to the hyperon stars discussed in sect. 5.1. We should stress however that this minimum mass (Keplerian configuration for f = 700 Hz) corresponds to the star composed entirely of non-strange, nucleon matter (P c < P NE ). The large value of M min is an indirect consequence of the EOS softening at high density and the M max > 2M requirement, which results in a large radius (i.e., small compactness) of the C 0 configuration.

Keplerian (mass-shedding) limit
Consider a static (f = 0) NS of gravitational mass M s and baryon mass M b . Then construct a sequence of rigidly 5 These results were obtained for a simplified quark EOS with constant sound velocity. It was shown however, that this form of P (ρ) dependence approximates well the EOSs obtained for more sophisticated, microscopic models of quark matter [108,104].  [110] with and without phase transition to quark matter. Quark phase is approximated by linear dependence P = a · (ρ − b) with a = 0.5c 2 and the density jump at the phase transition λ = 1.2. For details see [108,104].
rotating configurations with the same M b and increasing f . This sequence will terminate at the mass-shedding limit with the Keplerian frequency f K . At this limit the spin frequency of the NS is equal to the orbital frequency of a test particle on a circular orbit corresponding to the NS equator. For f > f K hydrostatic equilibria of a NS with baryon mass M b do not exist.
As we have seen in sect. 5, the EOS that predicts hyperonization and therefore a softening after the hyperon threshold at ρ 2ρ 0 , has to contain a sufficiently stiff nucleon segment ρ 0 < ρ < 2ρ 0 , in order to get M max > 2M in spite of the hyperon softening. This implies that the radii of the pre-hyperonic NS with ρ 0 < ρ c < 2ρ 0 , and in particular the radius of 1.4M NS with this EOS, R (H) 1.4 , are expected to be larger than those of "standard" NS with nucleon cores, R (N) 1.4 . A rough Newtonian argument about the centrifugal force ∝ R 2 eq , predicts therefore that for the same moderate M s mass: 1) the Keplerian frequency for EOS with hyperonization is smaller than for purely nucleonic EOS, 2) the difference R eq grows with increasing f (see fig. 11). It should be stressed that both properties are valid for M < 0.9 M (stat) max [11]. The results reviewed above can be corroborated quantitatively by theoretical considerations. For a non-rotating spherically symmetric NS of gravitational mass M and circumferential radius R the orbital frequency of a test particle in a circular orbit of radius r orb > R is [111] f orb = 1 2π According to the Birkhoff theorem the formula is the same as for a central point mass M . For a static, spherically symmetric NS of mass M we therefore get which is identical with the Newtonian formula for a massshedding limit for a spherically symmetric self-gravitating star. It has been shown, with a precision of a few percent, that the formula for f K , eq. (13) holds also for realistic NS at the Keplerian limit provided we replace R by the equatorial radius at the Keplerian limit R K [96]. The high precision comes as a surprise, because a NS at the Keplerian limit is strongly flattened and rapid rotation produces strong frame-dragging effects in the exterior space-time. The formula which is (surprisingly) so precise for NS, but not close to maximum allowable mass: M s < 0.9 M stat max , holds strictly for the relativistic Roche model with extreme central condensation of mass [112][113][114]. It should be stressed that validity of this model breaks down near M stat max . Moreover, the formula (14) is much less precise for strange quark stars built of self-bound quark matter which are characterized by a rather uniform density [114].
For practical applications, one can use the empirical formula of ref. [114], valid for NS with and without an exotic core: where M s and R s are mass and radius of a static configuration of the same baryon mass as the Keplerian configuration. It holds for 0.5M < M s < 0.9M stat max and implies K (M s ) for M < 1.6M . Note that the empirical prefactor 1.08 kHz is larger than 1.00 kHz which corresponds to the relativistic Roche model.

Maximum frequency of stable rotation
The empirical formula for the absolute maximum of f for stably rotating configurations with a given EOS, f max ( [6] and references therein), is of a different character than eq. (15). Namely, the formula for f max results from an approximate but quite precise correspondence between two extremal configurations: the static configuration with a maximum allowable mass, M stat max , R M stat max , and another extremal configuration, which is stably rotating with a maximum allowed frequency (called maximally rotating configuration). This extremal configuration is stable with respect to mass-shedding and stable with respect to axisymmetric perturbations (large filled dots in fig. 3), (16) The prefactor 1.22 kHz is significantly larger than 1.08 kHz in eq. (15). Moreover, in contrast to eq. (15), the formula for f max is valid for both NS and self-bound quark matter stars. However, it can be used to constrain only static configuration with the maximum allowable mass. Assuming M stat max = 2.0M and R M stat max = 10 km one gets f max = 1725 Hz.

Constraint from f obs max = 716 Hz
Figures 3, 8 and 10 readily illustrate that the minimum mass of a rotating NS, M f min , is sensitive both to f and to the EOS. This is to be contrasted with the minimum mass of a static NS, which is rather weakly dependent on the EOS (provided it is a unified EOS, so that the crust and liquid core are calculated using the same nuclear interaction model). We get M stat min 0.1M (see [6] and references therein).

Back-bending and stability of rotating stars
The appearance of a new phase at the center of a NS results always in a softening of the EOS. For the global properties of NS the consequence of this softening is a slower increase of the stellar mass and total baryon mass as central pressure increases. For significantly strong softening of the EOS and sufficiently large star mass this feature can be observed in the case of the evolution of an isolated pulsar as the so-called back-bending phenomenon -the epoch at which the angular momentum loss due to the evolutionary processes leads to the spin up of the star. As a rigidly rotating star of a fixed total baryon mass M b loses its angular momentum J, the central pressure and density increase. For a strong softening of the equation of state above a critical density it is possible (for some range of M b ) that for a slowing-down star the central density crosses this critical value and the core of a new, dense phase develops in the center. Then the star shrinks with a significant decrease of the moment of inertia I and this has to be compensated by the increase of rotational frequency Ω to fulfill the equality δJ = ΩδI + IδΩ. The Fig. 13. The back-bending phenomenon for an EOS softened by a mixed phase at high density. Shown are evolutionary tracks for an isolated NS losing its energy and angular momentum J (the total baryon mass is fixed along each curve, M b1 < M b2 < M b3 ). For a sufficiently large mass (M > M b2 ), spin-up due to angular momentum loss is observed.
back-bending phenomenon (associated in fig. 13 with an "S" shape of the J(f ) curve) was proposed in [116] as a signature of a phase transition to an exotic state of matter in the center of a spinning-down pulsar (for example the appearance of a mixed-phase core at the star center). It should be however mentioned that a similar behavior can be also caused by hyperonization of dense matter, provided that the EOS softening above the hyperon threshold is strong enough [117].
The back-bending phenomenon is not the only consequence of the softening of the dense matter EOS. For a significant softening the result could be even more spectacular -a region of configurations which are dynamically unstable appears. Increasing the softening of the EOS finally results in a non-monotonic behavior of J along an evolutionary track with a fixed M b . Configurations in a region where J increases with increasing ρ c are unstable with respect to small axisymmetric perturbations (fig. 14).
The discussion of the existence of the back-bending phenomenon or the stability of the rigidly rotating configurations can be performed by the analysis of the extrema of three basic, macroscopic parameters of the rotating star that are well defined in GR -M , M b and J as functions of any variable which parametrizes stationary rotating stellar configurations, for example central density ρ c , pressure P c or equatorial radius R eq . The change in stability of a rigidly rotating configuration (from stable to unstable or vice versa) corresponds to an extremum of the two parameters from the {M, M b , J} set, with the third parameter fixed [53]. According to the turning-point theorem, the stability of a rotating configuration can be stated by checking that The back-bending phenomenon corresponds to the existence of a region in which, due to some evolutionary processes, the rotational frequency f increases with a decreasing total angular momentum J. It can be written as or, equivalently, The condition for the back-bending phenomenon in eq. (19) is similar to the second relation in eq. (17) with J replaced by f . The back-bending manifests itself in the existence of a local minimum of the M b (x) curves plotted for fixed frequency f . However in order to decide if the configurations which are subject to back-bending are stable one has to analyze M b (x) relations for fixed total angular momentum J.
In fig. 15 rotating configurations exhibit back-bending for frequencies f > f 2 and baryon masses larger than  In fig. 16 the case of strong softening is presented. There exists a region of dynamically unstable configurations defined by the condition (∂M b /∂n c ) J < 0.
The decreasing parts of M b (n c ) relations, marked by dotted lines in fig. 16, which correspond to unstable configurations define the instability strip separating two branches of stable, rotating configurations -one, less compact, with maximum mass defined by the onset of instability due to the softening of EOS and the second one with maximum mass corresponding to the threshold for the collapse to black hole. The analysis of the large set of EOS with softening through a mixed phase or first-order phase transition at constant pressure indicates that the existence of these two families does not depend on the rotational frequency (however, the width of the instability strip does depend on f ). Equivalently the instability strip starts at non-rotating configurations and continues up to the Keplerian limit-these two families are disjoint in the M (R eq ) or M (P c ) plane.
The turning-point theorem is a sufficient condition. Recent numerical simulations by [118] study the onset of the dynamical instability for rotating stars and obtain it slightly below the maximum mass (at fixed angular momentum J). The relative difference of M is of the order of 10 −3 for rapidly rotating configurations; for most astrophysical purposes the turning-point theorem is therefore precise enough to locate the instability regions by means of stationary calculations.

First-order phase transition and instability
(20) However, if the EOS exhibits a first-order phase transition (e.g. the Maxwell construction and phase transition to a new phase at fixed pressure P NE , see sect. 4.1) the P (ρ) relation is discontinuous at P = P NE with a density jump ρ N → ρ E . Let us stress that in the present section we assume full thermodynamic equilibrium (no metastability).
In this case an increase of the central pressure due to evolutionary processes, resulting in P c > P NE , leads to the appearance of a small core of the dense phase at the center of a star with a density jump (ρ N → ρ E ) at the boundary separating the two phases. For a non-rotating star this density jump results in the non-continuous change of the derivatives of the main global parameters of a star (X = M, M b ) with respect to the central pressure [119]: where λ = ρ E /ρ N is the density jump and x N ≡ P NE /ρ N c 2 . The derivatives with subscript N and E are taken, respectively, at a pressure infinitesimally smaller and larger than P NE . The direct consequence of eq. (21) is the stability criterion for a stellar configurations with a small core of the new, dense phase. The stability condition is given by the inequality [120,121]: We define a "weak" first-order phase transition with a relatively small density jump for which the condition in eq. (22) is fulfilled. This case is presented in fig. 10 with a phase transition to quark matter with λ = 1.2. The softening of the EOS manifests itself as a sudden change of the slope of the M (R) curve as a new phase of matter appears at the center. For a "strong" first-order phase transition (λ > λ crit ) the derivatives dM/dP c , dM b /dP c change their sign at P c = P NE and stellar configurations with a small core of dense phase in the center are dynamically unstable. The oscillatory mode (radial) which is in this case unstable has no counterpart for one-phase configurations.
The main feature of this oscillations is the flow of a matter through the pulsating boundary between two phases -the frequency of this mode is proportional to √ 3 − 2λ + 3x N , which directly means instability for λ > λ crit and collapse into a new configuration with a sizeable E-core [122].

Dynamics
The existence of two disjoint families has important consequences for the evolutionary tracks of isolated and accreting NS. The loss of angular momentum moves an isolated NS leftwards along constant M b lines on the M b (R eq ) plane ( fig. 17). Once it reaches the instability strip, it collapses into a more compact counterpart with the same total angular momentum J (arrows in fig. 17). The dynamical properties of this minicollapse were studied using general-relativistic (GR) numerical codes [76,123]. The prefix mini is reflecting the fact that the changes of stellar parameters (radius, moment of inertia, frequency of rotation) associated with collapse under consideration are usually small. However, for a strong first-order phase transition the radius can shrink by ∼ 10% (see figs. 17,18). The distribution of a specific angular momentum in a collapsing star differs from those in the rigidly rotating initial and final configurations, because a differential rotation profile develops during collapse. However, the degree of differential rotation is found to be small and the assumption about the conservation of the total angular momentum of the collapsing star J seems to be justified. The energy release in the minicollapse turned out to depend weakly on J. However, it strongly depends on parameters of the EOS and of the phase transition itself. The formation of an E-core is followed by NS pulsations and the generation of gravitational waves, as described in more detail in the following subsections.   . 18. Total baryon mass M b vs. circumferential equatorial radius R eq for stationary non-rotating NS (thick lines) and NS rotating at a fixed angular momentum J, for two models of the phase transition weak and strong. Dotted segments correspond to unstable configurations. The minicollapses are marked by arrows. The arrows above the C 0 configuration correspond to a minicollapse from the metastable (overcompressed) phase N .

Metastability, rotation, and the total energy release
The metastability of the N phase at the center of evolving NS can be parametrized by the overpressure ΔP = P crit − P NE at which the E-phase nucleates in the Nphase triggering a minicollapse from the initial configuration with central pressure P c = P crit > P NE .
Mini-collapses for different ΔP are schematically presented in fig. 18 in the case of weak and strong phase transition and for rotating and non-rotating stars. Because for a strong phase transition two branches of stable configurations are separated, the minicollapse is possible even without a overpressure ΔP = 0 (thick arrows in fig. 18) with a total energy release (difference) of the order of ∼ 10 51 erg. For a weak first order phase transition the overcompression and the existence of a metastable core is the only cause of a minicollapse and without overpressure the energy release is ΔE(ΔP = 0) = 0. An overcompression of the order of 5-10% leads to a ΔE ∼ 10 50 erg for a weak phase transition and to an increase of ΔE by ∼ 10 51 erg for a strong one [67,68]. The dependence of the energy release on ΔP is very weakly affected by the rotation rate and can be approximated with a very high accuracy by the values obtained for non-rotating stars (see [67,68]).

Astrophysical signatures of a minicollapse of neutron star
A phase transition in a NS core, that induces a dynamical minicollapse of NS, is associated with heating of stellar interior, matter flow, and NS pulsations. It is important to remind, that the very essential dynamical character of minicollapse is based on the assumption that the conversion of the matter into an exotic phase is of a detonation type. We assume that the N −→ E conversion front moves at supersonic speed, driving a shock wave in the N-envelope.

Surface glowing, and a delayed re-brightening
The dynamical and thermal effects were studied in hydrodynamical simulations of spherically symmetric minicollapses induced by pion-condensation in hadronic matter in [124], where the references to previous works can be found. The Newtonian approximation was used, and thermal effects, such as heating and cooling, were considered. The initial temperature of the NS interior was assumed to be 10 8 K. For an assumed model of pion-condensation, the total energy release ΔE was typically ∼ 10 50 erg (notice that this quantity is strongly dependent on the phasetransition model). This energy was liberated in ∼ 0.4 ms, and split into heating of the exotic core due to a latent heat and matter compression, and into heating of the Nmatter outside the E-core due to the compression. The kinetic energy of the matter flow was mostly imparted into the NS pulsations. The exotic (inner) core was heated to ∼ 10 11 K, and the outer core and crust to some 10 10 K. We expect that the shock wave strongly heated the NS surface [125][126][127] albeit this was not modelled in [124] due to a too low spatial grid resolution. The NS-core was rapidly cooled by neutrino radiation, so that after few hours most energy generated by the minicollapse has been carried out by neutrinos. By that time NS pulsations have been damped. The heat content in the N-envelope, with cooling timescale orders of magnitude longer than for the E-core, was diffusing to the NS surface, leading to its delayed brightening. This brightening was associated with surface X-luminosity L X increasing during 30 yr after minicollapse. Then L X dropped by many orders of magnitude, on a much shorter timescale, because the heat content of the N-envelope had been exhausted.

A gamma-ray burst?
A minicollapse due to pion condensation was proposed to explain a famous very energetic burst of gamma rays of 5th March 1979, coming from supernova remnant N49 in the Large Magellanic Cloud [125,126]. The initial gammaray pulse had a very short rise time of < 0.4 ms, indicating a dynamical character of the burst mechanism. It was then followed by a pulsed decaying photon flux, suggesting NS rotation (8 s period was however a puzzle at that time). A shock wave sent from the collapsed pion-condensed core propagated towards the NS surface, heating it to high temperature and allowing for an energetic burst of gamma rays. The model was further elaborated in [127]. These authors used the minicollapse parameters following the models of [128] as far as the energetics of the burst was concerned. The minicollapse scenario for the source of the 5th March 1979 event was further studied in [82].
An essential progress in the gamma-ray astronomy, started in 1990s, led to a conclusion that the source of 5th March 1979 was not a typical gamma-ray burster. It was actually an exceptionally strong outburst from a softgamma ray repeater (SGR). SGRs belong to a subpopulation of magnetars, slowly rotating (rotation period of several seconds, consistent with 8 s period of pulsations in the tail of the 5th March 1979 burst) and highly magnetised (10 14 -10 15 G) NS. A SGR emits GRB at irregular intervals, from hours to years and longer. The bursts from SGR are powered by magnetic field, and they are triggered by crust-quakes and are associated with magnetic field annihilation and reconfiguration. Only seven SGR are known today, and they are obviously not related to minicollapse which occurs only once in a NS lifetime.
Can a minicollapse due to a phase transition in NS core produce a class of regular non-repetitive GRBs that occur at cosmological distances larger than 100 Mpc, thousands of which were detected since 1990s? Unfortunately, here the minicollapse and shock heating model faces a very basic problem, pointed out using detailed numerical simulations in [129]. Consider a hot layer of relativistic plasma created by the shock wave at NS surface. It is an element of a fireball which is a precursor of GRB. In order to produce a typical GRB at cosmological distance such shock-produced fireball has not only to be sufficiently energetic (E f.b. > 10 51 erg) but should also contain a not too large amount of baryons (M f.b. b < 10 −5 M ), so that the so-called Lorenz factor b c 2 > 100 (see, e.g., [130]). In other words, the initial fireball should contain a nearly pure mixture of photons, neutrinos, and e + e − pairs that make it opaque: the contribution from the kinetic energy of baryons should be negligibly small. Unfortunately, even for the most optimistic models of the shock-wave propagation (no neutrino losses, no photodisintegration of nuclei) the best one can get is E f.b. = 10 46 erg for Γ f.b. = 40, which is only 10 −5 of the required energy [129].

A burst of gravitational waves
Such astrophysical signature of a minicollapse was studied in [76,123], which give also references to the previous work. The phase transitions considered were associated with quark deconfinement and kaon condensation. GR hydrodynamics was used in the 3+1 formulation, but thermal effects were not included. Both studies were concentrated on dynamics of minicollapse of a rotating neutron star and on the potential importance of minicollapse as a detectable source of gravitational waves (GW). The degree of differential rotation due to a minicollapse was found to be small. Pulsations induced by a minicollapse, when coupled to rotation, break the axial symmetry. This opens the possibility of a GW radiation in a burst (10-100 ms long), which if occurred at 10 kpc, could be detectable by the current second-generation interferometric detectors, like the Advanced LIGO (back online since September 2015), the Advanced Virgo (which will resume operations in the middle of 2016) [76,123], and the planned Einstein Telescope, a third-generation underground detector [131].

Effect of the crust formation scenario on the M-R relation
It is usually assumed that the NS crust is composed of cold, catalyzed matter. This is a good approximation when the outer layers of a NS are formed at the birth in a stellar core collapse, when T > 10 10 K allows for nuclear equilibrium in dense matter. However, for a NS that passed through the long stage of accretion of matter from a stellar companion in a binary system (e.g. NS recycled to millisecond pulsars in low-mass X-ray binaries) the crust is formed from the accreted layers of matter and its composition is widely different from that of catalyzed matter [132,133]. In view of these two possible formation scenarios we have different EOS of NS crust, composed of catalyzed or accreted matter. The latter EOS is stiffer than the former, which results in a different thickness of the crust and different radii, R acc (M ) > R cat (M ) [132,134].
The formation of a fully accreted crust would take ∼ 10 7 yr for a mean accretion rate 10 −9 M yr −1 typical of low-mass X-ray binaries. For a 1.4M star and R cat = 12 km one gets R acc −R cat 100 m. The difference in the equatorial radius grows faster than quadratically with the rotation frequency and depends quite strongly on R cat (f = 0). For 716 Hz, 1.4M and R cat (f = 0) = 12 km the difference in equatorial radii is 140 m [134]. This means that the effect of the formation scenario is significantly smaller than the uncertainties in the measurements of R (sect. 2).

Discussion and conclusions
In order to unveil the structure of dense matter of density up to ten nuclear densities, we confront theoretical models with measured NS parameters. Theoretical models are legion. The discovery of two radio pulsars of 2M resulted in an essential progress by putting strong constraints on the EOS, but did not produce a satisfactory answer to our fundamental question: do massive NS contain exotic cores? Crucial for solving this problem are precise measurements of NS radii with known masses. As for today, the determinations of radii are neither sufficiently precise nor reliable to give a definite answer. Hopefully, the situation may change in the future, thanks to the progress in X-ray astronomy. The task is very challenging, and requires, e.g., knowledge of distances to NS with precision of the order of two percent. It is regrettable that the LOFT mission will not fly in the near future, because it was offering some very special opportunities for NS radii measurements [26]. We can only hope that other missions like NICER [24] and Athena [25] will be successful in this respect.
All NS rotate, and many of those which are good targets for radius measurement are millisecond pulsars rotating at more than 400 Hz. There is a rather strong interplay between the rotation and the NS EOS, and therefore for both principal and practical reasons it is advantageous and very often mandatory to use accurate GR formalism, with a consistent choice of space-time coordinates, to calculate 2D hydrostatic equilibrium of rotating NS. The slow-rotation approximation is not suitable for checking stability, and does not allow the correct description of processes where fulfilling conservation laws is crucial. Therefore we encourage to use public domains precise 2D codes such as LORENE/nrotstar or RNS 7 in the studies that will eventually determine the true EOS of NS.
The use of precise 2D codes is particularly important for studying phenomena associated with softening of the EOS due to the appearance of the exotic phase. In particular, we explained this using the example of the backbending phenomenon and the loss/gain of stability in the spin evolution of NS.
Since the discovery of 2M pulsars, there were numerous works showing the possibility of the existence of hyperonic matter models that yield M max > 2M . As we showed in the present review, these dense matter models had a rather stiff pre-hyperon segment of the EOS, resulting in rather large radii of stars with M ∼ 1.4M , usually 7 http://www.gravity.phys.uwm.edu/rns/ . R 1.4 > 12 km. We argue that this EOS feature is amplified by rotation, and leads also to a rather high minimum mass of rotating NS.
A fully consistent calculation of the NS radius should in principle be performed for a unified EOS, where the crust and the outer layer of the core are described using the same nuclear model. Besides, "gluing" two different EOS based on different nuclear models is, to a large extent, arbitrary and unphysical.
The scenario for the formation of the crust (accretion in a close binary system or cooling down after the NS birth in supernova explosion) has a small effect on the M -R relation. The accreted crust is thicker by ∼ 100 m at 1.4M . Rotation at 716 Hz increases this difference by some 50%.
The presence of a quark core in NS is not excluded by observations of 2M pulsars. However, if a quark core in such hybrid stars contains a sizeable fraction of the star's mass, quark matter has to be stiff enough (sound speed > 0.6 c) and deconfinement should occur at not too low density (2.5n 0 -3.5n 0 , and the density jump at the core edge should be small). Finally, the last purely hadronic star should have a rather large radius > 12 km and a mass not much higher than 1.5M [104,135]. Response to rotation is strong, with M min (716 Hz) > 1M .
For a strong first-order phase transition to quark matter (large density jump at the quark core edge) we obtain a family of very massive 2M hybrid stars separated from less compact massive hadron stars; in this way one finds massive configurations of hadron and hybrid stars of the same M 2M but with a different structure and radius (hybrid and hadron twins; the hadron twin has a significantly larger radius > 14 km than the quark-core twin [109]). We may expect that the M -R curve of these disjoint families will be quite sensitive to rotation.
Generally, phase transitions in NS cores can have a strong effect on the M -R relation for NS. Our discussion in the case of the first-order phase transitions studied in the past was actually very general. As we stressed, one has to separate two aspects of the phase transition N −→ E. A crucial parameter for the M -R relation for hydrostatic equilibrium configurations is the energy density ratio at the phase coexistence point P = P NE , λ = ρ E /ρ N . Consider first non-rotating NS. For λ < λ crit = 3 2 · (1 + P NE /ρ E c 2 ) equilibrium configurations form a continuous family with M stat min < M < M stat max , with M stat min 0.1M , and M stat max > 2M . This continuous character is conserved for configurations rotating rigidly at f , with M rot max increasing a few percent at 716 Hz while M stat min (716 Hz) ∼ 0.7-1M . For λ > λ crit the family of stable configurations of NS splits into two disjoint families, separated by a segment of unstable configurations (instability with respect to spherically symmetric perturbations) that cannot exist in Nature. It has been checked that this topological feature is conserved for rotating configurations, with the instability being induced by axially symmetric perturbations.
Transitions between two stable segments of the M -R curve are associated with NS minicollapse, with energy release weakly depending on the angular momentum of the collapsing configuration, but strongly depending on the degree of metastability of the core undergoing the phase transition.
We reviewed various aspects of a minicollapse due to the formation of an exotic core. The most spectacular astrophysical signatures are associated with a minicollapse induced by a N −→ E conversion of detonation type. We described the history of the famous gamma-ray burst on March 5, 1979 and its initial explanation by a minicollapse in the core of NS in N49 supernova remnant in Large Magellanic Cloud. The source of this burst turned out to be a soft-gamma repeater, a magnetar emitting repetitively gamma and X-ray flares powered by the huge magnetic field. Finally, a positive message from recent numerical simulations of minicollapses of rotating NS was that they could produce a burst of gravitational waves detectable in the Galaxy and its close vicinity by the network of Advanced Virgo and Advanced LIGO detectors and in the future by the planned Einstein Telescope underground interferometric detector.