A new algorithm for investigating strongly correlated systems using Hubbard model

In this work for an interacting electrons system that obeys the Hubbard model, a new quantum Monte Carlo algorithm introduces for calculation of average Green function. This algorithm is applied to investigate the effects of electrons onsite Coulomb repulsion on the band structure of a square lattice in both single-site approximations such as dynamical mean field theory (DMFT) and multi-site approximations such as effective medium supercell approximation (EMSCA). The advantages of our algorithm in comparison to the Hirsch-Fye algorithm and also the Blankenbecler, Scalapino, and Sugar (BSS) algorithm are the elimination of instabilities resulting from the Metropolis algorithm in the accepting and rejecting configurations, stability at low temperatures, the elimination of systematic errors resulting from the update of the Green's function in the quantum Monte Carlo process, and considering different probabilities for each possible configuration. Finally, by using our algorithm, it is possible to calculate the interacting three-dimensional system's band structure and the density of states that obey the Hubbard model. We have applied our algorithm to an interactive two-dimensional square lattice. As a result, phase transition boundaries can be easily recognized through calculated bands for different repulsions. Our results show that critical Coulomb repulsion values for Mott transition are u c =9.05t and u c =2.4t for DMFT and BEMSCA respectively. This means that DMFT significantly overestimates band splitting due to electrons' Coulomb repulsion. We found by starting at low repulsions and then increasing electrons' Coulomb repulsion, a partially flatted valance band around the center of the first Brillouin zone appears, but disappears at high repulsions.


Introduction
Studying the effects of the on-site Coulomb interaction on the two-dimensional square lattice is one of the interesting topics of condensed matter physics because the two-dimensional square lattice in the presence of the on-site Coulomb interaction can be a suitable candidate for high-temperature superconductors [1]. Hubbard model [2] describes the dynamics of interacting electrons in a lattice. For a specific on-site Coulomb interaction, Hubbard Hamiltonian is as follows. Methods such as dynamic mean field theory(DMFT) [17,18], the variational cluster approximation [19], the dynamic cluster approximation(DCA) [20], and cellular dynamic mean field theory [21] are used to calculate the phase diagrams of the Hubbard model, and the magnitude of the gap in the spectrum of fermions. Unfortunately, there is a wide range of u c values (u c is the minimal value of the on-site interaction at which the gap in the spectrum opens) calculated by various methods and approximations [4].
Finite-size quantum Monte Carlo method (QMC) using the Hirsch-Fye algorithm and Blankenbecler, Scalapino, and Sugar (BSS) algorithm, are the suitable methods for studying this type of Hamiltonian [22][23][24][25]. This model has been widely used to explain superconductivity of high critical temperature [26] based on local electrons' Coulomb repulsion in the CuO 2 plane [27][28][29][30]. The effective medium theory is used widely for these systems. The advantage of it concerning finite size exact diagonalization is that the medium effect on the system is considered. The BSS algorithm was introduced by Blankenbeckler, Scalapino, and Sugar [31] and has been used for various studies of lattice systems [32,33], and magnetic-impurity systems [34]. This algorithm is unstable at low temperatures. For more details on the BSS algorithm, refers to the reference [31]. The BSS algorithm is as follows.
1) Determine the desired initial configuration i s ( ) , which is usually assumed to be zero.
2) Green's function calculation in time slice 1 = through the following relation where R  is the relative weight of two configurations; ie

5)
Calling a random number in the interval [0,1] 6) If 'r' is greater than the random number, the configuration is accepted and the Green function of the new configuration must be updated via the following relation If 'r' is smaller than the random number, the previous configuration is accepted as the new configuration. 7) After completing the lattice points in the time slice , Green's function should be calculated in the next time slice through the relation 8) Go to step 3 9) After completing all time slices and all spatial points, the traveler returns to the first time slice and the new Green's function must be recalculated for the updated configurations through 10) After several sweeps in the space-time lattice, the calculation of physical quantities begins, and by using a number of random and independent configurations, averaging over these configurations begins.
One of the weaknesses of the BSS algorithm is that after passing time and many time using the relation (V), the rounding errors increase and the accuracy of the Green's function decreases. Hence the BSS algorithm becomes unstable at low temperatures [23].
A modified algorithm has been introduced by Hirsch and Fye (HF) [35] for magnetic-impurity systems that are stable at low temperatures. In the HF algorithm that proposed for magnetic impurity, the Green function is calculated only in the places where the interaction takes place (the place of impurities), through Dyson's equation. In the Hubbard's model interacting electrons are on all lattice points. Therefore, Green's function must be calculated through Dyson's equation for all lattice points, and equation (III) should be used to update the accepted Green function. The similarity of the two algorithms is using of relation (III) to update the new configuration Green function. In the HF algorithm, using equation (III) at low temperatures does not lead to instability because using of relation (III) is not the origin of the instability in the BSS algorithm. In the process of calculating the Green function in the next time-slice and successively using the equation (V), the rounding errors will quickly increase and cause the instability of the BSS algorithm at low temperatures. In the HF algorithm, the Green function is not calculated in the next time slice through equation (V), but it is directly recalculated. Therefore, the HF algorithm requires much more time than the BSS algorithm [23].

