EQUATION OF STATE AND NEUTRON STAR PROPERTIES CONSTRAINED BY NUCLEAR PHYSICS AND OBSERVATION

Neutron stars, apart from being systems for investigating such diverse topics as theories of gravity and the interstellar medium, are unique laboratories for studying matter at high densities. In neutron stars, matter ranges from nuclei embedded in a sea of electrons at low densities in the outer crust, to increasingly neutron-rich structures in the inner crust, to the extremely neutron-rich uniform matter in the outer core, and possibly exotic states of matter at high densities in the inner core (Haensel et al. 2006; Lattimer 2012). The theoretical understanding of nuclear matter and atomic nuclei over such a range of densities and isospin asymmetry is a current frontier in nuclear science. The equation of state (EOS) of dense matter is also a key ingredient in modeling neutron star and black hole formation and, in particular, the gravitational wave signal from mergers of binary neutron stars (Andersson et al. 2011; Bauswein & Janka 2012; Bauswein et al. 2012) and neutron star-black hole mergers (Lackey et al. 2012) is sensitive to it. Consequently, future searches with advanced LIGO and LISA are expected to provide information about the high-density EOS.

In the past, the EOSs of dense matter most commonly used in astrophysical simulations have been based on phenomenological nuclear interactions (for a review, see, for example, Heiselberg & Pandharipande 2000). Two-nucleon (NN) interactions are usually constructed to fit NN scattering data at low energies. In addition, three-nucleon (3N) interactions were introduced because, first, they exist on theoretical grounds and, second, NN interactions alone cannot reproduce the properties of nuclei and nuclear matter. The many-body problem based on those interactions is very challenging due to strongly repulsive forces at small relative distances, which lead to highly correlated wave functions and require nonperturbative many-body methods.

In recent years, the development of chiral effective field theory (EFT), following the pioneering work of Weinberg (1990, 1991), has provided the framework for a systematic expansion for nuclear forces at low momenta, where nucleons interact by pion exchanges and short-range contact interactions whose parameters can be fixed on the basis of two- and few-body observables (Epelbaum et al. 2009). Chiral EFT explains the hierarchy of two-, three-, and weaker higher-body forces and provides estimates of the theoretical uncertainties. In a recent Letter (Hebeler et al. 2010), we have shown that microscopic calculations based on chiral EFT interactions constrain the properties of neutron-rich matter up to nuclear saturation density to a high degree. On the basis of laboratory experiments and theory, our knowledge of the EOS at densities greater than one to two times the saturation density is limited. However, information may be obtained from measurements of neutron star masses. In particular, the recent discovery of a neutron star with a precisely determined mass of 1.97 ± 0.04 M_☉ (Demorest et al. 2010), the heaviest to date, is extremely important, as it rules out a large number of EOSs based on exotic degrees of freedom like hyperons or deconfined quarks. Such constituents soften the EOS at high densities and, without fine-tuning, are generally incompatible with a neutron star of such a large mass. Further information may be obtained, for example, from modeling X-ray bursts and quiescent low-mass X-ray binaries (Özel et al. 2010, 2012; Steiner et al. 2010, 2013; Güver & Özel 2013), but the neutron star properties deduced are more model-dependent than the direct mass constraint from a very heavy neutron star.

By extending our microscopic results for the EOS at low densities in a general way to higher densities, we showed in Hebeler et al. (2010) that it is possible to derive systematic constraints on the EOS and on the radii of neutron stars. The high-density extensions we used were only constrained by causality and by the heaviest observed neutron star at that time which has a mass of 1.65 M_☉. In this paper, we present details of these calculations and generalize and improve our approach in several ways: first, we require that the EOS is consistent with the observation of a 1.97 M_☉ neutron star (Demorest et al. 2010). To investigate the sensitivity to the possible future discovery of neutron stars of higher mass, we also consider a second case where the EOS supports a neutron star of a mass of 2.4 M_☉. Second, we generalize the microscopic neutron matter calculations. Our previous results were based on renormalization-group-evolved NN interactions plus the leading 3N interactions from chiral EFT. The renormalization-group evolution improves the many-body convergence (Bogner et al. 2010), but introduces uncertainties because the evolution was limited to NN forces. Here, we show that calculations based on unevolved chiral EFT interactions are in good agreement with the previous calculations (see Section 2). Third, we improve the way beta equilibrium is incorporated (see Section 3) and include explicitly the crust EOS below the crust-core transition density, which we calculate in Section 4. Finally, we generalize the piecewise polytropic extensions of the EOS to higher densities by allowing more density regions and also refine the step size of the variation of the polytropic parameters (see Section 5). We present our results for the nuclear EOS and the structure of neutron stars in Section 6. This shows that the constraints are robust and not significantly altered by the generalizations and improvements of the microscopic calculations and EOS extensions. For use in astrophysical simulations, we construct in Section 7 three representative EOSs consistent with the constraints from nuclear physics and observations and provide numerical data in Tables 5–7.

Our microscopic neutron matter calculations are based on chiral NN and 3N interactions. Neutron matter presents a unique system in chiral EFT because only the long-range two-pion-exchange parts of the leading 3N interactions contribute (Hebeler & Schwenk 2010). This is because three neutrons cannot interact via point-like S-wave interactions due to the Pauli principle. Moreover, the leading one-pion-exchange 3N interaction does not contribute in neutron systems because of the particular spin-momentum structure of this interaction.

The leading chiral 3N forces have been shown to give important contributions to the nuclear EOS and to properties of nuclei (Hammer et al. 2013). In particular, the saturation of symmetric nuclear matter is driven by 3N forces (Hebeler et al. 2011). While 3N contributions to the neutron matter EOS are smaller in size than for nuclear matter, they are still significant and crucial for predictions of observables like the nuclear symmetry energy (see Hebeler & Schwenk 2010 and Section 3). In this work, we include only the leading 3N forces, because it is presently possible to include only these beyond the Hartree–Fock level (Tews et al. 2013).

For calculations of the neutron matter EOS, the theoretical uncertainty is dominated by the uncertainties in the low-energy couplings c₁ and c₃, which determine the two-pion-exchange three-body interactions between neutrons, not by the many-body approximations (Hebeler & Schwenk 2010). These low-energy couplings relate πN, NN and 3N interactions. Their determination from πN scattering is within uncertainties consistent with the extraction from NN scattering, leading to c₁ = −(0.7–1.4) GeV⁻¹ and c₃ = −(3.2–5.7) GeV⁻¹ (Entem & Machleidt 2003; Epelbaum et al. 2005; Rentmeester et al. 2003; U.-G. Meißner 2006, private communication). Beyond the leading 3N forces, part of the subleading 3N forces are simple shifts of c₁ and c₃ (Epelbaum et al. 2009). Therefore, as in Hebeler et al. (2010), we take c₁ = −(0.7–1.4) GeV⁻¹ and c₃ = −(2.2–4.8) GeV⁻¹ for our calculations.

In general, nuclear forces depend on an intrinsic resolution scale Λ. Consequently, the nuclear Hamiltonian can be written in the form

$$\begin{equation} H(\Lambda) = T + V_{\rm NN}(\Lambda) + V_{\rm 3N}(\Lambda) + V_{\rm 4N}(\Lambda) + \cdots , \end{equation} \tag{ 1 }$$

where T denotes the kinetic energy, V_NN the NN interactions, V_3N the 3N interactions, etc. The renormalization group provides a powerful tool to systematically change the scale Λ, while low-energy observables are preserved (Bogner et al. 2010). The evolution to low momentum scales improves the many-body convergence due to a decoupling of low and high momenta in the Hamiltonian (Bogner et al. 2010). In general, the renormalization-group evolution changes all terms in Equation (1). The consistent evolution of 3N interactions in momentum space is a complex task and has been achieved only recently (Hebeler 2012).

We start from the chiral N³LO NN potential with Λ = 500 MeV of Entem & Machleidt (2003) and use the renormalization group to evolve this NN potential to low-momentum scales Λ = 1.8–2.8 fm⁻¹ (Bogner et al. 2007; Hebeler et al. 2007). At these scales, NN interactions derived from different initial potentials are very similar (Bogner et al. 2003). This universality can be attributed to common long-range pion physics and phase-shift equivalence. As a consequence, results of many-body calculations are rather insensitive to the particular choice of the initial NN interaction. Because the leading chiral 3N forces are of long-range character in neutron matter, they are expected to be to a good approximation invariant under the renormalization-group evolution for these cutoffs. Therefore, we use the leading chiral 3N forces determined by c₁ and c₃ also with low-momentum cutoffs.

