  X-ray astronomy provides information regarding the electromagnetic emission
of active galactic nuclei and X-ray binaries. These events provide details
regarding the astrophysical environment of black holes and stars, and help us
understand gamma-ray bursts. They produce estimates for the maximum mass of
neutron stars and eventually will contribute to the discovery of their equation
of state. Thus, it is crucial to study them in order to enhance the yield of
X-ray astronomy when combined with multimessenger astrophysics. An exact
solution of the field equations does not exist for rotating neutron stars.
There exist a variety of approximate solutions for compact objects that may
characterize relativistic stars. The most studied approximation is the
Hartle-Thorne metric that represents slowly-rotating compact objects, like
massive stars, white dwarfs and neutron stars. Recent investigations of photon
orbits and shadows of such metric revealed that it exhibits chaos close to
resonances. Here, we thoroughly investigate particle orbits around the
Hartle-Thorne spacetime up to second order in rotation. We perform an
exhaustive analysis of bound motion, by varying all parameters involved in the
system. We demonstrate that chaotic regions, known as Birkhoff islands, form
around resonances, where the ratio of the radial and polar frequency of
geodesics, known as the rotation number, is shared throughout the island. This
leads to the formation of plateaus in rotation curves during the most prominent
$2/3$ resonance, which confirms that generic geodesics are nonintegrable. We
measure their width and show how each parameter affects it. The
nonintegrability of Hartle-Thorne metric may affect quasiperiodic oscillations
of low-mass X-ray binaries, when chaos is taken into account, and might
potentially improve estimates of mass, angular momentum and multipole moments
of astrophysical compact objects.
