{'fraction_non_alphanumeric': 0.03535081322722805, 'fraction_numerical': 0.01984718303160965, 'mean_word_length': 4.290394722129039, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recent Progress in Neutron Star Theory

This review contains chapters discussing: Energy density fluctionals of nuclear matter, Many-body theory of nucleon matter, Hadronic and quark matter, Mixtures of phases in dense matter, Neutron star observations and predictions.


Introduction
Neutron stars are among the most fascinating bodies in our universe. They contain over a solar mass of matter within a radius of ∼ 10 km at densities of order 10 15 g/cc. They probe the properties of cold matter at extremely high densities, and have proven to be fantastic test bodies for theories of general relativity. In a broader perspective, neutron stars and heavy ion collisions provide access to the phase diagram of matter at extreme densities and temperatures, that is basic for understanding the very early Universe and several other astrophysical phenomena.
The discovery of the neutron by Chadwick in 1932 prompted Landau [1] to predict the existence of neutron stars. The first theoretical calculations of neutron stars were performed by Oppenheimer and Volkoff [2] in 1939 assuming that they are gravitationally bound states of neutron Fermi gas. The calculated stars had a maximum mass of ∼ 0.7 M ⊙ , central densities up to ∼ 6 × 10 15 g/cm 3 and radii ∼ 10 km. For comparison the density of nuclear matter inside a large nucleus like 208 Pb is ∼ 0.16 nucleons/fm 3 , i.e. ≃ 2.7 ×10 14 g/cm 3 [3]. Their predicted maximum mass was less than the Chandrasekhar mass limit of ∼ 1.4 M ⊙ for white dwarfs made up of iron group nuclei, and having densities up to ∼ 10 9 g/cm 3 [4]. The pressure to balance the gravitational attraction in white dwarfs and Oppenheimer-Volkoff neutron stars is supplied by degenerate electron and neutron Fermi gases respectively.
In 1934 Baade and Zwicky [5] suggested that neutron stars may be formed in supernovae in which the iron core of a massive star exceeds the Chandrasekhar limit and collapses. The large amount of energy released in the collapse blows away the rest of the star and the collapsed core may form a neutron star. For efficient production of neutron stars with this mechanism, the maximum mass of neutron stars should exceed 1.4 M ⊙ . In the 60's, using schematic models of nuclear forces, Tsuruta and Cameron [6] showed that they could increase the neutron star masses beyond 1.4 M ⊙ .
Bell and Hewish discovered radio pulsars in 1967, and they were soon identified as rotating neutron stars by Gold [7]. The subsequent detection of the Crab pulsar in the remnant of the Crab supernova, observed in China in 1054 A.D., confirmed the link to supernovae, and initiated the present efforts to better understand neutron stars.

A Brief Overview of Observations
Almost 1200 pulsars have been discovered by the turn of this millennium. In these stars the magnetic and rotational axes are misaligned, thus they emit dipole radiation in the form of radio waves that appear to pulse on and off like a lighthouse beacon as the pulsar beam sweeps across the Earth. The rotational energy loss due to dipole radiation isĖ where the moment of inertia for a typical neutron star is I ∼ 10 45 g cm 2 . Pulsars have magnetic fields B of ∼ 10 12 G, deduced from the observedΩ, and independently confirmed by cyclotron absorption lines found in X-ray spectra. Their periods, P = 2π/Ω, ranging from 1.5 ms to 8.5 s, are increasing with derivatives P ∼ 10 −12 − 10 −21 . The pulsar age is approximately given by P/2Ṗ [8]; most pulsars are old and slowly rotating with relatively small period derivatives, except for a few young pulsars, e.g., those found in the Crab and Vela nebulae. In 1969 the Crab and the Vela pulsars were observed to "glitch", i.e. to suddenly speedup with period changes ∆P/P of the order of 10 −8 and 10 −6 respectively [9]. In post-glitch relaxation most of the period increase ∆P decays. These pulsars have glitched several times since then. The glitches suggest that the neutron stars have a solid crust containing superfluid neutrons. The interesting structure of their crust has been recently reviewed [10], and we discuss it rather briefly in this report.
The first binary of two pulsars was found by Hulse and Taylor in 1973 and they could determine many of its parameters including both masses, orbital period and period derivative, orbital distance and inclination. General relativity could be tested to an unprecedented accuracy by measuring the inward spiralling of the neutron stars in the Hulse-Taylor binary PSR 1913+16 [11]. The periastron advance in PSR 1913+16 is 4.2 • per year as compared to 43" per century for Mercury, which originally was used by Einstein to test his theory of general relativity. Six double neutron star binaries are known so far, and neutron stars in all of them have masses in the range 1.36 ± 0.08M ⊙ [12]. They confirm that nuclear forces have a large effect on the structure of neutron stars and increase their maximum mass beyond 1.4 M ⊙ . Neutron stars are estimated to have a binding energy of ∼ 10% of their mass. Thus ∼ 1.5 M ⊙ of nuclei are needed to obtain a 1.35 M ⊙ star.
A distinct subclass of radio pulsars are millisecond pulsars with periods < ∼ 100 ms. The fastest pulsar known has a period of 1.56 ms [13]. The period derivatives of millisecond pulsars are very small corresponding to low magnetic fields ∼ 10 8 − 10 10 G. They are believed to be recycled pulsars, i.e. old pulsars that have been spun up by mass accretion whereby the magnetic fields have decayed. About 80% of the millisecond pulsars are in binaries whereas less than 1% of normal radio pulsars are in binaries. About 20 -almost half of the millisecond pulsarsare found in binaries where the companion is either a white dwarf or a neutron star.
With X-ray detectors on board satellites since the early 1970's about two hundred X-ray pulsars and bursters have been found of which the rotational period has been determined for about sixty. The X-ray pulsars and bursters [14] are believed to be neutron stars accreting matter from high (M > ∼ 10M ⊙ ) and low mass (M < ∼ 1.2M ⊙ ) companions respectively. The X-ray pulses are attributed to strong accretion on the magnetic poles emitting X-rays (as northern lights). The observed radiation is pulsed with the rotational frequency of the accreting star. X-ray bursts are thermonuclear explosions of accreted matter on the surface of neutron stars. After accumulating hydrogen on the surface for hours, pressure and temperature become sufficient to trigger a runaway thermonuclear explosion seen as an X-ray burst that lasts a few seconds [14]. Masses of these stars are less accurately measured than for binary pulsars. We mention recent mass determinations for the X-ray pulsar Vela X-1: M = 1.87 +0. 23 −0.17 M ⊙ [15], and the burster Cygnus X-2: M = 1.8 ± 0.4)M ⊙ [16]. They are larger than the typical 1.36 ± 0.08M ⊙ masses found in pulsars binaries, presumably due to accreted matter.
A subclass of half a dozen anomalous X-ray pulsars has been discovered. They are slowly rotating, P ∼ 10 sec, but rapidly slowing down. This requires huge magnetic fields of B ∼ 10 14 G and they have appropriately been named "magnetars" [17]. Four gamma ray repeaters discovered so far are also believed to be slowly rotating neutron stars. The magnetars and likely also the gamma ray repeaters reside inside supernova remnants.
Recently, quasi-periodic oscillations (QPO) have been found in 12 binaries of neutron stars with low mass companions. If the QPO originate from the innermost stable orbit [18,19] of the accreting matter, their observed values imply that the accreting neutron stars have masses up to ≃ 2.3M ⊙ . In this case the QPO's also constrain the radii of the accreting star.
Non-rotating and non-accreting neutron stars are virtually undetectable but the Hubble space telescope has observed one thermally radiating neutron star [20]. Its surface temperature is T ≃ 6 × 10 5 K≃ 50 eV and its distance is less than 120 pc from Earth. Circumstantial evidence indicate a distance of ∼ 80 pc which leads to a radius of 12-13 km for this star. In recent years much effort has been devoted to measuring pulsar temperatures, especially with the Einstein Observatory and ROSAT. Surface temperatures of a few pulsars have been measured, and upper limits have been set for many [21].
¿From the human point of view supernova explosions are rare in our and neighboring galaxies. The predicted rate is 1-3 per century in our galaxy and the most recent one was 1987A in LMC. No neutron star associated with this explosion has been detected; however, 19 neutrinos were detected on earth from 1987A [22], indicating the formation of a "proto-neutron star". It has been suggested by Bethe and Brown [23] that an upper limit to the mass of neutron stars can be obtained assuming that the remnant of SN 1987A collapsed into a black hole.
Astrophysicists expect a large abundance of ∼ 10 8 neutron stars in our galaxy. At least as many supernova explosions, responsible for all heavier elements present in our Universe today, have occurred. The scarcity of neutron stars in the solar neighborhood may be due to production of black holes or other remnants in supernovae, or due to a high initial velocity (asymmetric "kick") received during their birth in supernovae. Recently, many neutron stars have been found far away from their supernova remnants; and of the ∼ 1200 discovered radio pulsars only about ∼ 10 can be associated with the 220 known supernova remnants. Neutron stars thrown out of the galactic plane may be detected by gravitational microlensing experiments [24] designed to search for dark massive objects in the galactic halo.
The recent discovery of afterglow in Gamma Ray Bursters (GRB) allows determination of their very high redshifts (z ≥ 1). They imply that GRB occur at enormous distances. Evidence for beaming has been observed [25], and the estimated energy output is ∼ 10 53 ergs. Such enormous energies can be produced in neutron star mergers eventually forming black holes. From abundance of binary pulsars one can estimate the rate of neutron star mergers; it is compatible with the rate of GRB of approximately one per day. Another possible mechanism, is a special class of type Ic supernova (hypernovae) where cores collapse to black holes [26].
The future of neutron star observations looks bright as new windows are about to open. A new fleet of X-and Gamma-ray satellites have and will be launched. With upgraded ground based observatories and detectors for neutrinos and gravitational waves [27] our knowledge of neutron star properties will be greatly improved.

Theory of Neutron Star Matter
Neutron stars are made up of relatively cold, charge neutral matter with densities up to ∼ 7 times the equilibrium density ρ 0 = 0.16 nucleons/fm 3 of charged nuclear matter in nuclei. The matter density is > ρ 0 over most of the star, apart from the relatively thin crust [10]. The Fermi energy of neutron star matter is in excess of tens of MeV, and hence, at typical temperatures of < ∼ KeV, thermal effects are a minor perturbation on the gross structure of the star.
Matter at such densities has not yet been produced in the laboratory, its properties must be theoretically deduced from the available terrestrial data with guidance from observed neutron star properties. The quantities of interest are the phase and composition of cold catalyzed neutral dense matter, its energy density ǫ(ρ) and pressure P (ρ), where ρ denotes the baryon number density. The baryon number is conserved in all known interactions, therefore it is convenient to find the composition by minimizing the total energy E T (ρ) per baryon, including rest mass contributions. This gives: The equation of state (EOS) P (ǫ) is found by eliminating ρ from the above two.
The gravitational equilibrium of a nonrotating star is described by the Tolman-Oppenheimer-Volkoff (TOV) [4] Eq: where G is the gravitational constant, P (r) and ǫ(r) are the pressure and mass density at radius r in the star, and is the mass inside r. If we neglect the general relativistic corrections of order 1/c 2 the TOV Eq. reduces to the Newtonian hydrodynamic equation. The TOV Eq. can be easily integrated starting from the central density ǫ c ar r = 0 to find the density profile ǫ(r). At the radius R of the star P (R) = 0, and m(R) = M is the mass of the star as seen from outside. The stability of the star can be deduced from the M (ǫ c ) as discussed in [4], and the equations for rotating stars are given by [28]. The effect of rotation on the structure of most observed neutron stars seems to be rather small, however, it could be significant at periods less than a millisecond [29]. At densities < 2 × 10 −3 ρ 0 matter is believed to have the form of a lattice of nuclei in a relativistic degenerate electron gas [10], qualitatively similar to that of metals. The main focus of the theory reviewed here has been on determining the properties and EOS of matter in the density range 2 × 10 −3 ρ 0 < ρ < 10ρ 0 from terrestrial data. In the lower part of this range we expect to find nucleon matter (NM) composed of nucleons and electrons. In contrast to matter in nuclei, it has mostly neutrons with a small fraction of protons and equal number of electrons to maintain charge neutrality. The large Fermi energy, µ e ∼ 100 MeV, of the electron gas limits the fraction of protons in NM.
At higher densities there are several possibilities including condensation of negatively charged pions and kaons, occurrence of hyperons, and the transition from hadronic to quark matter. All these possibilities exploit the large electron Fermi energy of NM, therefore only one of these, if any, may occur and lower the µ e . In addition, neutron star matter can have interesting mixed phase regions in which the mixing phases are charged but the matter is overall neutral [10].
We begin with a review of NM, and later consider the more exotic possibilities. In the last sections the range of neutron star structures predicted by theory is presented along with a comparison with the observational data.

Energy-Density Functionals of Nucleon Matter
The simplest description of nuclei is obtained within the mean field approximation. It assumes, following the nuclear shell model, that nucleons occupy single particle orbitals in an average potential well produced by nuclear forces. The energy of the nucleus is assumed to be a functional of the orbitals occupied by the nucleons, and the orbitals are determined variationally as in the Hartree-Fock approximation. In reality the mean field approximation is not exactly valid for nuclei. The observed proton knockout reaction rates [30] indicate that the shell model orbitals are occupied with a probability of ∼ 70% in the simplest closed shell nuclei like 208 Pb. The differences between the real and the mean field wave-functions, due to correlations induced by nuclear forces, are subsumed in the energy functional as suggested by Kohn [31] in the context of atomic and molecular physics.
The energy density of hypothetical, uniform NM at zero temperature is the main term in the energy functionals. The nucleon orbitals in uniform matter are simple plane waves, and the ground state in mean field approximation is obtained by filling the proton and neutron states up to their Fermi momenta k F,N = (3π 2 ρ N ) 1/3 , where N = n, p for neutrons and protons. The energy density, denoted by E(ρ n , ρ p ), includes kinetic and strong interaction contributions, but excludes rest masses and the Coulomb interaction, which destabilizes uniform charged matter. The total density is denoted by ρ = ρ n + ρ p , the asymmetry of the matter is defined as β = (ρ n − ρ p )/ρ, and the energy per nucleon, E(ρ, β), is given by E(ρ, β)/ρ.
Analysis of nuclear properties with the liquid drop model [3] reveals that, in the absence of electromagnetic forces, the ground state of NM is symmetric (i.e. β = 0), has total equilibrium density ρ 0 = 0.16 ± 0.01 fm −3 , and binding energy E 0 = −16±0.5 MeV per nucleon. The symmetry energy E sym (ρ 0 ) = 34±6 MeV, is defined as 1 2 ∂ 2 E/∂β 2 at equilibrium. The NM energy E(ρ, β) can be expanded about its minimum value at β = 0 in powers of β 2 , assuming charge symmetry of nuclear forces. In variational [32] as well as Brueckner [33] theories the coefficients of terms with β n≥4 are estimated to be small, and E(ρ, β) ≈ (1 − β 2 )E(ρ, 0) + β 2 E(ρ, 1). In this approximation the symmetry energy is the difference between the energy of pure neutron matter and symmetric nuclear matter. The incompressibility K 0 = 240 ± 30 MeV [34] of symmetric nuclear matter is defined as K 0 = k 2 F ∂ 2 E/∂k 2 F at equilibrium. The energies of the collective breathing mode vibrations of nuclei are sensitive to K 0 ; however, in all stable nuclei the surface effects are significant. It is difficult to extract the density and β-dependence of the incompressibility, and the density-dependence of the symmetry energy from available nuclear data.
Analysis of elastic scattering of nucleons by nuclei shows that the nuclear mean field has a dependence on the energy of the moving nucleon [3]. Over a wide energy range this dependence is approximately linear, suggesting that nucleons in equilibrium nuclear matter have an effective mass m ⋆ ∼ 0.7m, where m is the free nucleon mass. This effective mass should not be identified with the Landau effective mass which describes the density of single particle states in a narrow energy interval about the Fermi energy [35,36]. The Landau m ⋆ in uniform matter is difficult to extract from nuclear data, since nucleons at the Fermi energy are strongly coupled to nuclear surface dynamics. Some of the phenomenological energy functionals are chosen to fit the observed nuclear level densities at the Fermi energy, while others fit the value of m ⋆ (ρ 0 , 0) obtained from the energy dependence of the optical model potential [37]. The nonrelativistic functionals based on Skyrme effective interactions [37] generally contain the following terms: Here τ N = 0.6k 2 F,N ρ N /m are the kinetic energy densities, and τ = τ n + τ p . The parameters x 1 to x 4 and α describe the ρ and β dependence of the volume integral of the static part of the effective interaction between nucleons in matter, while the x 5 and x 6 describe effective masses produced by the momentum dependence of the effective interaction. In principle the values of the seven parameters in a typical Skyrme functional are constrained by the empirically known values of ρ 0 , E 0 , E sym (ρ 0 ) and K and the choice made for m ⋆ (ρ 0 , 0). However, since the constraints are insufficient, there are many Skyrme models of the energy functional.
The simple form of the functional (Eq.5) chosen by most Skyrme models is convenient, but the real functional can be much more complex. The analytic form of the energy density predicted by realistic models of nuclear forces, as discussed in the next chapter, has been studied by Ravenhall [38]. A much more elaborate function of the type: is required to reproduce the predicted E(ρ n , ρ p ) up to β = 1 and ρ ∼ 1 fm −3 . This functional also explains the nuclear binding energies and the empirically known values for symmetric nuclear matter [39], however, it is unlikely that the values of all of its thirteen parameters can be obtained by fitting nuclear data. The energy of NM can be easily calculated from a covariant effective Lagrangian in the mean field approximation, as shown by Walecka [40], and in the past decade many properties of medium and heavy nuclei have been studied with this approach [41]. The effective Lagrangian used in the recent work has the form: Hear ψ, ω µ and ρ µ are respectively the nucleon and ω and ρ vector-meson fields. Overhead arrows are used to denote isospin vectors, and Ω µν = ∂ µ ω ν − ∂ ν ω µ etc. The effective scalar field σ is responsible for nuclear binding, and the σ 3 and σ 4 terms are necessary to obtain the empirical incompressibility of nuclear matter [42]. The isovector ρ field is required to obtain the empirical symmetry energy.
The observed values of the masses m, m ω and m ρ are used, and the coupling constants g ω , g ρ , g σ , g 2 and g 3 , as well as the mass m σ of the effective scalar field are adjusted to fit the nuclear data. The above Lagrangian, without the σ 3 and σ 4 terms but including pion fields and their coupling to the nucleon, is also used to model the two-nucleon interaction discussed in the next section. The relativistic mean field theory of nuclei is very elegant and often used to study properties of neutron star matter [43]. It has provided important insights into relativistic effects in nuclei and NM. However, the effective mean-field Lagrangian (Eq.7) is unlikely to have a simple physical meaning. The inverse masses of the vector and scalar fields correspond to lengths of ∼ 0.25 and 0.4 fm, which are much smaller than the unit radius r 0 = (4πρ/3) −1/3 ∼ 1.2 fm for equilibrium nuclear matter. The naive condition for the validity of the mean field approximation, that r 0 be much less than the inverse masses of the fields is totally violated in nuclei as well as in neutron stars. Pions are omitted from the effective Lagrangian because they do not contribute to the energy of matter in the mean field approximation. Their higher order contributions are subsumed in the effective scalar field. Therefore, the effective mean-field Lagrangian must be interpreted as a relativistic generalization of Kohn's energy functional. Eq. (7) assumes the simplest form necessary to fit nuclear data. A more general form, necessary to explain the properties of NM over the wide density-asymmetry range in neutron stars, can have additional fields, isovector scalar for example, density dependent coupling constants to take into account the changes in correlations with density, and field energies containing high powers of the fields, etc.
The energy of low density neutron matter is well determined by realistic models of two-nucleon interaction obtained by fitting the nucleon-nucleon (NN) scattering data. Different models and methods of calculation give very similar results up to ρ ∼ ρ 0 , beyond which three-nucleon interactions and relativistic effects, as well as computational difficulties may become appreciable. These energies thus provide a test of the ability of the Skyrme and relativistic mean field theories to find neutron matter properties by extrapolating data on nuclear binding energies, sizes, vibrations, etc. As shown in Fig.1, the neutron matter energies predicted by the various functionals are widely different, and not in agreement with the results of many-body calculations at ρ < ρ 0 . It thus appears likely that the simple forms of effective interactions or Lagrangians used in the present mean field theories are inadequate to predict the properties of neutron star matter by extrapolating the observed nuclear properties. Nevertheless, effective mean-field Lagrangians have been widely used in neutron star studies due to their simplicity [43,51].

Many-Body Theory of Nucleon Matter
Many properties of nuclei and nuclear matter can be understood from the Hamiltonian: which includes the kinetic energy, two-nucleon interactions denoted by lower case v and three-nucleon interactions by capital V . The ellipsis denote neglected four and higher body interactions. In this section we review the present status of this approach and its limitations.

Models of Two Nucleon Interaction
Our understanding of QCD has not yet progressed enough to predict the twonucleon interaction v ij ab initio. The long range part of v ij is known to be mediated by pions, the lightest of all the mesons, and it is denoted by the one pion exchange potential (OPEP) given by: The pion-nucleon coupling constant f 2 πN N /4π = 0.075 ± 0.002 [52], and the radial functions associated with the spin-spin and the tensor parts are: where x = m π r. The tensor operator, S ij = 3 σ i ·r ij σ j ·r ij − σ i · σ j , and since T π (r) ≫ Y π (r) in the important x < ∼ 1 region, the OPEP is dominated by its tensor part. The complete v ij is expressed as v π ij + v R ij , where v R ij contains all the other, heavy meson, multiple meson and quark exchange parts. It has to be obtained along with the short range cutoffs ξ Y (r) and ξ T (r) in the OPEP by fitting NN scattering data. In boson-exchange models the v R is approximated by a sum of attractive scalar meson, and repulsive vector meson exchange potentials, while other models use attractive two-pion exchange and repulsive core potentials. The OPEP contains a δ−function term omitted from Eq. (10). Due to the finite size of nucleons this term acquires a finite range, and is difficult to separate from v R .
In the early 1990's the Nijmegen group [48] carefully examined all the data on elastic NN scattering, at energies below the pion production threshold of ∼ 350 MeV, published between 1955 and 1992. They extracted 1787 protonproton and 2514 proton-neutron "reliable" data, and showed that these could determine all NN scattering phase shifts and mixing parameters quite accurately, thus claiming that the experimental information on elastic NN scattering is now complete. Additional measurements are being carried out at several laboratories including the Indiana University Cyclotron Facility, and the CELCIUS facility in Uppsala, Sweden, to test the accuracy of the claim and improve on the quality of the data base. Nevertheless the Nijmegen analysis has been a major step.
NN interaction models which fit the Nijmegen data base with a χ 2 /N data ∼ 1 are called "modern". These include the Nijmegen models [49] called Nijmegen I, II and Reid 93, the Argonne v 18 [53] (A18) and CD-Bonn [54]. In order to fit both the proton-proton and neutron-proton scattering data simultaneously and accurately, these models include a detailed description of the electromagnetic interactions and terms that violate the isospin symmetry of the strong interaction via the differences in the masses of the charged and neutral pions, etc.
These five models use different parameterizations of the v R ij , and the Nijmegen-I and CD-Bonn also include nonlocalities suggested by boson-exchange representations. Thus, like the older models, they make different predictions for the manybody systems. However, the differences in their predictions are much smaller than those between older models, and can be partly understood.
The interaction in the spin-isospin T, S = 0, 1 has the largest model dependence, which gets carried over to the estimates of the energy of symmetric nuclear matter (sec. 3.4). Fortunately the deuteron structure provides significant information on this interaction. Fig.2 shows the deuteron wave functions obtained with the five modern potentials [55]. The three potentials, Reid 93, Nijmegen II and A18, which are local in each NN partial wave, give essentially the same deuteron wave function. We expect that they will give rather similar matter properties. The Nijmegen I and CD-Bonn potentials have momentum dependent terms associated with heavy meson exchange. These two potentials give larger 3 S 1 wave functions at r < 0.8 fm, because they have softer repulsive cores than the local models; however, this effect is not very large as can be seen from Fig.2. At r > 0.8 fm only the CD-Bonn predictions differ from the rest. This is because CD-Bonn has a strongly nonlocal OPEP as suggested by pseudoscalar pion-nucleon coupling, which suppresses the 3 D 1 wave function.
These results indicate that the main difference between the preset models of v ij is from the assumed nonlocality in OPEP. However, relativistic field theories permit use of OPEP with different nonlocalities [56] related with the Dyson transformation. The three-nucleon interaction, V ijk , depends upon the choice of OPEP [57], and the final results obtained after including it should be independent of the choice. Therefore, if relativistic field theories can be used to describe pion exchange forces, one can use either the local OPEP in A18, or the nonlocal one in CD-Bonn. In this case it may be better to use the local representation because it is simpler, and more accurate many-body calculations can be carried out with it.
Realistic models of nuclear forces, the modern as well as the older, predict the existence of a dense toroidal inner core in the deuteron [58]. The density distribution, ρ(r ′ ), of the deuteron in spin projection M = 0 state, in the center of mass frame is shown in Fig.3. Here r ′ = ±r/2 are the nucleon coordinates in the center of mass frame, and r is the internucleon distance. This density distribution is symmetric under rotations about the Z ′ -axis, and the top part shows its cross section in the X ′ − Z ′ plane as predicted by the A18 model. The bottom part shows the toroidal shape of the equi-density surface for half maximum density. The density peaks on a ring of diameter of ∼ 1.0 fm inside this torus having thickness of ∼ 0.8 fm. The shape is produced by constructive (destructive) interference between the S-and the D-wave functions shown in Fig.2 along the X ′ (Z ′ ) axis of the deuteron in the M = 0 state.
The dominant, static part of the v ij in the deuteron, obtained by omitting the terms dependent on the angular momentum, is anisotropic due to the tensor part of the OPEP. It strongly depends upon the angle between the unit vector r ij and the spin directions. The expectation value of the interaction in the S = 1, M S = 0 state, 1 √ 2 | ↑↓ + ↓↑ , is shown in Fig.4 as a function of r ij for θ = 0 and θ = π/2. The interaction is attractive for θ = π/2 and repulsive for θ = 0, like that between two magnetic dipoles. The OPEP is also shown in Fig.4 by dashed lines. The NN interaction in all states except those with T, S = 1, 0 is dominated by the OPEP at r ij > 1 fm. The OPEP is weakest in the T, S = 1, 0 states. At small r ij the repulsive core of v R ij dominates in all states. Most of the M = 0 deuteron wave function has also M S = 0. Thus it is possible to understand the density distribution of the deuteron from the potential shown in Fig.4. The two peaks in the density shown in Fig.3 correspond to the twonucleons in the deuteron being ∼ 1 fm apart at θ ∼ π/2 where the potential has its minimum value of ∼ −200 MeV in the A18 model. Other models also have a deep minimum at this position. The smallness of the density along the Z-axis is due to the repulsive potential at θ = 0. Even though these features existed in the potential models of the sixties their experimental confirmation came via a series of measurements of the electromagnetic form factors of the deuteron up to momentum transfers of ∼ 8 fm −1 conducted since the mid eighties at SLAC and Bates. These and the more recent, high precision measurements [59,60,61] carried out at the Jefferson Lab are in good agreement with the predictions of the A18 model, and verify the predicted deuteron structure beyond r ∼ 0.7, or equivalently r ′ ≥ 0.35 fm.
The T, S = 1, 0 two-proton distribution functions, believed to be similar to the two-neutron distribution functions due to isospin symmetry, are predicted [58] to have a dip at r ∼ 0 and a peak at r ∼ 1 fm, due to the repulsive core and the minima of v N N respectively. The experimental information on these spherically symmetric distribution functions is less direct. It comes from sums of longitudinal response functions of light nuclei [62]. The sums observed in 3 He and 4 He are in fair agreement with theory, and show evidence of the predicted structure, however, the relativistic and other corrections to the observed sums are significant.
The observed deuteron form factors and, to a lesser extent, the sums of longitudinal response of light nuclei indicate that the modern two-nucleon potentials and the wave functions they predict have validity at internucleon distances larger than ∼ 0.7 fm. This may appear surprising because the rms charge radius of the proton is known to be 0.8 fm. However, the nucleons seem to have a small and dense core. The charge form factor of the proton, as well as the magnetic form factors of the proton and neutron are well approximated by the dipole (1 + q 2 /q 2 0 ) −2 with q 0 = 840 MeV/c. Inverting this form factor gives the proton charge density as ρ p ch (r) = 3.3 e −r/0.23f m fm −3 . The charge densities of two protons, one fm apart are shown in Fig.5. It should be noted that the charge densities shown can have corrections due to the neglect of relativistic effects in inverting the form factor at r < ∼ 1/m ∼ 0.21 fm, and that recent measurements of the proton form factor [63] show deviations from the dipole form at momenta > ∼ 4 fm −1 . They suggest that proton charge density flattens out at r < 0.3 fm. Nevertheless Fig.  5 indicates that nucleons a fm apart can retain their identities. The NN interaction includes the change in the energy due to their overlap, and it has minima near r ∼ 1 fm. In absence of the quantum kinetic energy term (−∇ 2 /2m) in the Hamiltonian (Eq.8) the deuteron will shrink to a ring of radius ∼ 0.5 fm, and the equilibrium density of nuclear matter will be ∼ 1 fm −3 . The density of matter in most neutron stars is less than that.

Models of Three Nucleon Interaction
All realistic models of v ij , the modern and the older, underbind the triton and other light nuclei and predict too high equilibrium density for symmetric nuclear matter. In both cases the deviation from experiment is not too large, particularly when compared with the expectation values of v ij . For example, the expectation value v i<j in 3 H is about ∼ −50 MeV, while the underbinding is by < 1 MeV. It is likely that these differences are due to three-nucleon interactions expected and predicted since the fifties [64].
High precision modeling of the two-nucleon interaction is possible because scattering cross sections can be easily and exactly calculated from the assumed v ij , and a complete set of ∼ 4000 cross sections has been measured. Such an approach is not practical at present for the three-nucleon interaction V ijk . In principle deuteron-nucleon scattering and reactions can be used to study V ijk in the isospin T = 1/2 state. However, this scattering is dominated by the two-nucleon interaction [65], and very high precision data is necessary to extract the effects of V ijk , as is being attempted in new experiments at the Indiana University cyclotron, focused on spin observables.
The large v i<j is cancelled to a large extent by the kinetic energy in nuclear binding energies. Thus V i<j<k in nuclei is expected to be of the order of 10 % of their binding energy. Hence we can construct realistic models of V ijk by fitting binding energies of light nuclei, which can now be calculated with an accuracy of the order of 1 % using Greens Function Monte Carlo (GFMC) methods [66,67] and the estimated equilibrium properties of nuclear matter. As will be discussed in sec. 3.7, it is not yet possible to calculate the EOS of symmetric nuclear matter with comparable accuracy. The models depend upon the v ij used in the Hamiltonian. This is inevitable, because unitary transformations make correlated changes in v ij and V ijk [57]. Only the combinations of v ij and V ijk in the Hamiltonian (8) are meaningful.
The information contained in nuclear binding energies and equilibrium properties of nuclear matter is limited. Therefore the realistic models of V ijk rely on theory to a much larger extent than the models of v ij , and contain very few parameters. The Urbana models of V ijk contain two isoscalar terms: The first term represents the Fujita-Miyazawa [64] two-pion exchange interaction: with strength denoted by A 2π . The functions T π (r ij ) and Y π (r ij ) in X ij , Eq.(10), are taken from the A18 model of v ij . This interaction is due to the pion exchanged by nucleons j and k being scattered by the nucleon i via the ∆ resonance in π-N scattering. In classical terms it is due to the polarization of the quark spins in nucleon i, due to the pion field of j (k), interacting with k (j); and it is similar to the three-body earth-moon-satellite gravitational interaction due to the polarization of the ocean water on earth by the moon's gravity.
The V R ijk is purely phenomenological, and has the form:  It was meant to represent the modification of the two-pion exchange part of v ij by other particles in matter, however, ∼ 40 % of it is due to relativistic effects discussed in the next subsection. The parameters of the present model U-IX, A 2π = −0.0293 MeV and U 0 = 0.0048 MeV, have been determined from exact GFMC calculations of 3 H and approximate variational calculations of the equilibrium density of nuclear matter with the A18 NN interaction [68].
The results of essentially exact GFMC calculations [66,67] with the A18 and U-IX interactions are shown in Table-1. We note that the better known pion exchange parts of these interactions give the largest contributions, but the contributions of the phenomenological parts, v R ij and V R ijk , are significant. The column ∆E expt gives the difference between the experimental energies and the calculated, while ∆E V M C is that between the variational Monte Carlo (VMC) upper bounds and the exact GFMC energies.
Table-1 shows that the U-IX interaction underbinds A = 8 nuclei, and since 8 He is more underbound than 8 Be, it misrepresents the isospin dependence of V ijk . The new Illinois models of V ijk [69] resolve this problem by including the leading three-pion exchange term, V 3π ijk , that is attractive in triplets having isospin T = 3/2, but has little effect on the T = 1/2 triplets in 3 H and 4 He. A much improved fit, with errors < 2%, to the observed energies is obtained as shown in Fig.6. The three parameters in Pieper's model IL-2R are the strengths of the V 2π ijk , V R ijk and V 3π ijk . Calculations of NM properties with this more accurate V ijk are in progress; the preliminary results are similar to those with U-IX since the V 3π is much weaker than the V 2π and V R . For example, the expectation values of V 2π ijk , , +19 (+14) and −2 (−5) MeV. In the following sections we will review the properties of neutron star matter calculated with A18 and U-IX interactions.
In principle there can be four-nucleon interactions (FNI) neglected in the Hamiltonian (Eq.8). It seems that they are very weak in nuclei. All the models of V ijk studied so far reproduce the energy of 4 He with an error < 0.5%, after fitting the observed energy of 3 H. Since this error is close to the accuracy of the 4 He calculation, there is no indication of FNI in that nucleus. The IL-2R model also gives the experimental energies of A = 8 nuclei within ∼ 1 MeV. In 8 Be, for example, the expectation values of v ij and V ijk are respectively −308 and −21 MeV respectively. By comparison with experiment we estimate that the possi-ble contribution of FNI in this nucleus is < 1 MeV. The V ijk presumably has additional smaller terms neglected in IL-2R, but it is difficult to determine their strengths from the nuclear spectra that can be calculated accurately from bare forces at present.

Relativistic Boost Interaction
In all models, the NN scattering data is reduced to the center of mass frame and fitted using phase shifts calculated from the NN interaction, v ij , in that frame. The v ij obtained by this procedure describes the interaction between nucleons having total momentum P ij = p i + p j = 0. In general, the interaction between particles depends on their momenta p i and p j . For example the electromagnetic interaction [70] between two particles of mass m and charge Q, contains a term with the factor −p i ·p j /2m is the relative momentum. The terms containing p ij are included in the momentumdependent parts of v ij , while the boost interaction δv(P ij ) contains parts dependent on the total P ij . Even though contributions of the boost interaction to the binding energy of SNM and 3 H were estimated by Coester and coworkers in the seventies and eighties [71,72], they were first included in studies of dense matter rather recently [73]. Walecka's relativistic mean field theory naturally contains the boost interactions [74].
Following the work of Krajcik and Foldy, Friar [75] obtained the following equation relating the boost interaction of order P 2 to the interaction in the center of mass frame: The general validity of this equation in relativistic mechanics and field theory was recently demonstrated [74]. Including boost interaction, the nonrelativistic Hamiltonian assumes the form: where the ellipsis denotes the three-body boost, and four and higher body interactions. This H * contains all terms quadratic in the particle velocities, and is therefore suitable for complete studies in the nonrelativistic limit. Studies of light nuclei using the VMC method [76,77] find that the contribution of the two-body boost interaction to the energy is repulsive, with a magnitude which is ∼ 37% of that of V R ijk in the UIX model. The boost interaction thus accounts for a significant part of the V R ijk in Hamiltonians which fit nuclear energies neglecting δv. The V * ijk in Eq.(16) has a V R of appropriately smaller strength than that in the V ijk in Hamiltonian H given by Eq. (8).
We should expect additional relativistic corrections to the Hamiltonian (16). However, when nonrelativistic potentials are fit to the experimental data, relativistic effects present in the data are automatically buried in these potentials. In order to study the magnitude of a chosen relativistic correction, such as that due to the approximation of the kinetic energy or the boost interaction, or the nonlocalities of OPEP, it is necessary to refit the same data set, fitted to obtain the nonrelativistic Hamiltonian, and then study the differences. Such comparisons [77] indicate that the relativistic corrections associated with kinetic energies and nonlocalities of OPEP are small, whereas the boost corrections are significant. This is not surprising since the boost interaction was totally omitted from the conventional nonrelativistic nuclear Hamiltonian (Eq.8).

Brueckner Calculations of Nucleon Matter
Calculating the properties of matter from the interaction v ij between pairs of its constituents is a well known problem in many-body theory. It is particularly challenging for NM due to the strong spin-isospin dependence of v ij . In the method developed by Brueckner, Bethe and Goldstone the perturbation expansion of the energy of NM is cast into a series ordered according to the number of independent hole lines (HL). This method has been used extensively since the sixties [78] to study symmetric nuclear matter (SNM) and nuclei. It was also used then to predict properties of pure neutron matter (PNM) [79] soon after the discovery of pulsars. The convergence of the expansion depends upon the choice of the single-particle energies in the assumed unperturbed Hamiltonian. For the hole states with momenta less than k F , they are chosen self consistently as suggested by Brueckner and Gammel [80] via the Brueckner-Hartree-Fock (BHF) procedure. The older calculations by Day [81] used kinetic energies for particle states. Since this leads to a discontinuity in the unperturbed single-particle energies at k F , Day's choice is called "discontinuous". In 1976 the Liège group [36] advocated the "continuous" choice by extending the definition of BHF single-particle energies to particle states. If all the higher order terms of the HL expansion are computed the final results should be independent of the choice.
Detailed calculations of SNM have been carried out with the older Argonne v 14 interaction by the Catania group [82] using both choices. The results of the lowest order 2-HL calculation depend significantly on the choice, however, those including 2+3-HL terms are almost the same for the two choices. They find that the 2-HL results with the continuous (discontinuous) choice are ∼ 15 % below (30 % above) those of the 2+3. Day's calculations [81], with the same interaction, include additional 4-HL terms, and give energies below the 2+3 Catania results, closer to the 2-HL continuous. For example, the energies of SNM at ρ = 0.28 fm −3 are −11.3, −16.1, −18.3 and −17.8±1.3 for the Catania 2-HL discontinuous, 2+3-HL, 2-HL continuous, and Day's calculations respectively, while the variational upper bound at this density is −16.2 ± 0.4 [81] with the methods described in the next subsection. The size of the error due to truncation of the expansion is estimated in Day's and variational calculations.
A significant advantage of Brueckner's method is that it can be easily applied to local as well as nonlocal modern interactions. The present calculations, called lowest-order BHF (LOBHF), include only the 2-HL terms and use the continuous choice [83]. Their results for PNM and SNM are shown in Fig.7; those for matter with intermediate values of proton fraction can be estimated from these by interpolating with β 2 as metioned in sec. 2. These results provide an estimate of the uncertainty in the predicted matter energy due to that in the NN interaction.
The LOBHF energies for neutron matter are essentially model independent up to ∼ 0.3 fm −3 (Fig.7); at higher densities they deviate partly for the following reason. All the models fit the NN scattering data up to 350 MeV lab energy, i.e. up to maximum relative momenta k = 2.05 fm −1 . The maximum value of the relative momentum of two hole states in matter is k F , and the density of PNM (SNM) at k F = 2.05 is 0.29 (0.58) fm −3 . In PNM at ρ > 0.29 the interactions are being used at relative momenta larger than in the fitted data. Baldo et al [84] have shown that the 3 P 2 phase shifts predicted by the five modern potentials vary from 8 (A18) to 19 (Nijmegen II) degrees at k = 3 fm −1 . The average value, k rms = 3/10 k F , is smaller than k F , and exceeds 2.05 fm −1 at much larger densities of 1.77 (3.54) fm −3 of PNM (SNM). Presumably, this helps to keep down the model dependence.
The LOBHF energies of SNM have a larger model dependence starting at lower densities (Fig.7). Here the main cause seems to be the assumed nonlocality of v ij discussed in sec. 3 v i<j . TNI, added to the Hamiltonian to obtain the observed properties, naturally depend upon the choice of v ij . For example, those to be added to CD Bonn have to have stronger repulsive parts, which dominate at large ρ. The combinations of v ij +V ijk , constrained with experimental data, will have smaller model dependence than seen in Fig.7.
Baldo, Bombaci and Burgio [85] have carried out LOBHF calculations with the older Paris and Argonne interactions including Urbana TNI. Like the modern models, these older models give too large ρ 0 (∼ 0.26 fm −3 ) and E 0 (∼ −18 MeV) in LOBHF without TNI. By averaging over the position of the third nucleon the TNI is expressed as a density dependent NN interaction to be added to the v ij . Relativistic effects are included in the LOBHF calculations via Dirac-Brueckner approximation suggested by Celenza and Shakin [86]. The calculations include contributions of the boost interactions as well as TNI [87], and many nucleon interactions generated via the Z-graphs containing anti-nucleon lines. However, they do not contain contributions of the Fujita-Miyazawa and other TNI due to internal structure of nucleons. Results have been reported by Brockmann and Machleidt [88] for the older Bonn meson exchange NN interaction models; those with the Bonn-A model come close to reproducing the empirical properties of SNM at ρ 0 .

Variational Calculations of Nucleon Matter
Variational calculations of NM with realistic interactions have been carried out since 1970 [89]. The present calculations [68,73] use variational wave-functions, Ψ v , consisting of a symmetrized product of pair correlation operators, F ij , operating on the Fermi gas wave-function. In PNM, the F ij include four terms generating spatial, σ i · σ j , tensor and spin orbit correlations. The SNM F ij have eight terms; the additional four have τ i · τ j factors.
This wave-function is clearly too simple to accurately describe the ground state of nuclear matter. Monte Carlo studies of few-body nuclei use additional threebody correlations induced by both v ij and V ijk , in the variational wave-function; they reduce the energy of 16 O by ∼ 1 MeV/nucleon. The results shown in Table.1 indicate that the VMC energies of A = 8 nuclei, obtained after including three-body correlation operators, are above the exact GFMC values by ∼ 1 MeV/nucleon. From these we estimate that the present Ψ v may underbind SNM by a few MeV. In contrast, the three-body correlations have a smaller effect on the energy of pure neutron drops [90]. The variational energy of a drop with eight neutrons, calculated with the simple Ψ v , is greater than the exact value by ∼ 0.5 MeV/nucleon. Thus the variational energies are relatively more accurate for PNM than for SNM. This is as expected, since SNM has much stronger tensor correlations. Despite the aforementioned shortcomings, the simple Ψ v having only pair correlation operators describes the gross features of the nuclear wave-function rather well. For example, the spin-isospin dependent two-nucleon distribution functions calculated in this approximation are close to the exact distribution functions [66]. The correlation operators F ij are determined from Euler-Lagrange equations that minimize the two-body cluster contribution of an interactionv ij − λ ij . The interactionv ij is related to the v ij via a parameter α meant to simulate the quenching of the spin-isospin interaction between particles i and j, via their interaction with other particles in matter. The operator λ ij simulates screening effects in matter; it is determined from the ranges d t of tensor correlations, and d c of all the other correlations. The Ψ v thus depend on three variational parameters: α, d c and d t , determined by minimizing the energy. Two additional parameters are used in Ref. [91] to further lower the variational upper bound by small amounts.
The energy expectation value is evaluated using cluster expansion. The onebody term is just the Fermi gas kinetic energy, and the large two-body (2B) term, analogous to the interaction energy in LOBHF, is calculated exactly. The most important of the many-body (MB) cluster contributions are summed using Fermi-hypernetted chain (FHNC) and single operator chain equations [68,73], and constrains are imposed to satisfy the fundamental identities of pair distribution functions. The kinetic energy can be calculated using different expressions related by integration by parts. If all MB contributions are calculated, these expressions yield the same result. However, they yield different results when only selected parts of the MB clusters are summed using chain equations. Studies of atomic helium liquids with FHNC summation methods find the exact result to be between the energies obtained using the Jackson-Feenberg (JF) and Pandharipande-Bethe (PB) expressions. The average of these two is used as the result with half the difference as an estimate of the error.
More general pair correlations can be calculated by separately minimizing the two-body cluster contribution to each partial wave, specified by l, S, J and the relative momentum k [92]. These correlations f (l, S, J, k) depend on all the quantum numbers, and yield a lower 2B energy than the F ij operator with the same α, d c and d t . The MB contributions cannot be easily calculated with the general f (l, S, J, k), however. The differences between optimum f (l, S, J, k) and F ij can be included via the second order two-particle, two-hole contribution, ∆E 2 , in correlated basis perturbation theory [35,93]. In recent calculations [68,73] the ∆E 2 is approximated by the difference δE 2B between the 2B cluster energies calculated using f (l, S, J, k) and F ij . However, the values of α, d c , d t are determined by minimizing the energy calculated from the F ij .
The variational [73] and LOBHF energies obtained with the A18 interaction are compared in Fig.8. At densities below 0.6 fm −3 there is fairly close agreement between them, however, we expect the true results to be a few MeV below the variational upper bound. At higher densities the SNM LOBHF energy is significantly below the variational bound. The convergence of the HL expansion is expected to deteriorate at higher densities, and it may be the cause of the large difference.

Neutral Pion Condensation
Variational calculations with the Hamiltonian (16) with A18+δv+UIX * interactions, indicate the occurrence of a phase transition in both PNM and SNM; the energies of the two phases are shown in Fig.9. The tensor correlations have a longer range in the higher density phase (HDP) than in the low density phase (LDP). Detailed studies of the pion fields in the two phases [68] indicate that the HDP has a large enhancement of virtual neutral pions with momenta ∼ 1.5 fm −1 , and therefore this transition is believed to be due to neutral pion condensation in matter. Although the effect of this type of transition on the EOS is relatively small, it can have important consequences for the cooling and evolution of neutron stars [21].
Since the pioneering work of Migdal [94] and of Sawyer and Scallapino [95], many investigators have used effective interactions to study the possibility of pion condensation in SNM and PNM. These efforts were recently reviewed by Kunihiro et al [96]. Neutral pion condensation occurs when PNM (SNM) becomes unstable towards the development of a spin (spin-isospin) density wave as discussed in sec. 4.2. The variational wave functions used in the Urbana calculations [68] are not adequate to describe the long range order expected with π 0condensation; better calculations may be possible with the quantum Monte Carlo method described in the next subsection. NM with spin or spin-isospin density wave naturally has a pion field of the same wavelength; it causes N ⇀ ↽ ∆ transitions, which help lower the energy of this phase. The effects of the ∆-resonance are absorbed into the interactions in the Hamiltonian (16). In particular, the Fujita-Miyazawa V 2π ijk has a large effect on π 0 condensation. Without it there is no condensation predicted in SNM, while that in PNM occurs at a higher density of ∼ 0.5 fm −3 . This transition was found by Wiringa, Fiks and Fabrocini to occur with the older Argonne v 14 for PNM but not for SNM [91], while it does not occur in either PNM or SNM with the Urbana v 14 interaction of 1981. In Migdal's approach [94], the transition of SNM to the pion condensed phase is inhibited by a positive, short-range σ i · σ j τ i · τ j NN interaction whose strength is represented by the Landau parameter g ′ . In the case of PNM the sum of the σ i · σ j and σ i · σ j τ i · τ j interactions occurs since τ i · τ j = 1. The contact part of OPEP gives a negative σ i · σ j τ i · τ j potential at small r, and thus favors pion condensation. The σ i · σ j τ i · τ j part of the modern A18 NN interaction does become negative at small r, however, it is positive in the older Urbana-Argonne models which do not fit the NN scattering data as well as the A18.
The energies of β-stable LDP and HDP phases have been calculated by interpolation between PNM and SNM [73]. The LDP phase with ρ = 0.204 fm −3 and proton fraction x p = 0.073 is found to be in equilibrium with the HDP at ρ = 0.237 and x p = 0.057. In between there will be mixed phase regions as discussed in sec. 5.

Quantum Monte Carlo Calculations
The variational energy of the ground state of SNM, calculated with the A18+δv+UIX * model is ∼ −12 MeV, against the empirical value of −16 MeV. As mentioned earlier, this difference is believed to be mostly due to the inadequacy of the present variational wave functions. Accurate calculations of the energies of many-body systems in which the interactions do not depend upon the spins of the particles have been carried out using quantum Monte Carlo (QMC) methods [97]. In these systems one can work with a wave function Φ(R) that depends only upon the positions of all the particles represented by the configuration vector R = r 1 , r 2 , ...r A . A VMC calculation is used to obtain a good approximation Φ v (R), and in the GFMC method one operates on Φ v (R) with the imaginary time evolution operator exp(−[H − E 0 ]τ ) to project out the exact ground state Φ 0 (R).
The main difficulty in applying QMC methods to nuclear problems is that nuclear forces change spins and isospins of the interacting nucleons, and thus nuclear wave functions contain superpositions of all possible spin-isospin states of A-nucleons. Their number, , increases very rapidly with A. For this reason exact QMC calculations have been carried out only for nuclei having up to 8 nucleons [67,69], and attempts to calculate A = 9 nuclei are in progress. Carlson [98] has also calculated the ground states of 14 neutrons in a periodic box. Calculations of pure neutron systems are simpler because there is only one isospin state, and all τ i · τ j can be replaced by the unit operator.
In addition, GFMC calculations of Fermion systems suffer from the "Fermion sign problem". The real wave functions of simple Fermi systems have nodal surfaces because the Φ(R) must equal −Φ(R ′ ) when the configurations R and R ′ are related by the exchange of a pair of particles. GFMC configurations which diffuse across nodal surfaces, as the system evolves in imaginary time, increase the variance of the calculated quantity, making unconstrained propagation impractical for large systems. In the fixed node method [99] for simple systems this growth of variance is eliminated by restricting the configurations to domains enclosed by the nodal surfaces of the variational wave function. Such calculations generally have an error due to imperfections in that structure, but it is much smaller than that in variational.
A similar problem comes in nuclear GFMC with the additional complexity due to nuclear wave functions having many spin-isospin components, each with a different nodal structure. The growth of the variance is tolerable for A ≤ 7, but when A ≥ 8 it is necessary to use constrained path methods [67] to control the variance. The constraint can be removed at large τ to test if it influenced the calculated energy significantly. It appears that calculations with ∼ 2% accuracy in the binding energy are possible in this way for systems having up to 14 neutrons.
Auxiliary field diffusion Monte Carlo (AFDMC) [100] seems to be the longsought breakthrough needed to eliminate the exponential (2 A ) growth of spin states in GFMC calculations of neutron matter. This method combines two major themes in QMC. Auxiliary-field methods are used in the shell model Monte Carlo calculations [101], and several condensed matter systems in which the continuous spatial degrees of freedom have been eliminated, while diffusion Monte Carlo is another name for GFMC.
In the approach developed by Schmidt and Fantoni [100], the spatial parts, (i.e. kinetic energy and spin-independent interactions), of the Hamiltonian are propagated as in GFMC and the spin-dependent interactions between neutrons are replaced by interactions of neutrons with auxiliary fields. Integrating over the auxiliary fields reproduces the original spin-dependent interaction. In addition, a constraint analogous to the fixed-node approximation in GFMC is introduced, by requiring that the real part of the overlap with a trial function remains positive.
More recently, Schmidt and Fantoni [102] have carried out calculations with a realistic Hamiltonian consisting of Argonne v ′ 8 NN interaction used in the GFMC calculations [66], and the UIX TNI. The v ′ 8 contains the main parts of A18, and the difference between the two is treated perturbatively. Results of calculations with 38 neutrons in a periodic box with finite size corrections have been obtained. They are ∼ 5% below the variational energies obtained with the methods described above. For example, the AFDMC and variational energies for Argonne The trial functions used to constrain the present AFDMC calculations are rather simple without any spin correlations. In contrast it is possible to use more accurate variational wave functions with spin correlations to constrain the GFMC calculations. Carlson [98] has compared AFDMC and GFMC results for 14 neutrons in a periodic box. At ρ = 0.15 fm −3 the GFMC energy (220 ± 1) is about 7 % below the AFDMC result of 236.4 ± 1.5 MeV. From these we conclude that the variational PNM energies given in the last section may be ∼ 12 % above the exact values for the A18 and UIX interactions. The error in SNM E(ρ) is probably twice as large. As mentioned in the conclusions (sect.7), an overestimation of the E(ρ) of neutron matter by 12 % has a rather small effect on the predictions of neutron star properties.
The AFDMC is more accurate than the present variational method, and it is also more versatile. For example, it can be used to study matter with long range spin-isospin order induced by π 0 -condensation discussed in sect 4.2.

Hadronic and Quark Matter
It is likely that at high densities more general form of matter containing hadrons besides the nucleons, called hadronic matter (HM), has lower energy. The possibilities are that it contains negatively charged mesons like pions or kaons, or other hyperons such as Σ − or Λ. Finally, it is expected that at a high enough density there will be a transition to quark matter (QM) in which the quarks are not clustered into nucleons or hadrons.
The interactions between hyperons and nucleons, and between kaons and nucleons are not as well known as those between nucleons, and the energy of quark matter is difficult to calculate realistically. Therefore the transition densities from NM to HM or QM are rather difficult to calculate reliably. We review the recent estimates.

Kaon Condensation
Kaon condensation in dense matter was suggested by Kaplan and Nelson [103], and has been discussed in many recent publications [104,105]. Due to the attraction between K − and nucleons the kaon energy decreases with increasing density, and eventually if it drops below the electron chemical potential in NM, a Bose condensate of K − will appear. The key quantities of interest are the electron and kaon chemical potentials in NM. The former is obtained from the β-equilibrium condition, µ e = µ n − µ p relating electron, neutron and proton chemical potentials. The µ e obtained from the A18+δv+UIX * interactions with variational calculations is shown in Fig.10.
In neutron matter at very low densities, the interparticle spacing is much larger than the range of the K − n interaction, and the kaon interacts many times with the same nucleon before it encounters and interacts with another nucleon. Therefore one can use the scattering length, a K − n , as the "effective" kaon-nucleon interaction. In this low density limit the kaon energy deviates from its rest mass by the Lenz potential, and is given by [106] : where m R = m K m n /(m K + m n ) is the kaon-neutron reduced mass. The scattering length extracted from data is (−0.37 − i0.57) fm; its imaginary part is due to the open Λ − π − channel in vacuum. In the density region of interest to kaon condensation the kaon energy is too small for the K − + n → Λ + π − reaction to occur. Using effective lagrangians based on chiral perturbation theory Brown et. al. [104] estimate a K − n to be −0.41 fm in absence of reaction channels. There are two corrections to the Lenz energy at small densities. Including these, the kaon energy ω obtains the form [105,106]: where ξρ is the inverse correlation length. The first correction factor is a relativistic effect obtained from the Klein-Gordon equation, while the second factor is from the theory developed by Ericson and Ericson for propagation of mesons in nuclear matter. The relativistic correction decreases the kaon energy, while the correlation correction increases it. As the density increases further, and the interparticle spacing becomes of the order of the range of the interaction, the kaon will simultaneously interact with two or more nucleons and the Lenz approximation will break down. In the high density limit the kaon energy deviates from its rest mass by the Hartree potential: where v K − n is the K − n interaction potential. As shown in Ref. [107], the Hartree potential is considerably less attractive than the Lenz potential and thus ω Hartree > ω Lenz . The transition from the Lenz to Hartree limits has been recently studied [106] with a variety of methods including exact calculations for simple cubic crystal model of neutron matter. For reasonable interaction range the transition begins at very low densities (< 0.1ρ 0 ), and the Hartree limit is essentially reached by ∼ 3ρ 0 . There are no relativistic corrections to the Hartree energy of kaons condensed in the state with zero momentum [106] provided the K − n interaction is dominated by the Weinberg-Tomozawa vector potential. The typical resent results for kaon energy in neutron star matter are shown in Fig.10. The top solid line is obtained with a Wigner-Seitz [106] calculation for pure neutron matter; it is exact in both the low and high density limits, and gives essentially the Hartree energy at ρ > 3ρ 0 . The next curve shows the estimated Hartree results for NM containing ∼ 15 % protons, below that is the ω Lenz for neutron matter, while the lowest curves uses K − N scattering amplitude calculated in matter including relativistic, correlation and proton fraction corrections [105]. The recent [73] estimate for the central density of a 2.0 M ⊙ star is ∼ 5ρ 0 , and the results shown in Fig.10 indicate that kaon condensation is unlikely at densities lower than that. The heaviest stars (M = 2.2M ⊙ ) made up of NM are predicted to have central densities of ∼ 7ρ 0 , and the possibility of kaon condensation in their cores can not be ruled out. However, presence of π − , or Σ − or quark drops will decrease µ e as discussed in the following subsections and may make kaon condensation unlikely in even the most massive stars.

Charged Pion Condensation
Negatively charged π − -mesons will condense in matter when their chemical potential becomes lower than µ e , as suggested by Migdal [94] and by Sawyer and Scallapino [95]. In the seventies and eighties this possibility was studied by many researchers. Their work has been reviewed in Ref. [96], and we discuss it rather briefly. Fig. 10 suggests that, in absence of interactions π − with zero momentum will condense in matter at a rather low density of ∼ 1.5ρ 0 when µ e exceeds their rest mass of 139 MeV. However, the π − n S-wave interaction is repulsive, and raises the energy of zero momentum pions sufficiently above the estimated µ e . The recently found [108], deeply bound π − -nucleus atomic states, are influenced by this repulsion. The K − n S-wave interaction, on the other hand is attractive, leading to the possibility of kaon condensation discussed in the last section.
The energy of π − , having momenta of the order of 2 fm −1 , ∼ 400 MeV without interactions, is reduced in matter by the π − n P-wave interaction due to the coupling of the π − to p-n and N -∆ particle-hole states. Their energies and couplings to π − are calculated with effective forces described with Landau parameters whose density dependence is not well established. With plausible values for the Landau parameters a second order transition with π − condensation is predicted at ∼ 2ρ 0 ; however it is not expected to have a large effect on the EOS [96].
The Japan group also predicts a first order transition to an interesting phase with both π − and π 0 condensation at a density of ∼ 4.5ρ 0 . It has a significant effect on the EOS. Matter in this phase has spin aligned layers as discussed earlier in sec. 3.6 and illustrated in Fig. 11. The condensed π 0 -mesons have momenta in the Z-direction, perpendicular to the layers, while that of π − is in the X-direction in the plane of the layers. This way, when a proton absorbs a π − form the condensate and becomes a neutron its spin direction also gets flipped maintaining the attractive interaction with the π 0 field. It is a challenge to calculate the energy of matter in this interesting phase from bare nuclear forces.
The Illinois group, working with bare nuclear forces, predicts the first order π 0 -condensation at a much lower density of ∼ 1.5ρ 0 (sec. 3.6). However, the total decrease in energy of matter at ∼ 5ρ 0 due to pion condensation, estimated by the two groups: ∼ 80 by the Japan and ∼ 60 MeV/nucleon by Illinois, are not too different. The recent EOS of the Illinois group [73] contains this energy gain.

Hyperonic Matter
The possibility of hyperons contributing to the ground state of dense matter has been considered since 1959 [6,89]. The negatively charged Σ − and ∆ − will occur when their chemical potential becomes less than µ n + µ e in matter. With A18+δv+UIX * (A18+δv) interactions, the µ n + µ e exceeds the rest mass (1197 MeV) of Σ − hyperon at density of ∼ 2.2 (2.6) ρ 0 . If the Σ − -nucleon interactions are negligible, they will occur in matter via the weak interaction e − +n → Σ − +ν e at these densities mainly on account of the large µ e ∼ 170 MeV. Similarly the neutral hyperons such as the Λ will occur when their chemical potentials become less than µ n . For the above interaction models the µ n exceeds the Λ rest mass (1116 MeV) at densities of 2.9 and 3.7 ρ 0 , while the predicted central densities of 1.4 M ⊙ stars are 3.4 and 5.1 ρ 0 .
Much less data exist on hyperon-nucleon (YN) interactions than on the NN, and therefore their models are less constrained. LOBHF calculations [109] using the older Nijmegen soft core YN [110] and either Paris or Argonne v 14 NN interactions show that the thresholds for Σ − and Λ to appear are not much moved by YN interactions. For example, with the Paris NN interaction the threshold densities with (without) YN interactions are 3.0 (2.7) and 3.6 (3.7) ρ 0 for Σ − and Λ respectively. The Nijmegen group [111,112] has recently constructed boson exchange interaction models based on YN and NN data base using SU(3) symmetry. LOBHF calculations [113] using these models for NN as well as YN interactions give 2.2 ρ 0 for the threshold density of Σ − , however for Λ it is pushed beyond 7 ρ 0 .
The above LOBHF calculations use only two-body interactions without boost corrections and TNI, and fail to explain the saturation density of SNM as mentioned in sec. 3.4. Both the BHF groups [109,113] find that the threshold densities are lowered after including TNI effects. For example, the Catania group predicts them to be ∼ 2.1 and 2.4 ρ 0 with models including TNI and also with Dirac-Brueckner calculations. However, Σ − N N and ΛN N three body forces should also be included along with the TNI for consistency. The binding energies of Λ-hypernuclei [114] suggest that the ΛN N interaction is as strong as the TNI, while there is no data on the Σ − N N .
Ignoring three-body forces and boost interactions in both nucleon and hyperon matter, the LOBHF calculations with the unified NN, YN and YY interaction model [113] indicate that admixtures of Σ − lower the NM energy by only ∼ 25 MeV per nucleon at 5 ρ 0 . At this density the A18+δv+UIX * model gives an energy of ∼ 200 for β-stable NM. If the energy gain due to Σ − , is not much changed by three-body interactions, we should expect ∼ 10 % effects on the EOS of neutron star matter at ρ > 2ρ 0 due to hyperons; however, if present they would lower the µ e and increase the proton fraction x p significantly.

Quark Matter
When matter is compressed to densities so high that the quark cores of nucleons overlap substantially, one expects the nucleons to merge and undergo a phase transition to quark matter (QM). The EOS of both HM and QM are necessary to calculate effects of this transition in neutron stars. The interface properties are also needed to study the important mixed phase regions.
At present, lattice QCD can only treat the case of zero baryon chemical potential and is therefore not useful for neutron stars. Lacking a full theory, the simple Bag model is used to estimate the QM EOS. In this model the QM energy has a volume term denoted by the bag constant B; it represents the difference in the energies of the vacua occupied by hadron and QM, and is responsible for the confinement of quarks within nucleons in nuclei. The term dominating at high densities is the energy of noninteracting u, d and s-quarks; it is calculated neglecting the mass of u, and d-quarks, and using a typical value of ∼ 150 MeV for that of the s-quark. Since the quark Fermi energies are much larger than their masses, the QM properties are not too sensitive to the chosen mass for the s-quark. There is no one-gluon exchange interaction energy between quarks of different flavor, while that between quarks of of the same flavor i is given by (2α/3π)E i per quark of flavor i [115]. Here E i is the average kinetic energy per quark, and α is the strong interaction coupling constant, assumed to have a value of ∼ 0.5. All higher order gluon-exchange interactions are neglected; their contribution is presumably subsumed in the bag constant whose value is poorly known. Two representative values are B = 122 [116] and 200 MeV fm −3 [117].
The equilibrium conditions for uniform QM containing u, d, s-quarks, electrons and muons are: The energy densities of charge neutral QM and NM are compared in Fig.12. In the interesting region of ρ ∼ 1 fm −3 the total energy density of quark matter is about 1200 MeV fm −3 , of which only 122 or 200 MeV fm −3 comes from the bag. Assuming the A18+δv+UIX * EOS for NM the first order phase transition to QM is indicated at ∼ 1 baryon/fm 3 . It has been recently suggested that at high densities QM may have color superconductivity resulting from non-perturbative attraction between quarks. In QM with only u and d-quarks this invariably leads to the possibility of a diquark condensate which breaks global color invariance [118]. The associated color gap is estimated to be of the order of 100 MeV. When the s-quarks are included there are many possible phases [119].

Mixtures of Phases in Dense Matter
The phases of matter considered in the past sections are uniform and locally charge neutral, whereas bulk matter needs only to be charge neutral on average. For example, iron metal has positively charged regions occupied by iron nuclei, and the space in between is negatively charged by electrons. Generally the ground state of matter can have a mixture of regions occupied by different phases with the constrain of overall charge neutrality.
Matter in the outer crust of neutron stars, at ρ ≤ 0.002ρ 0 , is believed to be like terrestrial matter made up of neutron rich nuclei in electron gas. The inner crust, on the other hand, is believed to have a mixture of regions with positively charged NM composed of neutrons and protons and PNM with only neutrons [10], both immersed in a nearly uniform electron gas. These mixtures occur in the density range of ∼ 0.002 to 0.6 ρ 0 , beyond which uniform NM is believed to be the ground state. More recently mixed phases of NM and QM [120,121], and condensates of pions, kaons [122], or hyperons in NM have also been considered.
The Coulomb energy of a single phase uniform matter, due to fluctuations in the electron and hadron or quark densities, is negligible; however that of matter with mixed phases is not. For example, matter at ρ = 0.3ρ 0 has drops of NM in PNM. Their size is determined by a competition between Coulomb and surface energies; large drops have too much Coulomb energy and small drops to much surface energy per nucleon. It is necessary to know the energy of the interface to predict the nature of the mixed phase region quantitatively. Those for the interface between NM and PNM have been calculated from energy-density functionals of NM [10]; while those for interfaces of NM and QM are not well estimated. In this section we discuss QM and NM mixed phases, to avoid duplicating the review by Pethick and Ravenhall [10], and because they can influence the mass limit of neutron stars.

Equilibrium Conditions for Coexistence of QM and NM
Neglecting surface and Coulomb effects, the equilibrium conditions for the coexistence of of QM and NM at zero temperature are that they have equal pressures, and it costs no energy to convert a neutron or a proton in NM into quarks in QM. The last condition amounts to The electron density, and hence the µ e , is assumed to be the same in QM and NM, therefore the β-equilibrium conditions are: And finally, the charge neutrality condition is given by: Here ρ ch QM,N M are the charge densities of QM and NM, ρ e,µ are the electron and muon densities, and f is the fraction of space filled by QM. A graphical representation of these equilibrium conditions is given in ref. [73].
As was first pointed out by Glendenning [120], the QM and NM mixed phases can span a wide density region. For example, with the A18+δv+UIX * NM EOS, the mixed phases begin at 4.6 (3.4) and end at 11.3 (9.1) ρ 0 for B = 200 (122) MeV fm −3 (Fig.12). At lower densities we have uniform NM and QM at higher. The QM fractional volume increases almost linearly from 0 to 1 in this transition range over which µ e → 0. The uniform QM has almost equal number of u, d, and s-quarks, and is essentially charge neutral without electrons. An important consequence of β-equilibrium is that the QM is negatively charged at the beginning of the transition, where µ e is large. By immersing negatively charged drops of QM in the positively charged NM we can remove some of the high energy electrons and increase the proton fraction in NM. Both of these lower the energy density of matter.

Structure of Mixed Phase Matter
Matter with mixed phases has additional structure due to interfaces dividing the regions occupied by the two phases. The surface and Coulomb energies associated with these interfaces, neglected in above section, raise the energy of the mixed phase matter as well as determine the topology and structural length scales. The surface energy can presumably be estimated from the surface tension, σ. For NM σ ≃ 1MeV/fm 2 whereas for QM it poorly known; typical values are in the range 10-100 MeV/fm 2 [121]).
Denoting the dimensionality of the structures by D (D = 3 for droplets and bubbles, D = 2 for rods and D = 1 for plates) the surface energies are generally [10] and, for D = 3 the Coulomb energy, E C = (3/5)Z 2 e 2 /R, where R the size of the structure and is the excess charge of the droplet compared with the surrounding medium. General equations for E C (R, D) are given in ref. [10]. Minimizing the energy density with respect to R we obtain the usual result that E S = 2E C at equilibrium. Minimizing with respect to the continuous dimensionality as well thus determines both R and D. For droplets (D = 3) the equilibrium radius is found to be: A droplet of symmetric NM in vacuum has a surface tension σ = 1 MeV·fm −2 for which (26) gives R ≃ 4 fm, which agrees with the fact that nuclei like 56 Fe are the most stable form of matter at low density. For QM droplets both the surface tension and charge densities are larger but the estimates of R are similar. At the beginning of the mixed phase region we expect that spherical droplets of QM, with R ∼ 5 fm, will form a BCC lattice in uniform NM. They would have baryon number of the order of few hundred, and a negative charge of similar magnitude. As density increases and f approaches 0.5, the drops would merge and form rods, which merge further on to form sheats. When f > 0.5 the NM sheats break up into rods, and then into drops and eventually disappear when f = 1. This scenario is similar to that in the inner crust; at low densities there are drops of NM occupying a small fraction of space. By ρ ∼ 0.6ρ 0 NM occupies all space via a similar set of mixed phases.
An other effect of the Coulomb and surface energies is that they decrease the density range covered by the mixed phase region. In particular, the lower density edge of this interesting region may be pushed up by almost ρ 0 if σ is in the 10 to 50 MeV fm −2 range [121]. The energy density of the mixed phase matter is also raised by a few MeV fm −3 in this case. Finally, if σ were to be large ( > ∼ 70 MeV fm −2 ) the mixed phases may not be energetically favorable, and there will be a simple first order phase transition from NM to QM with a density discontinuity. One should bear in mind that even if the droplet phase were favored energetically, it may not be realized in practice if the time required to nucleate QM drops is too long compared to pulsar ages.

