Generalized Hartree-Fock Theory for Interacting Fermions in Lattices: Numerical Methods

We present numerical methods to solve the Generalized Hartree-Fock theory for fermionic systems in lattices, both in thermal equilibrium and out of equilibrium. Specifically, we show how to determine the covariance matrix corresponding to the Fermionic Gaussian state that optimally approximates the quantum state of the fermions. The methods apply to relatively large systems, since their complexity only scales quadratically with the number of lattice sites. Moreover, they are specially suited to describe inhomogenous systems, as those typically found in recent experiments with atoms in optical lattices, at least in the weak interaction regime. As a benchmark, we have applied them to the two-dimensional Hubbard model on a 10x10 lattice with and without an external confinement.


I. INTRODUCTION
Recent experiments with cold atomic systems have attracted the interest of several scientific communities. Many physical phenomena, like Bose-Einstein condensation, superfluidity, or the formation of Cooper pairs, have been observed in seminal experiments and carefully analyzed and characterized [1,2]. In the last years, the possibility of loading the atoms in optical lattices offers a new avenue to the observation of many other intriguing phenomena [3].
There exists several theoretical tools in order to describe most of those phenomena. For the case of Bosons, the Gross-Pitaevskii equation provides us with the basic tool to analyze many static and dynamic properties for cold atomic gases [4]. Such an equation is valid whenever the occupation number of a particular mode is large, and it is equivalent to a mean-field theory. Typically, it can be solved numerically, offering very accurate predictions and explanations of all the observed phenomena in the weak interacting regime. Furthermore, it can also be used to describe lattice systems, at least as long as strong correlations are absent. In the case of Fermions, the basic tools at hand are Hartree-Fock and BCS theory [5]. The first one assumes that the state of the particles can be approximated in terms of a Slater determinant, and thus it also correspond to a mean-field theory. The latter assumes a specific form of the wavefunction, whereby pairing between Fermions is allowed. For both approaches, a time-dependent version has been constructed [6][7][8][9][10].
Hartree-Fock and BCS theories can be unified in terms of a more general framework, the generalized Hartree-Fock Theory (gHFT) [11]. The main idea is that both, the Slater-Determinant-and BCS-states are members of a larger family of states, the so-called Fermionic Gaussian states (FGS). Therefore, one may try to find the state within this family which best approaches the real state of the system under study. FGS are those states whose density operator can be expressed as an exponential of a quadratic function of canonical creation and annihilation operators. This property immediately implies that they fulfill Wick's theorem: that is, all correlation functions can be expressed in terms of the second moments of all creation and annihilation operators. Such moments are gathered in a matrix, the so-called covariance matrix (CM) which completely characterizes the state [12]. Thus, FGS can be very efficiently described, since the number of parameters only scales quadratically with the number of available fermionic modes and the expectation values of observables can be easily computed. For the case of lattices, the number of modes is automatically finite, and thus it should be an ideal playground for such a gHFT. In fact, Bach, Lieb and Solovej [11] were the first who applied such a theory to a Hubbard model in 2D, and where able to solve the corresponding equations exactly for the homogeneous case, both for the ground state and in thermal equilibrium.
gHFT can be easily applied to any lattice system, homogeneous or not. The resulting equations are, however, difficult to solve in general. In this paper we derive such equations for the CM, and introduce different methods to solve them. First of all, we derive the equation for the real time evolution of the CM. Second, we consider the state within the FGS which minimizes the total energy, and thus approximates the ground state. For that, we derive the evolution equation (in imaginary time) for the CM in the presence of a lattice Hamiltonian. We show that under such an equation, the CM converges to the one which minimizes 1 the energy within the family of FGS. Thus, the evolution equation provides us with a practical way of determining the approximation to the ground state. Finally, we derive the equation for the CM that minimizes the free energy, and apply a simple fixed-point iteration method to solve it.
Then we apply the methods to the spin-1/2 Hubbard model in a 10 × 10 lattice. First, we benchmark them with exact solutions of the gHFT (for the homogeneous case) [11], finding an extraordinary agreement both for attractive and repulsive interactions. Second, in order to test the validity of gHFT itself, we compare the results for positive interactions, where strongly correlated effects are expected to be most pronounced, with the Monte Carlo results of Refs. [13,14]. The results are obviously less precise than those obtained with such methods, other Monte Carlo approaches [15], or the recently developed algorithms based on PEPS and other tensor networks states [16][17][18][19][20][21]. Although we obtain a good qualitative agreement with [13,14] in many regimes, the FGS are, as expected, not able to capture all possible fermionic phases in the strong-correlation regimes, as they constitute a subclass of all possible fermionic states. For example, the finite temperature Mott phase cannot be approximated well as a Gaussian state, and is thus absent in the phase diagram.
At this point we would like to warn the reader that in this work we will mainly concentrate on the possibilities offered by the numerical approaches, and not focus on the physics of the Hubbard model. To be precise, our goal is a) to formulate the standard gHFT for fermionic systems with a two-body interaction in the lattice using the language of covariance matrices, b) to derive numerical methods to determine the ground and thermal states as well as the time evolution of these systems within gHFT and c) to benchmark our techniques by applying them to the two-dimensional Hubbard model. This paper is organized as follows: In Sec. II we give a precise formulation of the problem we aim to solve, and introduce the concepts necessary for the understanding of FGS and gHFT. This Section includes further a discussion why FGS are supposed to capture well the physical properties of some of the most relevant fermionic models. Next, we derive the methods to approximate the real time evolution, the ground and thermal states of fermionic systems within gHFT in Secs. III A-III C. Then, we benchmark our methods by applying them to the two-dimensional translationally invariant Hubbard model, both in the attractive and the repulsive regime, in Sec. IV. In addition, motivated by the current possibility of implementing the Hubbard model in more than one dimension with the help of optical lattices, we also investigate the 2d Hubbard model with an external confinement in Sec. V, and close with a study of dynamical processes, both for the translationally invariant and the trapped system, in Sec. VI.

