  The energy level statistics of the Hubbard model for $L \times L$ square
lattices (L=3,4,5,6) at low filling (four electrons) is studied numerically for
a wide range of the coupling strength. All known symmetries of the model
(space, spin and pseudospin symmetry) have been taken into account explicitly
from the beginning of the calculation by projecting into symmetry invariant
subspaces. The details of this group theoretical treatment are presented with
special attention to the nongeneric case of L=4, where a particular complicated
space group appears. For all the lattices studied, a significant amount of
levels within each symmetry invariant subspaces remains degenerated, but except
for L=4 the ground state is nondegenerate. We explain the remaining
degeneracies, which occur only for very specific interaction independent
states, and we disregard these states in the statistical spectral analysis. The
intricate structure of the Hubbard spectra necessitates a careful unfolding
procedure, which is thoroughly discussed. Finally, we present our results for
the level spacing distribution, the number variance $\Sigma^2$, and the
spectral rigidity $\Delta_3$, which essentially all are close to the
corresponding statistics for random matrices of the Gaussian ensemble
independent of the lattice size and the coupling strength. Even very small
coupling strengths approaching the integrable zero coupling limit lead to the
Gaussian ensemble statistics stressing the nonperturbative nature of the
Hubbard model.
