Spin-polarized Neutron Matter, the Maximum Mass of Neutron Stars, and GW170817

Neutron-star observations, such as the recent detection of two merging neutron stars in the gravitational-wave (GW), gamma-ray, and electromagnetic (EM) spectra (Abbott et al. 2017a, 2017b, 2017c, 2019), designated as GW170817, GRB 170817A, and AT 2017gfo, respectively, provide crucial constraints on the equation of state (EOS) of dense strongly interacting matter. The EOS is a key quantity for astrophysics and sensitively depends on strong interactions. Hence, it connects astrophysical observations to laboratory experiments at rare isotope beam facilities for the most neutron-rich extremes, e.g., at the Radioactive Beam Factory (RIBF), Japan, and the future Facility for Rare Isotope Beams (FRIB) and Facility for Antiproton and Ion Research (FAIR) in the US and Germany. While there is a wealth of theoretical models for the EOS of neutron-star matter (see Hebeler et al. 2015; Lattimer & Prakash 2016; Özel & Freire 2016 for reviews), for densities beyond nuclear saturation density, ${n}_{\mathrm{sat}}\approx 0.16\,{\mathrm{fm}}^{-3}$, these models can only be confronted with a limited set of experimental data, e.g., from heavy-ion collisions (Danielewicz et al. 2002).

Neutron stars are the densest objects in the universe and probe the EOS up to several times saturation density. Neutron-star observations are therefore an ideal source of additional information that complement experimental data and provide powerful constraints for the EOS at higher densities (Hebeler et al. 2013; Tews et al. 2018a). The structure of a neutron star is described by the mass–radius (M–R) relation, which is an important observational quantity and in a one-to-one correspondence with the EOS; the M–R relation follows from the EOS by solving the Tolman–Oppenheimer–Volkoff equations (Tolman 1939; Oppenheimer & Volkoff 1939). Measuring the M–R relation, and therefore the EOS, observationally is however extremely difficult.

On the one hand, neutron-star masses can be determined very precisely for some neutron stars in binaries (Lattimer 2012). For example, the precise measurement of two-solar-mass neutron stars (Demorest et al. 2010; Antoniadis et al. 2013; Fonseca et al. 2016) established a robust and strong constraint on the EOS of strongly interacting matter, and implied that the EOS has to be sufficiently stiff at high densities to support neutron stars in that mass regime. This constraint was recently tightened by the observation of a ${2.14}_{-0.09}^{+0.10}\,{M}_{\odot }$ neutron star (Cromartie et al. 2020). In contrast to masses, radii are extremely difficult to measure because of a limited number of suitable neutron stars and large systematic and statistical uncertainties (Steiner et al. 2010). Recently, the NICER (Gendreau et al. 2012) mission was able to measure simultaneously the mass and radius of PSR J0030+0451, with a radius ${12.7}_{-1.2}^{+1.1}$ km (Riley et al. 2019) and ${13.0}_{-1.1}^{+1.2}$ km (Miller et al. 2019) at the 68% confidence level. These results are consistent with our present understanding of the nuclear EOS (Raaijmakers et al. 2019a). Future observations by the NICER and, e.g., eXTP missions (Watts et al. 2019) will improve this with target radius uncertainties of 5%–10%, corresponding to 1 km or better (Watts et al. 2016).

In this situation, the recent observation of a neutron-star merger (Abbott et al. 2017a, 2019) by the LIGO-Virgo collaboration, as well as many follow-up observations, have provided exciting additional insights. While the GW signal GW170817 has been used to constrain the radius of a typical 1.4 M_⊙ neutron star to be below 13.6 km (see also Annala et al. 2018; Most et al. 2018; Tews et al. 2018b), additional information can be obtained from the observed EM kilonova. Extracted ejecta properties disfavor a prompt collapse to a black hole, favoring a hypermassive neutron star supported by differential rotation as an immediate product of the merger. This object then collapsed to a black hole on the timescale of a few 100 ms, because a longer-lived supramassive (supported against collapse by uniform rotation) or stable neutron star would have been able to deposit large amounts of rotational energy into the ejecta, leading to the formation of an energetic relativistic jet, which was not observed (Margalit & Metzger 2017; Shibata et al. 2017; Rezzolla et al. 2018).