For low-momentum cutoffs, the EOS of neutron matter and nuclear matter can be calculated with theoretical uncertainties in a perturbative expansion in which the Hartree–Fock approximation is the first-order term (Hebeler & Schwenk 2010; Hebeler et al. 2011). Our neutron matter results also include second-order corrections (for calculational details, see Hebeler & Schwenk 2010). For low-momentum cutoffs 1.8 fm⁻¹ ⩽ Λ ⩽ 2.8 fm⁻¹, we have checked that contributions from third-order particle–particle diagrams give only small contributions to the energy per particle, about 25 keV for Λ = 1.8 fm⁻¹ and 300 keV for Λ = 2.8 fm⁻¹ at the saturation density n₀ = 0.16 fm⁻³. In addition, we find that our second-order results are independent of the resolution scale to a very good approximation, and at the saturation density the maximal variation is about 1 MeV. We also find that the results are insensitive to the single-particle spectrum used. All these findings indicate that neutron matter is perturbative for low-momentum interactions, and show that the theoretical uncertainties of the many-body calculation are small (Hebeler & Schwenk 2010).

Figure 1 shows our results for the energy per particle of neutron matter up to the saturation density (using a Hartree–Fock spectrum). The blue band is based on chiral NN and 3N interactions, with a renormalization-group evolution for NN interactions. The width of the band is due mainly to the uncertainties of c₁ and c₃ in 3N forces (Hebeler & Schwenk 2010). For comparison, the dashed red lines present the range based on unevolved chiral NN interactions plus the same leading chiral 3N interactions. The remarkable agreement indicates that neutron matter is, to a good approximation, also perturbative for chiral NN interactions (for a detailed study, see Krüger et al. 2013). We also explicitly checked the size of the contributions at third-order in the many-body expansion. As expected, they are larger for unevolved chiral NN interactions due to the stronger coupling between low and high momenta in the Hamiltonian. However, they are still significantly smaller than the second-order contributions (Krüger et al. 2013).

Zoom In Zoom Out Reset image size

In Figure 2 we present the uncertainty band for the pressure of neutron matter based on chiral NN and 3N interactions (with a renormalization-group evolution for NN interactions). For comparison, we also show the Quantum Monte Carlo results of Gandolfi et al. (2012; shaded red band) based on the phenomenological Argonne v₁₈ NN potential plus 3N forces fitted to a symmetry energy of 32.0–35.1 MeV (see Table I of Gandolfi et al. 2012). The agreement of the results is remarkable, given that the Hamiltonian and the many-body methods are completely different.

Zoom In Zoom Out Reset image size

We extend the microscopic results for neutron matter to matter containing both neutrons and protons. To this end, we use for the energy per particle of asymmetric nuclear matter an expression that interpolates between the properties of symmetric nuclear matter and neutron matter. For we take the kinetic energy plus an expression for the interaction energy that is quadratic in the neutron excess 1 − 2x:

$$\begin{eqnarray} &&\epsilon (\bar{n},x) T_0 = \frac{3}{5} [ x^{5/3} + (1-x)^{5/3}] (2 \bar{n})^{2/3}- [(2 \alpha - 4 \alpha _L) \nonumber \\ &&\quad\times x (1 - x)+ \alpha _L] \bar{n} + [(2 \eta - 4 \eta _L) x (1 - x) + \eta _L] \bar{n}^{\gamma } \,, \quad\quad \end{eqnarray} \tag{ 2 }$$

where $\bar{n} = n/n_0$ and x = n_p/n denote the density in units of the saturation density and the proton fraction, respectively. T₀ = (3π²n₀/2)^2/3ℏ²/(2m) = 36.84 MeV is the Fermi energy of symmetric nuclear matter at the saturation density. Equation (2) does not include the contributions from rest masses. The corresponding result for the pressure P = n²∂/∂n is

$$\begin{eqnarray} P(\bar{n},x) T_0 &=& \frac{2}{5} \, n_0 [ x^{5/3} + (1 - x)^{5/3}] (2 \bar{n})^{5/3} \nonumber \\ \quad && - [ (2 \alpha - 4 \alpha _L ) x(1- x) + \alpha _L ] \bar{n}^2 \nonumber \\ \quad && + \gamma [ (2 \eta - 4 \eta _L) x (1- x) + \eta _L ] \bar{n}^{\gamma +1} . \end{eqnarray} \tag{ 3 }$$

The parameters α, η, α_L and η_L can be determined from the saturation properties of symmetric nuclear matter combined with the microscopic calculations for neutron matter of Section 2. For γ = 4/3 and empirical saturation properties of symmetric nuclear matter,

$$\begin{eqnarray} &&\epsilon (\bar{n}=1, x=1/2) = - B = - 16 \, {\rm MeV} \quad {\rm and} \nonumber\\ && P(\bar{n}=1, x=1/2) = 0 \,, \end{eqnarray} \tag{ 4 }$$

this results in α = 5.87, η = 3.81, and a reasonable incompressibility parameter

$$\begin{equation} K = 9 \left. \frac{\partial ^2 \epsilon (\bar{n},x)}{\partial \bar{n}^2} \right|_{\bar{n}=1,x=1/2} = 236\,{\rm MeV} \,. \end{equation} \tag{ 5 }$$

The precise value of K can be adjusted by modifying the exponent γ in Equations (2) and (3). However, as shown in Table 1, the predicted range for the symmetry energy and its density derivative depend very weakly on the particular choice of γ and the resulting value of K, because the leading density dependence of the symmetry energy is linear.

The parameters α_L and η_L are extracted from the calculated bands for the neutron matter energy and pressure of Figures 1 and 2. Our results are based on the blue bands (with renormalization-group evolution) unless stated otherwise. We have first checked that the parameterizations (2) and (3) provide excellent global fits for the energy and pressure up to a density n₁ ≈ 1.1 n₀. To determine α_L and η_L, we sample their values systematically and require that the resulting energy and pressure be within the uncertainty bands shown in Figures 1 and 2 for densities from 0.45 n₀ to 1.1 n₀. This leads to the allowed range for α_L and η_L shown in Figure 3, with correlated limits α_L = 1.18–1.59 and η_L = 0.64–1.11.

Zoom In Zoom Out Reset image size

The proton fraction x for matter in beta equilibrium is determined by minimizing, for a given nucleon density, the total energy per particle, Equation (2), plus the contributions from electrons and from the rest mass of the nucleons. This amounts to the condition that μ_n + m_nc² = μ_p + m_pc² + μ_e, where μ_n and μ_p are the neutron and proton chemical potentials without the rest mass contribution, or equivalently

$$\begin{equation} \frac{\partial \epsilon (\bar{n},x)}{\partial x} + \mu _e(\bar{n},x) - (m_n - m_p) c^2 = 0 \,. \end{equation} \tag{ 6 }$$

For an ultrarelativistic, degenerate electron gas, the chemical potential is given by $\mu _e(\bar{n},x) = \hbar c \, (3 \pi ^2 x n_0 \bar{n})^{1/3}$. The allowed ranges for α_L and η_L imply ranges for the proton fraction and the neutron and proton chemical potentials in beta equilibrium, which are given for the saturation density n₀ and for n₀/2 in Table 2. In the calculations we neglected the difference between the neutron and proton masses (1.3 MeV), which is small compared with μ_e ∼ 100 MeV. These ranges provide anchor points for other equations of state.

The parameterizations (2) and (3) also make it possible to reliably extract the symmetry energy S_v and its density derivative L,

$$\begin{equation} S_v = \frac{1}{8} \frac{\partial ^2 \epsilon (\bar{n},x)}{\partial x^2} \left |\right. _{\bar{n}=1, x=1/2} \; {\rm and} \; L = \frac{3}{8} \frac{\partial ^3 \epsilon (\bar{n},x)}{\partial \bar{n} \partial x^2} \left |\right._{\bar{n}=1, x=1/2} \,. \end{equation} \tag{ 7 }$$

The region for α_L and η_L translates into an allowed region for the symmetry energy S_v and the L parameter shown in Figure 4, after Lattimer & Lim (2013). In addition, we give in Table 1 the predicted ranges for S_v and L for different γ values, corresponding to different incompressibilities K = 210 MeV, 236 MeV and 260 MeV. The predicted range for γ = 4/3 nearly spans the ranges for the other γ values. This demonstrates that the extrapolation (2) is robust and that the theoretical uncertainty due to the choice of γ is very weak and clearly much smaller than the empirical bands shown in Figure 4.

