CONSTRAINING THE SYMMETRY PARAMETERS OF THE NUCLEAR INTERACTION

The standard liquid drop model (Myers & Swiatecki 1966) holds that the energy of a nucleus is approximately

$$\begin{eqnarray} E(Z,A)&=&A(-B+S_vI^2)+A^{2/3}(E_s-S_sI^2) \cr &&+{3\over 5}{e^2Z^2\over r_0A^{1/3}}+E_{{\rm shell}}(Z,A-Z)+E_{{\rm pairing}}(A), \quad \end{eqnarray} \tag{ 3 }$$

where B ≃ 16 MeV is the bulk binding energy (BE), per baryon, of symmetric matter, and E_s ≃ 19 MeV is the symmetric matter surface energy parameter. I = (A − 2Z)/A is the neutron excess, and r₀ = (4πn₀/3)^−1/3, where n₀ is the nuclear saturation density. The last three terms in Equation (3) represent the Coulomb, shell, and pairing energies, respectively.

The net symmetry energy of an isolated nucleus is then

$$\begin{equation} E_{D,i}=I_i^2 \left(S_vA_i-S_sA_i^{2/3} \right), \end{equation} \tag{ 4 }$$

so a linear correlation between S_v and S_s is expected from minimizing the differences between model predictions and experimentally measured symmetry energies, i.e., minimizing

$$\begin{equation} \chi ^2=\sum _i(E_{{\rm exp},i}-E_{D,i})^2/\left({\cal N}\sigma _D^2 \right), \end{equation} \tag{ 5 }$$

where ${\cal N}$ is the total number of nuclei and σ_D is a nominal error. A χ² contour one unit above the minimum value represents the 1σ confidence interval, which is an ellipse in this linear example.

The properties of the confidence ellipse are determined by the second derivatives of χ² at the minimum,

$$\begin{eqnarray*} [\chi _{vv}, \chi _{vs}, \chi _{ss}]&=&{2\over {\cal N}\sigma _D^2}\sum _iI_i^4 \left[A_i^2, -A_i^{5/3}, A_i^{4/3} \right]\cr &\simeq &[61.6, -10.7, 1.87],\nonumber \end{eqnarray*}$$

where χ_vs = ∂²χ²/∂S_v∂S_s, etc. The specific values follow from the set of 2336 nuclei with N and Z greater than 40 included in the Audi et al. (2003) table. The confidence ellipse in S_s–S_v space has orientation α_D = (1/2)tan ⁻¹|2χ_vs/(χ_vv − χ_ss)| ≃ 98 with respect to the S_s axis, with error widths $\sigma _{v,D}=\sqrt{(\chi ^{-1})_{vv}}\simeq 2.3\sigma _D$ and $\sigma _{s,D}=\sqrt{(\chi ^{-1})_{ss}}\simeq 13.2\sigma _D$, where (χ⁻¹) is the matrix inverse. The correlation coefficient is $r_D=\chi _{vs}/\sqrt{\chi _{vv}\chi _{ss}}\simeq 0.997$. In this simple example, the shape and orientation of the confidence interval depend only on A_i and I_i and not the location of the χ² minimum or other drop parameters.

Although the specific values of the best-fit parameters and the exact correlation are sensitive to the inclusion of terms representing shell, pairing, Coulomb diffusion, exchange and Wigner energies, deformations, and the neutron skin, among others, their effects are relatively minor. For example, Van Isacker & Dieperink (2006) showed that the Wigner term is important in identifying the best-fit values of the liquid drop parameters, but does not significantly affect the deduced correlation between S_v and S_s.

In practice, the liquid droplet model (Myers & Swiatecki 1969), which differs from the liquid drop model by accounting for varying neutron/proton ratios within the nucleus that produce neutron skins in neutron-rich nuclei, is an improved treatment. Its symmetry energy is

$$\begin{equation} E_{{\rm LD},i}=S_vI_i^2A_i \left(1+S_sA_i^{-1/3}/S_v \right)^{-1}, \end{equation} \tag{ 6 }$$

and therefore predicts a linear correlation between S_s/S_v and S_v rather than between S_s and S_v as in the drop model. The same methodology as for the liquid drop model can be used to determine the confidence interval in S_s/S_v–S_v space. In addition to A_i and I_i, the properties of the confidence interval now depend, to a degree, on the measured masses, E_{exp, i}, as well as the parameters of the droplet model.

The correlation between S_s/S_v and S_v emerging from nuclear mass fitting is roughly equivalent to a correlation between L and S_v because the surface symmetry parameter, S_s, depends on the behavior of S(n) within the nuclear interior. The general features can be appreciated using the plane-parallel approximation for the nuclear surface. At zero temperature, for densities less than the nuclear saturation density n_s, bulk matter separates into two phases of differing densities and proton fractions (a dense phase and a dilute phase) having equal pressures and neutron and proton chemical potentials. The surface thermodynamic potential per unit area, σ, also known as the surface tension, is the difference between the total thermodynamic potential and that of two uniform fluids each having the properties of matter far from the interface on both sides. It is straightforward to show that in the case appropriate for finite nuclei, in which the density in the dilute phase (z = +∞) vanishes,

$$\begin{equation} \sigma =2\int _{-\infty }^{+\infty }n[e(n,x)-(1-x)\mu _n -x\mu _p]dz, \end{equation} \tag{ 7 }$$

where μ_n and μ_p are the neutron and proton chemical potentials, respectively. To lowest order in the dense phase neutron excess I = 1 − 2x_−∞, the surface tension can be expanded:

$$\begin{equation} \sigma \simeq \sigma _0-(1-2x_{-\infty })^2\sigma _x, \end{equation} \tag{ 8 }$$

where σ₀ is the symmetric matter surface tension. In this notation, the surface energy parameters of the liquid drop and droplet models are $E_s=4\pi r_0^2\sigma _0$ and $S_s=4\pi r_0^2\sigma _x$.

One can show that to leading order (Lipparini & Stringari 1989; Danielewicz 2003; Steiner et al. 2005)

$$\begin{equation} {S_s\over S_v}={E_s\over 2}{\int _0^{1}u^{1/2}[S_v/S_2(u)-1]f(u)^{-1/2} du\over \int _0^1u^{1/2}f(u)^{1/2}du}, \end{equation} \tag{ 9 }$$

where u = n/n_s and f(u) = e(u, 1/2) + B. A simple analytical relation is found if we approximate S₂(u) ≃ S_v + L(u − 1)/3 and f(u) ≃ K₀(u − 1)²/18, where K₀ ≃ 240 MeV is the incompressibility parameter:

$$\begin{equation} {S_s\over S_v}\simeq {135E_s\over 2K_0} [1-\eta \tan ^{-1}\eta ^{-1}], \end{equation} \tag{ 10 }$$

where $\eta =\sqrt{3S_v/L-1}$. Assuming L/S_v = 1, 3/2, or 2, one finds, respectively, S_s/S_v ≃ 0.7, 1.2, or 1.8, which shows the basic trend although the values are somewhat less than given by realistic calculations.

This illustration is model dependent, partly because it is based on the liquid droplet model and the semi-infinite surface approximation, but mostly because of the highly simplified assumed dependence of S(n) on L and S_v. In actuality, S(n) depends on higher order terms such as K_sym and Q_sym. We have found, applying the finite-range nuclear interaction (Myers & Swiatecki 1990) to Thomas–Fermi semi-infinite surface calculations, that S_s/S_v can be approximated by

$$\begin{equation} {S_s\over S_v}\simeq 0.6461+{S_v\over 97.85\,{\rm MeV}}+ 0.4364{L\over S_v}+0.0873\left({L\over S_v}\right)^2 \end{equation} \tag{ 11 }$$

with a mean error of 0.6%. We found that a very similar relation emerges from Thomas–Fermi calculations using generalized Skyrme functionals. Very approximately, S_s/S_v increases linearly with L, so that one expects the confidence interval in L–S_v space to have a highly elongated elliptical shape with a small angle with respect to the L axis.