II. STATEMENT OF THE PROBLEM
In this Section we give a statement of the problem, summarize the properties of fermionic Gaussian states that are necessary for the understanding of our work (for more details see e.g. [12]) and argue why these states capture the properties of many relevant physical models.
In our work we aim at understanding the physical properties of fermionic systems in a lattice, like for instance where the a k denote fermionic mode operators obeying canonical anti-commutation relations (CAR), {a k , a † l } = δ kl and t kl , u klmn ∈ C. These models describe fermionic atoms in optical lattices [22], and one of the most prominent ones, the Hubbard model [23], is expected to describe high-Tc superconductivity [24]. We approach the problem in the Majorana picture and define c k = a † k +a k , where M is the number of modes and k = 1, . . . , M . These operators satisfy the CAR {c k , c l } = 2δ kl . The Hamiltonian H in this picture is given by where T kl , U klmn ∈ R. The CAR allow us to antisymmetrize T and U such that T T = −T while U is antisymmetric under the exchange of any of two adjacent indices.
Our goal is to derive ground and thermal states as well as the time-evolution of systems characterized by the Hamiltonian given in Eq. (1) within the family of FGS. These are states that can be represented as an exponential of a quadratic form in the Majorana operators, and thus they are fully characterized by their second moments collected in the real and anti-symmetric covariance matrix (CM) Γ kl = i 2 [c k , c l ] from which all higher correlations can be obtained via Wick's theorem (see e.g. [12]): with real and antisymmetric Hamiltonian matrix h. All Gaussian states remain Gaussian under the time evolution governed by a quadratic Hamiltonian, and the CM transforms according to Let us now explain why FGS provide a very powerful technique for the description of fermionic manybody systems. First, note that all pure states within this family can be brought into a standard form |Ψ = where |u k | 2 + |v k | 2 = 1 ∀k via a change of basis [25]. Thus, the pure Gaussian states include the BCS-states introduced in the theory of superconductivity, as well as all Hartree Fock states. Second, FGS also include all mixed states for which Wick's theorem (3) applies, e.g., thermal states of quadratic Hamiltonians. Thus, all together, we see that FGS include a variety of states which form the basis of several manybody phenomena.

III. INTERACTING FERMIONS IN GENERALIZED HARTREE FOCK THEORY
In the following Section we derive the theoretical framework that will allow us later on to simulate fermionic many-body problems, like the Hubbard model, for big system sizes. This is achieved by approximating the state of the fermions by a FGS, something which is implemented by means of Wick's theorem.

A. Real Time Evolution
We start the section deriving an evolution equation that allows us to describe the time-evolution of a system at zero and finite temperature evolving under the interaction Hamiltonian (1) within generalized Hartree-Fock theory. The key idea is to use the covariance matrix for the description of the dynamical evolution. In the Heisenberg picture, where the time evolution of the Majorana operators is given by c k (t) = e iH(T,U )t c k e −iH(T,U )t , the covariance matrix evolves according to As H(T, U ) involves terms quartic in the Majorana operators this evolution clearly takes us out of the Gaussian setting. We truncate this transformation imposing that Wick's theorem holds, i.e. we write With the help of the commutation relation we can calculate the contributions for the quadratic term H q = kl T kl c k c l and the interaction term This implies the following time evolution of the covariance matrix: where we have defined tr B [U Γ] ij ≡ kl U ijkl Γ lk . This equation can be formally integrated: where T denotes the time ordering operator. Note that due to the anti-symmetry ofh(Γ(t )) Eq. (9) Further, this approximation scheme does not only ensure that we always remain within the set of Gaussian states, but a short calculation shows that it also conserves energy and, for a number conserving Hamiltonian, the particle number. We start with the energy, which is given by To prove the conservation of the mean particle number N , as long as [H,N ] = 0, we write the particle number operatorN = k a † k a k = M 2 + i Hence, we have shown that applying Wick's theorem to the evolution equation of a system governed by a twobody interaction Hamiltonian leads to a consistent dynamical equation for the CM. The truncation of the evolution is equivalent to an evolution under an effective state-dependent quadratic Hamiltonian within the set of Gaussian states, as depicted in Fig. 1 a). We remark that this approach works for pure as well as mixed FGS. Note also that the remaining equation resembles in some way the Gross-Pitaevskii equation for Bosons since that one can also be interpreted as giving the evolution in terms of a Hamiltonian that depends on the state; however, the main difference is that in the case of Bosons such an equation is obtained for the first moments of the field operators. Note that in the case of Bosons one could also perform the same approximation for the first and second moments, obtaining a set of coupled equations which would include the condensate part as well as thermal particles.

