We explore several dominant eigenvalues of the spectra of Markov matrices governing the dynamics of models in the universality class of the two-dimensional Ising model. By means of a variational approximation, we determine autocorrelation times of progressively rapid relaxation modes. The approximation of one eigenstate, associated with the slowest mode, is employed in a variance-reducing Monte-Carlo method. The resulting correlation times, for which statistical errors exceed the systematic errors associated with the variational approximation, are used for a finite-size scaling analysis which corroborates universality of the dynamic critical exponent z for three distinct Ising models on the square lattice. Tentative, variational results for subdominant states strongly suggest that the amplitudes of the divergent time scales associated with different relaxation modes differ solely by metric factors, setting a single non-universal time scale for each model. A by-product of our analysis is a highly accurate confirmation of static universality.
