ISSN 0021-3640, JETP Letters, 2011, Vol. 94, No. 10, pp. 768–773. © Pleiades Publishing, Inc., 2011. Original Russian Text © I.S. Krivenko, A.N. Rubtsov, 2011, published in Pis’ma v Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki, 2011, Vol. 94, No. 10, pp. 832–837. Analysis of the Nature of the Peak Structure of Hubbard Subbands Using the Quantum Monte Carlo Method I. S. Krivenkoa, b and A. N. Rubtsova a b Faculty of Physics, Moscow State University, Moscow, 119992 Russia Centre de Physique Theorique, Ecole Polytechnique, 91128 Palaiseau Cedex, France Received October 18, 2011 It has been shown that the description of the details of the electronic spectra obtained by the combination of the dynamical mean field theory, quantum Monte Carlo methods, and the maximum entropy method can be significantly improved by changing the last method to optimal regularization of the analytic continuation of the Green’s function to the real frequency axis. Starting with the quantum Monte Carlo data, this method has reconstructed peaks in the structure of Hubbard subbands with a maximum error of 0.001 to 0.01. Owing to the universality of the quantum Monte Carlo method, by varying hybridization, it is possible to determine the features of the hybridized function that are responsible for the formation of the structure of Hubbard subbands. It has been shown that there is no direct relation between the peak structure of subbands and the central Kondo peak. The result indicates the charge nature of resonances responsible for the formation of the peak structure. DOI: 10.1134/S0021364011220073 A significant advance recently achieved in the physics of structures and materials with strong electron correlations is primarily due to using numerical methods and combined numerical and analytical approaches. First of all, this concerns the combination of the dynamic mean field theory (DMFT) [1] with the family of quantum Monte Carlo (QMC) methods [2], which makes it possible to quantitatively calculate the electronic properties of strongly correlated materials, in particular, using realistic multiorbital models [3]. In this scheme, stochastic QMC algorithms are used to solve the impurity problem that describes an atom hybridized with Gaussian environment. The hybridization function is calculated self-consistently. The main advantages of the QMC methods in this situation are universality (the solution can be determined for an arbitrary hybridization function in a fairly wide temperature range) and acceptable resource consumption. At the same time, modern QMC algorithms have a significant disadvantage because calculations are performed for imaginary time; i.e., they give the electron Green’s function at Matsubara frequencies. The closed self-consistency scheme for the DMFT can be formulated at Matsubara frequencies. However, the determination of the physical electronic spectrum requires a certain analytic continuation procedure, which is an ill-conditioned mathematical problem. This circumstance, together with a rather high noise level of the results characteristic of stochastic algorithms, leads to serious problems in analyzing electronic spectra. The base model of the DMFT is the single-band Hubbard model with half-filling with the Hamiltonian H = U ∑ j † t 1 1 n j↑ – - n j↓ – - – ------c jσ c j'σ . 2 2 N 〈 jj'〉 , σ ∑ (1) Here, the first term describes the Coulomb repulsion of electrons with opposite spins at a site, the second term describes hops of electrons between neighboring atoms, and N is the number of neighbors of an atom in the lattice. The factor of the second term is chosen such that the width of the bare conduction band remains unchanged with a change in the dimension of the system (and, correspondingly, the number of neighbors). In the DMFT, the problem of determining the electronic properties of lattice model (1) is reduced to the single-site impurity Anderson model with the action S = S at + ∑Δ † ω c ωσ c ωσ , (2) ω, σ where Sat is the action of a single isolated atom with the Coulomb repulsion U and hybridization Δ should be determined self-consistently. In the limit of the infinite number of neighbors on the Bethe lattice [1], correlations are local in space and the transition from lattice problem (1) to impurity problem (2) is exact. For 768 ANALYSIS OF THE NATURE OF THE PEAK STRUCTURE 769 this reason, the DMFT is free of approximations in the limit N ∞. Hybridization is given by the expression 2 Δω = t gω , (3) where gω is the Green’s function of impurity problem (2). The solution is sought iteratively; i.e., the Green’s function of action (2) at each iteration step is calculated with hybridization obtained at the preceding step. This scheme ensures good convergence throughout the entire parameter region. The efficiency and accuracy of the calculation are determined by a program (solver) for the numerical determination of the Green’s function of problem (2). The main features of the solution of the single∞ are well band Hubbard model in the limit N known. With an increase in the Coulomb interaction constant U at fairly low temperatures, the electronic spectrum of the system exhibits two wide Hubbard subbands near the energies of ±U/2. At U = (4.7–6.0)t, the system undergoes a first-order phase transition and a gap appears in the spectrum, which corresponds to the formation of the Mott insulator phase. However, the spectral density at the Fermi level in the entire Fermi-liquid phase up to the transition point follows from the Luttinger theorem and locality of the selfenergy function. This corresponds to the narrow Kondo peak in the spectrum near the Fermi level, whose width vanishes when approaching the transition point. Figure 1 shows the typical spectrum near the transition point. It is noteworthy that the described three-peak structure of the density of state is characteristic of the infinite-dimensional limit. For example, the central peak for two-dimensional systems is completely suppressed due to the spatial nonlocality of correlations. Despite the simplicity, fundamentality, and fairly long history of investigations, the description of the ∞ still single-band Hubbard model in the limit N includes open problems regarding the so-called peak structure of Hubbard subbands. Hubbard subbands contain narrow peaks at the edge of the emerging band gap. Peaks are present only in the metal phase at quite large U values and disappear completely at the transition to the Mott insulator phase. The existence of this structure was apparently mentioned in [4], and its presence was reliably shown using the D-DMRG [5], which is an almost exact method that allows the determination of the spectrum at zero temperature immediately on the real frequency axis (the spectrum obtained in that work is also shown in Fig. 1). The physical nature of these peaks is still unclear. The authors of [5] attributed the fine structure of Hubbard subbands to the interaction of quasiparticles of the central Kondo peak with certain collective excitations. This is indirectly confirmed by the fact that the spectral weights of the central peak and fine-structure peaks decrease identically when approaching the transition point [6]. However, any strict justification of the JETP LETTERS Vol. 94 No. 10 2011 Fig. 1. Densities of states of the single-band Hubbard model with half-filling on the infinite-dimensional Bethe lattice at U = 2.4 and t = 0.5. The solid line is the result obtained by the combination of the DMFT, CTQMC, and optimal regularization (temperature T = 0.01). The dashed line is the D-DMRG result at zero temperature obtained in [5]. The thin horizontal straight line at a value of 2/π indicates the a(0) value following from the Luttinger theorem. nature of the peak structure of subbands has not yet been presented. It is even unclear whether this structure is attributed to the spin or charge degrees of freedom of the multielectron system. The description of the fine structure of Hubbard subbands at nonzero temperature using QMC solvers faces difficulties that are mentioned at the beginning of the paper and associated with the ill-conditioning of the analytic continuation from Matsubara frequencies to real frequencies. The typical random error of QMC data at Matsubara frequencies is about 0.003 or worse. Improvement in the accuracy requires a significant increase in computational resources and imposes increased requirements on the software implementation of the QMC method. In this case, the commonly accepted procedure involving the maximum entropy method (MAXENT) [7] does not reproduce the discussed structure; i.e., the reconstructed spectrum includes only smoothed subbands. It is noteworthy that the described problems increase significantly beyond the Hubbard model and for multiorbital models of correlated materials; as a result, the details of the electronic spectrum far from the Fermi level can hardly be described in the DMFT + QMC + MAXENT approach. In this work, we demonstrate that the description of the details of the electronic spectrum can be significantly improved by changing the MAXENT method to the regularization optimal for minimizing the error of the result. Starting with the QMC data, the use of the proposed optimal regularization allows the reconstruction of peaks in the structure of Hubbard subbands with a maximum error of 0.001 to 0.01. Using 770 KRIVENKO, RUBTSOV the QMC method, which is a universal method reliably applicable for arbitrary hybridization, and varying hybridization, it is possible to determine the features of the hybridization function that are responsible for the formation of the structure of Hubbard subbands and, correspondingly, to clarify the physical nature of this structure. It is shown that any direct relation between the peak structure of subbands and the central Kondo peak is surprisingly absent. The result indicates the charge nature of resonances responsible for the formation of the peak structure. The analytic continuation problem is reduced to solving the linear integral equation +∞ a(⑀) - d⑀ = g ( ω ) ∫ ----------iω – ⑀ (4) –∞ for the spectral density a(⑀), where the variable ⑀ is on the real frequency axis. After the decomposition of the functions g(ω) and a(⑀) into a pair of orthogonal bases Ᏺω and Ᏺ⑀, respectively, the integral equation is transformed to the linear system Ma = g Ᏺ [ a;R ] = Ma – g 2 + ( a, Ra ) = min. (6) Here, R is a regularizing Hermitian matrix. The minimum is reached at ∫ We propose the following optimality condition. The vector g is known with error; i.e., g = g + δg, where g is the mathematical expectation and δg is a random variable that has zero mean value and is charˆ = δgδg † (here acterized by the covariation matrix K g and below, the overbar means the mathematical expectation to which the QMC calculation converges). Let a be the exact solution of system (5) with the exactly known right-hand side, i.e., M a = g . Averaging the standard deviation of a from a over all possible values of the random vector δg, we obtain a–a 2 † = Tr { XAX – 2XB } + Tr { xx }, † (5) with a known matrix M. This system formally consists of an infinite number of equations. However, it can be expectedly reduced to a finite system with the appropriate choice of the bases. The specific choice of bases is beyond the scope of this brief paper. Note that a particular form of the reasonably chosen bases affects only the number of equations that should be retained in system (5), but it does not affect the final result of the calculation. Ill-conditioning means that small variations of g lead to (exponentially) large changes in the solution. The starting point of the proposed method is the Tikhonov regularization of the ill-conditioned problem [8]. We seek the vector a that minimizes the Tikhonov functional 2 2 d a- d⑀. In contrast to this ----2 d⑀ approach, we calculate the regularizing functional, requiring its optimality (i.e., minimum total error) on a certain class of solutions and for a given form of the error. tional (a, Ra) ∝ † † † A ≡ M Maa M M + M K g M, † † B ≡ M Maa . † –1 a = XM g, where X ≡ ( M M + R ) . (7) If the matrix R is nonzero, the solution thus obtained for problem (6) contains a systematic error. However, the correctly chosen regularizing matrix makes it possible to significantly reduce the random error of the result, eliminating exponential instability with respect to the input data. As a result, the total error is much smaller than that in the solution of system (5). Matrix R is usually chosen under certain a priori assumptions on the solution properties. For example, the assumption of the smoothness of the solution would correspond to the popular regularizing func- (9) (10) The optimal matrix R should ensure the minimum of 2 x – x . This condition itself is not constructive, because the exact solution a is unknown in each particular case. However, we can seek the matrix R optimal for a given set (class) of possible solutions. Marking the averaging over this set by triangular brackets, we obtain the following condition for R: 2 〈 a – a 〉 = Tr ( X 〈 A〉 X – 2X 〈 B〉 ) † (11) + Tr 〈 aa 〉 = min . R Solving this variational problem and taking into account an additional condition δR = δR†, we arrive at the following equation for matrix X: † † (8) 〈 A〉 X + X 〈 A〉 = 〈 B〉 + 〈 B 〉 . (12) Physical reasons should be used to choose the class of possible solutions for averaging. The use of this a priori information additional to initial problem (5) is obligatory for any regularization schemes. However, note that detailed information on the class of solutions is not required. All necessary information on possible solutions is contained in the first and second moments of the a priori distribution of a, because only these quantities enter into Eq. (11) (taking into account Eqs. (9) and (10) for A and B). These moments can be calculated for a certain model class that reproduces the main properties of possible solutions. JETP LETTERS Vol. 94 No. 10 2011 ANALYSIS OF THE NATURE OF THE PEAK STRUCTURE As this model class, we take the following superposition of several Lorentzian peaks with the same width γ: J a(⑀) = Zj -, ∑ ------------------⑀ – Ω + iγ j=1 (13) j where the positions of the peaks Ωj and residues Zj at them are random variables with known statistical distributions. We assume that all residues Zj are distributed identically and independently of each other; 2 therefore, only two model parameters, 〈 Z 〉 and 〈 Z〉 , are necessary to completely describe Zj in the correlation function 〈 a ( ⑀ 1 )a* ( ⑀ 2 )〉 . The relation 〈 Z〉 = 1/J immediately follows from the sum rule for the spectral density Jj = 1 Z j = 1. Finally, we assume that the positions of poles Ωj are distributed independently according to a certain model distribution. One of the simple possibilities is the following Lorentzian distribution with zero mathematical expectation: ∑ 1 1 -, P ( Ω j ) = --------- -------------------------πΩ M 1 + ( Ω j /Ω M ) 2 (14) where ΩM is the estimated FWHM of the spectrum. We do not assume that the desired spectral density is the superposition of Lorentzian peaks with a given FWHM. Model spectrum (13) is used only to calculate the first and second moments of the a priori distribution of a. In the chosen basis Ᏺ⑀, these moments are a column vector and a square matrix, respectively, whose components were calculated through numerical integration. Then, matrix X is determined from Eqs. (12). To determine this matrix, it is simplest to perform a unitary transformation to the eigenbasis of the matrix 〈 A〉 in which the solution of Eq. (12) can be represented in the explicit form B 'ij + ( B *ji )' X 'ij = --------------------, λi + λj 1 ≤ i, j ≤ N, (15) where primes mark the matrix elements in the eigenbasis of the matrix 〈 A〉 and λi are the eigenvalues of the matrix 〈 A〉 . In this approach, it is sufficient to diagonalize the N × N matrix (∝N3 operations) and then to use Eq. (15) (∝N2 operations); this is the advantage of this approach compared to the direct solution of the system of N(N + 1)/2 equations (12) (∝N6 operations). It can be proved that resulting formula (15) certainly yields finite expressions for all matrix elements X 'ij , because the denominator in it is always positive. Nevertheless, the software implementation of the method faces difficulties, because the numerical solution is strongly sensitive to round-off errors appearing in the diagonalization process. In order to suppress JETP LETTERS Vol. 94 No. 10 2011 771 round-off errors, we used the MPFR program library, implementing arithmetic with numbers of arbitrary precision. A typical calculation involved operations with matrices of dimension N ≈ 40; elements of these matrices were calculated with 100 decimal digits. A further increase in the accuracy or in the number of equations did not affect the result. Using the obtained matrix X, the spectral density can be easily reconstructed by Eq. (7) from the Green’s function. In the case of processing a set of similar input data, one calculation of matrix R is sufficient in view of the linearity of the problem. An important remark should be made. The existing algorithms for the reconstruction of the spectral density are significantly nonlinear, because they include the requirement on the positivity of the spectral density and the satisfaction of the sum rule +∞ a (⑀)d⑀ = 1. These con- ∫ –∞ ditions (as well as any other exact relations) can be satisfied in our scheme if the regularizing functional R is determined from the known matrix X and conditional minimum (6) is sought on the set of accessible solutions instead of calculations by Eq. (7). Such additional constraints were not introduced in our calculations, the results of which will be presented below. In particular, due to this fact, the spectral function has nonphysical negative values that are small in absolute value in some regions. However, it is clear that, if a conditional minimum is sought on the set of positive solutions, the spectral weight in these regions will be zero and the remaining part of the plot changes insignificantly. At the same time, the small deviation to the regions of negative values indicates the general applicability of the scheme. An important advantage of the linear algorithm is also the possibility of reconstructing unnormalized and sign-alternating quantities such as off-diagonal components of the spectral density of states in multiorbital problems and the self-energy function. Below, we present the results of calculating the density of states in the Hubbard model on the infinitedimensional lattice. To solve the DMFT equations, we used the TRIQS (Toolbox for Research on Interacting Quantum Systems) [9] involving ALPS (Algorithms and Libraries for Physics Simulations) [10, 11]. The impurity problem was solved by the hybridization expansion method [12] with collection of the statistics in an orthogonal polynomial representation [13]. To reconstruct the spectral function, the model distribution with the parameters γ = 0.1 and ΩM = 1.7 was used. For simplicity, the covariation matrix is assumed to be diagonal with matrix elements of about 10–5. To test the correctness of the determination of this value (which includes stochastic noise of the QMC method and, possibly, other errors inherent in a particular numerical implementation), we compared the input data and “smoothed” Green’s function found by substituting the calculated spectral function into Eq. (4). 772 KRIVENKO, RUBTSOV a smoothed section (Fig. 2, long-dashed line). Using this modified spectral function, a new hybridized ˜ was calculated by Eqs. (4) and (3). For function Δ ˜ , the Green’s function g̃ impurity problem (2) with Δ Fig. 2. Densities of states of the single-band Hubbard model with half-filling on the infinite-dimensional Bethe lattice at U = 2.4 and t = 0.5. The short-dashed line is the result obtained by the combination of the DMFT, CTQMC, and optimal regularization (temperature T = 0.01). The long-dashed line is the density of state a(⑀) after the removal of the central peak. The solid line is the density of state ã (⑀). The maximum and typical deviations of the smoothed Green’s function from the input data are about 0.01 and 0.001, respectively. This confirms the correctness of the choice of the covariation matrix. We recall that the diagonal matrix elements of this matrix are square in the error of the input data. It was also verified that the choice of the parameters of the model distribution, as well as the covariation matrix, insignificantly affects the result. It was shown that the described optimal regularization method makes it possible to reconstruct the peak structure of Hubbard subbands. The typical reconstructed spectral function containing this structure is shown in Fig. 1. When U exceeds the critical value, a gap appears in the spectrum and the peak structure disappears. Note that the spectral function obtained in our work differs from the D-DMRG calculation: the spectral peaks are spread and the spectral weight near the pseudogap in our spectrum is much higher than that for D-DMRG. These differences are due not only to the regularization-induced smoothing of dependences, but also to the physical effects of nonzero temperature. Due to the possibility of studying the peak structure of Hubbard subbands in the results of the QMC calculations, the relation between the central Kondo peak and peaks in Hubbard subbands can be examined simply. For U = 4.8t (i.e., for the system with the pronounced peak structure presented in Fig. 1), we performed an additional QMC calculation for an artificially generated spectral function without the central peak. After reaching self-consistency, we plot the spectral function shown in Figs. 1 and 2 (short-dashed line). The central peak in this function was replaced by was calculated using the QMC method, and the modified density of state ã , which is shown by the solid line in Fig. 2, was reconstructed from this Green’s function. Note that the procedure of the removal of the central peak is generally accompanied by the loss of a small part of the spectral weight. The spectral weight ˜ by the correspondcan be recovered by multiplying Δ ing normalizing factor. However, as was verified, this normalization almost does not change the result for ã . It can be seen that the central peak is almost completely absent not only in the spectral function corre˜ , but also in the final sponding to the hybridization Δ spectrum ã . Thus, spin-flip effects (Kondo processes), responsible for the formation of the central quasiparticle peak, are completely removed from the system. Nevertheless, the peak structure of subbands holds, but the frequencies of the peaks are lower and, for this reason, the peaks are more pronounced. We conclude that the peak structure of subbands is not directly related to the quasiparticle peak at the Fermi level, contrary to earlier assumptions. Nevertheless, the nonzero density of state near the Fermi level is important for the formation of the peak structure of subbands, because this structure disappears when the gap is formed. The total number of states in a fairly wide region near the Fermi level (i.e., in a region about the width of the pseudogap), rather than the quasiparticle peak itself, is likely important. For the studied case of nonzero temperature, the spectral weight in the pseudogap region is significantly nonzero. For this reason, the removal of the central peak leads only to a comparatively small decrease in the number of these states and, correspondingly, to a decrease in the frequencies of the peak structure. In real materials, this situation is characteristic of plasma excitations whose frequencies are almost independent of Kondo processes. Thus, our result indicates that charge, rather than spin, excitations are responsible for the formation of the peak structure of Hubbard subbands. We are grateful to A.I. Likhtenshtein, M.I. Katsnel’son, S. Birman, and J. Mravlje for stimulating discussions. This work was supported by the Council of the President of the Russian Federation for Support of Young Scientists and Leading Scientific Schools (project no. NSh-4895.2010.2), the Dynasty Foundation, the Russian Foundation for Basic Research (project no. 11-02-01443-a), l’Agence Nationale de la Recherche (France, project CorrelMat), and Institut du Développement et des Ressources en Informatique Scientifique/Grand Equipement National de Calcul Intensif (France, project no. 20111393). JETP LETTERS Vol. 94 No. 10 2011 ANALYSIS OF THE NATURE OF THE PEAK STRUCTURE REFERENCES 1. A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996). 2. E. Gull, A. J. Millis, A. I. Lichtenstein, et al., Rev. Mod. Phys. 83, 349 (2011). 3. G. Kotliar, S. Y. Savrasov, K. Haule, et al., Rev. Mod. Phys. 78, 865 (2006). 4. X. Y. Zhang, M. J. Rozenberg, and G. Kotliar, Phys. Rev. Lett. 70, 1666 (1993). 5. M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 72, 113110 (2005). 6. M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 77, 075116 (2008). 7. M. Jarrell and J. E. Gubernatis, Phys. Rep. 269, 133 (1996). JETP LETTERS Vol. 94 No. 10 2011 773 8. A. N. Tikhonov and V. Ya. Arsenin, Solutions of IllPosed Problems, 3rd ed. (Nauka, Moscow, 1981; translation of the 1st ed.: Halsted, New York, 1977). 9. M. Ferrero and O. Parcollet, in preparation. 10. A. F. Albuquerque, F. Alet, P. Corboz, et al., J. Magn. Magn. Mater. 310, 1187 (2007). 11. B. Bauer, L. D. Carr, H. G. Evertz, et al., J. Stat. Mech. 2011, P05001 (2011). 12. P. Werner, A. Comanac, L. de’Medici, et al., Phys. Rev. Lett. 97, 076405 (2006). 13. L. Boehnke, H. Hafermann, M. Ferrero, et al., Phys. Rev. B 84, 075145 (2011). Translated by R. Tyapaev