Simulating multiphase flows remains challenging due to discontinuities across phases, nonlinear interactions, multiscale couplings, topological changes, and rich multiphysics effects. These flows critically depend on accurate representation of phase interfaces, motivating the development of interface-resolved simulation techniques. Phase field (diffuse interface) methods have emerged as a powerful approach for capturing interface dynamics. In particular, the conservative Allen–Cahn phase field model \cite{Chiu2011,Mirjalili2020} offers an attractive balance between accuracy and computational cost, exhibits strong parallel scalability, and is relatively straightforward to implement. In parallel, high-order spectral element methods achieve comparable accuracy with fewer degrees of freedom than low-order methods, exhibit excellent diffusive and dispersive properties, and remain efficiently parallelizable through matrix-free formulations. The smooth interface representation of phase field methods makes them well suited for such high-order discretizations. However, this combination remains largely unexplored for multiphase flows. \\ Neko \cite{Jansson2024} is a modern open-source high-order spectral element solver optimized for GPU architectures, and has demonstrated excellent performance for direct numerical simulations of single-phase turbulent flows. We extend Neko to multiphase flows by implementing the conservative Allen–Cahn equation and coupling it to the consistent Navier–Stokes equations \cite{Mirjalili2021}, aiming for a solver that can robustly and accurately simulate high-density ratio multiphase flows with surface tension effects. \\ A central aspect of this work is investigating the interaction between high-order spectral element discretization and the phase field method. Higher-order schemes offer improved accuracy but introduce challenges related to Gibbs-type oscillations near the interface that warrant systematic investigation. Using canonical tests, we quantify how polynomial degree affects interface resolution, boundedness, and the accuracy of calculating interfacial quantities. Finally, scaling tests on GPU clusters are carried out to verify that the multiphase extension preserves Neko's strong scaling characteristics, enabling efficient simulation of large-scale turbulent multiphase flows. \\