Exact and many-body perturbation solutions of the Hubbard model applied to linear chains

This study examines how the GW approximation, one of the techniques covered by Green's functions and on many-body approximations (GFMBA), fares compared to the treatment of the Hubbard model solved using an exact diagonalization (ED) approach. We show that, for small linear chains, the GW approximation corrects the usual mean-field (MF) approach by reducing the total energy as well as the magnetization from the MF approximation. The energy gap shows also a better agreement with ED, especially in even-number of atoms systems where no plateau is observed below the predicted phase transition as in MF approximation. In terms of density of states, the GW approximation induces quasi-particles and side satellites peaks via a splitting process of MF peaks. At the same time, GW slightly changes the localization (e.g., edges or center) of the states. We also extend to GW approximation the L\"owdin's symmetry dilemma and show that GW predicts a paramagnetic-antiferromagnetic phase transition at a higher Hubbard parameter than MF.


Introduction
The Hubbard model has been introduced more than 50 years ago 1 and is now a reference approach to describe correlations in various types of physical systems, ranging from electronic systems , 2,3 optical trapped and ultra-cold atoms, 4,5 and quantum simulation. 6 It is thus a versatile description of the correlated systems that is the subject of many fundamental studies, 7-9 even if it remains a simplified approach, containing only few parameters.
The Hubbard model is of particular importance to the description of magnetic properties of nano-systems. 3,10,11 Indeed, it includes the necessary interactions between electrons of opposite spins to predict the total spin as well as spin symmetry, which are key properties in the field of magnetic nano-materials and spintronics.
Despite its apparent simplicity, the Hubbard model is very difficult to solve. The exact treatment is performed by the so-called the exact diagonalization (ED) [12][13][14][15] that requires switching the point of view from the one-body basis usually related to an independent electron picture to a many-body basis. While the ED technique is conceptually well understood, its applicability is limited because of its heavy numerical cost, due to the fact that the size of the Hilbert space grows rapidly with the number of particles. For example, its recent application to nano-graphene has been limited to 16 carbon atoms only. 10 In the electronic domain, a mean-field (MF) approximation is often introduced to approximate the Hubbard model. 3,10,11,[16][17][18] This approximation allows to treat larger systems (up to thousands of atoms) because it can be resolved using a one-body basis set. However, the MF approach misses all the correlation effects introduced by the full Hubbard model. It is thus important to compare in details the MF results to ED or other advanced techniques (e.g. quantum Monte Carlo or configuration interaction) 10,17,19 to determine the validity of the MF approximation for different systems as a function of the model's parameter.
Alternatively, correlations can be re-introduced in the description of the electronic properties using a more tractable computational technique compared with ED, such as the Green's functions formalism, together with many-body approximations (GFMBA). The GW approximation is one among a number of other GFMA approaches, such as the second-order or T-matrix approximations 20,21 and it was first introduced by Hedin. 22 Recent reviews on GW can be found in Refs. [23,24]. The first one focuses mainly on the physical foundations of the GW approximation while the second highlights results of current implementations in quantum chemistry. Additional reviews and textbooks also provide good insight about the GW approximation. 20,[25][26][27] Although known for a long time and widely used in quantum chemistry and ab initio computations, it is not clear which part of the correlation many-body approximations include.
It is not clear either to what extend, for which properties and for which systems, GW improves the MF performance compared with ED for a given range of parameter values.
There has been some research devoted to this topic in the late 90's [28][29][30] but most of them focused on the density of states (DOS) and the total energy (E tot ), leaving unexplored the local density of states (LDOS) as well as the symmetry of the ground state. To address this shortcoming, the LDOS as well as density of the total states are investigated in this report.
A recent study 19 showed that the so-called Löwdin's symmetry dilemma 31 is not satisfactorily resolved using GFMBA, and more specifically in second-order Born approximation (SOA). Löwdin's symmetry dilemma was first described within the MF approximation where, in some cases, solutions with unphysical spatial and spin symmetries are more stable. More precisely, the authors of Ref. [19] highlight the case of one-dimensional (1D) single-orbital Hubbard chains with an even number of atoms for which the ED predicts a paramagnetic (PM) ground state for all values of the Hubbard parameters whereas MF predicts an antiferromagnetic (AFM) ground state above a critical Hubbard parameters value. If the PM symmetry is imposed in the MF algorithm, the total energy increases and deviates from the ED total energy. In other words, there is a phase transition predicted from PM to AFM in the MF approximation but not in the exact ED.
In this paper, we investigate the difference in energy between PM and AFM states calculated in MF and in GW for a wide range of parameter value. The symmetry of the predicted states and their difference in energies are of importance for potential technological applications such as spintronics. 11 As a key result of this study, we show that the GW correction extends the parameter range for which the ground states display the correct spin symmetry compare to the MF approximation. We also show that the results from SOA for linear chains 19 remain valid within the GW approximation, leading to a Löwdin's symmetry dilemma.
Analytical studies to aid in the interpretation of the reasons of success and failure of GW and other approximations have been published in the literature. 21,[32][33][34][35] These analytical results are mainly obtained for the Hubbard dimer (a chain with 2 atoms). Here, we also provide new insight about their generalization to larger systems.
The rest of the paper is organized as follows: we first introduce the Hubbard model as well as the different approximations and formalisms (Green's functions, MF, GW and ED). We then study the improvement of GW approximation on global properties of the system, such as total energy, energy gap, density matrix or magnetization before switching to frequency-dependent quantities (including density of states and scanning tunneling spectroscopy simulations). The last section also discusses phase transition and how GW fares with respect to the Löwdin's symmetry dilemma.

Hubbard model
In this paper, we focus on a single-orbital Hubbard model for fermions that contains two parameters:Ĥ where t is the hopping parameter, U is the interaction (or Hubbard) parameter,ĉ † iσ (resp. c iσ ) is a creation (resp. destruction) operator of an electron on atomic site i with spin σ, andn iσ =ĉ † iσĉ iσ is the density operator (of electron on atomic site i and with spin σ).
The signs under the summation symbol indicate that the sum runs over all pairs of nearest-neighbor atomic sites.
The first term of the Hubbard Hamiltonian in eq. 1 is the tight-binding Hamiltonian and can be written easily in the one-electron basis but the second term (Hubbard or interaction term) is the product of two density operators (consequently of four creation or destruction operators) such that it is called a two-body operator. It cannot be written in the one-electron basis, and a many-electron basis is required.

Mean-field approximation
The MF approximation consists of approximating the Hubbard Hamiltonian of eq. (1) with one-body operators only, so that it can be written in the one-electron basis. In the single orbital-per-atom approach used here, the one-electron basis is of size 2N at since it consists of basis vectors of localised electrons on atomic sites with a given spin (up or down), responsible for the factor 2 in the basis size. Doing so, the size of the basis is significantly reduced compared to the many-body case (see Exact diagonalization section).
In practice, the density operators (n iσ ) are decomposed in the mean value of the operator (n iσ ) and the deviation (n iσ − n iσ ) from this mean value:n iσ = n iσ + (n iσ − n iσ ). Then, products of density operators in eq. (1) are expanded and products of deviations are dropped as they are expected to be small, leading to the mean-field Hamiltonian: and total energy: where E l are the one-electron energies obtained by diagonalizing the Hamiltonian of eq. (2), f FD is the Fermi-Dirac and n i↑ , and n i↓ are obtained from the eigen-vectors of the Hamiltonian of eq. (2). More details about MF approximation can be found in Refs. [3,11,20].
We emphasize here that, since the Hamitlonian (2) depends on mean values of density operators, it has to be solved self-consistently starting from initial guesses for the mean values. Depending on the symmetry of the initial guess (e.g., PM or AFM), the self-consistent cycle can converge towards different (stable or metastable) states. The comparison of total energies associated with the states (eq. (3)) is needed to identify the predicted ground state with the lowest total energy. Usually, the MF ground state is found by performing the MF self-consistent algorithm several times with random initial states.

GW approximation
In order to account for a portion of the correlation, a Green's function based method to compute GW approximation on top of MF computations is used in this paper, following Refs. [ 36,37]. It is based on the Feynman diagram expansion of the Green's functions explained, for example, in Ref. [20]. In this section, we first give a brief reminder of the Green's function formalism. We then present the main steps of the algorithm for a numerically approach and, finally, we give the main expressions to retrieve the useful quantities such as LDOS, total energy and STS simulated images.
In the Green's function formalism, we consider retarded (R), advanced (A), greater (>), and lesser (<) components of several quantities (Green's function (G), self energies (Σ) and screened potential (W )). We stress that these components are functions of real time instead of complex time for the total quantities. For definitions, properties and links between the components, we refer to chap.5 of Ref. 20 The different properties and links between components require to consider equations involving ≶ components or equations involving The retarded component of the Green's function (G R (ω)) obeys Dyson's equation: where G R 0 (ω) is the non-interacting retarded Green's function and Σ R (ω) is the retarded self-energy. All quantities in eq. (4) are matrices in the one-electron basis (with dimension 2N × 2N , N being the number of sites).
The retarded and advanced non-interacting Green's functions in the GW approximation (G R 0 (ω)) are taken as the Green's functions of the MF solution, given by: where |jσ is the one-electron state containing one electron with spin σ on site j, iσ| is the hermitian transpose of |iσ , and η is a small real positive parameter.
The greater and lesser (> and < indices) components of the Green's functions are given by (the expressions remain valid for the exact and non-interacting Green's functions): withf F D (ω) = 1 − f F D (ω) and µ the chemical potential.
Dyson's equation (eq. (4)) provides the exact Green's function (G R (ω)) from the noninteracting Green's function (G R 0 (ω)) and from the exact self-energy Σ R (ω). But the exact self-energy is not known and approximations have to be made.
In the GW approximation, the self-energy is approximated as the product of the Green's function (G) and the screened potential (W): where • is the Hadamard product, or element-wise product between matrices.
G ≶ (t) in time domain is found by performing an inverse Fourier transform of G ≶ (ω). The screened potential W R (ω) obeys a Dyson-like equation: where W 0 is the interaction matrix containing Hubbard U term of Hamiltonian (1) and χ 0 (ω) is the non-interacting susceptibility found by inverse Fourier transform of: where Θ(t) is the Heaviside step function.
With this proviso, the greater and lesser components of the screened interaction of eq. (7) are found from the relations between greater and lesser, and retarded screened interactions: where f BE is the Bose-Einstein distribution.
From the GW Green's function, the spectral matrix A iσ,jσ is given by: The local density of states (LDOS, n iσ (ω)) is expressed using the diagonal elements of the spectral matrix and the total density of states (DOS, D(ω)) with the sum of diagonal elements of the spectral matrix (i.e. the trace): and The matrix elements of the density matrix (n) are given by the integrals of the spectral function multiplied by the Fermi-Dirac statistics: and the local densities are the diagonal elements of the density matrix: The total energy of the system is found using the Galitskii-Migdal formula: 19,20 where E kin is the kinetic energy and E GM is the interaction energy. They are given by: whereĤ 0 is the non-interacting Hamiltonian, i.e. the tight-binding part of Hubbard Hamiltonian or the first term of eq. (1) proportional to the t parameter.
The simulation of spatially-resolved dI/dV images from STS (scanning tunneling spec-troscopy) in the interval of energies [E 1 , E 2 ] can be obtained from the spectral matrix: 36,38 where z 0 is the tip's height of the simulated ST S and λ is a length parameter that parameterizes the spatial extension of localized orbitals. (x, y, z 0 ) = r is the location where the ST S is simulated and r i are the atomic positions.

