Combined Constraints on the Equation of State of Dense NeutronRich Matter from Terrestrial Nuclear Experiments and Observations of Neutron Stars
Abstract
Within the parameter space of equation of state (EOS) of dense neutronrich matter limited by existing constraints mainly from terrestrial nuclear experiments, we investigate how the neutron star maximum mass M, radius km and tidal deformability of canonical neutron stars all together constrain the EOS of dense neutronrich nucleonic matter. While the 3D parameter space of (curvature of nuclear symmetry energy), and (skewness of the symmetry energy and EOS of symmetric nuclear matter, respectively) are narrowed down significantly by the observational constraints, more data are needed to pin down the individual values of , and with quantified uncertainties. The largely controls the maximum mass of neutron stars. While the EOS with is sufficiently stiff to support neutron stars as massive as 2.37 M, to support the hyperthetical ones as massive as 2.74 M (composite mass of GW170817) requires to be larger than its currently known maximum value of about 400 MeV and beyond the causality limit. The upper limit on the tidal deformability of from the recent observation of GW170817 is found to provide upper limits on some EOS parameters consistent with but far less restrictive than the existing constraints of other observables studied.
1 Introduction
What is the nature of neutron stars and dense nuclear matter? To answer this question has been a longstanding and shared goal of both astrophysics and nuclear physics. The fundamental importance and broad impacts of answering this question have been well documented in the literature. It is a major scientific thrust for many major research facilities in astrophysics (National Academies, 2011) and nuclear physics (National Academies, 2012), such as, various advanced Xray satellites and earthbased large telescopes, the Neutron Star Interior Composition Explorer (NICER), various gravitational wave detectors and all advanced radioactive beam facilities being built around the world. In particular, to answer this question has been identified as a major goal in both the U.S. 2015 Long Range Plan for Nuclear Sciences (U.S. LRP, 2015) and the Nuclear Physics European Collaboration Committee (NuPECC) 2017 Long Range Plan (NuPECC LRP, 2017). However, answering this question accurately is very challenging. It is well known that properties of neutron stars are determined by the Equation of State (EOS) of neutronrich matter over a large density range from zero to about ten times normal nuclear matter density. Unfortunately, even within the simplest model considering particles only, the corresponding EOS of neutron star matter is still poorly known, not to mention possibly other particles and various phase transitions predicted to occur in the core of neutron stars.
To summarize what terrestrial experiments and nuclear theories have taught us about the EOS of neutronrich matter near the saturation density of symmetric nuclear matter (SNM), we first recall that the specific energy in asymmetric nucleonic matter (ANM) can be approximated parabolically in isospin asymmetry as
(1) 
where is the specific energy in SNM and is the symmetry energy. Around , the and can be Taylor expanded up to as
(2)  
(3) 
in terms of several EOS characteristic parameters: the incompressibility and skewness of SNM as well as the slope , curvature and skewness of the symmetry energy. While the above Taylor expansions are not used in our study for neutron stars, they provide asymptotic boundary conditions at for the parameterizations we shall use to describe the EOS of neutron star matter.
Generally speaking, the EOS of neutronrich nucleonic matter remains very uncertain mostly because of the poorly known nuclear symmetry energy especially at suprasaturation densities (Li et al., 2014). Nevertheless, thanks to the hard work of many people in both astrophysics and nuclear physics over many years, much progress has been made in constraining the EOS of neutronrich nucleonic matter. In particular, various analyses of terrestrial nuclear experiments and astrophysical observations have constrained the , and in reasonably small ranges around MeV (Shlomo et al., 2006; Piekarewicz, 2010), MeV and MeV. Especially worth noting, the quoted most probable values of and are based on surveys of 53 analyses of different kinds of terrestrial and astrophysical data available up to Oct. 2016 (Li & Han, 2013; Oertel et al., 2017). Moreover, extensive analyses of heavyion reactions at intermediate and relativistic energies, especially various forms of nucleon collective flow and kaon production, have provided a reasonably tight constraining band for the EOS of SNM up to about , see, e.g., ref. Danielewicz et al. (2002). However, the coefficients characterizing the highdensity behavior of neutronrich matter, such as the , and are only loosely known to be in the range of MeV, MeV, and MeV mostly based on analyses of terrestrial nuclear experiments and energy density functionals (Tews et al., 2017; Zhang et al., 2017), respectively.
While continuous efforts have been made to constrain the EOS of dense neutronrich matter using heavyion reactions which may involve rare isotopes with large neutron/proton ratios, presently no concensus has been reached yet from analyzing limited data available (Li, 2017). On the other hand, properties of neutron stars and events involving them, such as the mass, radius, moment of inertia, quadrupole deformation, pulsing frequency (Hessels et al., 2006), cooling curve (Yakovlev et al., 2001; Page et al., 2006), frequencies and damping times of various oscillating modes, spin parameter of pulsars, as well as the strain amplitude and phase evolution of gravitational waves from inspiraling neutron star binaries all depend significantly on the EOS of neutronrich matter, see, e.g., refs. (Lattimer & Prakash, 2016; Oertel et al., 2017; Özel & Freire, 2016; Watts et al., 2016; Grigorian et al., 2016; Newton et al., 2014) for recent reviews. While the observational data of neutron star properties are relatively limited so far, they also provide stringent constraints on the EOS and guide theories of dense neutronrich matter. In particular, the observed masses around 2M for the two most massive pulsars J16142230 (Demorest et al., 2010) and J0348+0432 (Antoniadis et al., 2013) restrict mostly the stiffness of dense SNM EOS, while the radii of neutron stars are known to be more sensitive to the symmetry energy around (Lattimer & Prakash, 2000, 2001). While much progress has been made in measuring the radii of neutron stars, because of the great difficulties involved especially in determining accurately the distance and modeling reliably the spectrum absorptions with different atmosphere models, the reported radii normally suffer from relatively large uncertainties. In fact, some radii extracted from different analyses and observations are still controversial, see, e.g., ref. (Miller & Lamb, 2016) for a recent review. Since we are not in a position to make any judgement on the reliability of any astrophysical observations, in our analysis here we use as an example the radius of M canonical neutron stars () in the range of km inferred from analyzing quiescent lowmass Xray binaries in ref. (Lattimer & Steiner, 2014). Moreover, we also use the upper limit of the dimensionless tidal deformation from the recent observation of GW170817 by the LIGO+Virgo Collaborations (Abbott et al., 2017; LIGO & Virgo, 2017).
The first discovery of a neutron star merger using multiple messengers further signifies and sets an excellent example of combining and crosschecking multiple probes of dense neutronrich matter on earth and in heaven. Given the aforementioned constraints on the EOS parameters of neutronrich matter mostly based on terrestrial nuclear laboratory experiments, here we study how the astrophysical observations of M, km and all together constrain the highdensity EOS in a way consistent naturally with the existing constraints from terrestrial nuclear experiments. For this purpose, fixing the , and at their most probably values mentioned earlier, we explore the intersections of constant surfaces with , km, km, and , respectively, in the 3dimensional (3D) parameter space of , and . The 3D parameter space allowed by all three observational constraints are identified. Moreover, in constructing the EOS of neutron star matter, the crustcore transition density and pressure have to be calculated consistently. While effects of the magnitude and slope of symmetry energy at on the crustcore transition have been extensively studied in the literature, effects of the curvature are less known. We therefore will also explore contours of the and in the versus plane (2D). The significant role of the is clearly revealed. Furthermore, within the currently known uncertainty ranges of and , by setting in Eqs. (2) and (3) as often done in the literature, we also explore how/what the same three astrophysical constraints may teach us about the EOS of dense neutronrich matter in the parameter plane. In both the 3D and 2D model frameworks, the upper limit of the tidal deformability from the recent observation of GW170817 is found to provide upper limits on some EOS parameters consistent with but less restrictive than the existing constraints on them. Overall, while combing exiting constraints from both terrestrial nuclear experiments and astrophysical observations allows us to limit significantly the EOS parameter space of highdensity neutronrich matter, data of more independent observables are needed to pin down the individual values of , and with quantified uncertainties.
This paper is organized as follows. The construction of the EOS of neutron star matter is presented in Section 2. The Section 3 is devoted to constraining the EOS of dense neutronrich nucleonic matter with astrophysical observations first in the 3D space of , and , then in the 2D plane of with . Finally, a summary is given in Section 4.
2 Modeling the equation of state of neutron star matter
The focus of this study is on the EOS of dense neutronrich matter. We shall thus adopt the NV EOS (Negele & Vautherin, 1973) as that for the inner crust and the BPS EOS (Baym et al., 1971) for the outer crust while focusing on the EOS of the dense core in a large 3D and 2D parameter space. In this section, we shall first investigate the crustcore transition density and pressure using a thermal dynamical approach. Our main goal is to examine effects of the symmetry energy, especially its curvature on the transition point. We will then discuss how we construct the EOS for the core. For the purposes of this work, it is sufficient to use the simplest model for nonrotating and chargeneutral neutron stars consisting of only particles at equilibrium. As we mentioned earlier, essentially all available manybody theories using various interactions have been used to predict both the EOS of SNM and symmetry energy from low to suprasaturation densities. Given their widely different predictions especially at suprasaturation densities, we try to extract parameters charactering the EOS of dense neutronrich matter from the astrophysical observations using as little as possible predictions of any particular manybody theory and/or interaction. However, we must ensure that the EOS parameters asymptotically become naturally the ones extracted from terrestrial nuclear experiments near .
To avoid misleading the reader, here we discuss in more detail the relationship between the Taylor expansions of Eqs. (1) to (3) used to characterize the EOS near the saturation density in both nuclear theory and in analyzing terrestrial experiments and the parameterizations that we shall use in modeling the EOS of neutron star matter. We also address the question why we can’t just adopt the widely used isospinindependent multiparameters polytropic EOSs, see, e.g., refs (Topper, 1964; Butterworth, 1976; Read et al., 2009), for the purposes of our work here. Near and , any nuclear energy density functional can be Taylor expanded using Eqs. (1) to (3). The relevant expansion coefficients are determined by the via for the symmetry energy, for the slope parameter of , and for the incompressibility of SNM and the curvature of , as well as and for the skewness of and , respectively. The Taylor expansions become progressively inaccurate for large densities and do not converge when . Therefore, for describing neutron star matter, we parameterize the and in exactly the same form as in Eqs. (1) to (3) but their coefficients are no longer given by the above expressions from any energy density functional . Instead, we treat them as unknown parameters to be extracted from astrophysical observations. More specifically, while the NV+BPS EOSs are used for the crust, starting from the crustcore transition density we parameterize the and as
(4) 
(5) 
using the parameters and . These parameterizations naturally become the Taylor expansions in the limit of . As parameterizations, mathematically they can be used at any density without the convergence issue associated with the Taylor expansions. Any parameterized EOS has to satisfy all the constraints (asymptotic boundary conditions) for the EOS near from terrestrial experiments. The above parameterizations facilitate comparisons with the terrestrial constraints around . When the above parameterizations are used at high densities for modeling the core EOS of neutron stars, the parameters in the above two equations are not necessarily to be the same as the Taylor expansion coefficients in Eqs. (1) to (3). While both mathematically and physically, the and are expected to be very close, the highorder coefficients (parameters) in the Taylor expansions (parameterizations) may be very different. Unfortunately, very little is known about these parameters from experiments/observations and the model predictions for the EOS of dense neutronrich matter diverge. Since the uncertainty ranges of the Taylor expansion coefficients , and are already so large, it is then meaningful to use their uncertainty ranges as the starting/reference ranges in our search for the highdensity EOS parameters and from studying properties of neutron stars.
With the above definitions and explanations, we now introduce/discuss a convention/simplification that we shall use in the following. It is consistent with the existing convention widely used in the literature for the Eq (1). The latter is often referred as the empirical parabolic law/approximation of the EOS of isospin asymmetric nuclear matter (Bombaci & Lombardo, 1991). Namely, it is well known in nuclear physics that the Eq (1) has the dual meanings of being a Taylor expansion in in the limit of on one hand, and on the other hand being a parameterization when it is used in very neutronrich matter where . Adopting this convention and be consistent with that used for the Eq (1), albeit confusing without the above explanations, it is really not necessary to write both the Taylor expansions of Eqs. (1) to (3) and the parameterizations of Eqs. (4) and (5) as they have the same form although different meanings for the coefficients/parameters involved. In the following discussions, when the Eqs. (1) to (3) are used in describing terrestrial nuclear EOS or predictions of nuclear energy density functional theories, they are Taylor expansions near and . On the other hand, when they are used in describing the EOS of neutron stars, they are parameterizations to be constrained by astrophysical observations. Then, all three equations in Eqs. (1) to (3) have the dual meanings in describing the EOS of isospin asymmetric matter encountered in both terrestrial experiments and neutron stars within broad ranges of and without the convergence problem.
In the sense that because there is no reliable theory for the EOS of neutronrich matter significantly above the saturation density one has to use parameterizations to describe the EOS of neutron star matter, the spirit of the parameterizations of Eqs. (4) and (5) is essentially the same as many other parameterizations used in the literature, see, e.g., refs. Gandolfi et al. (2009, 2012); Steiner et al. (2016). However, instead of requiring indirectly the involved parameters to reproduce the known properties of symmetric nuclear matter and the symmetry energy near through sometimes complicated equations, we use directly the coefficients of the Eqs. (2) and (3). These coefficients as the asymptotic boundary conditions of the EOS are themselves known near and either through experiments or converged predictions of many reliable models, see, e.g., refs. (Danielewicz et al., 2002; Shlomo et al., 2006; Chen et al., 2009; Steiner et al., 2010; Piekarewicz, 2010; Khan et al., 2012; Dutra et al., 2012, 2014; Li & Han, 2013; Colò et al., 2014; Cai & Chen, 2014; Zhang et al., 2017). In fact, parameterizations similar to the Eqs. (2) and (3), albeit often truncated at the first order in density for and second order for , i.e., using and only, have already been used successfully in studying various properties of neutron stars, see, e.g., refs.(Oyamatsu & Iida, 2007; Sotani et al., 2012). Also, a similar approach of describing approximately the EOS of dense neutronrich matter was recently proposed in ref. (Margueron et al., 2017a) and successfully used in studying properties of both neutron stars (Margueron et al., 2017b) and finite nuclei (Chatterjee et al., 2017) with Bayesian perspectives. With the above cautions, conventions and justifications, we shall use the parameterized EOS for the core of neutron stars as described in more detail in section 2.3.
The multiparameters polytropic EOSs, see, e.g., refs (Topper, 1964; Butterworth, 1976; Read et al., 2009) are widely used in modeling the core EOS of neutron stars at suprasaturation densities. Why do not we simply follow this popular approach? This is mainly because the polytropes used so far depend only on the nucleon number (energy) density but explicitly independent of the isospin asymmetry . To extract information about the highdensity symmetry energy from neutron stars, we need to know explicitly the underlying isospin composition of the pressure/density. For example, to our best knowledge, all Bayesian analyses done so far using masses and radii of neutron stars have adopted the polytropes. They infer only the total pressures at a few fiducial high densities without any information about the isospin content of the matter at those densities. Indeed, to invert mathematically the TOV equation one only needs to know the total pressure as a function of density regardless of its microphysics composition. However, to understand the entire structure, composition and cooling mechanism of neutron stars, it is necessary to know the isospin dependence of its EOS from the crust to the dense core consistently. Since the pressure has a term proportional to where the isospinasymmetry profile in neutron stars at equilibrium is determined uniquely by the as we shall discuss in more detail below, by directly parameterizing the and , separately, and then evaluate the resulting pressure as a function of both and , we shall obtain information about the highdensity underlying the pressuredensity relation . While directly parameterizing the do not facilitate the extraction of such information about the highdensity symmetry energy from astrophysical observations.
2.1 Sampling the density dependence of nuclear symmetry energy
First of all, we illustrate the broad variation of by varying independently the coefficients in Eq. (3) within their known uncertainty ranges. It is well known that the isospin asymmetry in neutron stars at a given baryon density is uniquely determined by the in Eq. (1) through the charge neutrality and equilibrium conditions as we shall recall more formally in section 2.3. Generally speaking, because of the term in the EOS, a higher value of will lead to a smaller at equilibrium. As quantitative examples, shown in Figure 1 are the and as functions of reduced density by varying only one coefficient each time while fixing all others: (a) , 50, 60, 70, and 80 MeV, (b) , 300, 200, 100, 0, and 100 MeV, and (c) , 0, 200, and 400 MeV. As their names indicate, the slope , curvature and skewness of symmetry energy play different roles and in order become increasingly more important at higher densities. Obviously, variations of them within their currently known uncertainty ranges allow us to sample very different behaviors of the and the corresponding .
It is worth noting that some combinations of the parameters lead to a decreasing that may even become negative at high densities. As summarized earlier in Szmaglinski et al. (2006) and reviewed very recently in Li et al. (2018), such kind of supersoft at high densities was predicted in a number of theoretical calculations using various interactions. At very high densities, when the shortrange repulsive tensor force due to the meson exchange makes the EOS of SNM increase faster with density than that of pure neutron matter where the tensor force is much weaker, the decreases or even becomes negative at high densities (Pandharipande et al., 1972; Wiringa et al., 1988; Li et al., 2008). To our best knowledge, such a seemingly unusual highdensity behavior of the is not excluded by neither any known fundamental physics principle nor experiments/observations so far. Possible consequences of such kind of symmetry energies are discussed in more detail in refs. (Szmaglinski et al., 2006; Li et al., 2018) and references therein. In fact, EOSs with a supersoft can still support massive neutron stars if the SNM parts of the EOSs are sufficiently stiff even without the help from the new light mesons proposed and/or possible modified strongfield gravity for massive objects, see, e.g., refs. (Krivoruchenko et al., 2009; Wen et al., 2009; Lin et al., 2014; Jiang et al., 2015). Interestingly, while not completely settled yet (Li, 2017), there are indeed some circumstantial evidences from intermediaterelativistic energy heavyion collisions indicating that the may become supersoft above about (Xiao et al., 2009). Currently, devoted efforts are being made by the intermediaterelativistic heavyion reaction community to pin down the highdensity behavior of nuclear symmetry energy, see, e.g., refs. (Li et al., 2014; Russotto et al., 2016; Xu et al., 2016; Tsang et al., 2017; Trautmann and Wolter, 2017).
As we shall discuss in detail next, the around the crustcore transition is mostly controlled by the and parameters when the is fixed at its most probable empirical value of 31.6 MeV. The variation of from 40 to 80 MeV and from 400 to 100 MeV allows us to sample the usual behavior of predicted by various nuclear manybody theories in the subsaturation density region and within their known empirical constraints at .
2.2 Imprints of the density dependence of nuclear symmetry energy on the crustcore transition in neutron stars
Although the crust of neutron stars contributes only a small fraction of the total mass and radius, it plays an important role in various astrophysical phenomena (Baym et al., 1971, 1971b; Pethick & Ravenhall, 1995; Link et al., 1999; Lattimer & Prakash, 2000, 2001; Steiner et al., 2005; Lattimer & Prakash, 2007; Chamel & Haensel, 2008; Sotani et al., 2012; Newton et al., 2013; Pons et al., 2013; Piekarewicz et al., 2014; Horowitz et al., 2015). Critical to many effects of the crust is the transition density between the inner crust and the outer core of neutron stars. Previous studies have found consistently that the transition density is very sensitive to the density dependence of nuclear symmetry energy. In particular, the role of the slope parameter has been extensively studied, see, e.g., refs. (Douchin & Haensel, 2000; Kubis, 2004, 2007; Lattimer & Prakash, 2007; Oyamatsu & Iida, 2007). Often the studies employ predictions of a particular nuclear manybody theory where the values of and are normally correlated. Here we shall first study the individual roles of the and in determining the corecrust transition properties, then contours of the transition density and pressure in the versus plane. Finally, effects of the correlation based on the systematics from analyzing over 500 nuclear energy density functions are examined.
In the present study, we employ the thermodynamical approach (Kubis, 2004, 2007; Lattimer & Prakash, 2007) to estimate the crustcore transition point where the uniform matter becomes unstable against being separated into a mixture of single nucleons and their clusters. The method is known to slightly overestimate the transition density compared to the dynamical approach (Xu et al., 2009; Ducoin et al., 2011; Providência et al., 2014) but sufficiently good for the purposes of this work. Specifically, the transition density is determined by the vanishing effective incompressibility of neutron star matter at equilibrium under the charge neutrality condition (Kubis, 2004, 2007; Lattimer & Prakash, 2007), i.e.,
(6) 
This approach has been used extensively in the literature to locate the transition point using various EOSs (see, e.g., refs. Xu et al., 2009; Ducoin et al., 2011; Providência et al., 2014; Routray et al., 2016). Enclosed in the bracket of the above expression for are the firstorder and secondorder derivatives of the symmetry energy, i.e., quantities directly determining the and . It is thus necessary and interesting to first explore separate roles of the latter on the transition density. Shown in Figure 2 are the transition density as functions of with different values of in the windowa and with different values of in the windowb, respectively. It is clearly seen that the changes more dramatically with the variation of than in their respective uncertainty ranges. This is mainly because the last two terms in the expression for largely cancel out, leaving the curvature of dominate. In addition, the value of is already relatively well constrained in a smaller range than the , making the variation of with look weaker.
Next, we examine the transition density and pressure by varying both the and within their uncertainty ranges continuously. Shown in the two windows of Figure 3 are contours of constant transition densities and pressures in the plane, respectively. For transition densities larger than fm, the required increases monotonically with , while different behaviors are observed for lower values of . The lowest transition density about fm appears around the boundary corner at MeV and MeV. Different from the contours of constant , for a fixed , the required always increases linearly with before reaching the boundary along the line . The latter is used as a limit in exploring properties of neutron stars in the EOS parameter space.
As indicated earlier, in using the Eqs. (2) and (3) we assume the coefficients are independent and intend to constrain them directly from observations using as little as possible predictions of any particular manybody theory. Thus, in determining the crustcore transition point for constructing the EOS of neutron stars, we freely vary the and within their uncertainty ranges specified earlier. Nevertheless, theoretically predicted values of the and are often correlated when model ingredients and/or interactions are varied. Thus, for a consistency check we also study how the predicted correlation between and may limit the transition point. As discussed by many people in the literature (see, e.g., refs. Farine et al., 1978; Ducoin et al., 2011; Providência et al., 2014; Pearson et al., 2014; Danielewicz & Lee, 2009; Chen et al., 2009; Vidaña et al., 2009; Ducoin et al., 2011; Mondal et al., 2017; Tews et al., 2017), based on the systematics of many predictions using various manybody theories and interactions, an approximately universal and linear correlation exists between the and . For example, using totally over 500 energy density functions including 263 Relativistic Mean Field (RMF) models, HartreeFok calculations using 240 Skyrme (Dutra et al., 2012, 2014) as well as some realistic and Gogny interactions, Mondal et al. (2017) found the following correlation
(7) 
Using essentially the same sets of energy density functionals but requiring fm, MeV, MeV, and MeV, Tews et al. (2017) rejected some of the energy functionals. They found the following correlation using the remaining 188 Skyrme and 73 RMF interactions
(8) 
where the and are error bars including 68.3% and 95.4% of the accepted EOSs, respectively. In Figure 3, the above two correlation functions are shown with the white and red dashed lines, separately. The 68.3% uncertainty range ( MeV) in Eq. (8) is indicated by the red dotted lines while that of Eq. (7) is too small to be shown here. While the parameterizations in Eqs. (7) and (8) are consistent, their mean values have slightly different slopes because of the different selection criteria used. Nevertheless, if the uncertainty range of Eq. (8) is enlarged from 68.3% to 95.4% of the accepted EOSs, the parameterization in Eq. (7) can then be fully covered by that in Eq. (8). Using the above two correlations, the transition density and pressure are then restricted to be about fm and MeV fm, respectively. These values are consistent with the crustcore transition properties often used in the literature. To this end, especially since some of the apparent correlations among EOS parameters from model calculations may be spurious (Margueron et al., 2017a), it is worth noting that while the correlation from the systematics of over 500 energy density functions is very useful for the consistency check, it is still necessary and important to determine the individual values of and from experiments/observations. In our following calculations, we thus use consistently the crustcore transition density and pressure by varying independently the and values within their respective uncertain ranges without using any of the about two correlation functions.
2.3 The core EOS of neutron stars
For completeness and the ease of our discussions in the following, we first recall here the formalism for calculating the EOS in the cores of neutron stars. The total energy density of charge neutral matter at equilibrium can be written as
(9) 
where and are the energy density of baryons and leptons, respectively. The can be calculated from
(10) 
where the specific energy of baryons is given in Eq. (1) and is the average rest mass of nucleons. The from the noninteracting Fermi gas model can be expressed as () (Oppenheimer & Volkoff, 1939)
(11) 
with
(12) 
and
(13) 
The chemical potential of particle can be calculated from
(14) 
The isospin asymmetry and relative particle fractions at different densities in neutron stars are obtained through the equilibrium condition and the charge neutrality condition The pressure of the system can be calculated numerically from
(15) 
The above expressions allow us to calculate the core EOS which is connected smoothly at the corecrust transition point to the NV EOS (Negele & Vautherin, 1973) for the inner crust followed by the BPS EOS (Baym et al., 1971) for the outer crust.
As mentioned earlier, the EOS parameters , and near are relatively well determined. To investigate how/what highdensity EOS parameters are constrained by the three astrophysical observations considered in this work, we construct the EOS of neutron star matter by varying the poorly known , and characterizing the EOS of dense neutronrich nucleonic matter. In principle, all coefficients used in Eqs. (2) and (3) should be varied simultaneously within a multivariant Bayesian inference (Steiner et al., 2010; Raithel et al., 2016; Margueron et al., 2017b). Such a study is in progress. In the present work, we shall perform traditional analyses first in the 3D parameter space spanned by , and while fixing all other parameters at their currently known most probable values. Equivalent to reparameterizing the EOS of SNM and with less parameters as often done in the literature, or expanding the Eqs. (2) and (3) only up to , we shall also explore the EOS in the 2D parameter space of and by setting while keeping other parameters at their known most probable values. The two cases studied here are similar in spirit to using different numbers of piecewise polytropes or parameters to model the EOS of dense neutronrich matter. Naturally, the values of the parameters involved may be different in the two cases, but they should all asymptotically approach the same existing constraints on them near .
Certainly, there have been continuous efforts in both astrophysics and nuclear physics to constrain the EOS parameters in both cases. For example, from the pressure of SNM constrained by nucleon collective flow data in relativistic heavyion collisions (Danielewicz et al., 2002), a constraint of MeV was obtained by Cai & Chen (2014). Combining it with the mass of neutron star PSR J0348+0432, they further narrowed it down to MeV. By analyzing Xray bursts, Steiner et al. (2010) extracted a value of or MeV assuming a photospheric radius of or , respectively. All of these constraints on overlap but have different uncertainty ranges. Similarly, the and have not been well constrained either by any experiments/observations so far. Nevertheless, the systematics of over 500 RMF and SHF energy density functionals indicates the following range: MeV (Dutra et al., 2012, 2014), MeV and MeV, respectively (see, e.g., Chen et al., 2009; Dutra et al., 2012, 2014; Colò et al., 2014; Zhang et al., 2017). Thus, we adopt these ranges for the parameters , and to be consistent with both existing experimental and theoretical findings.
Within the above uncertainty ranges of the EOS parameters, the pressure in neutron stars can be varied within the shaded band shown in the left window of Figure 4. Its upper and lower limit is obtained by using the parameter set of MeV, MeV, MeV and MeV and the set of MeV, MeV, MeV and MeV, respectively. The individual roles of these parameters are examined by varying them independently in the right window. As one expects, the variation of is most effective in modifying the pressure at suprasaturation densities.
3 Constraining the EOS of dense neutronrich matter with observed properties of neutron stars
With the EOSs prepared in the way described above, the mass (M)radius (R) relationship of neutron stars is obtained by solving the TolmanOppenheimerVolkov (TOV) equations (Tolman, 1934; Oppenheimer & Volkoff, 1939)
(16) 
(17) 
where is the gravitation constant, is the light speed and is the gravitational mass enclosed within a radius . The dimensionless tidal deformability is related to the Love number , neutron star mass M and radius R via
(18) 
The is determined by the EOS thorough a differential equation coupled to the TOV equation (Hinderer, 2008; Hinderer et al., 2010). More details about the formalism and code used in this work to calculate the can be found in, e.g., Fattoyev et al. (2013, 2014).
Within the 3D parameter space of , and and the 2D parameter space of and under the conditions discussed in the previous section, using the observational data of M, km and , we study how/what highdensity EOS parameters are constrained in the following two subsections, separately.
3.1 Observational constraints in the EOS parameter space
To be clear, we first emphasize again that in this case the following parameters are fixed at their currently known most probable values based on previous systematic surveys as discussed earlier: MeV, MeV and MeV. We explore constraints on the EOS in the 3D parameter space with the crustcore transition density consistently determined and the condition that . Technically, in exploring the 3D parameter space in three loops we change the to reproduce a specific observational data at given values of and which are varied independently. Shown in Figure 5 are the constant surfaces of neutron stars’ maximum mass of M (green), radius of km (magenta) and km (yellow) as well as the upper limit of the dimensionless tidal deformability (orange) of canonical neutron stars, respectively. For clarity, a causality surface limiting M M is not shown.
For ease of the following discussions, it is worth recalling first that, as shown in Eq. (1), the EOS of SNM , the symmetry energy and the isospin asymmetry at equilibrium are the three quantities together determining the total pressure in neutron stars. More specifically, the total pressure is proportional to . While the term is controlled by the parameter with a fixed , the and are determined by the and parameters when the is fixed. Moreover, the symmetry energy contribution to the total pressure is weighted by . When the is softer with smaller or negative and values, the system is more neutronrich as shown in Figure 1. In particular, for extremely small and values (e.g., MeV and MeV), the becomes negative and the reaches its maximum of 1 at high densities. Then, the necessary contribution to the pressure from the term will require a large value to support massive neutron stars. At the other extreme, however, when both the and are strongly positive (e.g., MeV and MeV at the bottom left corner), the symmetry energy is superstiff and the is very small as shown in Figure 1. The required to support massive neutron stars is then very small.
It is interesting to note several major features in this rather informationrich 3D plot summarizing very extensive calculations. Let us first focus on the constant surface of M. From the top right to the bottom left corner, the required first decreases quickly and then stays almost a constant with the increasing values of and from negative to positive. This feature is completely understandable based on the discussions in the previous paragraph. Namely, near the upper right corner, the symmetry energy is supersoft and the resulting is close to 1. The weight of the symmetry energy contribution to the pressure is significant while the value is small and may even be largely negative. To support neutron stars with the maximum mass of M, the value of has to be highly positive to compensate the small or overcome the possibly negative contribution from the symmetry energy. However, near the bottom left corner where the symmetry energy is superstiff, the resulting at equilibrium becomes so small such that the symmetry energy contribution to the pressure is strongly suppressed. The necessary value of is therefore small and the constant surface of M becomes rather flat at small values. More quantitatively, the required minimum value of is about MeV at MeV and MeV.
For a comparison and see more clearly the role of parameter in determining the masses of neutron stars, a constant surface of M (sky blue) corresponding to the composite mass of GW170817 is also shown. It is seen that the overall shapes of the constant surfaces of M and 2.74 M are rather similar. To support such a hypothetical neutron star would require a value above 400 MeV beyond its current uncertain range and the causality limit. While the fate of GW170817 is not completely determined yet, several analyses using observations of GW170817 combined with theories/simulations and the causality limit under some caveats have placed the upper bounds on neutron star masses in the range of M (Margalit & Metzger, 2017; Shibata et al., 2017; Rezzolla et al., 2018; Ruiz et al., 2018), see, e.g., ref. (Radice et al., 2018) for a very recent review. For instance, Margalit & Metzger (2017) used electromagnetic constraints on the remnant imposed by the kilonova observations after the merger and the gravitational wave information predicted a maximum mass of M with 90% confidence. To support neutron stars with such masses, the just needs to be slightly positive well within its current uncertain range. Thus, once the maximum mass is pinned down, it will put a stringent upper limit on the parameter.
In addition, it is also known that quadrupole deformations of neutron stars depend on the symmetry energy (Krastev et al., 2008a, b; Fattoyev et al., 2013, 2014, 2017; Krastev et al., 2018). Interestingly, Abbott et al. (2017) inferred that the dimensionless tidal deformability of GW170817 has an upper limit of at the 90% confidence level for the lowspin prior. However, it is seen in Figure 5 that the constant surface of (orange) locates far outside the constant surface of km. Thus, limits on the highdensity EOS parameters from the constraint alone are presently much looser than the radius constraint extracted from analyzing the Xray data. Nevertheless, the expected detection of gravitational waves from a large number of neutron star mergers has the potential to improve the situation.
Now, let us move to the EOS constraints from the radii of neutron stars. The radius of canonical neutron stars is known to depend strongly (weakly) on the nuclear symmetry energy (EOS of SNM) (Li & Steiner, 2006). It is thus not surprising that the two constant surfaces of radius at km and km for canonical neutron stars are essentially vertical in Figure 5, indicating a weak dependence on the as one expects. Indeed, they have significant dependences on both the and as indicated by the separation between the two constantradius surfaces. More quantitatively, the required values of in the constant surfaces of km and km decrease continuously with increasing and . The two surfaces can be approximately described by and , respectively. With a fixed value of , the nuclear pressure becomes stronger with increasing values of and . Thus, the constant surface of km is on the left (having stiffer symmetry energies) of the km surface.
As indicated by the arrows in Figure 5, the space enclosed by the three constant surfaces of M and km are the EOS parameter space allowed by the astrophysical observations of the maximum mass and radii of neutron stars. The space is constrained from the bottom by the maximum mass, its intersections with the two limits on the radii set the boundary lines restricting the and at small values around MeV. At higher values, the EOS parameter space is mainly bounded by the two radius constraints. To show the accepted EOS parameters more clearly, the allowed regions in the versus planes with , , , and MeV, respectively, are shown in Figure 6. The boundaries of these allowed regions in the planes are shown in Figure 7. The regions insides the lines are allowed for a given specified. The shadowed regions are excluded values of and for all values. More specifically, the boundary of the right shadowed region can be divided into two parts based on the slopes. The upper one with a negative slope is obtained by requiring the crustcore transition pressure to always stay positive, i.e, . It can be described approximately by the expression MeV. The lower one with a positive slope is obtained by the intersection line between the surfaces of M and km. It can be fitted by the expression MeV. This boundary sets an upper limit for at about MeV. The left shadowed region is excluded by the intersection line of km and MeV in Figure 5. It can be fitted by the expression MeV.
Overall, within the framework of our analyses, using the maximum mass of neutron stars as well as the upper and lower limits of the radii of canonical neutron stars, the three EOS parameters are only limited in a space shown in Figure 5 not completely closed with its boundaries partially given in Figure 7. Obviously, data of more independent observables from either or/both astrophysics and nuclear physics are needed to determine the individual values of , and .
3.2 Observational constraints in the EOS parameter plane
Within the currently known uncertainty ranges of and , one can parameterize the EOS with less parameters by setting in Eqs. (2) and (3) as often done in the literature. Then, the two most poorly known parameters are the and . The latter is known to be around MeV as mentioned earlier. In this 2D model framework, here we explore how/what the same three astrophysical constraints may teach us about the EOS of dense neutronrich matter.
Shown in Figure 8 are contours of constant maximum masses of neutron stars in the plane. The red, magenta, blue, and cerulean shadows represent regions where M, , km, and , respectively. In the white excluded region, the crustcore transition pressure . It is seen that the maximum mass of neutron stars increases with increasing as one expects. However, it is rather insensitive to the variation of in the region considered. Again and obviously, the tidal polarizability is much less restrictive than the radius constraint of km. It is also seen that the km constraint is covered by the constraint M. Consequently, the area bounded by the curves of M and km is the region allowed by the astrophysical observations considered. In this region, the maximum value of is about 52 MeV consistent with that found in the 3D analyses in the previous subsection. Also in this region, the upper limit of the maximum mass is about 2.37 M reached at MeV. Probably incidentally, the latter is also currently the most probable value of based on the surveys of 53 analyses of existing data (Li & Han, 2013; Oertel et al., 2017). Interestingly, the observed upper limit of the maximum mass is in good agreement with the findings by Fryer et al. (2015), Lawrence et al. (2015) and Alsing et al. (2017). They found that the upper limit of the maximum mass of neutron stars is between and M. However, it is necessary to caution here that we have taken in the 2D study. As we have discussed in the previous subsection, the value of affects significantly the maximum mass of neutron stars. Thus, without more precise knowledge about the parameter, the absolute maximum mass of neutron stars can not be pinned down. Based on our 2D model analyses here, neutron stars more massive than 2.37 M would require a positive value of .
4 Summary and outlook
In summary, within both the 2D and 3D EOS parameter spaces limited by the existing constraints from terrestrial nuclear experiments, we studied how the astrophysical observations of M, km and all together constrain the EOS parameters of dense neutronrich nucleonic matter. We also investigated effects of the curvature of nuclear symmetry energy on the crustcore transition in neutron stars. The consistently calculated transition density in the plane is used in constructing the EOS of neutron star matter from the surface to the core. The is found to affect significantly the crustcore transition density and pressure. Fixing the , and at their most probably values determined mainly by terrestrial nuclear experiments, we explored the intersections of constant surfaces with , km, km, and , respectively, in the 3D parameter space of , and . The 3D parameter space narrowed down significantly by the observational constraints is clearly identified. This helps guide nuclear theories for dense neutronrich matter and related studies in terrestrial experiments. However, to pin down the individual values of , and , data of additional independent observables from either astrophysical observations and/or laboratory experiments are needed. In particular, the skewness parameter of SNM largely controls the maximum mass of neutron stars. The 2D EOS with is found sufficiently stiff to support neutron stars as massive as 2.37 M, while to support a hypothetical neutron star as massive as 2.74 M (composite mass of GW170817) would require to be larger than about 400 MeV beyond its known limit. In both the 2D and 3D model frameworks considered in this work, the upper limit of the tidal deformability from the recent observation of GW170817 is found to provide upper limits on some EOS parameters consistent with but less restrictive than the existing constraints. In particular, its constraints on the symmetry energy parameters are far less restrictive than the observation of km from analyzing the Xray data.
While the analyses and results presented here are useful in their own rights, some aspects of our present work should be improved in future studies. In particular, more data and/or a better approach are necessary to determine precisely individual values of the highdensity EOS parameters with quantified uncertainties. In the era of gravitational wave astronomy accompanied by the planned new experiments using advanced radioactive beam facilities around the world, we are very hopeful that more data will come soon. Moreover, even with the limited data available, more quantitative information about the highdensity EOS parameters may be obtained by using a more robust statistical approach. For example, we fixed the parameters describing the EOS near the saturation density at their most probable values mostly from terrestrial experiments. By doing so we eliminated any possible correlation between these (fixed) parameters and the highdensity parameters , and that we want to determine from astrophysical observations. A Bayesian analysis where the prior distribution of parameters encodes what is already known and the new observational data refines such prior knowledge is particularly well suited for inferring the posterior probability density distributions of all the EOS parameters. Such a study using the masses and radii of 14 neutron stars extracted from Chandra XRay observations, the maximum masses of neutron stars from GW170817 and earlier observations of neutron stars as well as all the information from terrestrial experiments about the loworder EOS parameters will be carried out and reported elsewhere.
References
 Abbott et al. (2017) Abbott, B. P., et al. 2017, Phys. Rev. Lett., 119, 161101
 Alsing et al. (2017) Alsing, J., Silva, H. O., & Berti, E. 2017, arXiv:1709.07889
 Antoniadis et al. (2013) Antoniadis, J., et al. 2013, Science, 340, 448
 Baym et al. (1971) Baym, G., Pethick, C. J., & Sutherland, P. 1971, APJ, 170, 299
 Baym et al. (1971b) Baym, G., Bethe, H. A., & Pethick, C. J. 1971b, Nucl. Phys. A, 175, 225
 Bombaci & Lombardo (1991) Bombaci, I., & Lombardo, U. 1991, Phys. Rev. C44, 1892
 Butterworth (1976) Butterworth, E. M. 1976, APJ, 204, 561
 Cai & Chen (2014) Cai, B. J., & Chen, L. W. 2017, Nucl. Sci. Tech., 28, 185
 Chamel & Haensel (2008) Chamel, N., & Haensel, P. 2008, Living Rev. Relat., 11, 10
 Chatterjee et al. (2017) Chatterjee, D., Gulminelli, F., Raduta, Ad. R. & Margueron, J. 2017, Phys. Rev. C, 96, 065805.
 Chen et al. (2009) Chen, L. W., Cai, B. J., Ko, C. M., Li, B. A., Shen, C., & Xu, J. 2009, Phys. Rev. C, 80, 014322
 Colò et al. (2014) Colò, G., Garg, U., & Sagawa, H. 2014, Eur. Phys. J. A, 50, 26
 Danielewicz et al. (2002) Danielewicz, P., Lacey, R., & Lynch, W. G. 2002, Science, 298, 1592
 Danielewicz & Lee (2009) Danielewicz, P., & Lee, J. 2009, Nucl. Phys. A, 818, 36
 Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081
 Douchin & Haensel (2000) Douchin, F., & Haensel, P. 2000, Phys. Lett. B, 485, 107
 Ducoin et al. (2011) Ducoin, C., Margueron, J., Providência, C., & Vidaña, I. 2011, Phys. Rev. C, 83, 045810
 Dutra et al. (2012) Dutra, M., Loureno, O., Martins, J. S. S., Delfino, A., Stone, J. R., & Stevenson, P. D. 2012, Phys. Rev. C, 85, 035201
 Dutra et al. (2014) Dutra, M., Loureno, O., Avancini, S. S., Carlson, B. V., Delfino, A., Menezes, D. P., Providncia, C., Typel, S., & Stone, J. R. 2014, Phys. Rev. C, 90, 055203
 Fattoyev et al. (2013) Fattoyev, F. J., Carvajal, J., Newton, W. G., & Li, B. A. 2013, Phys. Rev. C, 87, 015806
 Fattoyev et al. (2014) Fattoyev, F. J., Newton, W. G., & Li, B. A. 2014, Eur. Phys. J. A, 50, 45
 Fattoyev et al. (2017) Fattoyev, F. J., Piekarewicz, J. & Horowitz, C.J., 2017, arXiv:1711.06615.
 Farine et al. (1978) Farine, M., Pearson, J. M., & Rouben, B. 1978, Nucl. Phys. A, 304, 317
 Fryer et al. (2015) Fryer, C. L., Belczynski, K., RamirezRuiz, E., Rosswog, S., Shen, G., & Steiner, A. W. 2015, APJ, 812, 24
 Gandolfi et al. (2009) Gandolfi, S., Illarionov, A. Yu., Schmidt, K. E., Pederiva, F., & Fantoni, S. 2009, Phys. Rev. C, 79, 054005
 Gandolfi et al. (2012) Gandolfi, S., Carlson, J., & Reddy, S. 2012, Phys. Rev. C, 85, 032801
 Grigorian et al. (2016) Grigorian, H., Voskresensky, D. N. & Blaschke, D. 2016, Eur. Phys. J. A, 52, 67
 Hessels et al. (2006) Hessels, J. W. T., et al. 2006, Science, 311, 1901
 Hinderer (2008) Hinderer, T. 2008, APJ., 677, 1216
 Hinderer et al. (2010) Hinderer, T., Lackey, B. D., Lang, R. N., & Read, J. S. 2010, Phys. Rev. D, 81, 123016
 Horowitz et al. (2015) Horowitz, C. J., et al. 2015, Phys. Rev. Lett., 114, 031102
 Jiang et al. (2015) Jiang, W.Z., Li, B.A. & Fattoyev, F.J. 2015, Eur. Phys. J. A, 51, 119
 Khan et al. (2012) Khan, E., Margueron, J., & Vidaña, I. 2012, Phys. Rev. Lett., 109, 092501
 Krastev et al. (2008a) Krastev, P. G., Li, B.A. and Worley, A. 2008, Phys. Lett. B, 668, 1
 Krastev et al. (2008b) Krastev, P. G., Li, B.A. and Worley, A. 2008, APJ, 676, 1170
 Krastev et al. (2018) Krastev, P. G., Li, B.A., 2018, arXiv:1801.04620
 Kubis (2004) Kubis, S. 2004, Phys. Rev. C, 70, 065804
 Kubis (2007) Kubis, S. 2007, Phys. Rev. C, 76, 025801
 Krivoruchenko et al. (2009) Krivoruchenko, M.I., Simkovic, F. & Faessler, A. 2009, Phys. Rev. D, 79, 125023
 Lattimer & Prakash (2000) Lattimer, J. M., & Prakash, M. 2000, Phys. Rep., 333, 121
 Lattimer & Prakash (2001) Lattimer, J. M., & Prakash, M. 2001, APJ, 550, 426
 Lattimer & Prakash (2007) Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109
 Lattimer & Steiner (2014) Lattimer, J. M., & Steiner, A. W. 2014, Eur. Phys. J. A, 50, 40
 Lattimer & Prakash (2016) Lattimer, J.M. & Prakash, M. 2016, Phys. Rep, 621, 127
 Lawrence et al. (2015) Lawrence, S., Tervala, J. G., Bedaque, P. F., & Miller, M. C. 2015, APJ, 808, 186
 Li & Steiner (2006) Li, B.A. & Steiner, A.W. 2006, Phys. Lett. B, 642, 436
 Li et al. (2008) Li, B.A., Chen, L.W. & Ko, C.M. 2008, Phys. Rep. 464, 113
 Li & Han (2013) Li, B. A., & Han, X. 2013, Phys. Lett. B, 727, 276.
 Li et al. (2014) Li, B. A., Ramos, À., Verde G., & Vidana, I. (Eds.) 2014, Eur. Phys. J. A, 50, 9
 Li (2017) Li, B. A. 2017, Nuclear Physics News, 27, 7
 Li et al. (2018) Li, B. A., Cai, B. J., Chen, L. W., & Xu, J. 2018, Progress in Particle and Nuclear Physics, 99, 29
 LIGO & Virgo (2017) LIGO Scientific Collaboration, Virgo Collaboration, et al. 2017, ApJL, 848, L12
 Link et al. (1999) Link, B., Epstein, R., & Lattimer, J. M. 1999, Phys. Rev. Lett., 83, 3362
 Lin et al. (2014) Lin, W.K., Li, B.A., Chen, L.W., Wen, D.H. & Xu, J. 2014, J. of Phys. G, 41, 075203
 Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017,APJ, 850, L19
 Margueron et al. (2017a) Margueron, J., Hoffmann Casali, R., Gulminelli, F., 2018, Phys. Rev. C 97, 025805
 Margueron et al. (2017b) Margueron, J., Hoffmann Casali, R., Gulminelli, F., 2018, Phys. Rev. C 97, 025806
 Mondal et al. (2017) Mondal, C., Agrawal, B. K., De, J. N., Samaddar, S. K., Centelles, M., & Viñas, X. 2017, Phys. Rev. C, 96, 021302(R)
 Miller & Lamb (2016) Miller, M.C. & Lamb, F. K. 2016, European Physical Journal A 52, 63.
 Negele & Vautherin (1973) Negele, J. W., & Vautherin, D. 1973, Nucl. Phys. A, 207, 298
 Newton et al. (2013) Newton, W. G., Gearheart, M., & Li, B. A. 2013, APJ Suppl. Ser., 204, 9
 Newton et al. (2014) Newton, W.G., et al., 2014, Euro Phys. J. A50, 41.
 Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007
 Oppenheimer & Volkoff (1939) Oppenheimer, J., & Volkoff, G. 1939, Phys. Rev., 55, 374
 Oyamatsu & Iida (2007) Oyamatsu, K., & Iida, K. 2007, Phys. Rev. C, 75, 015801
 Özel & Freire (2016) Özel, F. & Freire, P. 2016, Annual Reviews of Astronomy and Astrophysics, 54, 401
 Page et al. (2006) Page, D., Geppert, U., & Weber, F. 2006, Nucl. Phys. A, 777, 497
 Pandharipande et al. (1972) Pandharipande, V. R., et al. 1972, Phys. Lett. B, 39, 608
 Pearson et al. (2014) Pearson, J. M., Chamel, N., Fantina, A. F., Goriely, S. 2014, Eur. Phys. J. A, 50, 43
 Pethick & Ravenhall (1995) Pethick, C. J., & Ravenhall, D. G. 1995, Ann. Rev. Nucl. Part. Sci., 45, 429
 Piekarewicz (2010) Piekarewicz, J. 2010, J. Phys. G, 37, 064038
 Piekarewicz et al. (2014) Piekarewicz, J., Fattoyev, F. J., & Horowitz, C. J. 2014, Phys. Rev. C, 90, 015803
 Pons et al. (2013) Pons, J. A., Viganó, D., & Rea, N. 2013, Nat. Phys., 9, 431
 Providência et al. (2014) Providência et al., C. 2014, Eur. Phys. J. A, 50, 44
 Radice et al. (2018) Radice, D., Perego, A., Zappa, F., & S. Bernuzzi 2018, APJ, 852, L29
 Raithel et al. (2016) Raithel, C. A., Özel, F., & Psaltis, D. 2016, APJ, 831, 44; ibid, arXiv:1704.00737
 Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Phys. Rev. D, 79, 124032
 Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, APJ, 852, L25
 Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, PRD, 97, 021501(R)
 Routray et al. (2016) Routray, T. R., Viñas, X., Basu, D. N., Pattnaik, S. P., Centelles, M., Robledo, L. B., & Behera, B. 2016, J. Phys. G, 43, 105101
 Russotto et al. (2016) P. Russotto et al. (ASYEOS Collaboration) 2016, Phys. Rev. C, 94, 034608
 Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012
 Shlomo et al. (2006) Shlomo, S., Kolomietz, V. M., & Colò G. 2006, Eur. Phys. J. A, 30, 23
 Sotani et al. (2012) Sotani, H., Nakazato, K., Iida, K., & Oyamatsu, K. 2012, Phys. Rev. Lett., 108, 201101
 Steiner et al. (2005) Steiner, A. W., Prakash, M., Lattimer, J. M., & Ellis, P. J. 2005, Phys. Rep., 410, 325
 Steiner et al. (2010) Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, APJ, 722, 33
 Steiner et al. (2016) Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2016, Eur. Phys. J. A, 52,18
 Szmaglinski et al. (2006) Szmaglinski, A., Wojcik, W. & Kutschera, M. 2006, Acta Phys. Polon. B, 37, 277