Our algorithm a) Equation of motion
In the imaginary time τ, the equation of motion for electrons corresponding to the system Hamiltonian, Eq. (1), is given by, where I in the spin space is a 2 × 2 unit matrix, ij ij tt   = I and two particle Green function matrix The single particle equation of motion corresponding to Eq. (2) in the effective medium theory is given by in which self-energy matrix component The relation between interacting single particle Green function matrix G(τ ), average single particle Green function ()  G and two particle Green function G 2 (τ ) in the real and spin spaces that obtains from Eqs.
(2) and (5) is Eq. (6) can be reduce to a super cell by using effective medium supercell approximation (EMSCA) [34,35]. By applying EMSCA to the interacting system, Eq. (6) reduces to imp imp imp sc sc sc 2sc sc sc Equation (7) could be written as Equation (8) and imp imp sc sc sc 2sc where g sc (τ ) is called supercell cavity Green function with no interaction on its sites. Equations (9) and (10) imply that our calculations reduced to a supercell with sites I, J = {1, 2, …, Nc}. Supercell average of Eq. (10) leads to Eq. (9) where imp 2sc sc sc Another method for deriving Eqs. (9) and (10) is applying effective medium supercell approximation (EMSCA) to the interacting system total action, S [36]. Note that the EMSCA in the c N1 → limit recovers the DMFT and in the c NN → limit, the supercell Green's function and the supercell self-energy converts to exact average Green function and the self-energy of the whole effective lattice. Also, in this limit, the effective wave function becomes the Bloch wave function for the whole effective lattice.
To solve Eq. (10) we should find the two particle supercell impurity Green function By imaginary Fourier transform of Eq.(12), in the actual calculation the following imaginary time discretized are using [37]     (14) in which   s 0 e V is a diagonal matrix in the imaginary time space. This relationship is the basis of our algorithm, which after generating each configuration on the time slices, the Green's function corresponding to the desired configuration is calculated directly from relation (14). Average over all possible Ising fields' configurations of each site in the impurity super cell Green function imp sc ()   G gives us super cell real space average Green function sc () where, Z is the partition function , Deriving the relation between k-space and real space self energies and also the relation between k-space and real space average Green functions in the EMSCA is as follows. In the effective medium supercell approximation (EMSCA), electrons wave functions in EMSCA are periodic concerning supercell length where i r is the vector position of lattice points in the tight binding model. Eq. (17) should have Bloch in which qi U (r ) is the periodic part of wave function with respect to lattice vector. By inserting Eq. (18) into Eq. (17) we have where   Since n K i U (r ) are periodic with respect to the lattice vectors, i r , so we have n n 2 2 By substitution Eq. (21) into Eq. (20), we have Also, in the EMSCA, super cell self energy is periodic with respect to supercell lengths, which leads to ci i iq.N a e 1 = , hence Eq. (19) and N c different step function self energies in the FBZ, The average Green function in terms of average wave function and the non-interacting Green function is defined by n n Note that output self-energies of EMSCA are step functions in the FBZ that are discontinuous at small regions boundaries. To calculate the effective density of states, we need k-space continuous self-energies.
For this purpose, we use BEMSCA [42], n IJ n i(k K ).r sc sc Now we apply our method to an interacting two-dimensional square lattice and present obtained results. Figure 1(a) shows a two-dimensional lattice that is divided into similar supercells of 4 sites where each supercell has the original lattice symmetry, and Figure 1

