The direct simulation Monte Carlo (DSMC) method is applicable over a wide range of Knudsen numbers. However, the binary collision events make DSMC prohibitively expensive near the continuum regime. To address the computational cost issue in high-density regions, the particle-based Fokker-Planck (FP) method has been studied. In the FP method, particles evolve along independent stochastic paths, so the required cell size and time step do not need to resolve the collisional scale. While several monatomic and diatomic FP models have been proposed, the modeling of gas mixtures has received little attention so far. In this study, two new FP models are proposed to describe monatomic gas mixtures with an arbitrary number of constituents. One is the ellipsoidal statistical FP (ESFP) mixture model, which consistently evolves the relaxation of moments up to shear stress with the Boltzmann collision operator. The other is the ESFP+ mixture model, which includes a nonlinear drift coefficient to match the relaxation of moments up to heat flux with the Boltzmann collision operator. The numerical studies include the relaxation problem, Poiseuille flow, Couette flow, hypersonic flow around a vertical flat plate, and hypersonic flow around a cylinder. The results demonstrate that both the ESFP and ESFP+ models show good agreement with DSMC near equilibrium. However, the ESFP model fails to predict accurate shock structure at high Knudsen numbers, while the ESFP+ model better captures the shock structure.