Zoom In Zoom Out Reset image size

In Figure 4, we compare the S_v and L region predicted by our neutron matter results with values extracted from other data (Lattimer & Lim 2013). It is striking that the neutron matter results lead to the strongest constraints. These agree well with constraints obtained from energy-density functionals for nuclear masses (orange band; Kortelainen et al. 2010) and from the ²⁰⁸Pb dipole polarizability (yellow band; Tamii et al. 2011). In addition, there is good agreement with studies of the Sn neutron skin (light blue band; Chen et al. 2010), of isotope diffusion in heavy ion collisions (green band; Tsang et al. 2009), and of giant dipole resonances (red band; Trippa et al. 2008; Lattimer & Lim 2013). Moreover, there is very good agreement with an estimate obtained from modeling X-ray bursts and quiescent low-mass X-ray binaries (shaded region, labeled "Astrophysics"; Steiner et al. 2010). Remarkably, the constraints from these analyses have a common area of intersection (white area), which overlaps within uncertainties with the constraints from microscopic calculations of neutron matter. This suggests that quartic and high-order corrections are relatively small. In the future, we shall study them on the basis of chiral EFT. Based on the comparison in Figure 2, there is also very good agreement with the S_v–L correlation band obtained from the Quantum Monte Carlo results of Gandolfi et al. (2012; see Lattimer & Lim 2013).

The transition between the neutron star crust and a uniform state in the core is a first-order phase transition. However, the density jump across the boundary is small; consequently, a good estimate for the density at which the transition from crustal matter to uniform matter in the core takes place may be obtained by determining the conditions under which, on lowering the density, matter becomes unstable to formation of a small density modulation (Baym et al. 1971a; Pethick et al. 1995). The density at which this occurs provides a lower bound on the density of the uniform phase at which the transition occurs.

The electron screening length in the crust is large compared with typical nuclear separations, so it is a good first approximation to regard the electron density as remaining uniform when a proton density wave is formed. For matter to be stable to formation of a long wavelength density fluctuation, the energy density must increase when the density modulations are imposed subject to the condition that the total number of neutrons and the total number of protons remain constant. Our discussion is adapted from Baym et al. (1971a) with minor variations. One condition for stability when Coulomb and contributions to the energy from density gradients are included is

$$\begin{equation} v_0 + 2 (4\pi e^2\beta)^{1/2} - \beta k_{\rm FT}^2 > 0 \,. \end{equation} \tag{ 8 }$$

The quantity

$$\begin{equation} v_0 = \frac{\partial \mu _p}{\partial n_p} - \frac{(\partial \mu _p/ \partial n_n)^2}{(\partial \mu _n/\partial n_n)} \end{equation} \tag{ 9 }$$

is an effective proton–proton interaction when the Coulomb interaction is neglected. The first term is the contribution when the neutron density is constant and the second term is an induced interaction due to exchange of neutron density fluctuations. The effect of inhomogeneities in the density distribution is described by the quantity β, given by

$$\begin{equation} \beta =2 (Q_{pp}+2Q_{np}\zeta +Q_{nn}\zeta ^2) , \end{equation} \tag{ 10 }$$

and a more detailed discussion is given in the Appendix. Here Q_ij = B_ij/n₀ in the notation of Baym et al. (1971a) and

$$\begin{equation} \zeta =-\frac{\partial \mu _p/\partial n_n}{\partial \mu _n/\partial n_n} . \end{equation} \tag{ 11 }$$

The Thomas–Fermi wave number k_TF is given by

$$\begin{equation} k_{\rm TF}^2 = \frac{4}{\pi } \frac{e^2}{\hbar c} \, k_e^2 , \end{equation} \tag{ 12 }$$

with k_e = (3π²nx)^1/3. Generally, the second term in Equation (8), which is due to the nonuniformity of the electron density, is small compared with the first but we do not drop it. If Equation (8) and the condition ∂μ_n/∂n_n > 0 are satisfied, matter is stable to small density modulations. With decreasing density, for realistic equations of state, the first stability condition to be violated is Equation (8), and instability first sets in when this becomes an equality, which corresponds to Equation (9.18) of Baym et al. (1971a).

Figure 5 shows contour plots of the densities and pressures at which instability sets in for uniform matter in beta equilibrium, with Coulomb and density gradient contributions for Q = 75 MeV fm⁵ (right panel), and without these contributions, Q = 0 (left panel); see the Appendix. We find transition densities around $\bar{n} = 0.55\hbox{--}0.625$, and inclusion of the Coulomb and density gradient contributions lowers these by about 15% to $\bar{n} = 0.475\hbox{--}0.55$. These transition densities are somewhat smaller than those predicted by the FPS and FPS21 interactions in Pethick et al. (1995), and consequently neutron star models with the interactions used in this paper will have lower crustal masses and lower crustal moments of inertia.

Zoom In Zoom Out Reset image size

We investigate the structure of nonrotating neutron stars by solving the Tolman–Oppenheimer–Volkov equations. Since the central densities of neutron stars can significantly exceed the regime of our neutron matter calculations, we need to extend the EOS beyond this density. To this end, we employ a general polytropic extension, where the pressure of neutron star matter is piecewise given by P(ρ) = κ ρ^Γ, with mass density ρ = mn (Read et al. 2009; Hebeler et al. 2010). This is illustrated in Figure 6. At low densities we use the BPS crust EOS (Baym et al. 1971b; Negele & Vautherin 1973) up to the transition density ρ_crust. Figure 6 shows that for densities ρ₀/10 to ρ₀/2 the BPS crust EOS lies within the band predicted for the pressure of neutron star matter based on the neutron matter results and incorporating beta equilibrium following Section 3 (see also Hebeler et al. 2010). Therefore, our results are insensitive to the particular choice of ρ_crust within this region. In the following, we use ρ_crust = ρ₀/2, based on our results for the crust-core boundary given by Figure 5. The pressure from ρ_crust to ρ₁ = 1.1 ρ₀ is given by the band predicted by chiral EFT interactions. Our results are insensitive to the precise value of ρ₁ in the vicinity of saturation density ρ₀, so we have taken a conservative value for which the uncertainty band of the microscopic neutron matter calculations is reasonable.

Zoom In Zoom Out Reset image size

For the extension beyond ρ₁, we use three polytropes with exponents Γ₁, Γ₂ and Γ₃, which make it possible to vary the softness or stiffness of the EOS in the density regions 1: ρ₁ ⩽ ρ ⩽ ρ₁₂; 2: ρ₁₂ ⩽ ρ ⩽ ρ₂₃, and 3: ρ ⩾ ρ₂₃, respectively (see Figure 6). For densities just above ρ₁, the EOS is still rather well constrained by the neutron matter calculations. The band predicted for the pressure of neutron star matter at ρ₁ corresponds to values of Γ in the range 2.25–2.5. Therefore, we take a restricted range for the first polytropic exponent 1.0 ⩽ Γ₁ ⩽ 4.5. We vary the value of all Γ_i in steps of 0.5. At intermediate densities we allow for the possibility of a phase transition and take a broad range 0 ⩽ Γ₂ ⩽ 8. Finally, for densities beyond ρ₂₃ we allow for 0.5 ⩽ Γ₃ ⩽ 8. We exclude the value 0 for this density region in order to avoid artifacts connected with a phase transition up to arbitrarily high density. For the densities between polytropes, ρ₁₂ and ρ₂₃, we take 1.5 ρ₀ ⩽ ρ₁₂ < ρ₂₃ < ρ_max with a step size of ρ₀/2. We will show in the next section that the maximal density for our suite of EOSs of neutron stars is ρ_max ≈ 8.3 ρ₀.

The general polytropic extension leads to a very large number of EOSs, which cover all possible pressures in the gray region in Figure 7. Furthermore, we emphasize that this strategy is very general. It is based on a well defined uncertainty band at nuclear densities and does not rely on assumptions about the nature of the constituents of neutron star matter and their interactions at higher densities. The values of Γ_i and ρ_ij are limited by nuclear physics and observation. In the following, we will demonstrate that the recent observation of a 1.97 ± 0.04 M_☉ neutron star (Demorest et al. 2010) puts rather tight constraints on the EOS at high densities and the radii of neutron stars. We note that our results agree with those of our first study (Hebeler et al. 2010), which used only two polytropes. This shows that the general extension is robust and the conclusions would not change significantly if additional polytropes were introduced to characterize the pressure at high densities.

