  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.
