Finite element method is employed to solve the Maxwell’s equations in Cole-Cole dispersive media. First, the time-fractional derivative is transformed into an improper integral via a diffusive representation, and an equivalent local model is obtained by introducing auxiliary differential equations; in constructing the energy functional, correction terms are added to the classical energy to ensure a strictly monotonically decreasing energy functional. Second, the Raviart-Thomas-Nédélec mixed finite element method is adopted for spatial discretization, and the corresponding semi-discrete error is estimated. Finally, the improper integral is numerically approximated to yield a practical model, and a fully discrete scheme is derived by applying the Crank-Nicolson method for temporal discretization. Numerical experiments demonstrate that the proposed method exhibits energy stability for time-fractional derivatives of arbitrary order, and the numerical errors are in good agreement with the theoretical predictions.