Zoom In Zoom Out Reset image size

The piecewise polytropic extension described in the previous section is used to generate a very large number of equations of state that cover the pressure-density plane at higher densities. We solve the Tolman–Oppenheimer–Volkov equations for each of these EOSs and retain only those that fulfill the following two constraints.

• 1.  
  The speed of sound remains smaller than the speed of light for all densities relevant in neutron stars: $v_s(\rho) = \sqrt{{\it dP}/d{\mathcal {E}}} \, c \leqslant c$, where ${\mathcal {E}}$ is the energy density.
• 2.  
  The EOS supports a neutron star mass $M = \widehat{M}$, the mass of the heaviest neutron star observed.

We consider each EOS up to densities at which the maximal neutron star mass is reached or the EOS becomes acausal, whichever density is smaller. In Figure 7 we present the individual EOSs that fulfill both constraints for two cases: $\widehat{M} = 1.97 \, M_\odot$ (left panel), the mass of the heaviest known neutron star (Demorest et al. 2010), and $\widehat{M} = 2.4 \, M_\odot$ (right panel), an estimated mass of the black widow pulsar B1957+20 (van Kerkwijk et al. 2011). However, since the uncertainties of the latter determination are large, the $\widehat{M} = 2.4 \, M_\odot$ constraint should be considered as a hypothetical mass, which is used here to probe the sensitivity of our results to the constraint from observations.

The blue band at lower densities in Figure 7 represents the pressure predicted for matter in beta equilibrium as described in Sections 2 and 3. The lighter blue bands give the allowed EOS range, which is the envelope of the allowed polytropes at higher densities. The color of the individual lines indicates the maximal central density of the individual EOS (see figure caption). Clearly, the pressure accessed by the possible polytropic extensions (the gray area in Figure 7) is substantially reduced by causality and by the constraint from neutron star mass measurements. The higher the mass of the heaviest neutron star observed, the stronger the EOS band is constrained. In addition, we find that the maximal densities in neutron stars are limited: stiff EOSs with large polytropic exponents have smaller maximal densities (see yellow lines), which are strongly constrained by causality. Softer EOSs tend to have larger central densities. For M = 1.97 M_☉, we find central densities as high as ≈8.3 ρ₀, and for M = 2.4 M_☉, the densities reach only ≈5.8 ρ₀.

In Figure 8 we compare the EOS bands of Figure 7 with a representative set of EOSs used in the literature. This set contains EOSs calculated within different theoretical approaches and based on different degrees of freedom. For details and notation, we refer the reader to Lattimer & Prakash (2001). We find that a significant number of EOSs are not compatible with the lower density band based on chiral EFT interactions. In addition, at higher densities only very few EOSs, including the variational EOSs based on phenomenological nuclear potentials (Akmal et al. 1998), AP3 and AP4 in Figure 8, are located within the uncertainty bands over the entire density range.

Zoom In Zoom Out Reset image size

Figure 9 shows the uncertainty bands for the pressure as a function of energy density ${\mathcal {E}}$. These are the natural variables to study to what extent the causality constraint is responsible for the apparent softening of the EOSs at high densities. For comparison, we show the limit $P = {\mathcal {E}}$, represented by the dotted line. Furthermore, we generated a causal extension for each individual EOS for energy densities ${\mathcal {E}} > {\mathcal {E}}_{\rm limit}$ by choosing $P({\mathcal {E}}) = {\mathcal {E}} - {\mathcal {E}}_{\rm limit} + P({\mathcal {E}}_{\rm limit})$. This ensures that the energy density, pressure, and speed of sound are continuous at all densities, with speed of sound $v_s({\mathcal {E}}) = c$ for ${\mathcal {E}} > {\mathcal {E}}_{\rm limit}$. This can generate an EOS that has the speed of sound equal to the speed of light immediately after a phase transition at the start of the third polytrope. We exclude such a scenario. In this way we can extend all EOSs to arbitrarily high densities and probe the role of the causality constraint on the results. In Figure 9, the solid lines mark the allowed EOS range using the causal extension for $\widehat{M} = 1.97 \, M_{\odot }$, and the dashed lines mark the range for $\widehat{M} = 2.4 \, M_{\odot }$. The comparison to the blue bands shows that the causal extension changes only slightly the upper pressure limit and leads to somewhat higher maximal possible densities in a neutron star. We also observe that the causal extensions have a stronger impact on the lower limit of the pressure uncertainty band. For EOSs in this region, the speed of sound reaches the speed of light already for small neutron star masses. By employing causal extensions, the maximal neutron star mass for such EOSs can increase significantly and consequently more EOSs fulfill the mass constraint, which explains the significant difference between the blue bands and the regions within the black lines.

Zoom In Zoom Out Reset image size

Figure 10 shows the central densities of neutron stars with masses 1.4 M_☉, 1.97 M_☉, and 2.4 M_☉ based on the uncertainty bands of Figure 7. The results for the minimal and maximal central densities are given in Table 3. Since stiff EOSs along the upper limit of the pressure uncertainty band are not sensitive to the mass constraint, the minimal central densities are identical for the $\widehat{M}= 1.97 \, M_\odot$ and $\widehat{M} = 2.4 \, M_\odot$ cases.

Zoom In Zoom Out Reset image size

Table 3. Minimal and Maximal Central Densities ρ_c (in Units of the Saturation Density ρ₀) of the Neutron Stars Shown by the Points in Figure 10

	$\widehat{M} = 1.97 \, M_{\odot }$		$\widehat{M} = 2.4 \, M_{\odot }$
	Min	Max	Min	Max
ρc/ρ0 (1.4 M☉)	1.8	4.4	1.8	2.7
ρc/ρ0 (1.97 M☉)	2.0	7.6	2.0	3.4
ρc/ρ0 (2.4 M☉)			2.2	5.4

Download table as:  ASCIITypeset image

The uncertainty bands for the EOS directly translate into constraints for the radii of neutron stars. In Figure 11, we present the radius constraints obtained from the EOS bands of Figure 7, for the two cases $\widehat{M} = 1.97 \, M_{\odot }$ (left panel) and $\widehat{M} = 2.4 \, M_{\odot }$ (right panel). The mass–radius relationships for the individual EOSs were obtained by solving the Tolman–Oppenheimer–Volkov equations and from these an envelope was constructed. The blue regions in Figure 11 show the radius constraints based on the neutron matter results with renormalization-group-evolved interactions. In addition, we show results without renormalization-group evolution (see Figure 1) and for the allowed EOS range with causal extension (see Figure 9). The results with causal extension are also based on renormalization-group-evolved interactions.

Zoom In Zoom Out Reset image size

The predicted radius ranges are given in Table 4 for a 1.4 M_☉ neutron star and for the heaviest known neutron star, with M = 1.97 M_☉. The maximal radius is very robust and essentially independent of the details of the neutron matter calculation, the use of causal extensions, and the mass constraints. This can be understood intuitively. Very stiff EOSs, which determine the maximal radius constraint, lead also to large neutron star masses. Hence, the constraints $\widehat{M} = 1.97 \, M_{\odot }$ and $\widehat{M} = 2.4 \, M_{\odot }$ are always fulfilled. Furthermore, the central densities of such neutron stars are typically rather low (see Figure 10). For the stiffest EOSs we find the central density to be ρ_c ≈ 1.8 ρ₀ for an M = 1.4 M_☉ neutron star, which is typical for many observed neutron stars (see Table 3), a density region which is still rather well constrained by chiral EFT interactions (Hebeler & Schwenk 2010; Weise 2012). A 1.4 M_☉ neutron star with a radius significantly larger than R = 13.9 km would therefore be incompatible with constraints from chiral EFT interactions.

Table 4. Radius Constraints for a 1.4 M_☉ Neutron Star and for the Heaviest Known Neutron Star with M = 1.97 M_☉ Based on the Results of Figure 11 (See the Text for Details on the Column Labels)

	$\widehat{M} = 1.97 \, M_{\odot }$			$\widehat{M} = 2.4 \, M_{\odot }$
	Evolved	Unevolved	Causal	Evolved	Unevolved	Causal