Least-squares fitting of a liquid droplet model, together with the Thomas–Fermi result in Equation (11), leads to best-fit values of the symmetry parameters S_{v0, LD} = 29.1 MeV and L_{0, LD} = 53.9 MeV. Shell effects were modeled with linear and quadratic terms in the neutron and proton valence numbers (the smallest distances from their respective magic numbers) as in Van Isacker & Dieperink (2006). Pairing effects were modeled with a term a|A|^1/2 for odd–odd nuclei, −a|A|^1/2 for even–even nuclei, and zero for odd–even nuclei, where a ∼ 12 MeV was found from the fitting procedure. Conventional Wigner and Coulomb exchange energy terms were also employed, as was a further term representing the Coulomb diffuseness energy. The inclusion of these auxiliary terms, while improving the overall fit, which was χ_minσ_LD ≃ 1.26 MeV, did not significantly alter the inferred correlation between S_v and L, although the location of the χ² minimum does depend on them. The results scale with the assumed value of the nominal error, in this case σ_LD. We also found that the properties of the confidence interval were σ_{v, LD} = 6.2σ_LD, σ_{L, LD} = 80.0σ_LD, and α_LD = 153 (Figure 1). In spite of the simplified nature of this model, its predictions turn out to be robust when compared to those of more sophisticated microscopic models.

Zoom In Zoom Out Reset image size

The fitting of nuclear interaction parameters to BEs using microscopic nuclear models was undertaken by Farine et al. (1978). Using a number of different non-relativistic effective interactions, they performed a series of Hartree–Fock calculations and found a nearly linear correlation between L and S_v: L ≃ 10S_v − 290 MeV (Figure 1). Later, Oyamatsu & Iida (2003) showed that fitting non-relativistic and relativistic Hartree–Fock, as well as constrained variational Thomas–Fermi, nuclear models to mass data led to the linear correlation L ≃ 13.3S_v − 373.3 MeV (Figure 1). Neither publication explored the statistical significance of their results. The slope of the Oyamatsu & Iida (2003) correlation matches that of the liquid droplet model, but that of the Farine et al. (1978) correlation is discrepant.

Kortelainen et al. (2010), simultaneously fitting nuclear masses and charge radii with the UNEDF0 potential and the Hartree–Fock methodology, explored the significance of the symmetry energy correlation. They obtained S_{v0, UNEDF0} ≃ 30.5 MeV, L_{0, UNEDF0} ≃ 45.1 MeV, σ_{v, UNEDF0} ≃ 1.58σ_UNEDF0 MeV, σ_{L, UNEDF0} ≃ 20.0σ_UNEDF0 MeV, and α ≃ 152. Kortelainen et al. (2010) adopted the somewhat arbitrary value of σ_UNEDF0 = 2 MeV, resulting in the 68% and 95% confidence intervals shown in Figure 1. The rms error of the UNEDF0 fit was χσ_UNEDF0 = 1.45 MeV.

This result closely matches fits obtained with both liquid droplet and Thomas–Fermi nuclear models (Figure 1). The length of the confidence interval for the liquid droplet model has the same size if we choose σ_LD = 0.50 MeV. For the finite-range Thomas–Fermi (FRTF) fits, the parameters for the shell, pairing, Wigner, and Coulomb exchange energies were taken from the liquid droplet best-fit model rather than including them in the fitting procedure. We obtained χσ_FRTF ≃ 2.75 MeV, S_{v0, FRTF} ≃ 31.2 MeV, and L_{0, FRTF} ≃ 56.6 MeV. The length of the Thomas–Fermi finite-range confidence interval equals that of Kortelainen et al. (2010) if σ_FRTF = 3.1 MeV.

The correlation coefficient between S_v and L is approximately 0.97 for UNEDF0, less than the nearly unity value estimated from the liquid droplet or Thomas–Fermi models since the UNEDF0 model was also constrained by nuclear radii. It should be noted that the 95% confidence interval contains symmetry parameter values that may be somewhat unphysical from the point of view of the properties of neutron-rich matter near and above the saturation density. For example, negative values of L imply that pure neutron matter has negative pressures near n_s, implying stability for neutron clusters. The size of the confidence ellipse scales with the assumed value of σ_UNEDF0, so had a value σ_UNEDF0 = 1 MeV been chosen instead, the 95% confidence interval would have extended no lower than L ≃ 20 MeV.

Another measure that can be applied to estimate the sizes of confidence intervals is the rms error of the fit. For example, the liquid droplet 68% confidence interval displayed in Figure 1 corresponds to an rms error of approximately 0.013 MeV baryon⁻¹. In other words, applied to liquid drop dense matter equations of state, such as that of Lattimer & Swesty (1991), symmetry parameters chosen from within this confidence region would yield errors in neutron or proton chemical potentials of less than about 0.013 MeV (aside from those arising from neglected shell and pairing effects).

The UNEDF collaboration has published a second parameter set, UNEDF1 (Kortelainen et al. 2012), for a modified energy density functional suitable for studies of strongly elongated nuclei. In this study, BEs and charge radii, as well as excitation energies, of spherical nuclei, deformed nuclei, and fission isomers, were fit. The best-fit values for the symmetry parameters of UNEDF1 lie near the 1σ confidence ellipse of UNEDF0, as shown in Figure 1. However, this study did not discuss a correlation between S_v and L. Because this fit is heavily weighted by deformed nuclei, and because the fit includes excitation energies, it will not show the same trends as the UNEDF0 fit.

As a final demonstration that this correlation is robust, the symmetry parameters determined from nuclear mass-fitting with the finite-range droplet model (FRDM; Myers & Swiatecki 1990; Möller et al. 2012) have consistently produced values of S_v and L along the UNEDF0 correlation line (Figure 1). This occurs despite differences among nuclear models, including the algebraic forms of the nuclear interaction and treatments of shell, pairing, deformation, and Wigner contributions. Unfortunately, confidence intervals for the FRDM fits were not published.

The correlation between symmetry parameters established from nuclear mass fitting has previously been used to constrain symmetry energy variations in supernova simulations (Swesty et al. 1994), but, surprisingly, has rarely been combined with results derived from other experimental data. We show that doing so considerably enhances the constraints established from other experiments.

An example is the neutron skin thickness of neutron-rich nuclei. Neglecting Coulomb effects, the difference between the mean neutron and proton surfaces in the liquid droplet model (Myers & Swiatecki 1969) is

$$\begin{equation} t_{{\rm np}}={2r_o\over 3}{S_sI\over S_v+S_sA^{-1/3}}, \end{equation} \tag{ 12 }$$

which is, to lowest order, a function of S_s/S_v. Given that S_s/S_v is largely a function of L, it should thus form a nearly orthogonal constraint to that from nuclear masses. Values of t_np have been measured, typically with 30%–50% errors, for a few dozen nuclei. More frequently used is the neutron skin thickness r_np, the difference between the mean square neutron and proton radii, which is approximately equal to $\sqrt{3/5}t_{{\rm np}}$.

A recent study of Sn isotopes (Chen et al. 2010), where differential isotopic measurements with fixed Z potentially reduce errors, produced a correlation band in the S_v–L plane which is nearly orthogonal to mass-fit correlations (Figure 2). The results from the Sn isotopes can, approximately, be expressed as a range of values of the neutron skin thickness of ²⁰⁸Pb, $r_{{\rm np}}^{208}$, since numerous theoretical studies have shown them to be highly correlated (e.g., Reinhard & Nazarewicz 2010).

Zoom In Zoom Out Reset image size

A number of authors have demonstrated from Hartree–Fock and Thomas–Fermi studies that $r_{{\rm np}}^{208}$ varies linearly with L for a given value of S_v (see, for example, Furnstahl 2002; Roca-Maza et al. 2011). However, $r_{{\rm np}}^{208}$ can also be shown to vary linearly with S_v if L is fixed (Furnstahl 2002; Reinhard & Nazarewicz 2010; Chen et al. 2010). In general, $r_{{\rm np}}^{208}$ is a function of both S_v and L, and the high degree of correlation between S_v and L is responsible for these separate, nearly linear, correlations.

Chen et al. (2010) established a two-dimensional relation, $r_{{\rm np}}^{208}(S_v,L)$, using Skyrme Hartree–Fock calculations of ²⁰⁸Pb in which most Skyrme parameters were determined from nuclear matter properties, and S_v and L were then systematically varied. We can approximate their results with