Based on the EM observation and the previously discussed scenario, one can propose limits on the maximum mass of non-rotating neutron stars, M_max. In general, a larger M_max leads to a larger maximum mass for a rotating neutron star and, therefore, a longer lifetime of the merger remnant. The absence of a prompt collapse requires M_max to be sufficiently large to stabilize the hypermassive neutron star, while the absence of a longer-lived remnant forces M_max to be sufficiently small (Margalit & Metzger 2017). Based on this reasoning, several estimates on the upper limit of M_max were proposed. From the energy deposited into the kilonova ejecta, Margalit & Metzger (2017) concluded M_max to be bounded by M_max ≤ 2.17 M_⊙, similar to Shibata et al. (2017) who found M_max = 2.2 ± 0.05 M_⊙. This constraint was recently updated to M_max ≤ 2.3  M_⊙ (Shibata et al. 2019). Rezzolla et al. (2018) used empirical relations between M_max and the maximum masses of uniformly rotating or differentially rotating neutron stars to conclude ${M}_{\max }\leqslant {2.16}_{-0.15}^{+0.17}\,{M}_{\odot }$, and Ruiz et al. (2018) found M_max = 2.16–2.28 M_⊙.

From the theoretical side, the range of predicted M_max varies over a much wider range. Model EOSs for astrophysical simulations typically have M_max ∼ 2–2.5 M_⊙ (Lattimer 2012), while EOSs based on modern nuclear forces at nuclear densities and general extrapolations to higher densities usually allow for a much wider range of pressure in neutron stars (Hebeler et al. 2013; Krüger et al. 2013; Annala et al. 2018; Tews et al. 2018b; Greif et al. 2019), and can support extreme values for M_max, limited only by M_max ≲ 4.0 M_⊙ (M_max ≲ 2.9 M_⊙) when a nucleonic EOS is considered up to n_sat (2n_sat) in Tews et al. (2018b). Typically, approaches using general extensions find M_max ranging between 2.9–3.2 M_⊙ (Hebeler et al. 2013; Krüger et al. 2013; Annala et al. 2018; Greif et al. 2019), based on different sets of calculations and other physically motivated constraints. These ranges are consistent with earlier findings of M_max ≲ 2.9 M_⊙, where a nucleonic EOS up to 2n_sat was combined with the stiffest possible EOS at higher densities (Nauenberg & Chapline 1973; Rhoades & Ruffini 1974), or of M_max ≲ 4.0 M_⊙ in a similar approach, but considering a nucleonic EOS to a lower density (Kalogera & Baym 1996).

In this work, we propose a novel theoretical conjecture that limits M_max ≤ 2.6–2.9 M_⊙, and is therefore relevant for all extension schemes discussed above. In particular, we propose that a possible phase transition from unpolarized neutron-star matter to spin-polarized neutron matter (SPM) provides constraints on the properties of neutron-star matter and excludes areas with high pressure. Our novel conjecture is only based on theoretical calculations of SPM, does not take any observational constraints on the maximum mass from GW170817 into account, and hence, is complementary to the observational conclusions.

We have investigated several EOSs for SPM and found that such a phase transition drastically softens the EOS, and therefore limits neutron-star masses. This observation remains true for all possible unpolarized EOSs, even for those with regions of drastic stiffening beyond nuclear saturation density. Only for the stiffest possible SPM EOS do we find stars that experience a spin-polarized phase in their core. However, in these cases, the resulting EOSs do not considerably increase the maximum neutron-star mass. Furthermore, most of the resulting EOSs are ruled out based on radius constraints from GW170817 (Abbott et al. 2017a, 2019; Annala et al. 2018; Most et al. 2018; Tews et al. 2018b; Raaijmakers et al. 2019b). Hence, we do not find any neutron stars consistent with astrophysical observations that exhibit a considerable domain of SPM in their core, in agreement with Vidana et al. (2002). Therefore, the onset of such a phase transition mainly determines the end of the stable branch of the neutron-star M–R relation, limiting the range for ${M}_{\max }$. We show our main findings in Figure 1 for SPM calculated with various approaches that give very consistent results and limit M_max to be below 2.6–2.9 M_⊙. Our results exclude neutron stars that explore extreme pressures (gray areas in Figure 1).

Zoom In Zoom Out Reset image size

