Simultaneous Fitting of Neutron Star Structure and Cooling Data

Simultaneous Fitting of Neutron Star Structure and Cooling Data

Spencer Beloin    Sophia Han    Andrew W. Steiner    Khorgolkhuu Odbadrakh Department of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, USA Department of Physics and Astronomy, Ohio University, Athens, OH 45701, USA Department of Physics, University of California Berkeley, Berkeley, CA 94720, USA Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA Joint Institute for Computational Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA

Using a model for the equation of state and composition of dense matter and the magnitude of singlet proton superconductivity and triplet neutron superfluidity, we perform the first simultaneous fit of neutron star masses and radii determined from observations of quiescent low-mass X-ray binaries and luminosities and ages determined from observations of isolated neutron stars. We find that Vela is an outlier which strongly determines the values inferred for the superfluid/superconducting gaps, the neutron star radius, and the posterior mass distribution. We find, regardless of whether or not Vela is included in the analysis, that the threshold density for the direct Urca process lies between the central density of 1.7 and 2 solar mass neutron stars, but that 2 solar mass stars are unlikely to cool principally by the direct Urca process because of the suppression by neutron triplet superfluidity.

97.60.Jd, 95.30.Cq, 26.60.-c

Introduction: The determination of the equation of state (EOS) of dense matter at densities beyond the saturation density has been dominated by the analysis of photon-based observations which lead to determinations of neutron star masses and radii Cameron (1959); Lattimer and Prakash (2001); Steiner et al. (2010); Steiner et al. (2013a). The observation of gravitational waves from the binary neutron star merger GW170817 LIGO Scientific Collaboration and Virgo Collaboration (2017) by the LIGO-Virgo collaboration significantly changes this picture. Gravitational wave-based probes of neutron star structure have smaller systematic uncertainties than electromagnetic probes. Thus, depending on the mass distribution of progenitor masses of neutron stars which merge Fryer et al. (2015), gravitational wave observations can potentially drastically decrease the uncertainties in the EOS in the near future.

However, knowledge of the equation of state alone is insufficient to describe observed neutron star phenomenology. Neutron star cooling Yakovlev and Pethick (2004); Page et al. (2004) and pulsar glitches Piekarewicz et al. (2014) both depend strongly on the role played by superfluidity and superconductivity Page et al. (2014). On the other hand, neither neutron star mass and radius observations nor gravitational wave detections are likely to constrain superfluidity or superconductivity in the near future. Thus, a joint analysis of data that constrains the EOS and the nature of pairing in dense matter is required to refine our understanding of neutron star observations.

In addition, determination of the equation of state does not ensure knowledge of the composition of dense matter. This was demonstrated in Ref. Alford et al. (2005) where it was found that a mass-radius curve with neutrons and protons alone was nearly identical to that from a model containing deconfined quarks. Neutron star cooling, on the other hand, is sensitive to the neutron-to-proton ratio, since that controls the threshold for the direct Urca process Lattimer et al. (1991).

Before this work, technical limitations have forced quantitative analyses of neutron star masses and radii and neutron star thermal evolution to be almost entirely separate (for an exception, see Ref. Ho et al. (2015)). In this work, we perform the first combined quantitative analysis of mass and radius data and luminosity and age data for isolated neutron stars. In addition, we ensure that our model reproduces the structure of heavy nuclei. Since we allow the mass of each neutron star to vary, we can go beyond the “minimal cooling” model first envisaged in Ref. Page et al. (2004) by varying the EOS and including the direct Urca process.

Method: We presume that neutron stars contain only neutrons, protons, electrons, and muons, leaving the description of exotic hadrons or deconfined quarks to future work. We employ Bayesian inference, and all parameters below are taken to have uniform prior distributions. It is assumed that the data for each neutron star and all of the nuclear structure data are independent. The conditional probability (we will refer to this as the likelihood) will thus be a product of terms that we describe below.