B. Ground States
In principle, the ground state in gHFT can be found via a direct minimization min ρ Gaussian For the translationally invariant Hubbard model this problem could be reduced to an optimization over two parameters only [11]. Generically, this constrained quadratic minimization problem is a daunting task. We attack this problem from a different perspective. A wellknown approach for the determination of the ground state |φ 0 of a Hamiltonian H is imaginary time evolution. Due to the exponential growth of the state space with the number of modes this approach can be applied to small systems only. The idea is to apply the Gaussian approximation (in the form of Wick's theorem) to derive an evolution equation for the CM. In this way we get a simple quadratic scaling in the computational time with respect to the number of modes. Naively, one would think of obtaining the evolution equation in imaginary time by replacing t → it in the real-time evolution equation of the CM, Eq. (7). However, as it turns out, this is not the right road to take. Though, as we will show below, the following approach will lead us to the desired ground state: We start with an arbitrary pure Gaussian state ρ(0), and evolve it according to where the normalization of the density operator for any time, i.e. the denominator tr[e −2Ht ρ(0)], is crucial for obtaining the correct equations. Since we apply the nonquadratic operator H to the Gaussian state ρ, this procedure clearly takes us out of the setting of Gaussian states. One way around would consists of discretizing the evolution and finding for small time steps ∆t the best Gaussian approximation at each step. However, due to the truncation it is not clear that this procedure will converge. And even if we find a steady state it is not clear that this will be the ground state of the system. Another possible approach is to use the quadratic but state-dependent effective Hamiltonian H Q (Γ) that is derived from H(T, U ) for the imaginary time evolution. As H Q (Γ) is a quadratic operator we will stay in the space of Gaussian states. But as the Hamiltonian is state-dependent the outcome of this procedure is not clear. However, as we will show below, both approaches are equivalent and lead indeed to the desired solution. To be precise, we show that the following is equivalent: 1. Direct minimization of the energy Eq. (12) in generalized Hartree-Fock theory.
2. Imaginary time evolution of the state ρ with the full Hamiltonian H(T, U ) for small time steps ∆t followed by an approximation of ρ(t + ∆t) by a Gaussian state.
3. Imaginary time evolution of ρ with the quadratic but state-dependent Hamiltonian H Q (Γ).

Minimization of the energy
In order to obtain the generalized Hartree-Fock ground state we have to solve the minimization problem min ρ Gaussian It has been proven in Ref. [11] that the HF ground state is always pure, i.e. Γ 2 = −1. Using Lagrange multipliers, we arrive at the following necessary conditions for a local minimum (see Appendix A) 3 : These two equations are non-linear matrix equations and thus hard to solve, both analytically and numerically, for large systems. But we show next that these equations appear as the steady-state conditions of imaginary time evolution.

Imaginary time evolution
From Eq. (13) we see that the evolution of the density operator ρ under any Hamiltonian H in imaginary time is given byρ so that the covariance matrix evolves according to We show in Appendix B that both approaches for imaginary time evolution, number 2 and 3, lead to the same evolution equation of the covariance matrix: First, note that Eq. (19) ensures that a pure state remains pure under the evolution, since as all pure states fulfill Γ 2 = −1. Next, we show that the energy always decreases under the evolution: Here The evolution of Γ with the orthogonal O(t) has turned out to be especially fast and stable in numerical implementations compared to the direct integration of Eq. (19) via e.g. a Runge-Kutta solver, and we have used it for our numerical applications.

C. Thermal states
As we have seen in the last Section, the generalized Hartree-Fock ground state can be obtained via an imaginary time evolution. To learn about the finite temperature properties in gHFT we have to consider the approximation to the Gibbs state ρ ∼ e −βH , where β = 1 k B T is the inverse temperature. We recall that the Gibbs state minimizes the free energy, As in the case of the minimization of the energy, we reformulate the problem as an optimization problem using Lagrange multipliers. As we show in Appendix C, this approach leads to the following necessary conditions for a minimum of the free energy: where Now, assume for the moment that det(1+Γ 2 ) = 0. Then, according to Eq. (24), h F = 0 is the only possible solution. Hence, the CM has to be of the form Now, a CM of this form always fulfills Eqs. (23) and (24), and includes, in the limit β → ∞, the zero-temperature case, so that Eq. (26) is indeed a solution of Eqs. (23) and (24). Note that from the perspective of numerical implementation, we can solve Eq. (26) e.g. via a fixed-point iteration (see Fig. 1c).  27)) for an interaction strength u ∈ [−10, 10] at half-filling, i.e. µ = 0. The blue line shows the exact gHF ground state energy, Eexact, obtained from Ref. [11], while the triangles correspond to the energy obtained via imaginary time evolution, Eimag.