National Academies (2011)
The National Academies, New Worlds, New Horizons in Astronomy and Astrophysics, 2011,
https://www.nap.edu/catalog/12951/newworldsnewhorizonsinastronomyandastrophysics 
National Academies (2012)
The National Academies, Nuclear Physics: Exploring the Heart of Matter, 2012,
https://www.nap.edu/catalog/13438/nuclearphysicsexploringtheheartofmatter 
U.S. LRP (2015)
The 2015 U.S. Long Range Plan for Nuclear Science, Reaching for the Horizon,
https://science.energy.gov/~/media/np/nsac/pdf/2015LRP/2015_LRPNS_091815.pdf 
NuPECC LRP (2017)
The Nuclear Physics European Collaboration Committee (NuPECC) Long Range Plan 2017, Perspectives in Nuclear Physics,
http://www.esf.org/fileadmin/user_upload/esf/NupeccLRP2017.pdf.  Tews et al. (2017) Tews, I., Lattimer, J. M., Ohnishi, A., & Kolomeitsev, E. E. 2017, APJ, 848, 105
 Tsang et al. (2017) Tsang, M.B. et al. (SRIT Collaboration) 2017, Phys. Rev. C, 95, 044614
 Tolman (1934) Tolman, R. C. 1934, Proc. Natl. Acad. Sci. U.S.A., 20, 3
 Topper (1964) Topper, R. F. 1964, APJ, 140, 434
 Trautmann and Wolter (2017) Trautmann, W. & Wolter, H.H. 2017, arXiv:1712.03093
 Vidaña et al. (2009) Vidaña, I., Providência, C., Polls, A., & Rios, A. 2009, Phys. Rev. C, 80, 045806
 Watts et al. (2016) Watts, A.L. et al, 2016, Rev. Mod. Phys., 88, 021001
 Wen et al. (2009) Wen, D.H., Li, B.A., & Chen, L.W. 2009, Phys. Rev. Lett., 103, 211102
 Wiringa et al. (1988) Wiringa, R. B., Fiks, V., & Fabrocini, A. 1988, Phys. Rev. C, 38, 1010
 Xu et al. (2009) Xu, J, Chen, L. W., Li, B. A., & Ma, H. R. 2009, APJ, 697, 1549
 Xu et al. (2016) Xu, J. et al. (Transport Model Comparison Project) 2016, Phys. Rev. C, 93, 044609
 Xiao et al. (2009) Xiao, Z.G., Li, B.A., Chen, L. W., Yong, G. C. & Zhang, M. 2009, Phys. Rev. Lett., 102, 062502
 Yakovlev et al. (2001) Yakovlev, D. G., Kaminker, A. D., Gnedin, O. Y., & Haensel, P. 2001, Phys. Rep., 354, 1
 Zhang et al. (2017) Zhang, N. B., Cai, B. J., Li, B. A., Newton, W. G., & Xu, J. 2017, Nucl. Sci. Tech., 28, 181.