We start from the general EOS extension of Tews et al. (2018a, 2018b). This family of EOSs is constrained at nuclear densities by microscopic calculations using local chiral effective field theory (EFT) interactions and precise Quantum Monte Carlo (QMC) methods, see Gezerlis et al. (2013, 2014), Tews et al. (2016), and Lynn et al. (2016) for details. Chiral EFT provides a systematic theory for nuclear forces, based on the symmetries of quantum chromodynamics, in terms of nucleon and pion degrees of freedom (Epelbaum et al. 2009; Machleidt & Entem 2011; Hammer et al. 2013). It explicitly includes long-range pion-exchange interactions, and parameterizes short-range interactions by a general operator basis whose low-energy couplings are fit to nucleon–nucleon (NN) scattering data as well as few-body systems. Chiral EFT naturally provides three-nucleon (3N) interactions which have been found to be extremely important for calculations of nuclear matter, while four-body interactions have been found to be very small (Krüger et al. 2013; Tews et al. 2013). Due to their systematic organization, chiral interactions can be systematically improved and enable theoretical uncertainty estimates. Being a momentum expansion, chiral EFT is limited to low momenta as explored in atomic nuclei, but it allows to constrain the EOS at nuclear densities (Tews et al. 2013; Hebeler et al. 2015; Gandolfi et al. 2019). In addition, QMC methods are among the most precise methods to solve the many-body problem (Carlson et al. 2015). Using chiral interactions as input, QMC methods have been used to calculate nuclei and neutron matter with great success (Lynn et al. 2016, 2019), which shows that microscopic calculations can connect the physics of nuclei with the astrophysics of neutron stars.

For the results in this paper, we have used the auxiliary-field diffusion Monte Carlo (AFDMC) method (Schmidt & Fantoni 1999). Together with local chiral EFT interactions, this approach can be applied to neutron matter up to densities around 2n_sat (Tews et al. 2018b). In Tews et al. (2018a, 2018b), we demonstrated how to obtain the neutron-star EOS from these calculations. To extend these calculations to higher densities explored in the neutron-star core, we have used an extension in the speed of sound, c_S, which allows to model the most general family of EOSs consistent with our nuclear-density results (see also Greif et al. 2019). We stress that this general extension scheme is independent of a particular choice of degrees of freedom at higher densities. Instead, it explores all allowed density dependencies that are consistent with microscopic calculations at nuclear density, that are causal, and that lead to a stable neutron star, i.e., a star with monotonously growing pressure. As a consequence, this extension scheme covers all possible EOS models, e.g., models with phase transitions, sudden stiffening of the EOS etc.

Using chiral EFT constraints up to ${n}_{\mathrm{sat}}$ and the c_S extension at higher densities, we sample tens of thousands of EOSs within the allowed EOS range. For this ensemble of EOSs, we find $8.4\,\mathrm{km}\leqslant {R}_{1.4}\leqslant 15.2\,\mathrm{km}$ for the radius of a typical 1.4 M_⊙ neutron star, R_1.4, and M_max ≤ 4.0 M_⊙. This upper limit results from the stiffest nuclear EOS consistent with local chiral EFT constraints at nuclear densities and the stiffest possible causal EOS at higher densities, and reduces to 2.9 M_⊙ if nuclear-physics input is considered up to 2n_sat, as discussed before. We show the most general EOS band consistent with nuclear-physics constraints up to n_sat as gray bands in Figure 1.

For each of the EOSs within the EOS band, we construct a new EOS that includes a phase transition to SPM. We sketch our construction in Figure 2. To obtain the EOS for SPM, we use three different calculations: AFDMC calculations using local chiral interactions to next-to-next-to-leading order (N²LO), many-body perturbation theory (MBPT) calculations from Krüger et al. (2015) based on chiral EFT interactions to next-to-next-to-next-to-leading order (N³LO), as well as Brueckner–Hartree–Fock (BHF) calculations (Vidana et al. 2002) based on the phenomenological Reid 93 and Nijmegen (Stoks et al. 1994) potentials. We show the results for SPM for these different calculations in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using AFDMC, we performed calculations of SPM at leading order (LO), next-to-leading order (NLO), and for two different Hamiltonians at N²LO, which differ by the choice of the 3N parameterization, see Lynn et al. (2016) for details. For each Hamiltonian, we estimate the uncertainties from the order-by-order convergence (Epelbaum et al. 2015), giving the bands shown in Figure 3 (a). Our AFDMC results are in good agreement with the findings of Riz et al. (2018), which used the same method and interactions, but did not estimate the EFT uncertainties. In the following, we will only use the more conservative ${V}_{E{\mathbb{1}}}$ results (with higher energies for SPM, see the red band in Figure 3), but note that both 3N parameterizations overlap well. For the ${V}_{E{\mathbb{1}}}$ parameterization, we give the EOS up to 2n_sat in Table 1.

Table 1.  EOS of SPM from the AFDMC Calculation with the N²LO ${V}_{E{\mathbb{1}}}$ Parameterization of Figure 3(a) up to 2n_sat