IV. BENCHMARK: THE TWO-DIMENSIONAL TRANSLATIONALLY INVARIANT HUBBARD MODEL
In the following Section we apply our method to the two-dimensional translationally invariant Hubbard model with periodic boundary conditions: where x and y are points on a two-dimensional lattice, i, j denote nearest-neighbors and σ =↑, ↓ denotes the spin degree of freedom. We consider only the case where we have the same number of spin-up and spin-down particles, so that we can use the same chemical potential for the two species. For u < 0 (u > 0) the second term in H is an attractive (repulsive) on-site interaction between particles of opposite spin. In the following, we set t = 1 as the energy scale of the system. Note that the case of half-filling is characterized by µ = 0. The goal of this Section is two-fold: In the first part, Sec. IV A, we benchmark our numerical method by considering physical quantities for which an exact solution within gHFT is known. In the case of the translationally invariant Hubbard model it could be proven in Ref. [11] that the (free) energy for the (thermal) ground state can be found via a two-parameter optimization, and we will compare our numerical results with the numbers obtained from the optimization. After the demonstration of the power of our approach, we continue in the second part, Sec. IV B, with a comparison of gHFT itself to more powerful and and sophisticated methods, like Quantum Monte Carlo (QMC). In the following, all calculations are performed for a 10 × 10 lattice.

A. Comparison with the exact gHF solution
For the translationally invariant Hubbard model Eq. (27) the exact results for the energy of ground and thermal state within gHFT where presented in [11]. In Fig. 2 we compare the ground state energy for halffilling obtained via imaginary time evolution according to Eq. (21) (triangles) with the exact gHF solution (blue line) and find excellent agreement. The relative error in the energy is given by In Fig. 3 we show results away from half-filling. In this case, a closed solution for the energy could only be provided for the attractive Hubbard model in Ref. [11]. We have compared the exact gHF solution (dotted lines) with the results obtained via our approach (triangles and diamonds) for u = −2, −6, −10, proving that our numerical method also works well for a doped system.
Next, we consider the case of finite temperature, where the gHF Gibbs state is given by the implicit equation (26). Starting from the ground state we change the temperature in steps ∆β = 0.01 and compute the new thermal state via a fixed-point iteration. In Fig. 4 we present a comparison of our solutions (triangles and diamonds) for half-fiiling with the exact results (lines) for various values of u in a temperature range of β ∈ [0.4, 1.6], and find excellent agreement.
One of the most interesting questions concerning the Hubbard model is the appearance of superconductivity. As we have already pointed out in Sec. II, FGS describe among others superfluid and normal states, since every pure state of this family can be brought into the standard for various values of the interaction u at half-filling. We observe the occurrence of a phase transition from a paired phase to an unpaired normal state. The circles indicate the value of the transition temperature derived for gHFT in Ref. [11].
|0 via a change of basis. In this basis, the gap parameter ∆ = k u k v k is an order parameter for the superfluid phase. To consider superfluidity independent of the basis, it is more appropriate to consider the normalized pairing measure derived in [26]: Then, a positive value of P indicates that we are in a paired phase, while a normal state is characterized by P = 0. It has been shown in Ref. [11] that within gHFT the attractive Hubbard model exhibits a phase transition from a paired phase to an unpaired normal phase at a critical temperature β c for any finite lattice size, and the cor-responding transition temperature has been derived. In Fig. 5 we have depicted the pairing as a function of the inverse temperature for various values of the interaction u, and find a breakdown at a finite value β c that depends on the interaction strength. The values of β c that we obtain via our numerical method agree well with the results presented in Ref. [11], and that are depicted as circles in Fig. 5. Further, we give a list of the critical exponent γ defined via P = a(T c − T ) γ for various values of u in table I.  Finally, we give in Fig. 6 a phase diagram of the pairing for the ground state as a function of the chemical potential and the interaction strength. (Note that the value of the pairing is not unique in the case of halffilling due to the Gauge symmetries of the Hamiltonian (see Ref. [11]). We have restricted to represent the phase diagram for negative values of u, since we find that P = 0 for u positive.
In summary, we have demonstrated that our numerical methods are capable of obtaining the correct gHF ground and Gibbs state for the translationally invariant Hubbard model, both in the attractive as well as in the repulsive regime. Having demonstrated the validity of our approach, the next obvious question is now how gHFT itself compares to more powerful and sophisticated methods, like Quantum Monte Carlo (QMC).

B. Comparison with Quantum Monte Carlo
In this Section we compare the results obtained via our algorithm to numerical data from QMC simulations. To this end, we will concentrate on the repulsive Hubbard model at half-filling. The reason for this choice is twofold: First, as it has been shown in [27,28], the physics of the attractive Hubbard model is well-described by BCStheory, which is included in gHFT. This motivates to consider primarily the repulsive Hubbard model, and compare our results with results obtained from QMC. While QMC suffers from the notorious sign problem when applied to fermionic systems, it could be shown that in the case of a half-filled lattice the sign problem does not occur [29], and thus QMC becomes exact. As a consequence, the two-dimensional repulsive Hubbard model at half-filling has undergone an intense investigation in recent years. A partial list of these results is given in [30][31][32][33][34][35][36][37][38][39][40] and references therein. In this Section, we compare our results obtained in gHFT with the recent QMC  (28)) of the ground state as a function of the interaction strength and the chemical potential in the translationally invariant case. As a guide for the eye, we have added lines for the particle numbers 50, 100 (half filling) and 150.