$$\begin{equation} r_{{\rm np}}^{208}(S_v,L)\simeq a_0+a_1S_v+a_2L+a_3S_v^2+a_4S_vL+a_5L^2, \end{equation} \tag{ 13 }$$

where the coefficients are given in Table 1. The mean error of this fit is about 0.7% over the ranges of interest. We attempted to reproduce the results for $r_{{\rm np}}^{208}(S_v,L)$ with a Skyrme Hartree–Fock analysis, and found similar results to within about 0.01 fm over the relevant ranges of S_v and L. The coefficients of our fit are also given in Table 1. The fits are not identical due to differences in the assumed energy density functionals as well as errors in extracting values from Figure 7 of Chen et al. (2010). By combining the linear correlation between L and S_v from Kortelainen et al. (2010), namely L ≃ −348.45 MeV + 12.903S_v, with Equation (13), we are able to reproduce the approximate relations between $r_{{\rm np}}^{208}$ and L found by Furnstahl (2002), Chen et al. (2005), and Centelles et al. (2009).

Brown (2000) and Typel & Brown (2001) found an approximately linear correlation between r_np and the pressure (p_N(0.1)) of pure neutron matter at the density 0.1 fm⁻³,

$$\begin{equation} r_{{\rm np}}=(0.060\pm 0.015)+0.12(p_N(0.1)/{\rm MeV\; fm}^{-3})\;{\rm fm}. \end{equation} \tag{ 14 }$$

This correlation is phenomenological because matter in the neutron-rich nuclear surface is very different from pure neutron matter. p_N(0.1) cannot be completely specified by S_v and L, partly because of variations in the symmetric matter pressure and higher order terms in the expansion of S₂ with respect to density (see Section 6.5). Nevertheless, we have verified that Equation (13), with the parameters in Table 1, reproduces this correlation.

The limits on $r_{{\rm np}}^{208}$ deduced from the neutron skin thickness of Sn isotopes by Chen et al. (2010), $0.155\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.195\,{\rm fm}$, yield the band shown in Figure 2. The above results indicate that this band is not very model dependent.

A recent study (Krasznahorkay et al. 2013) of neutron skin thicknesses in Sn and Pb nuclei, deduced from measurements of antianalog giant dipole resonances (GDRs), yielded the constraint $0.16\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.21\,{\rm fm}$. Zenihiro et al. (2010), utilizing proton elastic scattering from ²⁰⁸Pb, calculated $0.148\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.265\,{\rm fm}$. Both agree with the analysis of Chen et al. (2010), but neither further constrains the symmetry parameters. Tsang et al. (2012) averaged results from various studies of photon emission during the decays of antiprotonic states of neutron-rich nuclei to deduce estimates for neutron skin thicknesses. They inferred a central value $r_{{\rm np}}^{208}=0.18$ fm, but with even larger uncertainties. First results (Abrahamyan et al. 2012) from the PREX experiment to directly measure the neutron skin of ²⁰⁸Pb were too imprecise to be definitive, yielding $0.15\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.49\,{\rm fm}$. The central value is discrepant with other results for neutron skins, but the uncertainties are still very large.

The hydrodynamical droplet model of Lipparini & Stringari (1989) showed that the GDR centroid energy is closely connected to liquid droplet parameters. In this model, the centroid energy is

$$\begin{equation} E_{-1}\simeq \sqrt{{3\hbar ^2(1+\kappa)\over m \langle r^2 \rangle }S_v\left(1+{5\over 3}{S_s\over S_v}A^{-1/3}\right)^{-1}}, \end{equation} \tag{ 15 }$$

where κ is an enhancement factor arising from the velocity dependence of the interaction and 〈r²〉 is the mean-square charge radius. The factor m〈r²〉/(1 + κ) does not significantly vary among interactions for a given isotopic chain, and its value therefore does not play a role in estimating a correlation between the symmetry parameters. A series of measurements of the GDR should result in a correlation between S_s/S_v and S_v that is similar to but slightly less steep (i.e., a greater value of α) than that established from nuclear masses. Equations (6) and (15) predict that the orientation angles (α) of the GDR and BE correlations in the S_s/S_v–S_v plane will have the ratio

$$\begin{equation} {\alpha _{{\rm GDR}}\over \alpha _{{\rm binding}}} \! \approx {5\over 3}\left(\! 1+2{S_{s0}\over S_{v0}}\langle A^{-1/3}\rangle \! \right)\!\left(\! 1+{8\over 3}{S_{s0}\over S_{v0}}\langle A^{-1/3}\rangle \!\right)^{-1}\simeq 1.45, \end{equation} \tag{ 16 }$$

where 〈〉 refers to an average over the measured nuclei. A similar ratio will then apply to the relative orientation angles in the L–S_v plane.

Using a variety of Skyrme forces to model the GDR in ²⁰⁸Pb, Trippa et al. (2008) concluded the measured centroid energy was best fit by those forces having a bulk symmetry energy S₂(0.1), evaluated at the density n = 0.1 fm⁻³, in the range 23.3 MeV <S₂(0.1) < 24.9 MeV. The relation between S₂(0.1) and S_v and L is obviously model dependent.

To explore this model dependence, we selected a number of functional forms for S₂:

$$\begin{equation} \begin{array}{ll}S_2(u)=(T_0/3)u^{2/3}(1+cu)+a_1u+b_1u^{4/3}, & \qquad ({\rm i})\\ S_2(u)=S_v+(L/3)(u-1)+(K_{{\rm sym}}/18)(u-1)^2 &\\ \qquad\quad +(Q_{{\rm sym}}/162)(u-1)^3+d_1(u-1)^4, & \qquad ({\rm ii})\\ S_2(u)=(T_0/3)u^{2/3}+[S_v-(T_0/3)]u^\gamma, & \qquad ({\rm iii}) \end{array} \end{equation} \tag{ 17 }$$

where T₀ = (ℏ²/2m_b)(3π²n_s/2)^2/3 and

$$\begin{eqnarray} a_1&=&4S_v-{2T_0\over 3}\left(1-{c\over 2}\right)-L, \!\quad b_1={T_0\over 3}(1-2c)+L-3S_v,\cr d_1&=&{L\over 3}-S_v-{K_{{\rm sym}}\over 18}+{Q_{{\rm sym}}\over 162}, \!\quad \gamma ={3L-2T_0\over 3(3S_v-T_0)}. \end{eqnarray} \tag{ 18 }$$

Form (i) of Equation (17) originates from Skyrme interactions. Here, c accounts for the nucleon effective mass in the kinetic energy contribution. We considered cases in which c = 0 or c = 0.4. Form (ii) is the polynomial density expansion with coefficients of Equation (2). The parameter d₁ ensures that S₂(0) = 0. To evaluate K_sym and Q_sym, which are not directly constrained by experiment, we considered two cases: Ducoin et al. (2011) argued them to be linearly correlated with L for Skyrme functionals fit to nuclear properties:

$$\begin{eqnarray} K_{{\rm sym}}&\simeq &3.33L-281\,{\rm MeV},\cr \ms Q_{{\rm sym}}&\simeq &-6.63L+765\,{\rm MeV}, \end{eqnarray} \tag{ 19 }$$

while Vidaña et al. (2009), using Q_sym = 0, found

$$\begin{equation} K_{{\rm sym}}\simeq 2.867L-260\,{\rm MeV}. \end{equation} \tag{ 20 }$$

This ensemble of four examples adequately describes the symmetry behavior of both Skyrme and relativistic field nuclear interactions. Finally, form (iii) has often been used in fitting results from heavy-ion collisions (Zhang et al. 2008).

The experimental constraints on S₂(0.1) using these variations of Equation (17) are displayed in Figure 3, from which we infer that the effective model dependence is moderate. For example, the terms involving Q_sym and d₁ have little influence on these results. The form Equation (17 iii) deviates at small values of S_v. We show in Section 6.5 that this parameterization is a poor representation of neutron matter at subnuclear densities. To be conservative, however, we have chosen to adopt the extreme ranges to represent the experimental GDR limits on S_v and L (Figure 2). The resulting band has a shallower slope, compared to that inferred from nuclear masses, as predicted.