Nuclei and the equation of state are modeled with a covariant field-theoretical Lagrangian from Ref. Steiner et al. (2005), based on earlier work in, e.g. Ref. Serot and Walecka (1986). In order to be as flexible as possible with regard to the description of high-density matter, the model in Ref. Steiner et al. (2005) (see Table 3) used 17 parameters. We intentionally employ a relatively large parameter space to ensure that models are not prematurely ruled out. We compute the structure of and in the Hartree approximation, and thus the first term in the likelihood is a Gaussian for the binding energy and charge radius of those two nuclei. We presume 2% uncertainties for the structure data, larger than the experimental accuracy, to allow for some systematics due to limitations in the Hartree approximation, similar to the method used in Refs. Steiner et al. (2005); Steiner et al. (2013b) (see also a similar calculation in Ref. Chen and Piekarewicz (2014)). Pairing is not included in the nuclear structure calculations. We constrain the nuclear incompressibility to be between 220 and 260 MeV Garg and Coló (2018). To be consistent with recent progress in nuclear theory, we restrict the symmetry energy at the nuclear saturation density, , to be less than 36 MeV and the slope of the symmetry energy at the nuclear saturation density, , to be less than 80 MeV Lattimer and Steiner (2014). A posterior correlation between and is automatically obtained from the fit to the nuclear structure data.

For the neutron star mass and radius data, we use the results obtained from the seven quiescent low-mass X-ray binary (QLMXB) neutron stars in the baseline data set from Ref. Steiner et al. (2018). The associated term in the likelihood function is constructed as in, e.g. Eq. 102 of Ref. Lattimer and Steiner (2014) (see also the more generic formalism developed in Ref. Steiner (2018)). A nuisance variable is required to parameterize the curve, and we use the neutron star mass for this purpose, giving 7 new parameters which are constrained to be larger than 1 solar mass. Each QLMXB (except for the neutron star in Cen) can have an atmosphere of either hydrogen or helium. We construct a new parameter, , for each neutron star ranging between 0 and 1. If is less than 2/3, then a hydrogen atmosphere is assumed, otherwise a helium atmosphere is assumed Steiner et al. (2018).

For the isolated neutron star cooling data, we use the luminosity and age data from Table I of Ref. Beloin et al. (2018), and the associated term in the likelihood function is based on Eq. 4 of Ref. Beloin et al. (2018). We do not include the neutron stars in the Cassiopeia A (“Cas A”) and HESS J1731347 supernova remnants. The thermal evolution of the neutron star in Cas A is still not well understood (see Refs. Page et al. (2011); Shternin et al. (2011); Bonanno et al. (2014); Elshamouty et al. (2013); Posselt et al. (2013)). XMMU J173203.3344518 (“J1732”) in HESS J1731347 may be described by either a carbon atmosphere or a hydrogen atmosphere Klochkov et al. (2015). Finally, the carbon-atmosphere stars were shown to not have as strong an impact on the analysis as the Vela pulsar (PSR B083345) does Beloin et al. (2018). Below we present results with and without Vela, due to its strong impact on the posteriors. Each of the 15 stars requires a mass parameter, and the 9 hydrogen-atmosphere stars may have a heavy- or light-element envelope, represented by the variable  Potekhin et al. (1997) varying between and . We assume spherically symmetric surface temperatures and that the magnetic field is small enough to be irrelevant for the cooling.

We use six parameters to describe proton superconductivity and neutron superfluidity as in Ref. Beloin et al. (2018). Each gap is described by a Gaussian function of the Fermi momenta, with height, centroid and width parameters denoted , and , respectively. To avoid double-counting, for both neutrons and protons, we constrain to be larger than the value of at the crust-core transition, and both and to be smaller than the value of in the core of the maximum-mass star. Note that we do not ensure that the superfluid or superconducting gaps are consistent with the Lagrangian used to describe the EOS of nucleonic matter. This kind of consistency is unnecessary, since the pairing interaction does not strongly impact the equation of state and our gaps will always be less than 1 MeV, which is relatively small compared to the nucleon Fermi energy.