Neutron Star Observations and Predictions
The gross structure of neutron stars has been predicted using very many EOS, phenomenological as well as based on realistic models of nuclear forces [43,51]. Of these we consider only those based on realistic models primarily because one can always find phenomenological energy density functionals or Lagrangians which reproduce their EOS.
Typical results for nonrotating stars with maximum mass, and with M = 1.4M ⊙ , obtained by recent calculations, are listed in Table 2. The results for A18 without boost correction δv are listed primarily for reference. This correction is unambiguous [74], and must be added to obtain reliable results. Those for A18+δv are also to be taken less seriously, because it gives too large value for ρ 0 . The TNI used with the Paris NN interaction [85] is of the Urbana form with parameters determined by reproducing the empirical SNM properties (see sec. 3.4). In A18+UIX and Paris+TNI models the δv is not considered explicitly; it is approximately subsumed in the TNI fitted to data. Of the three Bonn models, Bonn A comes closest to reproduce the empirical properties of SNM [88] with Dirac-Brueckner (DBHF) method. These calculations include the δv as well as many-body forces generated via Z-diagrams.
The A18+δv+UIX * , A18+UIX, Paris+TNI and Bonn-A (DBHF) models come close to reproducing the empirical ρ 0 ; the later two fit the SNM binding energy; while the former models fit binding energies of light nuclei via exact calculations, since the energy of SNM can not yet be calculated reliably. Nevertheless these four "realistic" models of NM give rather similar results which are not too different from those of the 1988 calculations of Wiringa, Fiks and Fabrocini [91] with the older Urbana-Argonne interactions now replaced with A18 and UIX.
The effect of the possible appearance of QM drops in high density matter has been studied with the A18+δv+UIX * model. The M m is reduced to 2.02 and 1.91 M ⊙ for bag-constant values B = 200 and 122 MeV Fm −3 respectively, while the predictions for 1.4 M ⊙ stars remain unchanged. Presence of either kaons or hyperons in dense matter is unlikely to have much of an effect on the 1.4 M ⊙ stars due to their low central density, while that on the mass limit is difficult to estimate quantitatively. For example, if kaons were to condense in matter at ρ = 5ρ 0 and limit ρ c to < 5 ρ 0 , the M m of the four realistic models will drop to ∼ 2.0, 2.3, 1.7 and 2.0 M ⊙ respectively; while if hyperons were to lower the energy of matter at ρ = 5ρ 0 by 25 MeV per baryon, the M m would be reduced by ∼ 0.2 M ⊙ .
The mass radius relation obtained with models based on the A18 interaction are shown in Fig.13. Results of A18 and A18+δv models are given primerily for comparison. As expected the harder EOS give larger M m and predict larger radii. The differences between the radii predicted by the realistic models is only ∼ 10 %.