Ground state properties for half-filling
We start with a comparison of our approach with the results presented in Ref. [13], where a half-filled system is considered. Since in the latter work the numerical results are presented for various lattice sizes, while our results are always given for a 10 × 10 lattice, we have to restrict to a qualitative comparison in certain cases. As the first physical quantity, we consider the momentum distribution n(k) that has been depicted for u = 2, 3, 4, 5, 6, 8 and very small, but non-zero temperature in Ref. [13]. In Fig. 8 we show our results for the momentum distribution n(k) for an interaction u = 1, . . . , 10 at half-filling in the entire Brioullin zone. Our results show an excellent qualitative agreement with Fig. 1a of Ref. [13]. As in Ref [13] we find that with increasing u, the momentum distribution fulfills n(0, 0) → 1 and n(π, 0) = 1/2 for all values of the interaction strength, and n(k) grows with increasing u for k ∈ [(π, π), (π, 0)]. Further, as in [13] we observe a broadening of the Fermi surface with growing interaction strength, i.e. the momentum distribution gets smeared out.
In the following, we consider two-particle properties of the system. To this end, we study as in Ref. [13] the following magnetic properties of the system: • The equal-time spin-spin correlation function C( y) = (n ( x+ y)↑ − n ( x+ y)↓ )(n x↑ − n x↓ ) .
In Fig. 9 we have depicted the local magnetic moment of the ground state, C(0, 0) = (n x↑ − n x↓ ) 2 , for an interaction strength u = 0, . . . , 10 at half-filling. The purple squares depict the results obtained via our algorithm, while the blue diamonds represent the QMC results of [13]. Singly occupied sites fulfill C(0, 0) = 1, while for empty or doubly occupied sites C(0, 0) = 0. In the non-interacting regime, all four states (full, empty, single occupied in the two spin-states) are equally probable, so that C(0, 0) = 1 2 . We find that the on-site repulsion suppresses doubly occupied, and for a half-filled lattice also the empty configuration. With increasing u the singly occupied configuration is enhanced, leading to a well-formed moment at each site.
Next, we study the equal-time spin-spin correlation function C( y) = (n ( x+ y)↑ − n ( x+ y)↓ )(n x↑ − n x↓ ) . This  quantity measures the extend to which two spins on site x and x + y are aligned. Note that the definition of C( y) is rotationally invariant. Our results are presented in Fig.  10. The correlations extend, already for small values of u, over the entire lattice, as it is predicted for a half-filled lattice. As one would expect intuitively, the antiferromagnetic correlations are enhanced with growing interaction strength. Our results show the same qualitative behavior as in [13], where C( y) has been calculated for a 24 × 24 square lattices. Finally, we consider the magnetic correlation function S(k) = l e i k· l C( l) in Fig. 11. A sharp peak occurs at momentum (π, π), emphasizing the antiferromagnetic correlations. (Compare with Fig. 8 of Ref. [13], where S(k) is depicted for u = 2). Further, we see that the magnetic correlations depend strongly on the interaction strength.