Rmin (1.4 M☉)	10.0	9.7	8.8	11.6	11.5	10.4
Rmax (1.4 M☉)	13.7	13.9	13.7	13.7	13.9	13.7
Rmin (1.97 M☉)	9.6	9.3	8.6	12.0	11.8	10.9
Rmax (1.97 M☉)	14.2	14.4	14.4	14.2	14.4	14.4

Note. The rows give the minimum and maximum radii.

Download table as:  ASCIITypeset image

The minimal radius is more sensitive to details of the EOS at higher densities and therefore less well constrained. The limits based on evolved and unevolved nuclear interactions agree well, except for very light neutron stars which are more sensitive to small differences in the low-density part of the EOS. In addition, the minimal radius strongly depends on the mass constraint. This implies that the lower limit of theoretical neutron star radii can be systematically improved by future observations of heavier neutron stars. Since the maximal mass of a neutron star based on soft EOSs can be increased by causal extensions (see the discussion of Figure 9), we find a reduction of the minimal radius for EOSs by about 1 km with causal extensions, as shown by the difference between the blue bands and the dotted lines. This is consistent with the lower radius limits of Lattimer (2012), giving R_min(1.4 M_☉) = 8.1 km and 9.1 km for $\widehat{M} = 1.97 \, M_{\odot }$ and $\widehat{M} = 2.4 \, M_{\odot }$, respectively.

The uncertainty bands for the EOS are an envelope of many individual EOSs. Therefore, the limits of the uncertainty bands are not given by one EOS in general (see Figure 7). In order to explore the EOS bands in astrophysical applications, we present three representative EOSs that probe the central region and the extremes of the uncertainty band for $\widehat{M} = 1.97 \, M_{\odot }$.

• 1.  
  A soft EOS (green dashed line in Figure 12), which agrees well with the lower limit of the EOS band for low and medium densities. This EOS excellently describes the minimal radius over the entire mass range (see the right panel of Figure 12). At higher densities, the pressure then reaches values above the lower limit of the band, as a consequence of the mass constraint. The polytropic parameters of this EOS are Γ₁ = 1.5, ρ₁₂ = 2.5 ρ₀, Γ₂ = 6.0, ρ₂₃ = 4.0 ρ₀, Γ₃ = 3.0, and ρ_max ≈ 7.0 ρ₀.
• 2.  
  An intermediate EOS (orange solid line in Figure 12), which is located in the central region of the EOS band and is also consistent with the $\widehat{M} = 2.4 \, M_{\odot }$ EOS band over the whole density range. The radius of neutron stars for this EOS is about 12 km for typical masses, and lies in the center of the predicted radius range. The polytropic parameters of this EOS are Γ₁ = 4.0, ρ₁₂ = 3.0 ρ₀, Γ₂ = 3.0, ρ₂₃ = 4.5 ρ₀, Γ₃ = 2.5, and ρ_max ≈ 5.4 ρ₀.
• 3.  
  A stiff EOS (red dash-dotted line in Figure 12), which closely follows the upper limit of the EOS band until the causal limit is reached. This EOS gives the largest neutron star radii consistent with the constraints. The polytropic parameters of this EOS are Γ₁ = 4.5, ρ₁₂ = 1.5 ρ₀, Γ₂ = 5.5, ρ₂₃ = 2.0 ρ₀, Γ₃ = 3.0, and ρ_max = 3.3 ρ₀.

Zoom In Zoom Out Reset image size

Detailed values of the pressure, energy density, and energy per particle for these three EOSs as a function of number density and mass density are listed in Tables 5 and 6, with the BPS crust EOS used at low densities in Table 7. In addition, Tables 5 and 6 include values for the neutron star radius and mass at that central density.

Table 5. Numerical Data for the Three Representative EOSs of Section 7 as a Function of Number Density n/n₀ or Mass Density ρ

n/n0	ρc2	Soft					Intermediate					Stiff
		P	${\mathcal {E}}$		R	M	P	${\mathcal {E}}$		R	M	P	${\mathcal {E}}$		R	M
0.5792	87.07	0.4470	87.90	8.920	3260	0.53	0.4470	87.90	8.920	3260	0.53	0.6960	87.99	9.937	77.04	0.10
0.7124	107.1	0.7162	108.2	10.06	248.8	0.10	0.7162	108.2	10.06	248.8	0.10	1.150	108.4	11.75	29.32	0.13
0.7861	118.1	0.9094	119.5	10.73	85.17	0.10	0.9094	119.5	10.73	85.17	0.10	1.473	119.7	12.82	23.22	0.15
0.8646	129.9	1.154	131.5	11.47	47.29	0.10	1.154	131.5	11.47	47.29	0.10	1.880	131.9	14.03	19.72	0.18
0.9483	142.5	1.464	144.4	12.31	32.05	0.12	1.464	144.4	12.31	32.05	0.12	2.392	144.8	15.38	17.50	0.20
1.0371	155.9	1.851	158.0	13.24	24.36	0.13	1.851	158.0	13.24	24.36	0.13	3.028	158.7	16.91	16.02	0.23
1.1	165.3	2.163	167.8	13.94	21.15	0.14	2.163	167.8	13.94	21.15	0.14	3.542	168.5	18.04	15.29	0.26
1.2	180.4	2.465	183.3	14.99	19.26	0.15	3.064	183.3	15.12	16.57	0.17	5.240	184.2	20.05	13.98	0.34
1.3	195.4	2.780	198.8	16.04	18.00	0.16	4.220	198.9	16.57	14.23	0.22	7.512	200.1	22.57	13.33	0.45
1.4	210.5	3.106	214.3	17.05	17.08	0.17	5.677	214.6	18.25	13.01	0.28	10.48	216.2	25.63	13.09	0.58
1.5	225.5	3.445	229.8	18.03	16.36	0.17	7.481	230.3	20.19	12.37	0.34	14.30	232.5	29.28	13.07	0.73
1.6	240.5	3.795	245.4	18.97	15.80	0.18	9.684	246.3	22.41	12.04	0.42	20.39	249.2	33.74	13.22	0.96
1.7	255.6	4.157	261.0	19.88	15.32	0.18	12.34	262.3	24.93	11.89	0.51	28.47	266.3	39.30	13.47	1.22
1.8	270.6	4.529	276.6	20.77	14.92	0.19	15.51	278.6	27.76	11.85	0.61	38.98	283.9	46.12	13.75	1.50
1.9	285.6	4.911	292.2	21.63	14.57	0.19	19.25	295.0	30.92	11.87	0.71	52.49	302.2	54.41	14.01	1.79
2.0	300.7	5.304	307.9	22.47	14.26	0.20	23.64	311.7	34.43	11.93	0.82	69.59	321.3	64.37	14.23	2.07
2.1	315.7	5.707	323.5	23.28	13.99	0.20	28.73	328.6	38.31	12.00	0.94	80.56	341.1	75.52	14.32	2.21
2.2	330.7	6.119	339.2	24.08	13.74	0.21	34.61	345.7	42.58	12.09	1.05	92.63	361.4	87.21	14.38	2.33
2.3	345.8	6.541	354.9	24.87	13.52	0.21	41.35	363.2	47.26	12.17	1.17	105.8	382.4	99.44	14.41	2.44
2.4	360.8	6.972	370.6	25.63	13.31	0.21	49.02	380.9	52.36	12.25	1.29	120.2	403.9	112.2	14.42	2.53
2.5	375.8	7.413	386.4	26.38	13.12	0.22	57.72	399.0	57.91	12.32	1.41	135.9	426.0	125.5	14.41	2.61
2.6	390.9	9.379	402.2	27.18	12.33	0.23	67.52	417.4	63.91	12.37	1.53	152.9	448.9	139.4	14.38	2.68
2.7	405.9	11.76	418.0	28.12	11.52	0.26	78.53	436.3	70.40	12.42	1.64	171.2	472.3	153.8	14.33	2.74
2.8	420.9	14.63	434.0	29.20	10.83	0.29	90.82	455.6	77.38	12.45	1.74	190.9	496.5	168.7	14.28	2.79
2.9	436.0	18.06	450.1	30.46	10.31	0.33	104.5	475.3	84.89	12.48	1.84	212.1	521.5	184.2	14.22	2.83
3.0	451.0	22.13	466.3	31.90	9.95	0.39	119.6	495.6	92.92	12.49	1.94	234.8	547.1	200.3	14.15	2.86
3.1	466.0	26.94	482.7	33.54	9.73	0.45	132.0	516.3	101.3	12.49	2.00	259.1	573.6	216.8	14.07	2.89
3.2	481.1	32.60	499.2	35.41	9.60	0.52	145.2	537.4	110.1	12.48	2.06	285.0	600.9	234.0	14.00	2.92
3.3	496.1	39.21	515.9	37.53	9.54	0.61	159.3	559.0	119.1	12.47	2.12	312.6	629.0	251.6	13.91	2.94