Zoom In Zoom Out Reset image size

The electric dipole polarizability, α_D, is related to the response of a nucleus to an externally applied electric field, and it is largely concentrated in the GDR. The approximate proportionality of α_D and $r_n^{208}$ was demonstrated macroscopically by Lipparini & Stringari (1989) and confirmed by microscopic calculations using energy-density functionals by Reinhard & Nazarewicz (2010). Tamii et al. (2011) measured the dipole polarizability of ²⁰⁸Pb to be α_D = (20.1 ± 0.6) fm³. Using the very high degree of correlation between α_D and $r_{{\rm np}}^{208}$ established by Reinhard & Nazarewicz (2010), this measurement of α_D is equivalent to the range $0.131\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.177\,{\rm fm}$. The analysis by Reinhard & Nazarewicz (2010) was based upon a non-relativistic density functional analysis. A systematic reanalysis (Piekarewicz et al. 2012), incorporating both relativistic and non-relativistic models calibrated to nuclear properties in an attempt to minimize model dependencies, revised this permitted range to $0.146\,{\rm fm}\lesssim r_{{\rm np}}^{208}\lesssim 0.190\,{\rm fm}$. This result is nearly the same as that established from neutron skin thicknesses of Sn isotopes (Figure 2).

Constraints from heavy ion collisions include those from isospin diffusion (Tsang et al. 2009), shown in Figure 2, and from multifragmentation in intermediate-energy heavy ion collisions (Shetty et al. 2007). The latter implied 40 MeV <L < 125 MeV, which is consistent with our other constraints, but does not further limit parameter space and is therefore not displayed in Figure 2 for clarity. Both limits have been established from specific types of nuclear interaction models, and may be subject to considerable model dependence.

Also not included in Figure 2 are results of a study of pygmy dipole resonances (Carbone et al. 2010) in ⁶⁸Ni and ¹³²Sn, which implied 31 <S_v/MeV < 33.6 and 49.1 <L/MeV < 80. Daoutidis & Goriely (2011) argued that theoretical and experimental uncertainties prevent them from being effective constraints. The results of Reinhard & Nazarewicz (2010) confirmed that pygmy resonance energies lack a significant correlation with L.

Danielewicz & Lee (2009) proposed isobaric analog states as constraints on symmetry parameters and concluded that S_v ≃ 33 MeV and S_s ≃ 96 MeV, which is equivalent to L ≈ 110 MeV. These values were derived from a fit including bulk and surface contributions to data for nuclei over a broad range of masses. However, we reanalyzed the fit including a curvature term (proportional to A^−2/3) and determined that S_v ∼ 30 MeV and S_s ∼ 45 MeV, which would be equivalent to L ≈ 50 MeV. Neither fit is statistically preferred. A more detailed analysis will be necessary to produce better constraints.

In summary, the various experimental results jointly confine the allowable ranges of symmetry parameters to a small region (Figure 2) bounded by 29.0 MeV <S_v < 32.7 MeV and 40.5 MeV <L < 61.9 MeV. The lower limit to L is set by the heavy-ion collision analysis, while the upper limit to L is determined by the dipole polarizability analysis.

The most direct connection between astrophysical observations and the nuclear symmetry energy concerns neutron star radii (R), which are highly correlated with neutron star pressures near n_s. Lattimer & Prakash (2001) found that the radii of neutron stars for masses near 1.4 M_☉ obey a power-law relation:

$$\begin{equation} R(M)=C(n,M)(p(n)/{\rm MeV\, fm}^{-3})^{1/4}, \end{equation} \tag{ 23 }$$

where R(M) is the radius of a star of mass M, p(n) is the pressure of neutron star matter at density n, and C(n, M) is a constant for a given density and mass. Lattimer & Prakash (2001) found that C(n_s, 1.4 M_☉) = 9.30 ± 0.60 km and C(2n_s, 1.4 M_☉) = 5.72 ± 0.25 km. These values were estimated by averaging the results of 31 disparate equations of state, several of which are now ruled out because of the maximum mass constraint imposed by PSR J1614−2230 (Demorest et al. 2010). Excluding those ruled-out equations of state, we find a revised value

$$\begin{equation} C(n_s,1.4\,M_\odot)=9.52\pm 0.49\,{\rm km}, \end{equation} \tag{ 24 }$$

with a standard deviation 20% smaller than previously obtained. For comparison, at 2n_s, we find the fractional errors in the revised constant are now smaller by nearly a factor of two: C(2n_s, 1.4 M_☉) = 5.68 ± 0.14 km. Although this relation is a tighter fit when applied for densities somewhat in excess of the nuclear saturation density (e.g., 1.5–2n_s), neutron star matter pressure cannot be accurately predicted using the symmetry parameters S_v and L alone at these densities. Hence, we utilize Equation (23), evaluated at n_s, for predictions.

The pressure of pure neutron matter at n_s in the quadratic approximation (Equation (1)) is p(n_s, x = 0) = Ln_s/3. However, neutron star matter is in β equilibrium with a finite proton fraction and has a smaller pressure by a factor

$$\begin{equation} f\equiv {p(n_s,x_{\beta s})\over p(n_s,0)}=1-x_{\beta s}\left[4(1-x_{\beta s})-{3S_v\over L}(1-2x_{\beta s})\right]. \end{equation} \tag{ 25 }$$

x_βs, the proton fraction in β equilibrium at the saturation density, is determined from S_v by the cubic equation

$$\begin{equation} 4S_v(1-2x_{\beta s})=\hbar c(3\pi ^2n_sx_{\beta s})^{1/3}. \end{equation} \tag{ 26 }$$

For the range of values within the overlap region indicated in Figure 2, 29.0 MeV < S_v < 32.7 MeV and 40.5 MeV < L < 61.9 MeV, these relations imply 0.035 < x_βs < 0.046 and 0.90 < f < 0.94.

To estimate the radius and its error for 1.4 M_☉ stars, we replaced the polygonal overlap region in Figure 2 with an error ellipse such that Monte Carlo sampling would result in 68% of the points lying within the polygon. The ellipse has a semimajor (semiminor) axis of 6.7 (0.55) MeV and an orientation angle of 86 from the vertical in the S_v–L plane. With these specifications, we find that, with Equation (24), R_1.4 = 11.9 ± 0.7 km to 1σ and R_1.4 = 11.9 ± 1.2 km to 90% confidence.

A study by Psonis et al. (2007), using the parameterized equation of state (EOS) of Prakash et al. (1988) and assuming S_v = 30 MeV, found a nearly linear correlation between R_1.4 and L, R_1.4 ≃ 8.702 + 0.070 (L/MeV) km, which agrees with Equations (23) and (24) to within their errors only for the interval 20 MeV < L < 60 MeV. The study of Gearheart et al. (2011) is only marginally consistent with predictions of Equations (23) and (24), perhaps because this study did not impose the neutron star maximum mass constraint.

Although more than 30 neutron star masses are now measured to a reasonable precision (Lattimer 2012), there are no precise measurements of neutron star radii. Among the most promising sources for determining radii are those stars having observed thermal radiation from their surfaces as well as certain X-ray bursts.

The interpretation of thermal emission has several difficulties, including uncertainties in distance, atmospheric composition, surface magnetic field, and interstellar absorption. Now-quiescent neutron stars in formerly accreting binaries in globular clusters (Heinke et al. 2006; Webb & Barret 2007; Guillot et al. 2011) are promising because it is believed that recent accretion has both suppressed the surface magnetic fields and rendered their atmospheres with a nearly pure hydrogen composition. In addition, there have been significant advances in estimating distances to globular clusters. If these sources were true blackbodies, measured fluxes and temperatures would directly yield their apparent angular radii, R_∞/D, where R_∞ = R(1 − 2β)^−1/2, β = GM/Rc², and M is the stellar mass. These sources imply a wide range of values for R_∞. A recent study found radii ranging from 8 km to 23 km (Guillot et al. 2013), using a low magnetic field, pure hydrogen atmosphere. Under the assumption that five quiescent sources would have a common radius, but possibly varying masses, motivated by the nearly vertical M–R trajectories found by Steiner et al. (2010, 2013; see below), they obtained a joint fit with $R=9.1^{+1.3}_{-1.5}$ km to 90% confidence, which is inconsistent with our results. However, this analysis has the theoretically unlikely result that, even with a common radius, the masses range from 0.86 M_☉ to 2.42 M_☉ (Guillot et al. 2013). Lattimer & Steiner (2013) suggest that uncertainties in the amount of interstellar absorption and atmosphere compositions contribute to the wide predicted ranges of mass and radius.

