The main aim of this contribution is to provide a mixed finite element for small strain elasto-viscoplastic material behavior based on the least-squares method. The L2-norm minimization of the residuals of the given first-order system of differential equations leads to a two-field functional with displacements and stresses as process variables. For the continuous approximation of the stresses, lowest-order Raviart–Thomas elements are used, whereas for the displacements, standard conforming elements are employed. It is shown that the non-linear least-squares functional provides an a posteriori error estimator, which establishes ellipticity of the proposed variational approach. Further on, details about the implementation of the least-squares mixed finite elements are given and some numerical examples are presented.