  The Neutron star Interior Composition Explorer (NICER) is currently observing
the x-ray pulse profiles emitted by hot spots on the surface of rotating
neutron stars allowing for an inference of their radii with unprecedented
precision. A critical ingredient in the pulse profile model is an analytical
formula for the oblate shape of the star. These formulas require a fitting over
a large ensemble of neutron star solutions, which cover a wide set of equations
of state, stellar compactnesses and rotational frequencies. However, this
procedure introduces a source of systematic error, as (i) the fits do not
describe perfectly the surface of all stars in the ensemble and (ii) neutron
stars are described by a single equation of state, whose influence on the
surface shape is averaged out during the fitting procedure. Here we perform a
first study of this systematic error, finding evidence that it is subdominant
relative to the statistical error in the radius inference by NICER. We also
find evidence that the formula currently used by NICER can be used in the
inference of the radii of rapidly rotating stars, outside of the formula's
domain of validity. Moreover, we employ an accurate enthalpy-based method to
locate the surface of numerical solutions of rapidly rotating neutron stars and
a new highly accurate formula to describe their surfaces. These results can be
used in applications that require an accurate description of oblate surfaces of
rapidly rotating neutron stars.