There are a few non-globular cluster sources having thermal emission, but only one, RXJ 1856−3754, has a precisely measured distance (Walter et al. 2010). Ho et al. (2007) modeled it with a carbon-rich surface with a thin layer of hydrogen. Including only distance uncertainties, they found M ≃ 1.33 ± 0.09 M_☉ and R ≃ 11.9 ± 0.8 km. However, no confirmation of this assumed atmosphere composition is available.

Another class of sources, Type I photospheric radius-expansion (PRE) X-ray bursts (Strohmayer & Bildsten 2004), are thought to be powerful enough to temporarily lift material from neutron star surfaces. They have observed maximum fluxes thought to be near their Eddington limit values, F_Edd = cGM/(κD²) (neglecting redshift), for which radiation pressure balances gravity (κ is the opacity). Combined with angular sizes determined from the tail of the burst, this permits simultaneous mass and radius determinations. For these sources, major uncertainties include distances as well as atmospheric properties. It is also uncertain if the location of the photosphere, when the maximum flux is observed, is close enough to the stellar radius to make redshift corrections to F_Edd significant. Estimated errors in M and R are of order 20% or more for each source.

Steiner et al. (2010) utilized M–R data from observations of an ensemble of four quiescent and four PRE burst sources to estimate the parameters for an EOS whose symmetry energy near the nuclear saturation density was parameterized by Equation (17 iii). Their Bayesian analysis found, to 1σ, 28.5 MeV ≲ S_v ≲ 34 MeV and 0.2 ≲ γ ≲ 0.4, which implies 41 MeV ≲ L ≲ 51 MeV. They also found, to 1σ for their baseline case, 11.5 km ≲ R_1.4 ≲ 12.2 km. Alternatively, utilizing their derived constraints for S_v and L, Equation (23) predicts 11.0 km ≲ R_1.4 ≲ 12.3 km, which is quite consistent. In contrast, studies of PRE bursts by Özel et al. (2010, 2012) deduced radii with 2σ uncertainties in the range 9–11 km. Also, Suleimanov et al. (2011) proposed that the short PRE bursts considered by Ozel et al. and Steiner et al. have systematic uncertainties, such as disc shadowing, that make them unsuitable for mass and radius determinations. They instead proposed the consideration of longer PRE bursts from which they find significantly larger radii, with estimates extending beyond 14 km, which necessarily imply larger L values (100 MeV or larger, according to Equation (23)), which are inconsistent with nuclear systematics. Collectively, these differences suggest that current models are subject to significant systematic uncertainties involving color correction factors, redshift corrections, and geometry.

A study by Steiner & Gandolfi (2012) combined neutron matter calculations and astrophysical observations to constrain L. They assumed the EOS in the vicinity of the nuclear saturation density was parameterized by an energy density functional with parameters taken from the neutron matter results of Gandolfi et al. (2012). They found, for neutron star models both with and without deconfined quark matter at high densities, that, to 1σ, 43 MeV ≲ L ≲ 52 MeV and 32 MeV ≲ S_v ≲ 34 MeV, and, to 2σ, 39 MeV ≲ L ≲ 54 MeV and 31 MeV ≲ S_v ≲ 34.5 MeV. To 1σ, they also inferred a neutron star radius range of 10.75–12.25 km for a 1.4 M_☉ neutron star, only about 0.35 km smaller, on average, than obtained by Steiner et al. (2010).

A more recent study (Steiner et al. 2013) considered additional sources, and in addition considered systematic variations to the fiducial analysis in which one or more sources were excluded from the Bayesian analysis. In extreme cases, all PRE X-ray burst sources or quiescent globular cluster sources were excluded. With these variations in input data, it was inferred that 43.3 MeV ≲ L ≲ 66.5 MeV (Figure 2), to 1σ, somewhat larger on average than that deduced by Steiner et al. (2010) and Steiner & Gandolfi (2012), but very close to the range found from experimental results. In agreement with the preceding studies, it was found that S_v is relatively poorly constrained. The 1σ radius range was found to be 10.6 km ≲ R ≲ 12.6 km. The average value is the same as determined by Steiner et al. (2010), but the range is broader due to the consideration of a greater variety of scenarios in an attempt to capture systematic uncertainties in models. It should be emphasized that neither Steiner et al. (2010) nor Steiner et al. (2013) assumed an a priori relation between the symmetry parameters, despite the fact that nuclear BEs imply that a strong correlation exists.

Finally, we note that recent analyses (Leahy et al. 2011) of pulsed light from rotating, accreting millisecond pulsars are consistent with a mass–radius relation similar to that deduced by Steiner et al. (2010).

The moment of inertia depends sensitively on the neutron star radius, scaling approximately as R². Hence, constraints to the neutron star radius will also constrain a neutron star's moment of inertia. The pulsar binary PSR 0737–3039 (Lyne et al. 2004) has detectable spin–orbit coupling, and the potential exists for a measurement of the most massive component's moment of inertia. Since the masses in this binary are known to a high degree of precision, this presents the possibility of measuring the radius of a neutron star. Lattimer & Schutz (2005) found that the moment of inertia, for neutron star EOSs capable of supporting at least 1.6 M_☉, can be well approximated by

$$\begin{equation} I=(0.237\pm 0.008)MR^2(1+2.84\beta +18.9\beta ^4), \end{equation} \tag{ 27 }$$

where β = GM/Rc² is the neutron star compactness parameter. For a 1.4 M_☉ neutron star, this relation, coupled with the radius constraints implied by the experimental limits to S_v and L, predicts I_1.4 = 71.0 ± 7.1 M_☉ km². Another approximation for the moment of inertia (Link et al. 1999), after correction of a typographical error, is

$$\begin{equation} I\simeq {2\over 7}MR^2\left(1-{33\over 28}\beta \right)^{-1}, \end{equation} \tag{ 28 }$$

which gives I_1.4 = 71.4 ± 7.2 M_☉ km², a nearly identical result.

Flanagan & Hinderer (2008) have recently pointed out that tidal effects are potentially measurable during the late stages of inspiral of a binary merger involving neutron stars. Gravitational waves emitted early during the detection phase (when the orbital frequency becomes larger than about 100 Hz) are nearly sinusoidal with gradually increasing frequency and amplitudes. Tidal forces induce quadrupole moments on the neutron stars, producing a very small correction which appears as an accumulated phase shift in the gravitational wave signal. This shift will be proportional to a weighted average of the induced quadrupole polarizabilities (λ_tid, also known as the Love number) for the individual stars. Since both neutron stars have the same EOS, this weighted average is relatively insensitive to the binary mass ratio (Hinderer et al. 2010). Since λ scales as R⁵, it would provide a sensitive measurement of the neutron star radius. However, the results of Hinderer et al. (2010) and Postnikov et al. (2010) show that for radii restricted to our favored range, λ_tid is small, making detection of the tidal signature by Advanced LIGO difficult (Hinderer et al. 2010). Much of the difficulty has to do with degeneracies between the spins of the neutron stars, their rotational quadrupole moments (Q_rot), and their tidal Love numbers λ_tid in their extraction from gravitational wave signals.

However, Yagi & Yunes (2013) pointed out that nearly equation-of-state (EOS)-independent relations exist among Q_rot, I, and λ_tid which can be used to break these degeneracies. In particular, the phenomenological relations Yagi & Yunes (2013) found relating I to λ_tid and Q_rot have interesting consequences. These are

$$\begin{eqnarray} \bar{I}\equiv {I/ M^3}&\simeq &\exp (1.47+0.817z+1.49z^2 \nonumber\\ && +0.287z^3-0.364z^4), \nonumber\\ \bar{Q}_{{\rm rot}}\equiv Q_{{\rm rot}}/M^5&\simeq &\exp (0.194+0.936z+0.474z^2 \nonumber\\ && -4.21z^3+1.23z^4), \end{eqnarray} \tag{ 29 }$$