C) Our Algorithm
Implementation of our algorithm is as follows: 1. An initial guess for self energies which is usually zero. hence the probabilities of the configurations are not same. In our algorithm, to avoid this problem, after producing each configuration, the probability of its occurrence is also calculated. In other words, in calculating the average Green's function, each configuration has its own probability.
In both BSS and HF algorithms, the Green's function of the new configuration is calculated directly from the update of the Green's function of the previous configuration. In addition to the dependence of the Green's function on the new and old configurations, due to the approximations in the relationship between the two new and old Green's functions, calculation errors are created, and these errors increase after calculating several Green's functions. To avoid these errors, in our algorithm, after generating each configuration, the Green's function of that configuration is calculated separately. The advantage of our algorithm is that the effects of the randomness of the selection of configurations and the increase of calculation errors due to the calculation of several Green's functions are removed, and in this way, the same results are obtained in several different executions. Considering all the configurations, it may increase program execution time, but despite advanced shared memory computers and full use of all processor capacity, this problem can be solved by parallel programming.
In parallel programming, all the capacity of processors that are connected to shared memory is used.

c) Renormalized band structure calculation
After finishing the execution of the main code, we can calculate the three-dimensional band structure and the density of states of the interacting system from the resulting self-energy. Having continuous self-energy in FBZ, we can calculate the effective Green's function as follows, Thus, we can calculate density of states obtained from effective Green's function as follows, For an exact effective medium system, which includes all lattice points, the imaginary part of the selfenergy tends to zero, so the real part of the self-energy renormalizes the band structure as follows.
In general, the self-energy of the effective medium is not exactly calculated. In the single-site dynamical mean field approximation (DMFT), where the self-energy in the FBZ is independent of the wave vector, the imaginary part of the self-energy is non-zero. Also, in the multi-site EMSCA approximation, where the self-energy is a step function, the imaginary part of the self-energy is not zero. Therefore, it is not possible to directly use the eq. (39) to calculate the renormalized band structure. In other words, the imaginary part of self-energy plays an important role in renormalizing the band structure. Therefore, our method to calculate the renormalized band structure is as follows.
For each wave vector, any energy that maximizes n N(k,E) N(nk,E) =  is the band structure in that wave vector [36].
where n,k E is the renormalized band structure. Using this renormalized band structure, the real density of states can be calculated as follows.
Our method to calculation of the renormalized band structure is applicable not only to the output of the BEMSCA method, sc (k;E i )   +  , which is continuous in the FBZ, but also to the output of the DMFT method, sc (E i )   +  , which is independent of the wave vector. In this work, using our method, the 3d renormalized band structure is calculated for the 2d square lattice, and through these 3d bands, the phase transition can be clearly observed.

