  Strongly correlated electron systems are generally described by tight binding
lattice Hamiltonians with strong local (on site) interactions, the most popular
being the Hubbard model. Although the half filled Hubbard model can be
simulated by Monte Carlo(MC), physically more interesting cases beyond half
filling are plagued by the sign problem. One therefore should resort to other
methods. It was demonstrated recently that a systematic truncation of the set
of Dyson-Schwinger equations for correlators of the Hubbard, supplemented by a
\textquotedblleft covariant" calculation of correlators leads to a convergent
series of approximants. The covariance preserves all the Ward identities among
correlators describing various condensed matter probes. While first order
(classical), second (Hartree-Fock or gaussian) and third (Cubic) covariant
approximation were worked out, the fourth (quartic) seems too complicated to be
effectively calculable in fermionic systems. It turns out that the complexity
of the quartic calculation\ in local interaction models,is manageable
computationally. The quartic (Bethe - Salpeter type) approximation is
especially important in 1D and 2D models in which the symmetry broken state
does not exists (the Mermin - Wagner theorem), although strong fluctuations
dominate the physics at strong coupling. Unlike the lower order approximations,
it respects the Mermin - Wagner theorem. The scheme is tested and exemplified
on the single band 1D and 2D Hubbard model.