where $z=(\ln \bar{\lambda }_{{\rm tid}})/10$ and $\bar{\lambda }_{{\rm tid}}=\lambda _{{\rm tid}}/M^5$. It is claimed that these relations are accurate to better than a percent for all realistic EOSs and for neutron star masses greater than about 0.5 M_☉, although they studied only four EOSs. We have verified the general validity of this relation, although we find it holds to slightly lesser accuracy than claimed. We find this relation holds to 2% accuracy for all the hadronic EOSs described in Lattimer & Prakash (2001) for stars greater than 1 M_☉, except for PS, WFF3, and WFF4. For the latter EOSs, none of which satisfy the minimum maximum mass constraint, the error increases beyond 2%, but only for M < 1.4 M_☉.

These equations imply that, for a 1.4 M_☉ star, $\bar{\lambda }_{{\rm tid}}\simeq 72\pm 32$ and $\bar{Q}_{{\rm rot}}\simeq 3.13\pm 0.45$.

These relatively EOS-independent relations can be combined with Equation (27), which itself is accurate to about 3% for the same mass range and the same set of EOSs, excluding those that cannot satisfy the minimum maximum mass constraint imposed by PSR 0737–3039. This results in relations relating $\bar{\lambda }_{{\rm tid}}$ and $\bar{Q}_{{\rm rot}}$ to the compactness β alone. We determined inverse relationships for Equation (29) for this purpose:

$$\begin{eqnarray} \bar{\lambda }_{{\rm tid}}&\simeq &\exp (-4.24+5.96u-3.01u^2+0.722u^3-0.0655u^4),\cr \ms \bar{Q}_{{\rm rot}}&\simeq &\exp (-1.16-0.700u+1.83u^2-0.593u^3+0.062u^4 ), \nonumber\\ \end{eqnarray} \tag{ 30 }$$

where

$$\begin{equation} u=\ln (\bar{I})\simeq \ln [0.237\beta ^{-2}(1+2.84\beta +18.9\beta ^4)]. \end{equation} \tag{ 31 }$$

We confirmed that these relations are valid to within about 3% for masses larger than 1 M_☉ and β ≲ 0.3 for EOSs that satisfy the minimum maximum mass constraint. Maselli et al. (2013) point out that a measurement of λ_tid to even 60% accuracy can be used to determine β to 10% accuracy, because λ_tid∝R⁵. This would translate into a similar accuracy for the radius assuming the masses can be determined to at least this accuracy. However, we note that a simpler expression for $\bar{\lambda }_{{\rm tid}}(\beta)$ obtained by Maselli et al. (2013) is considerably less accurate than Equation (30) when applied to the same suite of EOSs and for the same conditions on mass and compactness.

One of the more accurate measurements that can be determined from the neutrino signal from a galactic gravitational collapse supernova is the neutron star BE, assuming the distance to the supernova can also be found. Lattimer & Prakash (2001) determined a fit to the BE of a neutron star, for EOSs which support at least 1.65 M_☉:

$$\begin{equation} {{\rm BE}\over M} \simeq (0.60\pm 0.05)\beta \left(1-{\beta \over 2}\right)^{-1}. \end{equation} \tag{ 32 }$$

Using the radius constraints implied by the experimental limits to S_v and L implies BE $\simeq 279^{+35}_{-25}$ B for a 1.4 M_☉ neutron star, where 1 B is 10⁵¹ erg.

Black hole–neutron star mergers are among the most likely progenitors of short gamma-ray bursts (Lee & Ramirez-Ruiz 2007), sources of gravitational radiation (see, for example, Abadie et al. 2010; Somiya 2012), and sites for the astrophysical r-process (Lattimer & Schramm 1976; Lattimer et al. 1977). Gamma-ray bursts and r-process production are most likely in the case the neutron star is tidally disrupted, and radioactive decay in unbound neutron-rich ejecta could produce an isotropic optical (Roberts et al. 2011; Metzger & Berger 2012) or infrared (Kasen et al. 2013) transient. The potential for mass ejection depends on whether or not tidal effects are strong enough for the neutron star to be disrupted before reaching the innermost stable circular orbit (R_ISCO) of the black hole. In this case, some disrupted matter remains outside the black hole for times up to seconds in an accretion disk, a tidal tail, or in unbound matter. Primarily, the outcome depends on the mass ratio of the binary (q = M_BH/M_NS), the spin of the black hole (χ = a/M_BH), and the neutron star compactness (β). Low mass ratios, high spins, and large neutron star radii are favorable for large disk masses. Recently, Foucart (2012) has performed a study of numerical simulations in which these parameters were varied and formulated a two-parameter model that predicts the disk mass remaining outside the black hole a few milliseconds following a merger. The observed energy of a short gamma-ray burst (Giacomazzo et al. 2013) and, possibly, the mass of unbound matter, may be proportional to the disk mass, so observations of gamma-ray bursts can be compared to models of population synthesis of compact binaries in support of this hypothesis.

The fit of Foucart (2012) predicts that the disk mass fraction M_t/M_NS resulting from a merger varies with q, β, and χ according to

$$\begin{equation} {M_t\over M_{{\rm NS}}}={M^b_{{\rm NS}}\over M_{{\rm NS}}}\left[a_1(3q)^{1/3}(1-2\beta)-a_2\beta q{R_{{\rm ISCO}}c^2\over GM_{{\rm BH}}}\right], \end{equation} \tag{ 33 }$$

where $M^b_{{\rm NS}}$ is the baryon mass of the neutron star and a₁ = 0.288 ± 0.011 and a₂ = 0.148 ± 0.007 are parameters determined by a fit to merger simulations. The radius R_ISCO of the innermost stable circular orbit depends on χ according to

$$\begin{equation} \chi ={1\over 3}\sqrt{R_{{\rm ISCO}} c^2\over GM_{{\rm BH}}}\left[4-\sqrt{3{R_{{\rm ISCO}} c^2\over GM_{{\rm BH}}}-2}\right], \end{equation} \tag{ 34 }$$

valid for either prograde or retrograde orbits. We also note that Equation (32) implies that

$$\begin{equation} {M^b_{{\rm NS}}\over M_{{\rm NS}}}\simeq 1+{3\over 5}{\beta \over 1-\beta /2}. \end{equation} \tag{ 35 }$$

Foucart (2012) argues that the major effects of general relativity are taken into account by the factor 1 − 2β and the dependence of R_ISCO on M_BH and χ.

Equation (33) can be expressed in the form χ = χ(q, β, M_t, M_NS). Foucart (2012) presented results as a series of figures illustrating the ratio of disk mass to $M^b_{{\rm NS}}$ as a function of the three variables q, β, and χ. This four-dimensional function can be reduced to essentially two dimensions by employing experimental restrictions to the nuclear symmetry energy which confine neutron star radii to small ranges. Given that the range of L values deduced from experiments is nearly the same as that inferred from astrophysical modeling, we use the baseline results from Steiner et al. (2010, 2013) that 11 < R < 12.5 km is valid not only for 1.4 M_☉ neutron stars, but also for stars in the range 1.2–1.8 M_☉. We further restrict the likely values of M_NS to this range. With these constraints, 0.14 ≲ β ≲ 0.24 and Equation (33) therefore implies a small range of χ values for particular choices of q and M_t. As a result, the series of diagrams shown in Foucart (2012) can be distilled into a single picture (Figure 4). We find that the resulting width of χ bands corresponding to given M_t values is minimized if the abscissa is chosen to be qβ² (justified below). Nevertheless, due to the restricted range of β, the abscissa is still dominated by its q dependence. Generally, the bottom (top) of each χ band in Figure 4 corresponds to M_NS = 1.2(1.8)M_☉ and R = 12.5(11) km. The region to the right of the M_t = 0 band corresponds to the situation that tidal disruption occurs within the ISCO so that no disk is able to form and no short gamma-ray burst will be observed. It is readily seen that the disk mass increases rapidly with black hole spin, and that maximally spinning holes of virtually any stellar mass will form disks after merging. The disk mass also increases with decreasing q or β. This behavior is easily understood if Equation (33) is rewritten as