Exact diagonalization
In the exact diagonalization approach, we deal with the exact Hubbard Hamiltonian of eq. (1). This implies a treatment in the many-electron basis rather than the simpler oneelectron basis and, as a result, the rapidly growing size of the Hilbert space. With only one orbital per atom, the size of the many-electron basis is given by C 2Nat where N at is the number of atomic sites and N el is the number of electrons in the system, in strong contrast to the linear scaling of the single-electron basis.
We follow the method outlined in Refs. [12][13][14][15] It consists in labelling the state with one integer I, bijectively linked to two other integers I ↑ and I ↓ by the relations: and where // represents the integer division.
Writing I ↑ (resp. I ↓ ) in binary notation gives the space configuration of the states in the spin up (resp. down) sector. Organizing the Hilbert space in this way allows one to treat only integers and to easily find the actions of creation, destruction, and density operators on each state using simple standard binary operations (bin flip, bin counting,. . . ). We explicitly show the correspondence between basis states and I, I ↑ and I ↓ for the two-site Hubbard system at half-filling in table 1. Table 1: Illustration of the indexing scheme used to build many-electron basis states for the two-site Hubbard system at half-filling. I ↑ follows from the binary notation of spin up occupation of the states from right to left: the first state for example has one spin up electron on site 0 (right) and one spin up electron on site 1 (left) such that I ↑ = 2 0 + 2 1 = 3. The same reasoning for spin down electrons leads to determination of I ↓ and I follows from the formula (20).
In the ED framework, the retarded Green's function describing the system with N el electrons involves the Hamiltonian operators acting in the Hilbert space containing states representing N el + 1 and N el − 1 electrons (Ĥ N el +1 andĤ N el −1 ). Its definition for a general system of N el electrons in the frequency domain is: 39 where Ψ N el In order to evaluate the retarded Green's function (22), three exact diagonalizations are required: the first one is done in the N el sector to find the ground state Ψ N el 0 and its energy Then an ED is completed in the sectors N el ± 1 so that matrix elements of the type From the evaluation of the retarded Green's function (22), physical quantities such as the LDOS, DOS, density matrix, and the dI/dV STS simulated maps can be computed (eqs. (13), (14), (15), and (19)). The total energy of the system, E N el 0 , is obtained using a single ED.