Results and discussion
Now, to see the advantages of our algorithm, we apply this algorithm to investigate the band structure and the density of states of the square lattice with on-site Coulomb interaction. Now to highlight inter-site correlation effects in a 2d square lattice with on-site Coulomb repulsion, we compare calculated renormalized band structure and density of states of an electrons interacting system that obey Hubbard model in the single-site (DMFT) and multi-site N C =4 beyond effective medium supercell approximations by using our algorithm for solving impurity Green function with parameter T/t=0.28. The 2d square lattice system with on-site Coulomb repulsion is interesting because it can be a suitable candidate for explaining high-temperature superconductors mechanism [1]. The main advantage of our method concerning the previous works is not only preserving multi-site scattering effects but also calculating renormalizing band structure to obtain the renormalized density of states.
To observe the effect of Coulomb interaction on the imaginary time Green function, the on-site Green function in the single-site approximations (DMFT) and the multi-site approximations (EMSCA) in different Coulomb repulsions are shown in Figure 2. In the multi-site approximation, due to the consideration of inter-site effects, the Green's function goes to zero at a lower Coulomb interaction intensity, which indicates that in the multi-site approximation, the phase transition occurs at lower Coulomb repulsion than the single-site approximation. In the single-site approximation (DMFT), the self-energy is completely local in the real space, in other words, in the wave vector space, it is independent of the wave vector. In the multi-site approximation(EMSCA) , the self-energy is discontinuous in the form of step functions in the Brillouin zone, c 1 N { (K , E i ),..., (K , E i )}  +   +  .Therefore, using the BEMSCA approximation and using equation (34) , it is possible to calculate the continuous self-energy in the first Brillouin zone, (k, E i )  +  . In Figure 3, the real part and the imaginary part of the self-energy are drawn for the single-site approximation (DMFT) in u=9.1t, and for the multi-site (BEMSCA) in u=2.4t. To see the effects of repulsion increasing on the inter-site correlations, the density of states in the DMFT and multi-site BEMSCA are compared for different interaction repulsions strengths from weak to strong repulsions. From Figure 4(a), in the weak interaction regime (u=1.0t), the density of states of both approximations is almost the same. But as the interaction strength increases, these two approximations show different results, for example in the u=4.0t, the single-site approximation predicts a metallic phase, but BEMSCA predicts a Mott insulating phase (Figure 4(b)). This result is the superiority of the multi-site correlation effects concerning single-site approximation. As the interaction strength increases further, both predict a Mott insulating phase, but the gap between conduction and valence bands is larger for the case of multi-site approximation (Figure 4 (c)). In the weak interaction (u=1.0t), the density of states of both approximations is almost the same. But as the interaction increases, for example at the u=4.0t, the single-site approximation predicts a metallic phase, but the multi-site approximation predicts the Mott insulating phase. As the interaction strength increases further, for example for u=10.0t, both predict the Mott insulating phase, but the gap between the conduction band and the valence band is larger in the multi-site approximation.
Through calculated renormalized band structure, the phase transitions can be observed and the phase transition boundaries can be identified. In Figure 5, the band structure and density of states of twodimensional square lattice in DMFT (N c =1) are plotted. The critical interaction value for N c =1 is u c =9.05t, in other words, for repulsions smaller than this critical value, the system is in the metallic phase, but the phase transition occurs at u c =9.05t, and for u>u c , the valance and conduction bands are separated, and the system is in the Mott insulation phase. From Figure 5(b), at u=u c =9.05t, the energy gap between the partially flatted valance band and the conduction band is zero. The central peak in the density of states at the Fermi energy is due to this valance flattening. Our calculated bands and dos for Nc=4 are shown in Figure 6. The critical interaction value for Nc=4 is u c =2.4t, It is significantly smaller than the single-site approximation, which indicates the importance of considering inter-site correlations. In the multi-site approximation, the flattening of the valance and conduction bands in u=u c =2.4t can also be seen. From Figure 6(b), the height of the central peak is higher in the multi-site approximation, this could be explained by increasing the flat area of the valance band due to multi-site effects. Finally, the phase transitions in the calculated renormalized band structures were more accurately observed. Figure 6 : shows the calculated renormalized band structure and dos in Nc =4 beyond supercell approximation for (a) u=1.0t , in this case, the system is in a metallic phase.(b) for u=u c ==2.4t , that is a semi-metal.(c) for u=4 .0t, the system is in the Mott insulation phase.
To show the flattening of the valance band, its contour plot is drawn in Figure (7) in the multisite approximation for u=uc=2.4t, and the band structure is plotted in the paths that pass through the center of the first Brillouin zone. The flattening of the valence band can be seen through band structure in the determined paths at the Fermi level.

Figure 7 :
Shows the contour plot of the band structure in the multi-site approximation in u=uc=2.4t. The above diagram shows the band structure in the specified vertical path that passes through the center of the first Brillouin zone, and also the right diagram shows the band structure in the specified horizontal path that passes through the center of the first Brillouin zone.
The gap (Δ) between the bands according to the strength of the Coulomb interaction is shown in Figure 8. Through this figure, it is clear that the phase transition in the multi-site approximation occurs at a lower Coulomb interaction than the single-site approximation, which is due to the contribution of the inter-site scattering diagrams.

Conclusion
In this work, we introduced a new algorithm to investigate the electrons Coulomb repulsion effects in the Hubbard model. Advantage of our algorithm with respect to the HF and BSS algorithms are: the elimination of instabilities resulting from the Metropolis algorithm in accepting and rejecting configurations, stability at low temperatures, the elimination of systematic errors resulting from the update of the Green's function in the quantum Monte Carlo process, and considering different probabilities for each possible configuration. Our algorithm can be used in single-site and multi-site approximations. Finally, we applied our algorithm to a square lattice with on-site Coulomb interaction. The effect of on-site Coulomb interaction on the density of states and band structure of a two-dimensional square lattice has been investigated in both single-site approximation (DMFT) and multi-sites beyond effective medium supercell approximation (BEMSCA). Our results show that the density of states obtained from DMFT and BEMSCA are similar in the weak interactions, but at high repulsions strengths, they are different due to inter-site correlations contributions that are considered in the multi-site approximation where the density of states splits faster in the multi-site approximation. Another important result is observing partial flat valance bands at Fermi energy for special repulsion strengths. Our calculated bands show that a metallic to an insulator phase transition occurs at u c =9.05t and u c =2.4t for N c =1 and N c =4 beyond effective medium supercell approximation respectively. The difference between u c is due to the consideration of inter-site correlation in the multi-site approximation, which is ignored in the single-site approximation.