$$\begin{equation} r \!\equiv \! {R_{{\rm ISCO}}c^2\over GM_{{\rm BH}}} \!=\! {a_1(q\beta ^2)^{1/3}\beta ^{1/3}(1-2\beta)-(M_t/R)\left(M_{{\rm NS}}/M_{{\rm NS}}^b \right)\over a_2q\beta ^2}. \end{equation} \tag{ 36 }$$

We note that 0.32 ≲ β^1/3(1 − 2β) ≲ 0.37 and $0.86\lesssim M_{{\rm NS}}/M_{{\rm NS}}^b\lesssim 0.92$ are weak functions of β in the relevant range, so that r is largely a function of M_t and qβ² alone. Also recall that r is a decreasing function of χ = a/M_BH. Clearly, r decreases and a/M_BH increases with increasing q²β for small to moderate values of M_t. However, for large enough values of M_t, of order 0.5 M_☉ or greater, the numerator of Equation (36) vanishes for small values of q²β so that χ ≈ 1 for all values of the abscissa.

Zoom In Zoom Out Reset image size

A further consequence of Equation (33) is that the mass of the disk and, presumably, the unbound matter will be sensitive to the neutron star compactness and expected to decrease with increasing β.

The most energetic observed short gamma-ray bursts imply disk masses of order 0.5 M_☉ (Giacomazzo et al. 2013), which can only be accommodated for black hole–neutron star systems with χ ≃ 1. Alternatively, Giacomazzo et al. (2013) suggest that binary neutron star mergers could more likely produce large disk masses. They estimated from general relativistic simulations using idealized equations of state that

$$\begin{equation} M_t\approx 3(1-q)(1+q-M_{{\rm tot}}/M_{{\rm max}})\simeq 3(1-q^2)(1-M_1/M_{{\rm max}}), \end{equation} \tag{ 37 }$$

where q < 1 is the binary mass ratio, M_tot is the total gravitational mass of the two neutron stars, and M_max is the maximum mass of an isolated neutron star. This formula predicts that a disk mass of order 0.5 M_☉ could be produced if q is significantly different from 1 and if the larger neutron star is considerably less massive than the maximum mass. However, Equation (37) is suspect because it is based on simulations employing idealized fluid EOSs, which have unphysical M–R relations; in addition, the dependence on neutron star compactness was not considered. It can be expected, in particular, that the disk mass should decrease with compactness as was already noted for black hole–neutron star mergers above. In fact, this inverse behavior is observed for the mass of unbound matter in the extensive set of general relativistic simulations employing realistic EOSs produced by Bauswein et al. (2013).

In addition to neutron star radii and moments of inertia, several properties of the neutron star crust are also impacted by the nuclear symmetry energy. The most important effect is the location of the core–crust boundary, where a phase transition separates heterogeneous matter containing nuclei and a neutron fluid, and homogeneous, uniform nuclear matter. It has been demonstrated that the thickness of the crust, ΔR, defined as the distance between neutron drip density and the core–crust transition density, depends primarily on the stellar radius and the neutron chemical potential at the core–crust boundary (Lattimer & Prakash 2007):

$$\begin{equation} \Delta R/R\simeq ({\cal H}-1)[{\cal H}(1-2\beta)^{-1}-1]^{-1}. \end{equation} \tag{ 38 }$$

In this expression, ${\cal H}=\exp [2(\mu _{nt}-\mu _{n0})/m_bc^2]$, where μ_nt and μ_n0 ≃ −9 MeV are the neutron chemical potential at the core–crust transition and at the stellar surface, respectively. The crust thickness, together with its transport properties, determines the thermal relaxation time of the crust, which is important in the cooling of young neutron stars (Yakovlev & Pethick 2004), such as the recently detected rapid cooling of Cassiopeia A (Heinke & Ho 2010; Page et al. 2011) and in the cooling of soft X-ray transients (Page & Reddy 2012).

A related quantity is the fraction of the star's moment of inertia residing in the crust, ΔI/I, which depends primarily on the stellar mass and radius and the pressure at the core–crust boundary, and scales as p_tR⁴M⁻² (Link et al. 1999). The standard model for pulsar glitches holds that they are due to the neutron superfluid in the star's crust. In this case, the observed glitch rates and magnitudes for the Vela pulsar imply that

$$\begin{equation} \Delta I/I\gtrsim 0.016, \end{equation} \tag{ 39 }$$

and similar or smaller values are found for other glitching pulsars (Andersson et al. 2012).

An approximate method of determining the location of the core–crust boundary is by considering the stability of a homogeneous fluid of nucleons in beta equilibrium. Baym et al. (1971) showed that small density fluctuations in an otherwise uniform fluid lead to instability when

$$\begin{equation} \mu _{{\rm pp}}-\mu _{{\rm np}}^2\mu _{{\rm nn}}^{-1}+4e\sqrt{\pi \eta }-4\eta \alpha \left(9\pi n_p^2 \right)^{1/3}=0, \end{equation} \tag{ 40 }$$

where μ_ij = ∂μ_i/∂n_j and μ_i and n_i are the chemical potential and number density, respectively, of neutrons, protons or electrons. α ≃ 1/137 is the fine-structure constant. β is determined by density gradient terms in the nuclear Hamiltonian and can be approximated by η = D[1 − 4μ_np/μ_nn + (μ_np/μ_nn)²], where D ≃ 100 MeV fm⁵ (Hebeler et al. 2013). It has been assumed that electrons are relativistic and degenerate. The chemical potentials and derivatives in Equation (40) are equivalent to

$$\begin{eqnarray} \mu _n&=&{\partial (ne)\over \partial n}-x{\partial e\over \partial x},\quad \mu _p={\partial (ne)\over \partial n}+(1-x){\partial e\over \partial x},\cr \mu _{{\rm nn}}\mu _{{\rm pp}}-\mu _{{\rm np}}^2&=&{1\over n}{\partial ^2(ne)\over \partial n^2}{\partial ^2e\over \partial x^2}-\left({\partial ^2e\over \partial n\partial x}\right)^2,\cr \mu _{{\rm np}}&=&{\partial ^2(en)\over \partial n^2}+(1-2x){\partial ^2e\over \partial n\partial x}-{x(1-x)\over n}{\partial ^2e\over \partial x^2},\cr \mu _{{\rm nn}}&=&{\partial ^2(en)\over \partial n^2}-2x{\partial ^2e\over \partial n\partial x}+{x^2\over n}{\partial ^2e\over \partial x^2}. \end{eqnarray} \tag{ 41 }$$

Obviously, the functional form of the energy per baryon e(n, x) will influence the core–crust transition density u_t = n_t/n_s determined by Equation (40). We assume a quadratic symmetry energy with S₂(u) given by Equation (17), and two different functional forms for the symmetric matter energy,

$$\begin{equation} \begin{array}{ll}e(u,1/2)=(3T_0/5)u^{2/3}[1+cu]+a_0u+b_0u^{4/3}, &\qquad (i) \\ e(u,1/2) =-B+(K/18)(u-1)^2+k(u-1)^3 &\nonumber\\ \qquad\qquad +d_0(u-1)^4, &\qquad (ii) \end{array} \end{equation} \tag{ 42 }$$

where

$$\begin{eqnarray} a_0 &=& -\left[4B+{6T_0\over 5}\left(1-{c\over 2}\right)\right],\quad b_0={3T_0\over 5}(1-2c)+3B,\nonumber\\ d_0 &=& B-{K\over 18}+k. \end{eqnarray} \tag{ 43 }$$

In the above, we assume B = 16 MeV is the bulk BE at n_s and K = 240 MeV is the incompressibility. The coefficient d₀ ensures the energy per particle vanishes at zero density, e(0, 1/2) = 0. With these two forms, and the five associated with Equation (17), we have 10 different expressions for modeling nuclear matter. Excluding variations utilizing Equation (17 iii), we found that the energy and pressure predicted for pure neutron matter can be adequately approximated, as shown in Figure 5, for the density range of interest for determining the core–crust transition densities. In β-equilibrium, at these densities, pure neutron matter is very similar to neutron star matter since proton fractions are of order 1%.

Zoom In Zoom Out Reset image size