In order to account for the fact that our Lagrangian may not fully describe nucleons in the neutron star core, we add a two-parameter polytrope to our EOS following Ref. Han and Steiner (2017), which begins at twice the nuclear saturation density and either softens or stiffens the EOS depending on the sign of the coefficient and the magnitude of the exponent. We automatically reject any models that are acausal (because of the additional polytrope), have a maximum mass less than 2 , or for which any neutron stars have a mass larger than the maximum mass of the EOS.

Our model has a total of 64 parameters, and traditional Bayesian inference applied to this problem requires the solution to the TOV equations, the nuclear structure calculations, and several cooling curves (to handle the variation in mass and envelope composition) at each point. To make it computationally tractable, we do not fully compute every point. Instead, we populate the parameter space with a library of exact calculations, and then use an emulator (Gaussian for the likelihood and inverse-distance-weighted for the output quantities) to perform the final inference.

Results: The posterior distributions for the neutron and proton density distributions are given in panel (a) of Fig. 1. Regardless of whether or not Vela is included, the charge radii is within 2% of the experimental value of 5.5 fm. On the other hand, our 95% credible range for the neutron skin thickness of is 0.15 to 0.20 fm (0.143 to 0.157 fm) without (with) Vela. When Vela is included, the value of (and thus the value of the skin thickness) is pushed towards smaller values in order to decrease the specific heat and increase the cooling (this is clear also in the luminosity-age curves below). Note that when is smaller, the isospin asymmetry in the central region increases in order to ensure that the total neutron and proton number is fixed Steiner et al. (2005).

Figure 1: All panels except (d) show results with and without Vela. For display purposes, simulations with and without Vela are fixed to have the same overall normalization. Panel (a): 95% credible intervals for the neutron and proton density profiles for Pb. Panel (b): 95% credible intervals for the mass-radius curves. Panel (c): histograms for the maximum value of the proton singlet and neutron triplet gaps. Panel (d): histograms for the other two gap parameters for neutrons and protons for the analysis excluding Vela. Panel (e): histograms for the central proton fraction of 1.4, 1.7, or 2.0  neutron stars. Panel (f): Posterior mass distributions for the seven QLMXB neutron stars or the 14 isolated cooling neutron stars (or 15 stars when Vela is included).

The distribution for mass-radius curves is plotted in panel (b) of Fig. 1. We find the radius of a 1.4 solar mass star to be between 12.4 and 13.1 km (12.35 to 12.65 km) without (with) Vela. When Vela is included, the smaller value of leads to smaller neutron star radii. These ranges are consistent with that obtained from QLMXBs alone (11.0 to 14.3 km) Steiner et al. (2018), and also with the radius inferred by the GW170817 merger (9.1 to 12.8 km) LIGO Scientific Collaboration and Virgo Collaboration (2018). There is some scatter in the radius limits due to small statistics.

The peak values of critical temperatures for the proton and neutron superfluid are presented in panel (c) of Fig. 1. As in Ref. Beloin et al. (2018), we find significantly different results depending on whether or not Vela was included: the gaps are significantly larger when Vela is included, in order to maximize the neutrino emissivity from the breaking and formation of Cooper pairs Flowers et al. (1976); Leinson and Perez (2006); Leinson and Pérez (2006); Kolomeitsev and Voskresensky (2008); Steiner and Reddy (2009); Page et al. (2009). Panel (d) shows the centroid and width parameters for the gaps without Vela. We find that is large with or without Vela, which effectively prevents the direct Urca process from occurring in the core. This is consistent with the observational data, as no confirmed isolated neutron stars have very low luminosities (see Figs. 2 and 3 below).

Panel (e) in Fig. 1 displays histograms for the proton fraction in the core of a 1.4, 1.7, and 2.0  neutron stars along with a vertical line at 11%, the threshold for non-superfluid neutrons and protons to participate in the direct Urca process. We find it unlikely that 1.7  neutron stars have proton fractions larger than 11%, and very likely that all 2.0  stars do. However, the neutron superfluid quenches the direct Urca process, even in 2.0  stars.

