Extracting Nuclear Symmetry Energies at High Densities from Observations of Neutron Stars and Gravitational Waves

By numerically inverting the Tolman-Oppenheimer-Volkov (TOV) equation using an explicitly isospin-dependent parametric Equation of State (EOS) of dense neutron-rich nucleonic matter, a restricted EOS parameter space is established using observational constraints on the radius, maximum mass, tidal polarizability and causality condition of neutron stars (NSs). The constraining band obtained for the pressure as a function of energy (baryon) density is in good agreement with that extracted recently by the LIGO+Virgo Collaborations from their improved analyses of the NS tidal polarizability in GW170817. Rather robust upper and lower boundaries on nuclear symmetry energies are extracted from the observational constraints up to about twice the saturation density $\rho_0$ of nuclear matter. More quantitatively, the symmetry energy at $2\rho_0$ is constrained to $E_{\rm{sym}}(2\rho_0)=46.9\pm10.1$ MeV excluding many existing theoretical predictions scattered between $E_{\rm{sym}}(2\rho_0)=15$ and 100 MeV. Moreover, by studying variations of the causality surface where the speed of sound equals that of light at central densities of the most massive neutron stars within the restricted EOS parameter space, the absolutely maximum mass of neutron stars is found to be 2.40 M$_{\odot}$ approximately independent of the EOSs used. This limiting mass is consistent with findings of several recent analyses and numerical general relativity simulations about the maximum mass of the possible super-massive remanent produced in the immediate aftermath of GW170817.