Choosing values for S_v and L within the overlap region in Figure 2, as previously described, the above energy expressions can be used to find the core–crust transition density, pressure, and neutron chemical potential, u_t, p_t, and ${\cal H}_t$, respectively. For each EOS model, u_t (p_t and ${\cal H}_t$) decreases (increase) with L or S_v with narrow spreads. When all possible variations are considered, the values of these quantities, and in particular p_t, are poorly constrained. However, limiting S₂ to the forms (i) and (ii) of Equation (17), supported by neutron matter studies, results in relatively tight constraints. Monte Carlo sampling to 90% confidence limits the transition density to 0.49 ≲ u_t ≲ 0.60, the pressure to 0.40 MeV fm⁻³ ≲ p_t ≲ 0.60 MeV fm⁻³, and ${\cal H}$ to $1.043\lesssim {\cal H}\lesssim 1.053$. This represents a considerable improvement over the analysis of Lattimer & Prakash (2001) which considered interactions leading to 0.2 MeV fm⁻³ ≲ p_t ≲ 0.65 MeV fm⁻³.

In comparison, Hebeler et al. (2013), working directly from the properties of pure neutron matter, found 0.48 ≲ u_t ≲ 0.56 and 0.41 MeV fm⁻³ ≲ p_t ≲ 0.47 MeV fm⁻³. Both quantities are within the boundaries of our determinations, consistent with the relative S_v and L values implied by neutron matter studies relative to experimental constraints (Figure 2). Hebeler et al. (2013) also found that u_t decreased and p_t increased linearly with increasing L or S_v values. Conversely, Xu et al. (2009) found 0.25 ≲ u_t ≲ 0.40 and 0.01 MeV fm⁻³ ≲ p_t ≲ 0.26 MeV fm⁻³ using the MDI interaction, both below and outside our estimated bounds.

The crust thickness of a 1.4 M_☉ star, using Equation (38), the experimental constraints to S_v and L, forms (i) and (ii) of Equation (17), and Monte Carlo sampling, can consequently be limited to ΔR/R = 0.082 ± 0.008, or $\Delta R=0.96^{+0.17}_{-0.13}$ km, with 1σ errors.

The fraction of the moment of inertia of the neutron star crust depends primarily on the mass and radius of the star as well as p_t and, to a lesser extent, u_t or h_t. Lattimer & Prakash (2007) showed that this fraction is well approximated by

$$\begin{equation} {\Delta I\over I}\simeq {8\pi p_tR^4\over 3GM^2}\left[{MR^2\over I}-2\beta \right]e^{-4.8\Delta R/R}. \end{equation} \tag{ 44 }$$

Applying the experimental constraints to S_v and L, and using Equation (27) together with Monte Carlo sampling, we find $\Delta I/I=0.037^{+0.012}_{-0.007}$ for a 1.4 M_☉ star. For comparison, the approximation suggested by Link et al. (1999), in which the exponential factor in Equation (44) is replaced by

$$\begin{equation} \left[1+{2p_t\over n_sm_bu_t}\left({Rc^2\over GM}-2\right)\left({Rc^2\over GM}+7\right)\right]^{-1}, \end{equation} \tag{ 45 }$$

and I by Equation (28), yields [MR²/I − 2β] = (7/2)[1 − 7β/4] and $\Delta I/I=0.042^{+0.016}_{-0.008}$ for a 1.4 M_☉ star.

The inferred values of ΔI/I from observations, as discussed above, appear to be consistent with this estimate. However, this explanation has been called into doubt (Andersson et al. 2012; Chamel 2012) due to entrainment of superfluid neutrons in the crust, which raises their effective masses to several times their bare mass. As a consequence, the inferred lower limit to ΔI/I would have to be larger by a similar factor, so that

$$\begin{equation} \Delta I/I\gtrsim 0.07 \end{equation} \tag{ 46 }$$

would be required to explain glitches. Only very low mass neutron stars, with M ≲ 1.2 M_☉, could satisfy this criterion, calling the standard model for glitches into question.

Link (2013) has pointed out, however, that I in Equations (39) and (46) does not necessarily represent the total moment of inertia of the star. Although it is usually assumed that the crust couples to the entire core of the star, this need not be the case. Link (2013) argues that the timescale for the outer core to decouple from the crust is of the order of weeks to years, and could correspond to the observed post-glitch relaxation timescale, which may simply represent the dynamical recovery of the outer core. If, for example, the crust couples to a fraction of the core, the inferred value of ΔI/I is decreased by this fraction, perhaps a factor of two or more. If this is the case, observations would regain consistency with theoretical expectations for 1.4 M_☉ and, perhaps, more massive stars.

The restrictions on the density dependence of the symmetry energy from experiments and neutron matter theory expectations not only establish important constraints on neutron star radii, but, coupled with the discovery of a nearly 2 M_☉ neutron star (Demorest et al. 2010), considerably reduce uncertainties in predictions for the properties of high-density matter. The results of Steiner et al. (2010, 2013) suggest that even at 6n_s the pressure can be estimated to within a 1σ error of about 35%. The maximum mass constraint has also been used to limit the maximum strangeness (e.g., hyperon) fraction in a neutron star interior: Weissenborn et al. (2012) estimate this to be less than 20% even if the maximum mass of a purely hadronic star is about 3.2 M_☉, which is the largest mass allowed by causality if a 1.4 M_☉ star has R ⩽ 12 km (Lattimer 2012). The strangeness content would be smaller for less stiff hadronic EOSs and smaller maximum masses. This restriction in the strangeness content of neutron stars is consistent with studies using the chiral quark–meson coupling (Katayama et al. 2012) and relativistic mean-field (Bednarek et al. 2012) models.

Furthermore, many studies rule out even the possibility of hybrid stars (those with a mixed phase of quarks and hadrons; e.g., Dexheimer et al. 2013). The study of Endo (2011) found no stable hybrid stars with maximum masses greater than 1.8 M_☉. In general, the phase boundary can be modeled with either a Maxwell or Gibbs construction, but if the surface tension of quark matter is sufficiently large, Gibbs transitions are not allowed. Many studies indicate that if the phase transition is of the Maxwell type, the onset of quarks immediately destabilizes the star so that the existence of a hybrid star is not possible.

However, some studies do permit hybrid stars that can satisfy the minimum maximum mass constraint. Lenzi & Lugones (2012) find stable hybrid configurations using the Nambu–Jona–Lasinio (NJL) model with the relatively stiff (i.e., large values of L) GM1 and TM1 nucleon forces. However, when hyperons are included, the ranges of parameter space allowing stable hybrid stars are severely restricted, and the maximum mass can barely exceed 2 M_☉. However, they did not investigate stability using nucleon forces with values of L in the range constrained by experiments, which lowers the permitted maximum masses. Weissenborn et al. (2011) have shown that the existence of hybrid quark stars with a mixed phase and having a sufficiently large maximum mass is strongly dependent upon the assumed hadronic EOS, and is limited to small, specific ranges of quark matter parameters. They found a marked sensitivity to the stiffness of the EOS as determined by the radii of moderate (1.4 M_☉) stars, so that less stiff equations of state lead to smaller permitted parameter ranges for stable hybrid stars. Nevertheless, Weissenborn et al. (2011) only considered hadronic interactions with relatively large values of L (>100 MeV). Alford et al. (2013) also explored the effects of phase transitions in hybrid stars, but also included the dependence on L in their study. They demonstrated that L values in the experimentally preferred range severely restrict models capable of satisfying the observed minimum maximum mass constraint. Their models are parameterized by the transition density n_crit, the magnitude of the phase transition, Δ(ne) = μΔn, and the sound speed $c_s^2$, which is assumed constant for n > n_crit. For example, with Δn = 0, n_t = 2n_s, and $c_s^2=c^2$, L ≃ 40 MeV could still result in maximum masses of order 2.2 M_☉. However, increasing Δ(ne) to 0.2(ne)_crit and n_crit to 4n_s reduces the maximum mass to less than 2 M_☉. And if $c_s^2=1/3$, the value expected asymptotically for quark matter, no value of Δ(ne) can satisfy the maximum mass constraint for n_crit ≳ 1.5n_s.