$n\,(\,{\mathrm{fm}}^{-3})$	E/N (MeV)	$\epsilon \,(\,\mathrm{MeV}\,{\mathrm{fm}}^{-3})$	$P\,(\,\mathrm{MeV}\,{\mathrm{fm}}^{-3})$
0.04	21.31(7)	38.41(1)	0.55(1)
0.08	33.72(11)	77.81(1)	1.83(5)
0.12	44.69(83)	118.0(1)	3.8(4)
0.16	55.1(2.5)	159.0(4)	6.5(1.4)
0.20	65.5(5.3)	200.8(1.1)	10.1(3.4)
0.24	75.8(9.4)	243.4(2.3)	14.5(6.9)
0.28	85.8(15.1)	286.8(4.2)	19.9(12.4)
0.32	95.5(22.5)	331.0(7.1)	26.3(20.3)

Download table as:  ASCIITypeset image

As indicated by the uncertainty bands in Figure 3, chiral interactions become less reliable with increasing density. This is especially true for local interactions employed in AFDMC, because they suffer from sizable local regulator artifacts (Dyhdalo et al. 2016; Tews et al. 2016; Huth et al. 2017). To estimate their impact, we also explore MBPT calculations with nonlocal chiral interactions (Krüger et al. 2015), which do not suffer from these additional regulator artifacts. We show the total MBPT uncertainty band from Krüger et al. (2015), which was obtained by studying several chiral N³LO interactions as well as variations of the cutoff and the 3N couplings. While the uncertainty band is not based on a systematic order-by-order study, the described uncertainty estimation is very reasonable at N³LO, and consistent with the AFDMC calculation, but considerably smaller at larger densities. Finally, we compare to results from BHF calculations based on the Nijmegen and Reid93 phenomenological interactions (Vidana et al. 2002). These calculations do not provide uncertainties but describe NN scattering data with high precision and are in good agreement with the chiral EFT results.

