  The central speed of sound (SS) measures the stiffness of the Equation of
State (EOS) of superdense neutron star (NS) matter. Its variations with density
and radial coordinate in NSs in conventional analyses often suffer from
uncertainties of the specific nuclear EOSs used. Using the central SS and NS
mass/radius scaling obtained from solving perturbatively the scaled
Tolman-Oppenheimer-Volkoff (TOV) equations, we study the variations of SS,
trace anomaly and several closely related properties of NSs in an EOS-model
independent manner. We find that the SS increases with the reduced central
pressure $\widehat{P}_{\rm{c}}\equiv P_{\rm{c}}/\varepsilon_{\rm{c}}$ (scaled
by the central energy density $\varepsilon_{\rm{c}}$), and the conformal bound
for SS tends to break down for NSs with masses higher than about
1.9$M_{\odot}$. The ratio $P/\varepsilon$ is upper bounded as
$P/\varepsilon\lesssim0.374$ around the centers of stable NSs. We demonstrate
that it is an intrinsic property of strong-field gravity and is more relevant
than the perturbative QCD bound on it. While a sharp phase transition at high
densities characterized by a sudden vanishing of SS in cores of massive NSs are
basically excluded, the probability for a continuous crossover signaled by a
peaked radial profile of SS is found to be enhanced as $\widehat{P}_{\rm{c}}$
decreases, implying it likely happens near the centers of massive NSs.
Moreover, a new and more stringent causality boundary as
$R_{\max}/\rm{km}\gtrsim 4.73M_{\rm{NS}}^{\max}/M_{\odot}+1.14$ for NS M-R
curve is found to be excellently consistent with observational data on NS
masses and radii. Furthermore, new constraints on the ultimate energy density
and pressure allowed in NSs before collapsing into black holes are obtained and
compared with earlier predictions in the literature.