The Mass Limit
The observed mass of Hulse-Taylor pulsar B1913+16 of 1.4411 ± 0.00035 [12] shows that M m > 1.44M ⊙ . All the radio pulsars in known neutron star, and neutron star-white dwarf binaries have masses with lower limits less than 1.44M ⊙ . The X-ray pulsar Vela X-1, which orbits a supergiant, however is consistently estimated to have a larger mass of ∼ 1.9. The motion of this star is perturbed from being pure Keplerian, presumably by tidal forces exerted by the neutron star, and its present mass estimate, 1.87 +0.23 −0.17 [15], indicates that M m > 1.7M ⊙ at 95 % confidance level. Finally, if the QPO's indeed originate from the innermost stable orbit [18,19], then M m > 2M ⊙ . These mass limits are compatible with predictions of realistic NM models.
On the other hand there is no evidence that SN 1987A produced a neutron star. Its observed luminosity is now well below the 10 38 ergs/s Eddington limit, suggesting that no neutron star was produced in this supernovae [123]. If we assume that the total mass, M T ot , of the collapsed core plus the matter that fell back on to the core after the explosion, went into a black hole, then neutron star M m must be less than ∼ 0.9M T ot . The factor 0.9 takes into account the ∼ 10 % gravitational binding energy of the neutron star. Bethe and Brown [23] estimate M T ot ∼ 1.73M ⊙ using supernovae calculations by Wilson and Mayle, and conclude that M m < 1.56M ⊙ . Uncertainties in these arguments have been discussed by Zampieri et. al. [124]. If the conclusion is found to be valid, then there must be other explanations for the Vela X-1 observations and the QPO, and the NM prediction for the M m is too large.

