  Astrophysical charge-free black holes are known to satisfy no-hair relations
through which all multipole moments can be specified in terms of just their
mass and spin angular momentum. We here investigate the possible existence of
no-hair-like relations among multipole moments for neutron stars and quark
stars that are independent of their equation of state. We calculate the
multipole moments of these stars up to hexadecapole order by constructing
uniformly-rotating and unmagnetized stellar solutions to the Einstein
equations. For slowly-rotating stars, we construct stellar solutions to quartic
order in spin in a slow-rotation expansion, while for rapidly-rotating stars,
we solve the Einstein equations numerically with the LORENE and RNS codes. We
find that the multipole moments extracted from these numerical solutions are
consistent with each other. We confirm that the current-dipole is related to
the mass-quadrupole in an approximately equation of state independent fashion,
which does not break for rapidly rotating neutron stars or quark stars. We
further find that the current-octupole and the mass-hexadecapole moments are
related to the mass-quadrupole in an approximately equation of state
independent way to $\sim 10%$, worsening in the hexadecapole case. All of our
findings are in good agreement with previous work that considered stellar
solutions to leading-order in a weak-field expansion. The quartic in spin,
slowly-rotating solutions found here allow us to estimate the systematic errors
in the measurement of the neutron star's mass and radius with future X-ray
observations, such as NICER and LOFT. We find that the effect of these
quartic-in-spin terms on the quadrupole and hexadecapole moments and stellar
eccentricity may dominate the error budget for very rapidly-rotating neutron
stars. The new universal relations found here should help to reduce such
systematic errors.