To extrapolate the chiral EFT results to higher densities, we fit the simple functional form (Gandolfi et al. 2009)

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\mathrm{pol}}}{N}(n)=a\cdot {\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{\alpha }+b\cdot {\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{\beta }\end{eqnarray} \tag{ 1 }$$

to the results. In particular, in the case of AFDMC, we fit this functional to the result of our calculation as well as to the upper and lower bounds of the uncertainty bands up to n_sat. For MBPT, we fit the functional to the bounds of the uncertainty band. In the case of the AFDMC results, we test the quality of the extrapolation by comparing it to the data points between n_sat and 2n_sat and find that the fit provides a reliable extrapolation to these higher densities. For the BHF results, we fit this functional to the individual results over the whole density range. Finally, in case the EOS for spin-polarized matter becomes acausal, we replace it by a causal EOS with c_s = c.

As shown in Figure 2 for given neutron-star and SPM EOSs, unpolarized matter is energetically favorable at low-energy densities but becomes less favorable than spin-polarized matter at higher densities for sufficiently stiff EOSs. The reason is that interactions in SPM tend to be weak and results are close to the free Fermi gas (Krüger et al. 2015), while interactions in unpolarized matter are much stronger and become increasingly repulsive. We then identify the phase transition between unpolarized matter and SPM by a Maxwell construction, i.e., by matching pressure and Gibbs energy or chemical potential, and construct a new EOS using the unpolarized EOS below and the polarized EOS above the phase transition (blue dashed line in Figure 2).

We emphasize that this neglects corrections to the spin-polarized matter EOS from protons (or other particle species), which are expected at the level of ≈10% given typical proton fractions.

To explore the sensitivity to the exact construction of the phase transition, we have explored additional EOSs where the phase transition is smeared out, similar to a Gibbs construction. Such a transition would appear due to the formation of a mixed phase if protons and electrons were included (Glendenning 1992; Heiselberg et al. 1993). Instead of enforcing ${P}_{\mathrm{tr}}^{\mathrm{polarized}}\,={P}_{\mathrm{tr}}^{\mathrm{unpolarized}}={P}_{\mathrm{pt}}$, and connecting the EOSs by a segment with c_s = 0, we construct EOSs with the unpolarized phase up to P_pt−ΔP and the polarized phase after P_pt + ΔP. These EOS segments are then connected by a smooth interpolation. The resulting M_max as a function of ΔP is shown in the inset of Figure 2 for the given EOS. Smearing out the phase transition lowers M_max because the EOS gets softened earlier. The maximum M_max is therefore found for a strong first-order phase transition resulting from the Maxwell construction. This is in good agreement with similar findings for phase transitions to quark matter of Bhattacharyya et al. (2010) and Wu & Shen (2019) and holds for all EOSs in our sample. Therefore, our conclusions are robust with respect to the properties of the transition.

We have repeated this construction for all unpolarized neutron-star EOSs in our original band and the different EOSs for SPM: the AFDMC result, its upper and lower bounds, the upper and lower MBPT bounds, as well as the two phenomenological EOSs. This leads to the results of Figure 1, where the hatched areas are given by the uncertainty bands of the SPM calculations in the corresponding panels of Figure 3 (or by the two different Hamiltonians in the BHF case). For the AFDMC calculations, we find $1.75\,{M}_{\odot }\leqslant {M}_{\max }\leqslant 2.93\,{M}_{\odot }$, with the centroid being at M_max = 2.29 M_⊙. For the MBPT calculations, we find 1.84M_⊙ ≤ M_max ≤ 2.59 M_⊙, and for the BHF calculations, M_max ≤ 2.61 M_⊙ or M_max ≤ 2.44 M_⊙. Except in the stiffest possible case, for the upper AFDMC bound, these findings are in very good agreement with each other. However, for this stiffest case the uncertainty is increased by regulator artifacts and most likely overestimated. We find that the predicted M_max is in very good agreement with inferences from the EM counterpart of GW170817 of Margalit & Metzger (2017), Shibata et al. (2017), Rezzolla et al. (2018), which are shown in Figure 1 as horizontal lines or bands.

Because a phase transition to SPM softens the EOS drastically, we find that it is unlikely that a neutron-star with a spin-polarized core exists in nature. Typically, we find the mass of the SPM domain to be ≤0.005 M_⊙, largely a result of numerical discretization artifacts. Only the stiffest possible spin-polarized EOSs can stabilize any star with spin-polarized matter in their core. In this case, we find that the mass of spin-polarized domain in the neutron-star core is ≤0.19 M_⊙. However, most of the resulting EOSs lead to neutron stars with R_1.4 ≥ 13.6 km, and are therefore ruled out by GW170817 (Abbott et al. 2017a, 2019; Annala et al. 2018; Most et al. 2018; Tews et al. 2018b). If we were to exclude these EOSs from consideration, we find that the mass of spin-polarized matter in the core is ≤0.02 M_⊙.

We have ignored the effects of magnetic fields, which could impact the EOS of spin-polarized matter if there is a net magnetization. We have also ignored a gradual polarization of neutron matter which, however, softens the EOS sooner, leading to a lower M_max. Hence, the case we investigated presents an upper limit on M_max due to a transition to SPM in the core. Finally, a similar effect on the EOS might occur due to a phase transition to deconfined quark matter; see e.g., category A in Alford et al. (2013). However, since quark-matter properties cannot be predicted from first principles in the density range of interest, in this case no strong constraint on M_max can be obtained.

We have investigated the impact of a phase transition from neutron-star matter to SPM and found that such a phase transition limits the pressure in the neutron-star core. Combining information from AFDMC calculations, as well as previous MBPT (Krüger et al. 2015) and BHF (Vidana et al. 2002) calculations limits M_max of neutron stars to lie below 2.6– 2.9 M_⊙, depending on the microscopic nuclear forces used, while significantly larger M_max of 3–4 M_⊙ could be reached without these constraints. These limits can be improved if the uncertainty in SPM calculations is reduced. The lower M_max, and in particular the result for the AFDMC calculations without uncertainty estimates, M_max = 2.29 M_⊙, are in very good agreement with recent constraints from the EM counterpart of GW170817 from Margalit & Metzger (2017), Shibata et al. (2017), and Rezzolla et al. (2018). Finally, we find that stars with a large spin-polarized domain in their core are ruled out by the radius constraint from GW170817.

We are grateful to C.J. Pethick for helpful discussions and J. Carlson and C. Wellenhofer for comments on the manuscript. This work was supported by the US Department of Energy, Office of Science, Office of Nuclear Physics, under Contract DE-AC52-06NA25396, the Los Alamos National Laboratory (LANL) LDRD program, the NUCLEI SciDAC program, and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 279384907—SFB 1245. Computational resources have been provided by Los Alamos Open Supercomputing via the Institutional Computing (IC) program, by the National Energy Research Scientific Computing Center (NERSC), which is supported by the U.S. Department of Energy, Office of Science, under contract DE-AC02-05CH11231, and by the Jülich Supercomputing Center.