Panel (f) in Fig. 1 shows the mass posterior for the QLMXB and isolated neutron star (INS) data sets with Vela (dashed and dot-dashed lines) and without Vela (dotted and solid lines). Including Vela decreases the radius for the INSs, as shown above. Since the QLMXB observations tend to fix , decreasing tends to increase . When Vela is not included, probability distributions for both QLMXBs and INSs peak near 1.4 , similar to that observed in NS-NS and NS-WD binaries with accurate mass measurements Lattimer (2012).

In Fig. 2, we show the luminosity-age cooling curves excluding Vela. The bands refer to the 95% credible regions for the luminosity as a function of age. Note that while the gravitational mass is presumed to be constant over the time period shown, the amount of light elements in the envelope, parameterized by , may decrease over time as nuclear reactions fuse these light elements into heavier ones. Cooling curves with Vela included are plotted in Fig. 3. The uncertainty ranges are much smaller, demonstrating that a much smaller set of parameters is able to accommodate Vela’s small luminosity given its age.

Figure 2: The 95% credible intervals for the luminosity-age curves for different neutron star masses and different amounts of light elements in the neutron star envelope (parameterized by ), along with the luminosity and age values obtained from observations which were used in this work. Vela (labeled “0833”) is not included in the model results but the data point is included in the figure for comparison.
Figure 3: The 95% credible intervals for the luminosity-age curves for different neutron star masses and different amounts of light elements in the neutron star envelope (parameterized by ), along with the luminosity and age values obtained from observations which were used in this work. Vela (labeled “0833”) is included in the model results.

There are several other recent works describing the cooling of isolated neutron stars which analyze a smaller data set, and employ a smaller set of nucleon pairing models with a smaller set of EOS models, but often take into account exotic degrees of freedom. Ref. Raduta et al. (2018) studies stars containing hyperons and finds a significant contribution from hyperonic direct Urca processes. Ref. Negreiros et al. (2018) also studies stars containing hyperons and finds small neutron star radii and large masses are preferred. Ref. Potekhin, A. Y. and Chabrier, G. (2018) focuses on magnetic fields, which we do not include, and is thus able to include magnetars in their analysis. Ref. Ofengeim and Yakovlev (2017) contains a helpful analytical description of isolated neutron star cooling. Ref. Beznogov et al. (2018) uses a method similar to ours to analyze the cooling of J1732, and include the possibility of axion cooling. They apply a single EOS and find large proton and small neutron gaps.

Discussion: Field-theoretical Lagrangians give systematically different results than Skyrme models, in particular they tend to lead to larger neutron skin thicknesses and larger neutron star radii Steiner et al. (2005). However, Skyrme models are also likely to become acausal at densities below the central density of the maximum mass star Dutra et al. (2012). As also found for mass and radius data alone in Refs. Steiner et al. (2013a); Steiner et al. (2016), we do expect our results to be sensitive to our prior choices, particularly with regard to the equation of state.

We find that the most likely models selected by this inference are unlikely to be compatible with the cooling of the coldest transiently-accreting sources, such as SAX J1808.43658 Heinke et al. (2007); Han and Steiner (2017). The suppression of the direct Urca emissivity due to pairing found here cannot explain the low temperatures observed. This suggests that the best way forward is to simultaneously fit both isolated and accreting neutron stars, in addition to the information from neutron star mass and radius determinations. Initial progress on simultaneous fitting of these two data sets has been made in Refs. Beznogov and Yakovlev (2015a, b) using a handful of EOSs.

We thank Dany Page for publicly releasing his NSCool neutron star cooling code. We also thank Craig Heinke and Joonas Nättilä for helpful discussions. This work was supported by NSF grant PHY 1554876 and by the U.S. DOE Office of Nuclear Physics. S.H. is also supported by Chandra grant TM8-19002X, NSF grant PHY-1630782, and the Heising-Simons Foundation, Grant 2017-228. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) allocation PHY170048 supported by NSF grant number ACI-1548562 and computational resources from the University of Tennessee and Oak Ridge National Laboratory’s Joint Institute for Computational Sciences.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description