Notes. The mass density ρ times c², pressure P, and energy density ${\mathcal {E}}$ are given in MeV fm⁻³. The energy per nucleon is given in MeV. We also list the neutron star radius R in km and mass M in units of M_☉ at the central density ρ.

Download table as:  ASCIITypeset image

Table 6. Same as Table 5

n/n0	ρc2	Soft					Intermediate
		P	${\mathcal {E}}$		R	M	P	${\mathcal {E}}$		R	M
3.4	511.1	46.90	532.8	39.92	9.54	0.69	174.2	581.0	128.3	12.45	2.17
3.5	526.2	55.81	550.0	42.61	9.57	0.79	190.0	603.4	137.9	12.42	2.21
3.6	541.2	66.09	567.5	45.62	9.63	0.89	206.8	626.3	147.7	12.39	2.25
3.7	556.2	77.90	585.2	48.99	9.70	1.00	224.5	649.7	157.9	12.35	2.28
3.8	571.3	91.41	603.3	52.74	9.77	1.11	243.2	673.6	168.2	12.31	2.32
3.9	586.3	106.8	621.8	56.92	9.85	1.21	262.9	698.0	178.9	12.26	2.35
4.0	601.3	124.3	640.7	61.54	9.93	1.32	283.7	722.9	189.9	12.22	2.37
4.1	616.4	133.9	659.9	66.46	9.96	1.38	305.5	748.3	201.1	12.17	2.39
4.2	631.4	143.9	679.4	71.50	9.99	1.43	328.4	774.3	212.6	12.12	2.41
4.3	646.4	154.5	699.2	76.66	10.01	1.47	352.4	800.8	224.3	12.06	2.43
4.4	661.5	165.5	719.1	81.94	10.03	1.55	377.6	827.9	236.4	12.01	2.44
4.5	676.5	177.0	739.4	87.33	10.05	1.56	403.9	855.6	248.7	11.96	2.46
4.6	691.5	189.1	759.9	92.87	10.06	1.59	426.7	883.8	261.3	11.91	2.47
4.7	706.6	201.7	780.6	98.52	10.06	1.63	450.3	912.6	273.9	11.86	2.47
4.8	721.6	214.9	801.7	104.3	10.07	1.66	474.7	941.8	286.7	11.82	2.48
4.9	736.6	228.6	823.0	110.2	10.06	1.69	499.8	971.6	299.7	11.77	2.49
5.0	751.7	242.8	844.6	116.2	10.06	1.72	525.6	1002	312.8	11.72	2.49
5.1	766.7	257.7	866.5	122.3	10.05	1.75	552.3	1033	326.0	11.67	2.49
5.2	781.7	273.2	888.7	128.6	10.05	1.77	579.8	1064	339.3	11.63	2.49
5.3	796.8	289.2	911.2	135.0	10.03	1.80	608.1	1096	352.8	11.58	2.50
5.4	811.8	305.9	934.0	141.5	10.02	1.82	637.2	1128	366.4	11.53	2.50
5.5	826.8	323.2	957.1	148.1	10.00	1.84
5.6	841.9	341.2	980.6	154.8	9.99	1.86
5.7	856.9	359.8	1004	161.7	9.97	1.88
5.8	871.9	379.1	1028	168.7	9.95	1.89
5.9	886.9	399.0	1052	175.8	9.93	1.91
6.0	902.0	419.7	1077	183.0	9.90	1.92
6.1	917.0	441.1	1102	190.3	9.88	1.94
6.2	932.1	463.1	1128	197.8	9.85	1.95
6.3	947.1	485.8	1154	205.4	9.83	1.96
6.4	962.1	509.3	1180	213.1	9.80	1.97
6.5	977.1	533.6	1206	220.9	9.77	1.98
6.6	992.2	558.6	1233	228.9	9.75	1.99
6.7	1007	584.4	1261	237.0	9.72	1.99
6.8	1022	611.0	1289	245.2	9.69	2.00
6.9	1037	638.4	1317	253.5	9.66	2.01
7.0	1052	666.5	1345	261.9	9.63	2.01

Note. The maximal density for the stiff EOS is reached already in Table 5.

Download table as:  ASCIITypeset image

Table 7. Numerical Data of the BPS EOS for the Outer Crust (Baym et al. 1971b; Negele & Vautherin 1973)

n/n0	ρc2	P	${\mathcal {E}}$		n/n0	ρc2	P	${\mathcal {E}}$
2.956 × 10−14	4.444 × 10−12	6.303 × 10−25	4.385 × 10−12		1.972 × 10−4	2.964 × 10−2	3.713 × 10−5	2.931 × 10−2
2.975 × 10−14	4.472 × 10−12	6.303 × 10−24	4.407 × 10−12		2.482 × 10−4	3.732 × 10−2	5.048 × 10−5	3.692 × 10−2
3.069 × 10−14	4.613 × 10−12	6.303 × 10−23	4.547 × 10−12		3.125 × 10−4	4.698 × 10−2	6.865 × 10−5	4.649 × 10−2
4.369 × 10−14	6.568 × 10−12	7.551 × 10−22	6.472 × 10−12		3.934 × 10−4	5.914 × 10−2	9.330 × 10−5	5.852 × 10−2
6.187 × 10−14	9.302 × 10−12	8.737 × 10−21	9.150 × 10−12		4.953 × 10−4	7.445 × 10−2	1.269 × 10−4	7.376 × 10−2
1.700 × 10−13	2.556 × 10−11	1.061 × 10−19	2.516 × 10−11		6.235 × 10−4	9.373 × 10−2	1.621 × 10−4	9.284 × 10−2
7.937 × 10−13	1.193 × 10−10	3.632 × 10−18	1.183 × 10−10		6.906 × 10−4	1.038 × 10−1	1.805 × 10−4	1.029 × 10−1
4.331 × 10−12	6.511 × 10−10	1.186 × 10−16	6.416 × 10−10		7.850 × 10−4	1.180 × 10−1	2.053 × 10−4	1.169 × 10−1
3.934 × 10−11	5.915 × 10−9	6.081 × 10−15	5.825 × 10−9		9.881 × 10−4	1.485 × 10−1	2.791 × 10−4	1.473 × 10−1
9.881 × 10−11	1.485 × 10−8	3.100 × 10−14	1.463 × 10−8		1.244 × 10−3	1.870 × 10−1	3.630 × 10−4	1.855 × 10−1
2.483 × 10−10	3.732 × 10−8	1.517 × 10−13	3.675 × 10−8		1.608 × 10−3	2.417 × 10−1	4.871 × 10−4	2.398 × 10−1
6.235 × 10−10	9.373 × 10−8	7.183 × 10−13	9.228 × 10−8		1.669 × 10−3	2.509 × 10−1	4.924 × 10−4	2.488 × 10−1
1.566 × 10−9	2.355 × 10−7	3.286 × 10−12	2.319 × 10−7		1.954 × 10−3	2.937 × 10−1	5.212 × 10−4	2.917 × 10−1
3.934 × 10−9	5.914 × 10−7	1.447 × 10−11	5.825 × 10−7		2.469 × 10−3	3.712 × 10−1	5.678 × 10−4	3.688 × 10−1
9.881 × 10−9	1.485 × 10−6	6.088 × 10−11	1.463 × 10−6		2.974 × 10−3	4.471 × 10−1	6.135 × 10−4	4.443 × 10−1
2.483 × 10−8	3.732 × 10−6	2.441 × 10−10	3.676 × 10−6		3.633 × 10−3	5.461 × 10−1	6.759 × 10−4	5.427 × 10−1
3.125 × 10−8	4.698 × 10−6	3.282 × 10−10	4.627 × 10−6		4.464 × 10−3	6.711 × 10−1	7.601 × 10−4	6.673 × 10−1
6.235 × 10−8	9.373 × 10−6	8.955 × 10−10	9.233 × 10−6		5.491 × 10−3	8.255 × 10−1	8.731 × 10−4	8.207 × 10−1
1.244 × 10−7	1.870 × 10−5	2.392 × 10−9	1.842 × 10−5		6.250 × 10−3	9.396 × 10−1	5.944 × 10−4	9.390 × 10−1
2.482 × 10−7	3.732 × 10−5	6.278 × 10−9	3.676 × 10−5		2.500 × 10−2	3.758	5.799 × 10−3	3.762
4.953 × 10−7	7.445 × 10−5	1.625 × 10−8	7.337 × 10−5		5.000 × 10−2	7.517	1.166 × 10−2	7.530
9.881 × 10−7	1.485 × 10−4	4.166 × 10−8	1.464 × 10−4		7.500 × 10−2	11.27	2.085 × 10−2	11.30
1.244 × 10−6	1.870 × 10−4	5.453 × 10−8	1.843 × 10−4		1.000 × 10−1	15.03	3.216 × 10−2	15.08
1.972 × 10−6	2.964 × 10−4	1.017 × 10−7	2.922 × 10−4		1.250 × 10−1	18.79	4.515 × 10−2	18.86
3.125 × 10−6	4.698 × 10−4	1.890 × 10−7	4.631 × 10−4		1.500 × 10−1	22.55	5.961 × 10−2	22.64
3.934 × 10−6	5.914 × 10−4	2.577 × 10−7	5.830 × 10−4		1.750 × 10−1	26.31	7.544 × 10−2	26.42
4.952 × 10−6	7.445 × 10−4	3.143 × 10−7	7.342 × 10−4		2.000 × 10−1	30.07	9.260 × 10−2	30.21
6.235 × 10−6	9.373 × 10−4	4.281 × 10−7	9.245 × 10−4		2.250 × 10−1	33.82	1.110 × 10−1	34.00
9.881 × 10−6	1.485 × 10−3	7.938 × 10−7	1.465 × 10−3		2.500 × 10−1	37.58	1.308 × 10−1	37.79
1.566 × 10−5	2.355 × 10−3	1.470 × 10−6	2.323 × 10−3		2.750 × 10−1	41.34	1.518 × 10−1	41.58
2.482 × 10−5	3.732 × 10−3	2.722 × 10−6	3.683 × 10−3		3.000 × 10−1	45.10	1.742 × 10−1	45.38
3.125 × 10−5	4.698 × 10−3	3.533 × 10−6	4.637 × 10−3		3.250 × 10−1	48.86	1.980 × 10−1	49.18
3.934 × 10−5	5.914 × 10−3	4.807 × 10−6	5.836 × 10−3		3.500 × 10−1	52.62	2.231 × 10−1	52.98
4.952 × 10−5	7.445 × 10−3	6.540 × 10−6	7.353 × 10−3		3.750 × 10−1	56.37	2.497 × 10−1	56.78
6.235 × 10−5	9.373 × 10−3	8.893 × 10−6	9.256 × 10−3		4.000 × 10−1	60.13	2.777 × 10−1	60.58
7.850 × 10−5	1.180 × 10−2	1.209 × 10−5	1.166 × 10−2		4.250 × 10−1	63.89	3.073 × 10−1	64.38
9.881 × 10−5	1.485 × 10−2	1.562 × 10−5	1.468 × 10−2		4.500 × 10−1	67.65	3.385 × 10−1	68.19
1.244 × 10−4	1.870 × 10−2	2.124 × 10−5	1.848 × 10−2		4.750 × 10−1	71.41	3.713 × 10−1	72.00
1.566 × 10−4	2.355 × 10−2	2.888 × 10−5	2.328 × 10−2		5.000 × 10−1	75.17	4.054 × 10−1	75.81

