Summary: The resistive magneto-hydrodynamic (MHD) model describes the behavior of a strongly ionized plasma in the presence of external electric and magnetic fields. For problems involving the plasma with an open boundary to very low pressure/vacuum regions, the continuum assumption is no longer valid in the entire domain of interest. For example, this is seen in a plasma accelerator where a dense plasma expands into a vacuum background. A common practice to deal with this issue is to assign small values of density and pressure to the vacuum regions and proceed to solve the continuum based MHD equations throughout the domain. We show that this approach fails to produce solutions consistent with the physics of the problem and can give rise to unacceptable artifacts such as spurious shocks. We develop a plasma-vacuum interface tracking approach to mitigate this problem. The plasma-vacuum interface is tracked in a physically consistent manner and the MHD equations are solved only in the regions that contain the plasma. The interface tracking is achieved using a face-flux formulation derived from the theoretical solution to a 1D free expansion problem. Coupled with a threshold based approach, the interface tracking is implemented for both explicit and implicit time stepping frameworks on generalized unstructured grids. Simulations of magnetized thermal plasma jets expanding into a vacuum background indicate plume profiles devoid of the unphysical shock obtained using the small background density approach. In the context of resistive MHD simulations, the interface tracking approach overcomes the numerical stiffness induced by specifying a large background resistivity in the vacuum regions, resulting in significant wall-clock time gains.