Ground state properties of a doped system
Next, we consider the case of a doped system, i.e. µ = 0. This instance has been studied e.g. in Ref. [14]. There, the local density n x was calculated as a function of the chemical potential. The results obtained via our approach (see Fig. 12) show a qualitative agreement with Fig. 2a of Ref. [14], where a density of n(µ) = 1 is found for µ ≈ −1, . . . , 1. We find that for µ = 0 we have half-filling, as expected, while an increase in the chemical potential results finally in doubly occupied sites.
We include here a remark on the limitations of gHFT. It is believed that the doped Hubbard model with repulsive interaction is a highly correlated state, exhibiting pairing for some value of the doping. However, for u ∈ [0, 10] and µ ∈ [−5, 5] we only find an unpaired phase. This is in agreement with Thm. 2.11 of Ref. [11], where it has been proven that for any positive definite potential U the gHF ground state is unpaired, independent of the doping.

Thermal state properties for half-filling
In this Subsection we consider some finite temperature properties of the repulsive Hubbard model, and concentrate again on the case of half-filling. We have seen in the last Section that gHFT provides us with a powerful tool to study the magnetic properties of fermions with repulsive interactions in a lattice. Especially, we have seen that the ground state of the system shows antiferromagnetic correlations, already for small values of the interaction strength. These are expected to decrease with increasing temperature, since the thermal fluctuations destroy the antiferromagnetic order. We find indeed this behavior, which is depicted in Fig. 13. There, we have plotted the order parameter A(d) = n x↑ n x+d↓ (29) as a function of the of the inverse temperature for u = 10.
We have also derived the critical exponent α given by A(T ) ∼ (T − T c ) α for lattice sizes 8 × 8 and 10 × 10, leading to α ≈ 1.3 − 1.4. The antiferromagnetic phase is further characterized by the fact that the local fluctuations in the particle number vanish, i.e. we are in a Mott phase. This phase is predicted to last for even higher temperatures than the antiferromagnetic phase ("Mott-gap"), since additional energy is needed to delocalize the fermions [41]. We investigate this behavior in Fig. 14. There, the solid purple line shows the value of A 1,2 = A(1) (c.f. Eq. (29)), the antiferromagnetic correlations of nearest-neighbors, while the blue dotted line is the order parameter for the Mott phase, O M = (n x ) 2 − n x 2 . We see that independent of the interaction strength the Mott phase as well as the antiferromagnetic phase break down at the same temperature, i.e. we do not find a Mott-gap. This observation can be explained by the fact that the FGS cannot represent the Mott phase, as we have already stated in the introduction.
The fact that there exists no Gaussian state corresponding to the Mott phase except the state of the antiferromagentic phase also leads to an incorrect transition temperature from the antiferromagnetic to the normal phase. For the exact solution of the Hubbard model the antiferomagnetic correlations are predicted to be destroyed at a temperature T c ∼ t 2 /u which is the energy scale of the superexchange energy. However, in our approach the antiferromagnetic correlations appear as soon as we find one particle per site, leading to a transition temperature T c ∼ u (cf. Fig 14).

V. THE HUBBARD MODEL IN AN EXTERNAL TRAP
While the translationally invariant Hubbard model allows for an elegant analytic solution within generalized Hartree Fock theory, experimental setups break this symmetry. Since the particles have to be spatially confined, an external trapping potential is needed. In many setups, this trapping is realized via a harmonic confinement. The Hamiltonian of the system is then altered by a position- dependent chemical potential, i.e.
where for In the following we investigate how the external confinement alters the physical properties of the attractive Hubbard model (with open boundary conditions) at zero and finite temperature. To gain first insights into the behavior of the system, we have depicted in Fig. 15 the density at the center of the trap as a function of the interaction strength and the chemical potential for three different trapping potentials (Fig. 15 b-d), as well as for the translationally invariant case (Fig. 15 a). We see that the filling in the center of the trap is only altered slightly be the external potential. In the following, we call a system half-filled if there is one particle at the center of the trap. We find that the condition for half-filling does not depend much on the trapping potential, and we can use the results from Fig. 15 to adjust the chemical potential in order to obtain half-filling in the following.