Note. The units are as in Table 5.

Download table as:  ASCIITypeset image

Recent advances in nuclear theory combined with new astrophysical observations have systematically tightened the constraints on the EOS of neutron-rich matter over a wide range of densities. In this paper, we have presented constraints for the EOS and for the structure of neutron stars based on microscopic neutron matter calculations with chiral EFT interactions. At nuclear saturation density and below, the uncertainties are dominated by uncertainties in 3N forces. In the future, these uncertainties can be reduced by improved constraints on the low-energy couplings entering 3N forces and by systematic improvements in chiral EFT that incorporate higher-order many-body forces (Tews et al. 2013; Krüger et al. 2013) and Δ degrees of freedom explicitly. In addition, reduced uncertainties are expected when renormalization-group-evolved 3N forces (Hebeler 2012; Hebeler & Furnstahl 2013) and nonperturbative many-body calculations are employed.

The uncertainties of the neutron matter calculations directly translate into constraints for the symmetry energy S_v and its density derivative L, which is related to the pressure of neutron matter. On the assumptions that the interaction energy depends quadratically on the neutron excess, this leads to S_v = 29.7–33.5 MeV and L = 32.4–57.0 MeV. These ranges provides very tight constraints and are consistent within uncertainties with different empirical extractions of these parameters (see Figure 4). In addition, the neutron matter calculations predict the crust-core transition density to be in the range (0.475–0.55) n₀ when Coulomb and density gradient corrections are included.

We extended the EOS to higher densities by employing a general piecewise polytropic extrapolation, as illustrated in Figure 6. This led to a very large number of EOSs, which cover a wide pressure range, including the possibility of a soft phase transition region. From the individual EOSs we selected those that remain causal at all densities relevant in neutron stars and are able to support a neutron star with mass 1.97 M_☉, the heaviest known neutron star (Demorest et al. 2010). We also considered a second case where the EOS supports a neutron star with mass 2.4 M_☉. Note that the constraints are in accord with our previous work (Hebeler et al. 2010), where we used only two polytropes. Combined with the microscopic neutron matter calculations, this provides tight constraints on the nuclear EOS at sub- and supranuclear densities, and rules out many model EOSs, as demonstrated by Figure 8.

Based on the allowed EOS band, we predicted ranges for the radii and central densities of neutron stars (see Figures 10 and 11 as well as Tables 3 and 4). For the constraint $\widehat{M} = 1.97 \, M_{\odot }$, the radius of a 1.4 M_☉ neutron star is predicted to be 9.7–13.9 km (based on renormalization-group-evolved or unevolved chiral EFT interactions) with central densities up to 4.4 ρ₀. Note that the maximum radius is determined by the causality constraint, so that the discovery of a heavier neutron star only affects the minimum radius. If a 2.4 M_☉ neutron star were discovered, this would imply a radius range 11.5–13.9 km for a 1.4 M_☉ neutron star, with central densities up to only 2.7 ρ₀.

For use in astrophysical simulations, we have constructed three representative EOSs (soft, intermediate, and stiff), consistent with the constraints from nuclear physics and observation. In addition to the nuclear physics improvements mentioned above and the measurement of heavier neutron stars, information on the radii of neutron stars will significantly tighten the EOS band at high densities. Currently, observational radius limits are often conflicting. Studies of five photospheric radius expansion X-ray bursts by Steiner et al. (2010, 2013) are in agreement with the present analysis, but other analyses (Özel et al. 2010) of two of them (EXO1745−248 and 4U1820−30) indicate radii that are too small to be compatible with our results. On the other hand, a study of the burster 4U1724−307 (Suleimanov et al. 2011) yields a radius too large to be compatible. Using phase modeling of pulsar profiles, Bogdanov (2013) found a compatible 3σ lower limit to the radius of pulsar PSR J0437−4715 of 11.1 km, assuming the 1.56 M_☉ lower limit for its mass found by Verbiest et al. (2008). However, Hambaryan et al. (2011) found R = (10–12)(M/M_☉) km for the pulsar RBS 1223 using the same technique, which is compatible only if M < 1.4 M_☉. Also, using pulse profile modeling, Leahy et al. (2011) found radii for the sources SAX J1808.4−3658, XTE J1807−294, and XTE J1814-338 which were compatible, but that XTE J1824−338 did not, being too large. A recent analysis (Guillot et al. 2013) of five quiescent neutron stars in low-mass X-ray binaries in globular clusters, in which it was assumed that the stars all had the same radius, determined the radius to be $R = 9.1^{+1.3}_{-1.5} \, {\rm km}$ to 90% confidence, which is marginally consistent with our results. Additional radius measurements from the LOFT (Feroci et al. 2012) and NICER (Gendreau et al. 2012) X-ray missions are therefore eagerly anticipated.