Magnetic order parameters
We have already mentioned that the magnetic symmetry (or order) requires a good description of the correlation. We define here two order parameters that will allow us to characterize the phase transition between the PM, FM and AFM phases. The local magnetic moments (M i ) are defined as the difference between spin-up and spin-down local mean densities and the total spin (S) is the sum of the local magnetic moments: 3 and Finally, the staggered magnetization (M ) is defined as the sum of absolute values of the magnetic moments: We distinguish the cases of even and odd number of atoms in the chain because the total spin (see eq. (24)) of the ground state is predicted to be either 0 or 1/2 for even and odd cases respectively, in agreement with Lieb's theorem. 3,40 In the even case, this implies that the ground state cannot be FM but must be PM or AFM. In the odd case, the PM phase is forbidden due to Lieb's theorem.
Total energies computed by the three methods (MF, GW, and ED) as a function of the U parameter are shown in figure 1. For even numbers of atoms, we obtain results in agreement with those reported in Ref. [19] for both MF and ED approaches (Ref [19] does not report results for odd number of atoms). Random initial occupations have been used in the case of the MF approximation such that it has to be compared with unrestricted-spin Hartree-Fock ("usHF") method employed in Ref. [19]. For small values of U , MF and GW are in good agreement with the exact ED and display a linear increase of the total energy with U . 19 GW approximation lowers the total energy from MF and agrees very well with ED until U = 3t, thereby doubling the range of very good agreement compared to MF.
The energy gap (energy interval between the highest peak in the DOS below the Fermi level and the lowest peak above the Fermi level) of the same linear chains are shown in figure 2. For an odd number of sites ( Fig. 2 (a) and (b))), the energy gap is zero at U = 0 and grows with increasing U . This growth is overestimated by MF. MF deviates significantly from ED above U t where the GW approximation follows the ED curve more closely until U 2t − 2.5t. Then it deviates in turn with ED, lowering the energy gap from MF but sill overestimating it when compared with ED. For even number of atoms (( Fig. 2 (c) and (d))), the energy gap is non-zero for U = 0. It increases with U for ED but remains constant until U = 1.5t and starts increasing around U = 2t in MF. As shown in, 19 this sudden increase in MF is due to a phase transition as the ground state for small values of U is PM whereas it is AFM for larger values. This transition will be discussed later on considering the staggered magnetization as an order parameter. According to Ref. [19], the increase of energy gap in ED is due to an increase in correlation, that is absent in MF. In MF, the correlation is  We note however that the GW approximation results considered in this section exhibit the same phase transition as in MF approximation since the latter is the starting point of former.
This phase transition is not present in the ED and is not directly visible when looking at the total energies. We discuss the change in the phase transition in a following section.
The analysis of density matrices allows to compare local densities of the different states as well as non-diagonal coupling between atomic sites. The density matrices (see eq. (15)) for ED and the PM and AFM states of MF and GW approximations are represented graphically in figure 4 for U = 2t and in figure 5 for U = 5t. The diagonal elements are local electronic densities such that a uniform diagonal corresponds to a PM state whereas an AFM is characterized by diagonal elements that exhibit high values for spin up and spin down electrons, alternatively. We first observe that the PM states computed in MF and GW have a similar appearance compared to those obtained by ED (see e.g. figure 4 a), b) and c)). This is not the case for the AFM states (compare for example figure 4 d) and e) with c)).
For U = 2t, the density matrix of the PM state in MF ( figure 4 (a)) presents a strong negative coupling between sites 0 and sites 3, i.e. the edges of the chain. This coupling is reduced in GW ( figure 4 (b)), in better agreement with ED ( figure 4 (c)). Then, the coupling between neighboring sites (site 0 -site 1 and site 2 -site 3) is also reduced in absolute value from MF to GW, still better reproducing the ED result. Considering AFM states (figures 4 (d) and (e)), the anti-ferromagnetic character is decreased between MF and GW as the contrast between alternating diagonal elements is decreased. The negative couplings between sites 1 and 3 for spin up and between sites 0 and 2 for spin down are reduced from MF to GW, but is still non-zero in GW contrary to ED. Increasing the U parameter to U = 5t (see figure 5) leads to a similar GW correction: the negative coupling between sites 0 and 3 and the positive coupling between adjacent sites are reduced. However, even if GW leads to a more accurate description of density matrices, the differences with ED density matrix are greater than for U = 2t. This is an expected feature since the GW approximation is a perturbative method in terms of the U parameter.
For the AFM states, the alternating behavior in the diagonal of the density matrix is also reduced from MF to GW, and some other minor changes occur in non-diagonal elements but both MF and GW density matrices are far away from the ED one, that shows a PM state.
The correction from MF to GW is less important than for U = 2t.