Temperatures, Cooling and Radii
Neutron stars are born with interior temperatures of the order 10 12 K, but cool rapidly via neutrino emission to temperatures of the order 10 10 K within minutes and < ∼ 10 6 K in 10 5 yr. Spectra observed in X-ray or UV bands for nearby pulsars have in some cases black-body components from which surface temperatures of order T ∼ 10 6 K are extracted for pulsars of age 10 3 − 10 6 years. It is, however, unclear how much of the observed radiation is due to pulsar phenomena, to a synchrotron-emitting nebula or to the neutron star itself. In other cases upper limits have been set from the absence of X-rays. The surface temperatures are compatible with predictions from standard modified URCA cooling processes [125] n + n → n + p + e − +ν e , n + p + e − → n + n + ν e .
Faster cooling processes as direct URCA or due to quark matter, kaon or hyperon condensates generally lead to considerably lower temperatures [126]. To be consistent with observed surface temperatures the exotic coolant can only exist in a minor portion of the neutron star or it is superfluid whereby cooling is suppressed by factors of ∝ exp(−∆/T ), where ∆ is the pairing gap. The Hubble Space Telescope (HST) has observed one thermally radiating neutron star RX J185635-3754 with surface temperature T ≃ 6 × 10 5 K≃ 50 eV [20].
Its distance is less than 120 pc from Earth and should soon be determined more accurately by HST parallax measurements. Circumstancial evidence indicate a distance of ∼80 pc [51] which leads to a black-body radius of ∼ 12 − 13 km from its luminosity and temperature. Such radii would agree well with predictions of realistic NM EOS (Fig.13) for M ≃ 1 − 2M ⊙ .

