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 these configurations in order to enhance the yield of X-ray astronomy when combined with multimessenger gravitational-wave astrophysics and black hole shadows. Nevertheless, 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.