Progress in Constraining Nuclear Symmetry Energy Using Neutron Star Observables Since GW170817

New observational data of neutron stars since GW170817 have helped improve our knowledge about nuclear symmetry energy especially at high densities. We have learned particularly: (1) The slope parameter $L$ of nuclear symmetry energy at saturation density $\rho_0$ of nuclear matter from 24 new analyses is about $L\approx 57.7\pm 19$ MeV at 68\% confidence level consistent with its fiducial value, (2) The curvature $K_{\rm{sym}}$ from 16 new analyses is about $K_{\rm{sym}}\approx -107\pm 88$ MeV, (3) The magnitude of nuclear symmetry energy at $2\rho_0$, i.e. $E_{\rm{sym}}(2\rho_0)\approx 51\pm 13$ MeV at 68\% confidence level, has been extracted from 9 new analyses of neutron star observables consistent with results from earlier analyses of heavy-ion reactions and the latest predictions of the state-of-the-art nuclear many-body theories, (4) while the available data from canonical neutron stars do not provide tight constraints on nuclear symmetry energy at densities above about $2\rho_0$, the lower radius boundary $R_{2.01}=12.2$ km from NICER's very recent observation of PSR J0740+6620 of mass $2.08\pm 0.07$ $M_{\odot}$ and radius $R=12.2-16.3$ km at 68\% confidence level sets a tight lower limit for nuclear symmetry energy at densities above $2\rho_0$, (5) Bayesian inferences of nuclear symmetry energy using models encapsulating a first-order hadron-quark phase transition from observables of canonical neutron stars indicate that the phase transition shift appreciably both the $L$ and $K_{\rm{sym}}$ to higher values but with larger uncertaintie , (6) The high-density behavior of nuclear symmetry energy affects significantly the minimum frequency necessary to rotationally support GW190814's secondary component of mass (2.50-2.67) $M_{\odot}$ as the fastest and most massive pulsar discovered so far.