This work is a further development and the culmination along our research line of interface tracking in two dimensions [Zhang and Liu, J. Comput. Phys., 27 (2008), pp. 4063--4088; Zhang, SIAM J. Numer. Anal., 51 (2013), pp. 2822--2850; Zhang and Fogelson, SIAM J. Sci. Comput., 36 (2014), pp. A2369--A2400; Zhang and Fogelson, SIAM J. Numer. Anal., 54 (2016), pp. 530--560]. Previously, we have proposed an analytic framework for explicit interface tracking and a fourth-order interface tracking method. In this paper, we continue our approach to propose a new method, called the cubic MARS method, that features (1) a solid theoretical foundation consisting of well developed concepts from distinct branches of mathematics, (2) a representation of the interface with cubic splines, (3) a novel algorithm for Boolean operations on linear/curvilinear polygons, and (4) fourth-, sixth-, and eighth-order accuracy via a refined control of the relation between the Eulerian grid size $h$ and a Lagrangian length scale $h_L$ of interface markers. Volume conservation of the cubic MARS method also has convergence rates at 4, 6, and 8. A simple analysis and numerical results from several benchmark problems show that the cubic MARS method is superior to existing interface tracking methods in terms of accuracy and efficiency. In particular, for the single-vortex test and the deformation test, errors of an eighth-order cubic MARS method are close to machine precision on a 128-by-128 Eulerian grid! Currently the new method does not tackle topological changes of the interface, but its theoretical foundation provides a solid footing for future treatments of such phenomena.