DOS, LDOS and simulated ST S maps
The previous section focused on global properties such as total energy, total staggered magnetization, energy gap, and density matrix. We now turn to frequency-dependent properties.
These quantities are also of interest as they can be interpreted as single-particle features (energy, density,. . . ) in the MF approximation and can thus offer further insight into the effects of many-body interactions and the missing correlation in MF by comparison with ED results.
They are also measurable properties or can be used to compute measurable quantities (e.g., ST S maps, optical absorption, plasmons,. . . 36,37 ) We start by investigating the effect of GW on a two-site system so that we can make the connection with previously published analytical studies 21,32,33 and then extend our study to larger systems as those in previous sections. We note that in the cited analytical studies, the so-called "GW approximation" is actually the "G 0 W 0 approximation" in the sense that no self-consistency procedure was introduced and the GW results are identified to the results obtained after the first step of computation.
All results shown in this section are obtained at an intermediate value of the Hubbard parameter (U = 2t) and all the Fermi levels have been aligned to 0 in the study of DOS and LDOS for better visualization. The DOS for the two-site system in MF, GW, and ED are shown in figure 7. The DOS in MF consists only in two doubly degenerated peaks at ω = ±t that can be attributed to single-particle states with bonding (resp. anti-bonding) symmetry for the state below (resp. above) the Fermi level. Bonding states (resp. anti-bonding states) are characterized by a non-zero (resp. zero) intensity at the mid-point between the atoms. This is consistent with ED results. In GW, these main peaks are shifted away from the Fermi level at ω ±1.21t, increasing the gap and resulting in a value closer to the result of ED (ω ±1.24t). However, these main peaks do not integrate to 1 in GW and ED and cannot be attributed to single-particle states but to quasi-particle (QP) peaks, as explained in Ref. [ 23]. The symmetry of the QP states is preserverd and corresponds to ED (see figure 8 (a) for the QP peak below the Fermi level).
A closer look at the GW approximation DOS reveals smaller satellite peaks at ω ±2.69t for GW and ω ±3.24t for ED, surrounding the QP ones. 23 This is another evidence of the improvement from MF to GW although the peaks are not exactly at the correct energies. However, the GW approximation predicts many satellites (of very small intensities) whereas only one on each side of the Fermi level is obtained in ED. As shown analytically in Ref. [34], the satellite below (resp. above) the Fermi level have an anti-bonding (resp. bonding) character, that is the opposite symmetry than the nearest QP peak. This feature is reproduced well by the principal satellites introduced by GW as shown in figure 8 (b).
The DOS for N at = 3, 4, 5 and 6 chains are shown in figure 9. As previously in the two-site system, we observe the energy correction of main peaks and the appearance of We now focus on the simulated ST S maps for N at = 3. The one to one comparison of the STS maps is needed to correctly assign the DOS peaks to states with a given symmetry for each approaches. The lower state in MF (E 1.67, see the * symbol) is compared to less intense peaks below 2t in GW and ED (see fig. 10 a)). For the second lower energy state in MF, the GW approximation predicts two near peaks that better reproduce the STS map of ED greater peak (at E ±1.41t) when integrated together (see the + symbol). The initial state in MF at E ±1.67t thus splits into a satellite that is further away from the Fermi level and a QP peak that is closer to the Fermi level. The satellite is more localized at the ends whereas the QP peak is more centrally located than the initial MF state.  . Red curves correspond to MF approximation, green curves to GW approximation and blue curves to ED. All Fermi energies have been aligned to 0 for better visualization. +, * and # symbols are used to mark the peaks as references for figure 10. Figure 10: (a) Simulated STS maps for the particle peak of MF at E −1.66t, the satellite at E −2.25t for GW and the peak at E −2.09t for ED, for the N at = 3 Hubbard chain at half-filling. (b) Simulated STS maps for the particle peak of MF at E −1.33t, the two close peaks at E −1.37t and E −1.46t for GW and the large peak at E −1.41t for ED, for the N at = 3 Hubbard chain at half-filling. (c) Simulated STS maps for the QP peaks just bellow the Fermi level, for the N at = 4 Hubbard chain at half-filling. MF peak is at E −0.75t, GW peak is at E −0.80t and the ED peak at E −0.81t. For each method and each subplot, the spin up and spin down channels are shown separately and the sum of the channels is labelled "Total". The intensities have been independently normalized. The size of the scanning area is (N at * d i−a ) × d i−a in each case with d i−a the inter-atomic distance. +, * and # symbols are used as references for peak identification (see figure 9).
We observe that the trend of splitting into a more edge-localized state further away from the Fermi level and a more central-localized state closer to the Fermi level seems to be a general feature of the GW approximation in linear chains.
Another important feature of the GW approximation is that, in chains with an even number of atoms, the states tend to be more spin-symmetric than in the MF approximation.
We illustrate this by looking at figure 10 (c) that shows the simulated dI/dV maps for the first peaks under the Fermi level for the N at = 4 chain. The dI/dV map is more spinsymmetric in GW than in MF, correcting in the right direction since ED dI/dV map is totally spin-symmetric.
Symmetry dilemma and phase stability around the transition Löwdin's symmetry dilemma has been observed in GFMBA using the SOA and a high value of U (U = 4t), well above the phase transition . 19 We discuss here the case of a more complex self-energy (GW) with a focus near the phase transition in MF of even systems (no phase transition are predicted in odd systems).
First, we examine total energies of the N at = 4 chain for U = 4t, considering both PM and AFM states for MF and GW and the ED ground state (PM) (Figure 11). We observe that for both approximations (MF and GW), the energy of the AFM state is lower than the energy of PM state, even though PM is the symmetry of the exact ED ground state. This is the illustration of Löwdin's symmetry dilemma: breaking the spin-symmetry in MF and GW results in a lower energy that is closer to the exact one. This is an evidence that Löwdin's symmetry dilemma might also occur in the GW approximation, generalizing the observation from Ref. [19] for SOA. We note however that the gap between PM and AFM energies is reduced from MF to GW (from 0.96t to 0.70t), such that the GW approximation still destabilizes the PM in comparison with MF. Figure 11: Energies for the N at = 4 chain at U = 4t. The energies are shown for the PM (squares) and AFM (crosses) states in the MF (red) and GW (green) as well as for the ED (blue circle) ground state, which is PM. The differences between the exact result and the energies of the ground states from MF and GW as well as the differences between PM and AFM from the same approximation are depicted in red (for MF) and green (for GW) double arrows.
Looking closer at the transition, we observe the difference between the PM state and the AFM state. The transition occurs in MF around U = 1.8t. Figure 12 a) figure 12 b) where the AFM energy in GW is more uniformly reduced from MF energies (see figure 12 b)).

