The implementation of the ellipsoidal statistical Bhatnagar-Gross-Krook (ESBGK) method in the open-source particle code PICLas is extended for multi-species modeling of polyatomic molecules, including internal energies with multiple vibrational degrees of freedom. For this, the models of Mathiaud et al. [1], Pfeiffer [2] and Brull [3,4] are combined. In order to determine the transport coefficients of the gas mixture, Wilke's mixing rules and collision integrals are compared. The implementation is verified with simulation test cases of a supersonic Couette flow as well as a hypersonic flow around a 70° blunted cone. The solutions of the ESBGK method are compared to the Direct Simulation Monte Carlo (DSMC) method to assess the accuracy, where overall good agreement is achieved. In general, collision integrals should be preferred for the determination of the transport coefficients, since the results using Wilke's mixing rules show larger deviations in particular for large mass ratios. Considering the computational efficiency of the methods, a considerable reduction in computational time is possible with the ESBGK model.