Since its development in 2003, the Partially Averaged Navier-Stokes equations (PANS) method has been widely used in numerous scientific and engineering applications, e.g., from flows past bluff bodies, cars, and trains to the increasingly important problem of ocean and climate prediction. PANS is a bridging method and, as such, can operate at any range of resolved scales – from fully resolved to fully modeled. It aims to only resolve the turbulent scales not amenable to modeling, representing the remaining through a closure. This feature is responsible for its efficiency and unleashes the concept of accuracy, one of the main drivers for PANS' success. The PANS framework was initially developed for constant density and incompressible closures, which neglects density fluctuations and turbulent kinetic energy generation by buoyancy mechanics. However, this set of features is crucial for variable density flows, such as material mixing and high-speed problems. This limitation has motivated the author to extend the PANS framework to variable-density flow. The present work summarizes the progress and challenges of PANS simulating variable density flows. It starts by describing the framework to derive PANS closures of variable density. Afterward, the performance of PANS is evaluated in the simulation of the Taylor Green Vortex, Rayleigh-Raylor mixing problem, and a shock-driven mixing flow. The results show that the PANS method can predict these transitional flows accurately, at a fraction of other Scale-Resolviong Simulations, such as DNS, LES, and ILES. The importance of understanding the predicted flow physics is demonstrated. The presentation concludes by discussing the PANS progress over the last 21 years and future challenges.