Conclusion
We have investigated the fundamental properties of linear Hubbard 1D chains in the MF and GW approximations and in ED. We highlighted the improvements from MF to GW as well as some limitations of the GW approximation as a perturbation theory. We report a systematic study of several interest quantities of the systems as a function of the U parameter of the Hubbard model.
Our study showed how GW improves global quantities, such as the total energy, the energy gap, the staggered magnetization, and the density matrix. The improvement from MF to GW is globally significant in the small to intermediate regime of the U parameter.
In the large U regime, the GW does not correct the MF approximation in a very effective manner. We stress the surprising behavior of GW density matrices as U increases from intermediate to large values. Indeed, whereas some quantities, such as the total energies, evolve in the same way in GW than in ED, the negative coupling between sites 0 and sites 3 for N at = 4 increases (in absolute value) from intermediate to high U values in GW, that is opposite to the ED evolution.
We also studied frequency-dependent properties, such as the DOS and the simulated STS maps. To the best of our knowledge, we report the first comparison of STS maps or local DOS between GW and ED. We showed that our numerical results for the two-site system are in line with previous analytical results and extended the study to larger systems. We identified different trends to be general in different lengths of linear chains, such as the energy shifting of QP peaks and the splitting of single-electron peaks (in MF) to QP and satellites in GW. In this case, we also observed that the resulting peaks closer to the Fermi level are more central-located than the initial state and the peaks moved further from the Fermi level are more located at the edges.
It is noteworthy that for applications in several electronic systems such as carbon-based or graphene nano-materials, typical values of the U parameter reported in the literature are in the range of values where GW is found to significantly improve upon the MF results. For example, Refs. 3,11 reports values from 0.9t to 1.3t and Ref. 10 uses a value of 2t.
Finally, we considered Löwdin's symmetry dilemma that was recently investigated in GFMBA and more precisely with the SOA. We showed that the symmetry dilemma is also present in the GW approximation. We also focused on the region of U slightly above the phase transition in MF and revealed the particular fact that the GW approximation predicts a phase transition at higher U than the MF approximation. The GW phase transition arises in an abrupt manner: in the investigated system, the energy of the PM state in GW undergoes a strong increase, making it crossing the energy of the AFM state. The prediction of the correct spin-state (ground state) is crucial for application, e.g. in spintronics Despite its wide use in quantum chemistry and ab initio computations, what part of the correlations the GW approximation represents remains unclear. It is indeed difficult to predict which part of the correct correlations is introduced in a given system. This is an advantage of dealing with it inside the Hubbard model and small systems: we are able to perform calculations implementing all the correlations (ED). An open question is thus "is the GW behavior observed in small linear Hubbard chains of the same kind as in larger and higher dimensional systems and for other models (e.g., ab initio)?" As a possible direction, let us stress the unknown physical origin of the sudden increase observed in ∆E leading to specific magnetic orders. The changes in the predicted phase transitions and density properties in higher dimensional systems, such as 2D systems, might also be investigated in further work.