A. Attractive Hubbard model in a trap
We start with a closer investigation of the effects of an external trapping potential on the physics of the attractive Hubbard model. In Fig. 16 we present the density distribution for u = −5 at half-filling for four different trapping potentials V t = 0.01, 0.05, 0.1, 0.25. We see that for V t = 0.01 the density distribution is close to that of a free system with open boundary conditions, so that boundary effects are expected in this case. Hence, we will consider only V t = 0.05, 0.1, 0.25 from now on, and compare the results to the translationally invariant model with periodic boundary conditions. Now we investigate how the phase diagram of the pairing changes under the influence of an external confinement. In Fig.17 we have depicted how the pairing per particle as a function of the interaction strength u and the chemical potential µ changes under the influence of an external trap. We have added lines for 50, 100 and 150 particles as a guide to the eye.
To gain further insight into the influence of the trap, we have depicted in Fig. 18 the pairing per particle for a half-filled trap, i.e. a setting where we find one particle at the center of the trap, and compare with the translationally invariant case. Note that for a translationally invariant system an additional Gauge symmetry emerges at half-filling [11], and the underlying Gauge transformation does not leave the pairing invariant. However, this Gauge symmetry is broken for an inhomogeneous system, and we can only perform a qualitative comparison of the translationally invariant and trapped system. In this respect, we observe that the pairing decreases with decreasing interaction strength, both in the translationally invariant and in the trapped system, as expected. Further, we see that a very strong confinement can lead to a slight decrease in the pairing per particle. This effect can be understood in the following way: Since we fix the number of particles at the center of the trap to one, and the external trap acts like a position-dependent chemical potential, a strong trap will significantly reduce the number of particles at the boundary of the lattice, and the paired phase is suppressed. We continue our investigations on the influence of the external trapping potential by analyzing the temperature-dependence of the pairing in an external trap. In Fig. 19 we consider interactions u = −10 and −5 and a temperature regime of β ∈ [0, 2]. First, we observe that the external trapping alters only slightly the temperature-dependence of the pairing. Further, the value of β c , where a phase transition from a paired to a normal state occurs is not altered by the external confinement, and is identical to the translationally invariant case. At this point we would also like to mention that the Gauge freedom of the ground state is lifted with increasing temperature, and a unique gHF ground state is expected if the temperature becomes high enough [11].
With the results depicted in Fig. 19 we can now conclude that for high enough temperature the pairing as well as the transition from a paired to an unpaired phase in a trap are also approximated well quantitatively by a translationally invariant system. Finally, we address another question related to superconductivity, namely the symmetry of the pair wave function, i.e. the question if we have an s-or p-wave superconductor. Every pure Gaussian state has a wave function of the form |Ψ = exp x, y φ σ,σ ( x, y)a † x,σ a † y,σ |0 . Since we consider a translationally invariant system, we have φ σ,σ ( x, y) = φ σ,σ ( r), where r = x − y. In Fig. 20, we have depicted the pair wave function in momentum space, |φ ↑,↓ (k 1 , k 2 )| 2 for a translationally invariant and trapped system for u = −6 and half-filling. We find for all trapping potentials a spherically symmetric wave function in agreement with s-wave superconductivity.
In summary, we find that a weak external confinement alters only slightly the physics of ground and thermal state of the attractive Hubbard model, so that experiments that require only a weak confinement will be well suited to gain further understanding of the 2d Hubbard model.

VI. DYNAMICAL PROPERTIES OF THE HUBBARD MODEL
Current experimental setups with ultra cold gases in optical lattices have opened the door to an investigation of the dynamical behavior of fermionic lattice systems. To this end, it has become possible to study processes like the BEC-BCS-crossover or dynamical quenches in the laboratory. Unfortunately, existing numerical techniques, like QMC, are not capable of describing the physics of time-evolution, so that numerical benchmark results for these processes are still missing so far.
In this Section we apply the machinery of timeevolution within gHFT derived in Sec. III A to dynamical processes as they can be realized by current experimental setups. In the following we consider both, the case of a translationally invariant system with periodic boundary conditions as well as the case of an external trapping potential.
A. Dynamical process of the translationally invariant Hubbard model In this Section we want to learn more about dynamical processes of the translationally invariant Hubbard model with periodic boundary conditions. To this end we have investigated a linear ramp of the interaction strength u for various ramping times T f , addressing the question if we can achieve an adiabatic evolution of the system. As an example, we have started from the ground state for u = −7 and have increased u linearly with ramping times T f = 1, 10, 100 and compare the energy during the evolution with the ground state energy for various u (black triangles). The evolution is adiabatic for all negative u when T f = 100, while adiabaticity is no longer given for u positive. This problem cannot be resolved by longer ramping times (a). Further, starting from u > 0 and ramping to u < 0 leads to the same observation: As soon as u = 0 is reached the evolution of the system does not follow the ground state energy any more (b).
to the final value u = 7. We have used ramping times T f = 1, 10, 100, 200, 400 and have compared the energy of the evolving system, E(u, T f ), with its ground state energy, E 0 (u), for a given u.
As a result, we observe first that E(u, T f ) has the same form for T f = 100, 200, 400, so that we have only depicted the evolutions for T f = 1, 10, 100 in Fig. 21. We find that for T f = 100 the evolution is adiabatic for all negative u. However, from u = 0 onwards, the energy starts to deviate more and more from the true ground state energy. An increase in ramping time does not help to resolve this problem (Fig. 21a). In Fig. 21b we have depicted a ramp from u = 2 to u = −2, and also find that adiabaticity is destroyed when we cross the value u = 0.
We can gain an understanding of this phenomena by looking at Fig. 18: We see that for u = 0 the ground state undergoes a phase transition from a paired to an unpaired phase. However, as we have depicted in Fig. 22, the dynamical evolution of the system does not show this transition: When we start from u = −2, where we are in a paired phase, the system remains paired even when we ramp to the repulsive side (solid lines). On the other hand, the unpaired system from which we start for u = 2 remains unpaired even when we ramp to positive values of u (dashed lines).