Glitches and Superfluidity
Sudden spin jumps, called glitches, superimposed upon otherwise gradual spin down have been observed in most of the younger isolated pulsars [127]. Since their discovery the Crab and Vela pulsars have each produced about a dozen glitches with period changes ∆P/P of the order of 10 −8 and 10 −6 respectively. In post-glitch relaxation most of the period increase ∆P decays.
Many mechanisms have been proposed to explain the glitches [128]. The most plausible of these attributes glitches to the angular momentum stored in the rotating superfluid neutrons in the inner crust [129,10]. The magnetic torque slows down the crust and most of the star except for these superfluid neutrons. Their angular momentum is stored in vortices pinned to nuclei in the inner core, until an instability occurs that leads to vortex depinning and sudden angular momentum transfer to the crust, leading to the glitch. At subnuclear densities in the crust, 1 S 0 pairing between neutrons leads to gaps of order ∼1 MeV [128]. In NM at ρ > ρ 0 this pairing gap vanishes, but 3 P 2 pairing of neutrons and 1 S 0 pairing of protons may occur [128,84].
The size of the glitches sets a lower limit on the moment of inertia of the superfluid in the inner crust which in turn sets a lower limit on the neutron star radius for a given mass [127]. Assuming that the mass of Vela pulsar is 1.4M ⊙ , a conservative limit on its radius is R ≥ 9 km; it is compatible with predictions of most EOS.

