We present an improved PAM (iPAM) method as the first fourth-order interface tracking method whose convergence rates are independent of C 1 (derivative) discontinuities of the interface. As an improved version of the polygonal area mapping (PAM) method [Q. Zhang and P. L.-F. Liu, J. Comput. Phys., 227 (2008), pp. 4063-4088], the accuracy of the iPAM method is achieved via (i) augmenting the abstract data structure of PAM to faithfully represent multiple components of material regions within a single cell, (ii) removing restrictive assumptions of PAM, (iii) adjusting the volume of represented cell material regions via polygon ear removal, and (iv) maintaining a relation (h L = r h h α ) between the Eulerian grid size h and the Lagrangian length scale that measures the distance between adjacent interface markers. Unlike volume-of-fluid methods and level-set methods, merging and separation of the tracked material is not automatically handled by the iPAM method. Instead, flexible subgrid resolutions are provided so that algorithmic behaviors of interface merging and separation can be determined from the specific physics of the problem under study. The iPAM method is a product of interdisciplinary research of different fields such as ordinary differential equations, computational geometry, and general topology. Its superior accuracy, efficiency, and versatility over previous interface tracking methods are demonstrated by a set of benchmark tests. In particular, for the vortex shear test on the 128 × 128 grid, the fourth-order iPAM method could be hundreds of times more efficient than state-of-the-art volume-of-fluid methods.