  The study of ground-state properties of the Fermi-Hubbard model is a
long-lasting task in the research of strongly correlated systems. Owing to the
exponentially growing complexity of the system, a quantitative analysis usually
demands high computational cost and is restricted to small samples, especially
in two or higher dimensions. Here, we introduce a variational method in the
frame of fermionic Gaussian states, and obtain the ground states of one- and
two-dimensional attractive Hubbard models via imaginary-time evolution. We
calculate the total energy and benchmark the results in a wide range of
interaction strength and filling factor with those obtained via exact two-body
results, the density matrix renormalization group based on matrix product
states (MPS), and projector Quantum Monte Carlo (QMC) method. For both 1D and
2D cases, the Gaussian variational method presents accurate results for total
energy with a maximum systematic error ~4% in the intermediate interaction
region. The accuracy of these results has negligible dependence on the system
size. We further calculate the double occupancy and find excellent agreement
with MPS and QMC, as well as the experimental results of cold quantum gases in
optical lattices. The results suggest that the Gaussian pairing state is a good
approximation to the ground states of attractive Hubbard model, in particular
in the strong and weak coupling limits. Moreover, we generalize the method to
the attractive Hubbard model with a finite spin-polarization, which can be
mapped to the repulsive interaction case via particle-hole transformation, and
obtain accurate results of ground state energy and double occupancy. Our work
demonstrates the ability of the Gaussian variational method to extract ground
state properties of strongly correlated many-body systems with negligible
computational cost, especially of large size and in higher dimensions.