We thank T. Krüger and I. Tews for helpful discussions. This work was supported in part by the NSF under grant No. PHY-1002478, the U.S. Department of Energy under grants DE-FG02-87ER40317 and DE-SC0008533 (SciDAC-3 NUCLEI project), the ESF AstroSim and CompStar networks, the Helmholtz Alliance Program of the Helmholtz Association, contract HA216/EMMI "Extremes of Density and Temperature: Cosmic Matter in the Laboratory," the DFG through Grant SFB 634, and by the ERC grant No. 307986 STRONGINT.

Here we discuss details of the density-gradient contributions to the instability condition and also express the instability condition in terms of the variables n and x, as was done in Lattimer & Prakash (2007), rather than n_n and n_p. The contribution of density gradients to the energy density is expressed in the form (Baym et al. 1971a)

$$\begin{equation} {\mathcal {E}}_{\rm grad} = Q_{pp}(\nabla n_p)^2 + 2Q_{pn}(\nabla n_p)\cdot (\nabla n_n) + Q_{nn}(\nabla n_n)^2 , \end{equation} \tag{ A1 }$$

where the coefficients Q_ij are in general functions of the neutron and proton densities, but will here be treated as constants. We follow Baym et al. (1971a) and take Q_np = Q_pn = 2Q_nn = 2Q_pp, in which case the quantity β of Section 4 is given by

$$\begin{equation} \beta = 2 Q_{nn} (1+4\zeta +\zeta ^2) . \end{equation} \tag{ A2 }$$

Indeed, realistic Skyrme density functionals such as Ska and SkM* satisfy this approximate relationship among the Q_ij values (Lattimer & Prakash 2007). It would be interesting to determine the Q_ij from modern energy density functionals (Erler et al. 2012a, 2012b).

The quantity Q may be determined either from the surface energy of symmetric nuclear matter or from the surface thickness of symmetric nuclei. For symmetric nuclear matter, the energy density is given by

$$\begin{equation} {\mathcal {E}}_{\rm bulk} + {\mathcal {E}}_{\rm grad} = n (\epsilon _{\rm s}(n) + m c^2) + Q(\nabla n)^2 , \end{equation} \tag{ A3 }$$

where $n [\epsilon _{\rm s} (n) + m c^2] = {\mathcal {E}}(n,x=1/2)$, Q = 3Q_nn/2 = 3Q_np/4 and m = (m_n + m_p)/2. In the Thomas–Fermi approximation and for symmetric nuclear matter the optimum density profile for a plane surface is determined by minimizing the quantity ∫dz [n(_s(n) − μ₀) + Q(n')²] for a surface lying in the xy-plane (Bethe 1968). Here μ₀ = −B is the chemical potential of symmetric nuclear matter at the saturation density and the prime denotes a derivative with respect to z. This leads to the condition

$$\begin{equation} n^\prime = -\sqrt{\frac{n[\epsilon _{\rm s}(n) - \mu _0]}{Q}} = -\sqrt{\frac{n_0T_0}{Q}} \, f(n/n_0) , \end{equation} \tag{ A4 }$$

when the nuclear matter is assumed to occupy the region with negative z. This equation defines f(n/n₀), and, on integration, determines n(z). Note that f(n/n₀) → 0 for both n → 0 and n → n₀. The surface energy per unit area, the surface tension σ₀, is then given by

$$\begin{equation} \sigma _0 = 2 I_\sigma \, \sqrt{QT_0n_0^3} , \end{equation} \tag{ A5 }$$

where

$$\begin{equation} I_\sigma = \int _0^1 f(u) \, du \approx 0.1696 , \end{equation} \tag{ A6 }$$

for the bulk energy per particle given by Equation (2) with γ = 4/3, α = 5.87, and η = 3.81 (see details after Equation (3)). The 90%–10% surface thickness can be expressed as

$$\begin{equation} t_{90\hbox{--}10}=I_t \, \sqrt{\frac{Qn_0}{T_0}} , \end{equation} \tag{ A7 }$$

where

$$\begin{equation} I_t= \int _{0.1}^{0.9}\frac{du}{f(u)} \approx 4.776 . \end{equation} \tag{ A8 }$$

For Q = 75 MeV fm⁵, one finds σ₀ ≈ 1.15 MeV fm⁻² and t_90–10 ≈ 2.7 fm, both reasonable values although perhaps 5% too large for the surface thickness. Note that Hartree–Fock calculations of surface energies can differ by about 10% from the Thomas–Fermi approximation (Ravenhall et al. 1972). However, we have no freedom to fit both observables with the adopted functional for the bulk energy per particle. With this value for Q, we find Q_np = 100 MeV fm⁵ and Q_pp = Q_nn = 50 MeV fm⁵, which are the values used in the calculations in Section 4.

We now comment briefly on the stability conditions when expressed in terms of the variables n and x (Lattimer & Prakash 2007), rather than n_n and n_p. To take care of the constraints on the numbers of neutrons and protons, it is convenient to work with the thermodynamic potential

$$\begin{equation} \Xi = {\mathcal {E}} - \mu _n^0 n_n - \mu _p^0 n_p \,, \end{equation} \tag{ A9 }$$

where $\mu _n^0$ and $\mu _p^0$ are the chemical potentials in the initial state. With this choice, the first order terms in an expansion of Ξ in powers of deviations of the densities from those in the initial state vanish. With n_n = n(1 − x) and n_p = nx, it follows from Equation (A9) that

$$\begin{equation} \delta ^2\Xi = \frac{1}{2} \frac{\partial ^2\cal E}{\partial n^2} \, (\delta n)^2\ + \left(\frac{\partial ^2\cal E}{\partial n\partial x} +\mu _n^0-\mu _p^0 \right) \, \delta n \delta x + \frac{1}{2} \frac{\partial ^2\cal E}{\partial x^2} \, (\delta x)^2 \,, \end{equation} \tag{ A10 }$$

where the term with the chemical potentials comes from the nonlinear dependence of the neutron and proton densities on n and x. Since the baryon pressure is given by $P = n^2 \partial ({\cal E}/n)/\partial n$ and the difference of the neutron and proton chemical potentials by $n \, (\mu _p-\mu _n)=\partial {\cal E}/\partial x$, it follows that

$$\begin{eqnarray} \frac{\partial ^2{\cal E}}{\partial n^2} &=& \frac{1}{n} \frac{\partial P}{\partial n} \,, \quad \frac{\partial ^2{\cal E}}{\partial x\partial n} = \frac{1}{n} \frac{\partial P}{\partial x} - \mu _n+\mu _p \,, \quad {\rm and} \nonumber\\ \frac{1}{2} \frac{\partial ^2\cal E}{\partial x^2} &=& n \frac{\partial (\mu _p-\mu _n)}{\partial x} \,. \end{eqnarray} \tag{ A11 }$$

Therefore, since in the coefficients we may put $\mu _n^0=\mu _n$ and $\mu _p^0=\mu _p$, the quadratic form, Equation (A10), may be rewritten as

$$\begin{equation} \delta ^2\Xi = \frac{1}{2n} \frac{\partial P}{\partial n} \, (\delta n)^2 + \frac{1}{n} \frac{\partial P}{\partial x} \, \delta n \delta x + \frac{n}{2} \frac{\partial (\mu _p-\mu _n)}{\partial x} \, (\delta x)^2 \,. \end{equation} \tag{ A12 }$$

The conditions for the quadratic form to be positive definite are that the diagonal terms be positive,

$$\begin{equation} \frac{\partial P}{\partial n}>0 \quad {\rm and} \quad \frac{\partial (\mu _p-\mu _n)}{\partial x}>0 \,, \end{equation} \tag{ A13 }$$

and that the determinant of the quadratic form be positive (Lattimer & Prakash 2007),

$$\begin{equation} \frac{\partial P}{\partial n} \, \frac{\partial (\mu _p-\mu _n)}{\partial x} - \left(\frac{\partial P}{ \partial x}\right)^2>0 \,. \end{equation} \tag{ A14 }$$

If the condition (A14) and one of the conditions (A13) are satisfied, the other condition (A13) is satisfied automatically. If Equation (A13) is satisfied, the inequality (A14) may be rewritten as

$$\begin{equation} \frac{\partial ^2 {\cal E}}{\partial n^2} - \left(\frac{\partial P}{\partial x}\right)^2 \frac{\partial x}{\partial (\mu _p-\mu _n)}>0 \,. \end{equation} \tag{ A15 }$$

The first term may be regarded as a direct interaction between density fluctuations, without changes in the composition, while the second represents an induced interaction between density fluctuations due to changes in the composition. For the system to be stable, the total interaction consisting of direct and induced contributions must be positive.