Conclusions
Since the discovery of pulsars a significant effort has been devoted to accurately calculate properties of dense NM from realistic models of nuclear forces. Exact calculations of NM are still out or reach, however the new AFDMC methods (sect. 3.7) may eventually succeed. The present variational upper bounds seem to be above the true energies by ∼ 12 %. Such an error does not have serious consequences on the predicted properties of neutron stars. For example, an EOS obtained by reducing the variational energies, without rest mass terms, by 12 % reduces the maximum mass of A18+δv+UIX * model by 2.3 % to 2.14 M ⊙ , and the radius of 1.4 M ⊙ star by 2.9 % to 11.2 km. Larger uncertainties stem from the fact that the double π 0 and π − condensation scenario illustrated in Fig.11 has not yet been calculated with realistic interactions, though it appears unlikely that it will influence the NM EOS by much more than 10 %.
Local models of two-nucleon interaction seem to be now converging. The predictions based on the 1988 calculations with Argonne 14 interaction are not too different from those of the 1998 calculations with the more accurate A18. It also seems likely that the local models give a fairly accurate description of two-nucleon interaction. A concern is that the present models of TNI are based on fits to a rather limited set of data, and are not as precise as the NN-interaction models. However, addition of the UIX * TNI to the A18+δv increases the maximum mass by ∼ 20 % and R(1.4) by 13 % (Table 2). These changes may be important but they are not very large.
The present models of kaon-nucleon and hyperon-nucleon interactions are based on very limited data, and we have none on K − N N and Σ − N N three-body forces. These could have significant effect on the threshold densities for kaons and hyperons to appear in dense matter. Hopefully advances in QCD and quark-models will provide a more rigorous framework to describe these interactions, and calculate properties of quark matter. The bag model estimates of QM EOS may have significant corrections at densities of interest in neutron stars.
¿From present observations there seem to be three possible scenarios for the limiting mass of neutron stars. If QPO's are indeed due to accretion from the innermost stable orbit, then the NM predictions of M m ∼ 2.2M ⊙ are reasonable, and strange baryons and quark drops do not soften the EOS of matter at ρ < 7ρ 0 significantly. If the Vela X-1 mass measurement is correct, but QPO's have some other origin, then M m could be ∼ 1.8 M ⊙ , indicating some softening of the NM EOS. However, if the present interpretation of QPO's and Vela X-1 mass measurements are both faulty, and M m is as small as 1.56 M m as estimated from the absence of a neutron star in SN 1987A, then a significant softening of the NM EOS by phase transitions is indicated. Further observations will hopefully clear this situation.
Phase transitions such as NM to QM, can soften the EOS significantly. Fortunately these can have a measurable effect on the spin down of a rapidly rotating star in favorable cases, as has been recently pointed out [130,131]. Consider the case of a rapidly rotating star whose central density is close to a first order phase transition. As the star slows and the central pressure increases due to decrease of the centrifugal force, the core matter will change its phase and become more dense at a critical angular velocity Ω c . This decreases the moment of inertia, which assumes the characteristic form: for Ω < Ω c . Here, c 1 and c 2 are small parameters proportional to the density difference between the two phases, and c 2 = 0 for Ω > Ω c . In order to make contact with observation, the temporal behavior of angular velocities must be considered. The pulsars slow down at a rate given by the loss of rotational energy, believed to be given by: d( 1 2 IΩ 2 )/dt ∝ −Ω n+1 , where n = 3 for dipole radiation, Eq. (1) and n = 5 for gravitational radiation. With the moment of inertia given by Eq. (28) the angular velocity can be calculated. The corresponding braking index, n(Ω) =ΩΩ/Ω 2 , depends on the second derivative of the moment of inertia, I ′′ = dI/d 2 Ω. Using Eq. (28) we obtain: which exhibits a characteristic (Ω c − Ω) −1/2 singularity as Ω approaches Ω c from below. Observations of the braking index of a rapidly rotating, new born pulsar would be very interesting. All realistic NM EOS predict that the radius of neutron stars with a mass of 1 to 1.5 M ⊙ is ∼ 11 to 12 km. Future high resolution Chandra and XMM space observatories will hopefully be able to measure black-body spectra and detect gravitationally redshifted spectral lines from several stars. Such observations will help determine masses, radii and temperatures uniquely if the distance of the star is known. It is important to know the radius of a 1.4 M ⊙ star, because that would test the EOS in the ρ < ∼ 3ρ 0 region in which large modifications of NM EOS are not expected on the basis of our present, naive estimates of kaon-nucleon and Σ − -nucleon interactions.    Figure 5: Charge densities of two protons located one fm apart at Z = ±0.5, obtained by inverting the dipole approximation to proton charge form factor. The sharp peaks at Z = ±0.5 are unphysical, they will be rounded off by relativistic corrections and improved data on proton form factor.    [73] and LOBHF (LOB) method [83]. The true results for this interaction are expected to be few MeV below the VCS.   Figure 13: Neutron star gravitational mass,M (R), in solar masses vs. radius in kilometers for the four models described in the text. Full curves are for β-stable matter and dotted ones for pure neutron matter. The dashed curve, FPS, is from [50]. Fig. from [73].