Symmetry energy of super-dense neutron-rich matter from integrating barotropic pressures in neutron stars and heavy-ion reactions

Within the minimum model of neutron stars (NS) consisting of neutrons, protons and electrons, a new approach is proposed for inferring the symmetry energy of super-dense neutron-rich nucleonic matter above twice the saturation density $\rho_0$ of nuclear matter directly from integrating iteratively barotropic pressures in both neutron stars and heavy-ion reactions. Simultaneously, the proton fraction of NSs at $\beta$ equilibrium is extracted as a function of baryon density from the same procedure. An application of this approach using the NS pressure from GW170817 and the pressure in cold symmetric nuclear matter (SNM) extracted earlier by analyzing nuclear collective flow data in relativistic heavy-ion collisions provides a useful constraining band for the symmetry energy above $2\rho_0$.

Motivation: A longstanding and shared goal of many astrophysical observations and terrestrial nuclear experiments is to understand the Equation of State (EOS) of dense neutron-rich nuclear matter. Indeed, much efforts have been devoted to extracting information about the EOS in both fields using various observatories and facilities during the last few decades. On one hand, by analyzing X-rays from individual neutron stars (NSs) and/or gravitational waves from their mergers, the pressure P NS as a function of energy density ǫ or baryon density ρ has often been extracted. However, generally no information about the composition of the source (e.g. the proton fraction in NSs) is extracted simultaneously. This is partially because normally isospin-independent polytropic EOSs are used in the analyses and a barotropic pressure is sufficient for solving the Tolman-Oppenheimer-Volkov (TOV) equations [1,2]. For example, the blue band shown in Fig.1 is the P NS (ρ) as a function of baryon density at 90% confidence level in canonical NSs at β equilibrium reported by the LIGO & VIRGO Collaborations from their analyses of the GW170817 event [3]. On the other hand, systematic analyses of terrestrial nuclear experiments over the last few decades have set a reasonably tight bound on the EOS of cold, symmetric nuclear matter (SNM) from sub-saturation densities to about four times the saturation density ρ 0 of nuclear matter [4,5,6,7]. As an example, shown in Fig.1 with the red band is the pressure P SNM (ρ) in SNM extracted from transport model analyses of nuclear collective flow [4] in heavy-ion collisions at beam energies from about 0.4 to 10 GeV/nucleon. Since there are confusions and incorrect statements in the literature about what EOS information can be extracted from heavy-ion reactions, we emphasize here that while heavy-ion reactions create hot matter, the underlying EOS of cold matter used as a basic input of transport model simulations of these reactions can be reliably extracted from comparing the simulations with experimental observables, see, e.g., refs. [4,8,9,10,11,12,13]. Of course, there are still many uncertainties involved in extracting the cold EOS from heavy-ion reactions. However, to our best knowledge, currently they are not more than those involved in extracting the EOS from gravitational waves emitted by merging NSs. In fact, both the P NS (ρ) and P SNM (ρ) shown in Fig. 1 have been used extensively in many different ways to constrain various aspects of nuclear many-body theories and interactions in the literature, see, e.g., refs. [14,15,16]. While both pressures still have rather large error bands, to our best knowledge, they represent the state of the art of the fields.   [3]. The red band is the pressure in cold symmetric nuclear matter extracted from analyzing nuclear collective flow [4] in heavy-ion collisions with beam energies from about 0.4 to 10 GeV/nucleon.
It has been broadly recognized that the high-density behavior of nuclear symmetry energy E sym (ρ) is the most uncertain aspect of the EOS of dense neutron-rich nucleonic matter besides possible phase transitions [12,17,18,19,20,21,22,23,24,25,26,27,28,29]. Indeed, the pressure P NS (ρ) from GW170817 has already been used to constrain some features of the E sym (ρ) within several approaches, see, e.g. refs. [31,32,33]. Since information about the composition of NSs is rarely extracted from analyses of astrophysical data so far, one often assumes that NSs are made of pure neutron matter (PNM) in many studies for various purposes. Such an extreme assumption implies that the E sym (ρ) is zero in NSs at β equilibrium. It is intrinsically inconsistent with the stated goal of extracting the underlying symmetry energy which determines uniquely the proton fraction x p (ρ) in NSs at βequilibrium (see, e.g., Eq. (6)).
Moreover, information about the density dependence of x p (ρ) has many important ramifications on properties of NSs especially the cooling mechanisms of protoneutron stars. So, is there any straightforward and self-consistent way to infer the high-density E sym (ρ) directly from the P NS (ρ) and P SNM (ρ) shown in where it is still observationally/experimentally largely unknown. Moreover, we extract simultaneously the proton fraction x p (ρ) in the same density range in NSs at β equilibrium.
Approach: It is well known that the EOS of cold asymmetric nucleonic matter (ANM) of isospin asymmetry δ and density ρ can be approximated as [34] in terms of the energy per nucleon E SNM (ρ) ≡ E(ρ, δ = 0) in SNM and the nuclear symmetry energy E sym (ρ). The corresponding pressure in ANM can then be written as is the pressure in SNM while the pressure in PNM P PNM (ρ) ≡ P (ρ, δ = 1) can be written as Therefore, if both the P SNM (ρ) and P PNM (ρ) are known from terrestrial nuclear experiments and/or astrophysical observations, the above equation can then be easily inverted to obtain the symmetry energy E sym (ρ) via where ρ i is a reference density below which the E sym (ρ) is known. To apply the Eq. (4), one needs the density dependence of both the pressure P SNM (ρ) in SNM and the pressure P PNM (ρ) in PNM. While the P SNM (ρ) has been extracted explicitly from analyzing heavy-ion reactions using mostly transport models where the EOS of cold SNM is a basic input [4,8], the pressure P PNM (ρ) in PNM is not the same as the NS pressure P NS (ρ) extracted directly from studying properties of NSs unless one explicitly assumes that NSs are completely made of just neutrons. In fact, it is well known that the proton fraction x p (ρ) at a given density in NSs at β equilibrium is determined uniquely by the symmetry energy E sym (ρ). Thus, the pressure extracted from analyzing observations of NSs without knowing explicitly their compositions should not be simply used as the P PNM (ρ) in extracting the underlying density dependence of nuclear symmetry energy E sym (ρ) in whatever approaches one may choose.
In the minimal model of charge neutral neutron stars consisting of neutrons, protons and electrons (npe matter) at β equilibrium, the pressure is given by where ρ e = 1 2 (1−δ)ρ and µ e = µ n −µ p = 4δE sym (ρ) are, respectively, the density and chemical potential of electrons. The value of the isospin asymmetry δ (or the corresponding proton fraction x p = (1−δ)/2) at β equilibrium is determined completely by the symmetry energy through the the chemical equilibrium and charge neutrality conditions. The resulting x p (ρ) can be written as [35] x We emphasize that the NS pressure P NS (ρ, δ) becomes barotropic, i.e., depending only on the density, once the density profile of the proton fraction x p (ρ) (and the corresponding δ(ρ)) at β equilibrium is determined. In the following discussions about how to find the right x p (ρ) (or δ(ρ)) through an iteration procedure, we keep using the notation P NS (ρ, δ) instead of the final P NS (ρ) for the model NS pressure.
To see how one may use progressively more accurately the available pressure P NS (ρ) extracted from observations without any NS composition information in evaluating the E sym (ρ), it is useful to first recast the NS pressure in Eq. (5) in terms of the P PNM (ρ), P SNM (ρ) and x p (ρ) as After expanding the second term (1 − 2x p ) 2 and using the definition of Eq. (3) for P PNM (ρ), the above equation can be rewritten as Thus, to calculate the E sym (ρ) using the Eq. (4) one may use as input Obviously, at the PNM limit, x p (ρ) = 0, one has P PNM (ρ) = P NS (ρ). The last two terms proportional to x p (ρ) and x 2 p (ρ) in Eq. (9) are progressively smaller quantities as x p (ρ) is always smaller than 0.5. Considering the coupled Eqs. (4), (6) and (9) together, it is seen that they constitute a closed iteration process for calculating directly the E sym (ρ) and the corresponding x p (ρ) simultaneously using the pressure P NS (ρ) from analyzing NS observations and the P SNM (ρ) extracted from heavy-ion reactions. Namely, first setting x p (ρ) = 0 (lowest accuracy) and inserting the resulting relation P PNM (ρ) = P NS (ρ) in Eq. (4), we calculate the symmetry energy E sym (ρ) at the zeroth order in x p . The resulting E sym (ρ) is then used to calculate a new density profile x p (ρ) and the corresponding δ(ρ) at β equilibrium using Eq. (6). Then, using the Eq. (9) one can calculate a new P PNM (ρ) to be used in the Eq. (4) to recalculate the symmetry energy E sym (ρ). The latter is then put back to the Eq. (6) to repeate the process until both the E sym (ρ) and x p (ρ) do not change anymore, namely, until the input and output E sym (ρ) on the two sides of Eq. (4) become and stay the same.
Inputs for the first application: As the first application of the procedure detailed above, we evaluate the upper and lower limits of the E sym (ρ) and x p (ρ) using the two pressures shown in Fig. 1. Since the Eq. (4) for calculating the E sym (ρ) involves the difference of P NS (ρ) and P SNM (ρ), we calculate the upper (lower) bounds of the E sym (ρ) and x p (ρ) by using the upper (lower) limit of P NS (ρ) from GW170817 but the lower (upper) limit of P SNM (ρ) from heavyion collisions. Naturally, this leads to the most conservative upper and lower limits of the E sym (ρ) and x p (ρ). We note that there are some ambiguities in assigning a statistical significance to the obtained boundaries because the P NS (ρ) from GW170817 is at 90% confidence level while the P SNM (ρ) from heavy-ion reactions is given in terms of its absolute (100%) upper and lower boundaries.
In addition, the resulting bands are expected to be wide because of the large uncertainties of the P NS (ρ) and P SNM (ρ) shown in Fig. 1. Nevertheless, we feel it is reasonable to consider the results as at 90% confidence level. Moreover, as we shall show next, the results are still very useful compared to not only the extreme PNM model for NSs but also the E sym (ρ) extracted from the recent Bayesian analyses using the currently available radius data of canonical NSs [38].
To infer the E sym (ρ) at high densities using Eq. (4), we need the value of E sym (ρ i ) at the reference density ρ i where the integration starts. In this work, since the P SNM (ρ) from analyzing nuclear collective flow data are only given in the range of 2ρ 0 to 4.5ρ 0 [4] although there are reports of P SNM (ρ) at lower densities but with different uncertainties from analyzing kaon productions in heavy-ion collision [44], we start the integration at 2ρ 0 as our interest here is mostly providing a guide for the E sym (ρ) at super-high densities. Moreover, we take the central value of E sym (2ρ 0 )=55.7 MeV from the ASY-EOS E sym (ρ) band (shown in red in Fig. 4) from analyzing relative flows and yield ratios of light mirror nuclei in heavy-ion collisions at SIS/GSI [7,45,46]. studied the E sym (ρ) and x p (ρ) along their upper and lower boundaries as functions of density and iteration number from the beginning to 300 steps. As expected, they both saturate quickly especially at relatively low densities and along the lower boundary. As examples, shown in Fig. 2 are the E sym (ρ) and x p (ρ) as functions of the step number at 2.5ρ 0 and 3.5ρ 0 along the two boundaries. It is seen that the x p (ρ) becomes stable after only about 5 steps for most cases while the E sym (ρ) takes a few more steps especially at high densities along the upper boundary. Generally speaking, because of the E sym (ρ)δ 2 term in the EOS of ANM, a higher E sym (ρ) requires a lower δ (or higher x p ) simply from the energetics while their exact relation is given by Eq. (6) as we discussed earlier. The x p reaches its upper limit of 0.5 quickly when the E sym (ρ) increases with density along the upper boundary. As a cubic root of the Eq. (6), the x p does not increases further once it has reached its upper limit even as the E sym (ρ) is still increasing, namely the NS becomes SNM with electrons when the E sym (ρ) becomes super-high. On the contrary, when the E sym (ρ) becomes super-low or even negative along the lower boundary, the NS becomes PNM.
These features are all expected. As indicated by the parabolic approximation of ANM EOS of Eq. (1), the nuclear symmetry energy is essentially the energy cost of converting all protons in SNM into neutrons to transform the system into PNM. When the E sym (ρ) is too high, this will not happen because it is not energetically favorible. Instead, the system becomes SNM with equal numbers of protons, neutrons and electrons even when we artificially initialized the system as a PNM. On the other hand, if the E sym (ρ) is low, then it is easier to convert more protons into neutrons. When the E sym (ρ) is super-low or negative, the NS quickly becomes PNM. These features from our numerical calculations clearly demonstrate that our iteration procedure worked successfully as it is supposed to.  The stabilized x p (ρ) values along both the upper and lower boundaries taken at the step 300 are shown as functions of density in Fig. 3. Below about 3.5ρ 0 , the x p (ρ) band is useful for better understanding properties of NSs. For example, there are at least about 10% protons at 2.5ρ 0 in NSs at β equilibrium. As we mentioned earlier, the band is expected to be wide given the uncertainties of the input pressures we used. Nevertheless, in our opinion, even getting the limited NS composition information from combining the barotropic pressures from both neutron stars and heavy-ion reactions within the novel approach proposed here represents a significant progress in the field. At higher densities, however, the proton fraction inferred is between zero and 0.5, namely, the data we used does not constrain at all the proton fraction of NS matter at densities above about 3.5ρ 0 . Thus, narrowing down the error bands of the pressures extracted in both astrophysical observations and terrestrial nuclear experiments in the ultra-dense region will be very useful.
Now we turn to the inferred nuclear symmetry energy in the super-dense region. For comparisons and completeness, it is first necessary to comment on the current status of constraining the E sym (ρ) at supra-saturation densities. As discussed extensively in the literature, see, e.g., refs. [17,18,12,19,20,21,22,23,24,25,26,27,28] and ref. [36] for the latest review, the E sym (ρ) has been relatively well constrained up to about 1.3ρ 0 using various terrestrial laboratory data, while in the region of 1.3ρ 0 − 2.0ρ 0 there are still some outstanding uncertainties and controversies [7,36,37]. On the other hand, the latest Bayesian analysis of the radii of canonical NSs indicates that the symmetry energy at twice the saturation density E sym (2ρ 0 ) =39.2 +12.1 −8.2 MeV at 68% confidence level [38]. This value is in very good agreement with findings of several other studies of neutron star radii, tidal deformability and maximum masses using different approaches within their statistical errors [30,31,32,39,40,41,42,43]. Interestingly, as shown in Fig.4, the value of E sym (2ρ 0 ) extracted from analyzing astrophysical data is also in reasonably good agreement with the ASY-EOS result (shown as the red bands in Fig. 4) [45,46]. However, above twice the saturation density, even the trend of the E sym (ρ) is not so clear. But it is encouraging to note that it has been found consistently in several very recent studies [32,38,43] that the very recently reported mass M = 2.17 +0. 11 −0.10 M ⊙ of PSR J0740+6620 [47] can tighten significantly the lower boundary of E sym (ρ) above about 2.5ρ 0 . However, this newly observed NS maximum mass remains to be confirmed and its statistical errors are still relatively large. Therefore, reliable observational/experimental information about the E sym (ρ) at densities above 2ρ 0 is presently very lacking. The density dependence of nuclear symmetry energy at suprasaturation densities from the present work in comparison with that from the Bayesian inference of ref. [38] using the radius data of canonical neutron stars. Fig.4 is a comparison of the E sym (ρ) from the present work with that from our recent Bayesian inference [38] from the radius data of canonical NSs. The ASY-EOS band (red) is also shown. As discussed earlier, our integration of Eq. (4) started from 2ρ 0 using the central value of the ASY-EOS band E sym (ρ)(2ρ 0 ) = 55.7 MeV. It is interesting to see that while the two constraining bands from this work and our recent Bayesian analyses are both quite broad at densities above 2ρ 0 , they overlap. If both bands are equally trustable, their overlapping area provides a reasonably tight constraint on the symmetry energy of super-dense neutron-rich nuclear matter above 2ρ 0 . With respect to the current status of the field mentioned above, this constraining band for the E sym (ρ) above 2ρ 0 is a welcomed step forward. Testing and further reducing it using future observational/experimental data and/or other approaches will certainly be very useful.

Shown in
Summary: A new approach is proposed for inferring the proton fraction and symmetry energy of super-dense neutron-rich nucleonic matter directly from the barotropic pressures P NS (ρ) in NSs and P SNM (ρ) for SNM from heavy-ion reactions. We have demonstrated that the approach is very efficient and viable.
An application of this approach using the NS pressure from GW170817 and the SNM pressure from relativistic heavy-ion collisions provides a useful constraining band for the symmetry energy above twice the saturation density of nuclear matter.