  We propose a framework to study the properties of the Lefschetz thimbles
decomposition for lattice fermion models approaching the thermodynamic limit.
The proposed set of algorithms includes the Schur complement solver and the
exact computation of the derivatives of the fermion determinant. It allows us
to solve the gradient flow (GF) equations taking into account the fermion
determinant exactly, with high performance. We can find both real and complex
saddle points and describe the structure of the Lefschetz thimbles
decomposition for large enough lattices to extrapolate the results to the
thermodynamic limit. The algorithms are described for a general lattice fermion
model, with emphasis on two types of lattice discretizations for relativistic
fermions (staggered and Wilson), as well as on interacting tight-binding models
for condensed matter systems. We apply these algorithms to the Hubbard model on
a hexagonal lattice, dealing with lattice volumes as large as 12x12 with 256
steps in Euclidean time, in order to capture the properties of the thimbles
decomposition as the thermodynamic, low-temperature, and continuum limits are
approached. The complexity of the thimbles decomposition appears to be
dependent on the form of the Hubbard-Stratonovich (HS) transformation. In
particular, the evidence is provided for the existence of an optimal regime for
the hexagonal lattice Hubbard model, with a reduced number of thimbles being
important in the overall sum. We have performed quantum Monte Carlo (QMC)
simulations using GF to deform the integration contour into the complex plane
and demonstrated the agreement with exact diagonalization on small volumes (8
sites in space) if optimal setup for HS transformation is used. The residual
sign problem is compared with the state-of-the-art BSS-QMC. We show that the
average sign can be kept substantially higher using the Lefschetz thimbles
approach.