B. Dynamical process for the Hubbard model in a trap
Finally, we investigate dynamical properties of the Hubbard model in a trap. Here, we want to investigate how the attractive Hubbard model reacts to a change of the trapping potential V t . To this end, we start with the ground state of the system in a shallow trap of strength V t = 0.1 and then squeeze the trap, by changing the trapping potential linearly in time, to a final trapping potential of V t = 0.25. As an example, we have considered an interaction strength of u = −6 and we have set the chemical potential to half-filling (here: µ = 3) (I). When we have reached the final value V t = 0.25 we reverse the process and now decrease the trapping potential linearly to its initial value of V t = 0.1 (II).
As the first question, we study if our algorithm allows for an adiabatic evolution, i.e. the energy of the system for process (I) and (II) is the same. To this end we have investigated ramping times of T f = 100 and T f = 200, and find a reversible process in both cases (see Fig. 23).
Next, we are interested in the evolution of the pairing under the change of the trapping potential. The result is depicted in Fig. 24. We find that the process gets reversible when we go to longer ramping times, and that the evolution of the pairing gets smoother. Further, we can understand the decrease (increase) of the pairing strength with an increasing (decreasing) trapping potential by looking at Fig. 6: The change of the trapping potential can be understood as a change of the chemical potential µ for fixed value of u. . We see that the process gets reversible when we go to longer ramping times, and that the evolution of the pairing gets smoother.

VII. CONCLUSION
In this work we have formulated the gHFT for fermionic lattice systems with a two-body interaction in terms of the covariance matrix, and have derived numerical methods to find the time-evolution, the ground and thermal states in this framework. We have benchmarked our methods with the two-dimensional translationally invariant Hubbard model, and have found that our numerical methods work very well to solve the general Hartree-Fock theory. Our algorithms run fast, stable and can be applied to inhomogeneous systems, as well as to large system sizes, since all algorithms are formulated in terms of the CM, leading to a quadratic scaling in the system size. However, we also found certain limitations of our approach that are due to the fact that gHFT is intrinsically not capable of describing certain strongly correlated phases, like the Mott phase at finite temperature.
In this respect, we believe that the numerical methods developed in this work represent a powerful tool to simulate fermionic systems in optical lattices, at least in the weakly interacting regime. Further, we have demonstrated that our algorithms can be easily applied to inhomogeneous systems, as they are realized, e.g., by the Hubbard model with an external trapping potential.
Concerning the question if properties of the system in the thermodynamic limit can be obtained via finite-size scaling the results presented in Fig. 7 and table I suggest that already a moderate system size of about 100 sites is enough to get a good estimate of the properties of the attractive Hubbard model in the thermodynamic limit. In addition, our numerics done for a lattice size up to 16 × 16 show that energies can be obtained with only a slight increase in the running time, while correlation functions need some more effort, but can still be obtained in a reasonable amount of time Thus, in a future work, we aim at applying our approach to some particular experimental setup, which will require an extension of our current algorithms to much larger system sizes, and constitute a further benchmark of our methods.
In the following, we will show how the CM evolves under an evolution in imaginary time with the quadratic, state-dependent Hamiltonian H Q = i ijh ij c i c j , i.e., we have to calculate (c.f. Eq. (18) and we obtainΓ = −4(ΓhΓ +h).

Evolution with the interaction Hamiltonian
In the following, we will show how the CM evolves under an evolution in imaginary time with the interaction Hamiltonian H(T, U ) = i ij T ij c i c j + ijkl U ijkl c i c j c k c l , i.e.
we have to evaluate (c.f. Eq. (18) We split the Hamiltonian H(T, U ) in two terms, H q = i kl T kl c k c l and H I = klmn U klmn c k c l c m c n , and perform the calculation for both terms independently. First, the calculation for H q is the same as for H Q , so that we can deduce from Eq. (B1) We treat the two terms independently, starting with i,j,m,n =k,l U ijmn Γ ij Γ mn : Next, we calculate i,j,m,n =k,l U ijmn Γ ik Γ jm Γ nl : with the help of Lagrange multipliers. An expression for the energy has already been given in Eq. (14). The entropy S = −tr[ρ ln ρ] is given by S = M ln 2 − 1 2 tr [(1 + iΓ) ln(1 + iΓ)] [43], so that we have to solve the following optimization problem: Then, using Eqs. (C1) and (C2) we obtain the following necessary conditions for a minimal entropy: [h F (Γ), Γ] = 0, (C5) h F (Γ)(1 + Γ 2 ) = 0. (C6)