{'fraction_non_alphanumeric': 0.03366037150502191, 'fraction_numerical': 0.01415441889324078, 'mean_word_length': 4.380346338410182, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient treatment of two-particle vertices in dynamical mean-field theory

We present an efficient and numerically stable algorithm for calculation of two-particle response functions within the dynamical mean-field theory. The technique is based on inferring the high frequency asymptotic behavior of the irreducible vertex function from the local dynamical susceptibility. The algorithm is tested on several examples. In all cases rapid convergence of the vertex function towards its asymptotic form is observed.


I. INTRODUCTION
Electronic correlations in materials have been one of the central topics of the condensed matter physics throughout its history encompassing topics such as high-temperature superconductivity, colossal magnetoresistance or heavy fermion physics. The introduction of dynamical mean-field theory (DMFT) in early 1990's 1-3 marked a big step forward in the theory of correlated electrons in providing an approximate, but non-perturbative, computational method with several exact limits. Moreover, numerical DMFT calculations are feasible also for multi-band Hamiltonians necessary for the description of real materials.
In the past twenty years DMFT was applied, at first, to models 3 and, later, to real materials 4,5 . Naturally, most of the studies focused on single-particle quantities, featured explicitly in the DMFT equations, which can be compared to photoemission spectra and which can provide information about the metal-insulator transitions. Also local two-particle correlation functions, static as well as dynamical, can be computed with little additional effort in most DMFT implementations, providing information about the local response to applied fields. Computation of the non-local response is a more difficult venture, however, with great potential gains. With the dynamical susceptibilities available it is possible to compare to the experimental data from inelastic neutron, x-ray or electron loss spectroscopies. Perhaps even more interesting possibility opens with the static susceptibilities. Monitoring their divergencies as a function of temperature and the reciprocal lattice vector allows investigation of the second order phase transitions and an unbiased determination of the order parameters. So far only calculations for simple models, such as the single-band Hubbard model 6,7 ,the periodic Anderson model [8][9][10] , or Holstein model 11 , have been reported. Similar calculations were also performed with cluster 12 or diagrammatic 13 extentions of DMFT.
In this paper, we present a scheme for computation of the static two-particle response functions in the multiband Hubbard model within DMFT. The key development is splitting the particle-hole irreducible vertex into low frequency (LF) and high frequency (HF) parts, and expressing the HF part in terms of the local dynamical susceptibilities. This allows reformulation of the Bethe-Salpeter equation in terms of the LF quantities plus corrections, thus reducing the numerical cost of the calculations and improving their stability. The procedure is applied to the single-band Hubbard model at and away from the half filling and a two-band bilayer model. In all the studied cases we find a rapid convergence of the vertex function towards its asymptotic form, which leaves the size of the LF problem rather small and manageable also for multi-band systems.
The paper is organized as follows. After a general introduction, the two-particle formalism is reviewed in section II, followed by the discussion of the numerical implementation. In section III applications to simple model systems are reported. Discussion of the asymptotic behavior of the irreducible vertex and the blocked form of the irreducible Bethe-Saltpeter equation is left to Appendices A and B.

II. THEORY
Our starting point is the multi-band Hubbard Hamiltonian where t ij R−R ′ is the hopping amplitude between orbital j on site R ′ and orbital i on site R, c † Ri (c R ′ j ) are the corresponding creation (annihilation) operators, and U ij,kl is the anti-symmetrized local Coulomb interaction. Throughout the text we do not distinguish between spin and orbital degrees of freedom. Summation over repeated orbital indices is assumed.
The developments and calculations reported here fall within the framework of dynamical mean-field theory. Detailed discussion of the formalism and basic applications of DMFT can be found in Ref. 3. The main feature of DMFT is that the irreducible vertex functions (see below for details) are build only from local propagators and thus can be obtained from an effective quantum impurity problem. Evaluation of the two-particle response functions can be viewed as a postprocessing of the solution to the DMFT equations, which self-consistently determine the fermionic bath for the effective impurity problem.

A. Linear response formalism
We review briefly the DMFT linear response formalism. For details the reader is referred to Ref. 3. We are interested in the response of a system controlled by (1) to static fields which couple to the spin, charge or a more general orbital density described by the susceptibility Here q is a vector from the first Brillouin zone and the site indices were dropped in the local averages. In the DMFT approximation the susceptibility (2) can be obtained from whereχ ij,kl (q; ω m , ω n ) is the solution of coupled integral equations 3,6 The local and corresponding 'lattice' quantities are distinguished by the presence of parameter q. The quantities is the equations are functions of the discrete fermionic Matsubara frequencies, which are at temperature T given by ω n = (2n − 1)πT (n is an integer). The equations have the form of the irreducible Bethe-Salpeter (IBS) equation, with Γ mn,pq (ω 3 , ω 4 ) being the particle-hole irreducible vertex at zero energy transfer, in the diagrammatic expansion of which the mn pair of external vertices cannot be separated from pq pair by cutting two single-particle propagators. In the DMFT approximation only the local diagrams contribute to Γ mn,pq (ω 3 , ω 4 ) 14 , which allows simplifying the IBS equation for the full two-particle correlation function to equation (4) in terms of the reduced quantityχ ij,kl (q; ω 1 , ω 2 ) and k-integrated particle-hole bubblẽ G ik (k; ω 1 ) is the single-particle propagator obtained from the Dyson equation where µ is the chemical potential, h k is the Fourier transform of the single-particle part of the Hamiltonian (1) and Σ(ω n ) is the local self-energy. Equation (5) is the IBS equation for the impurity, relating the local two-particle correlation functioñ to the local particle-hole bubblẽ

Asymptotic behavior of the vertex function
Separation of low and high frequency parts of a given problem is common in theoretical physics. In the DMFT practice it is widely used in calculation of the self-energy, where the numerical solution of the Dyson equation at low-frequencies is augmented by a high-frequency asymptotic expansion obtained from the moments of the spectral function 16 . In Appendix A we start by considering the HF asymptotic expansion for the self-energy, which serves a precursor for the HF expansion of the vertex function Γ ij,kl (ω 1 , ω 2 ).
Unlike the self-energy, where a separate equation exists for each frequency and so the LF and HF parts are decoupled, in case of equations (4) and (5) we are dealing with matrices in the Matsubara frequencies and have to solve coupled equations for the HF and LF blocks. The key ingredient of our procedure is replacing the vertex function outside the LF block by its asymptotic form where χ ph (ν) and χ pp (ν) are the local dynamical susceptibilities in the particle-hole and particle-particle channel, respectively, as functions of bosonic Matsubara frequency ν, and ω c is the cut-off frequency defining the LF block.
Derivation of expression (10)is given in Appendix A. Viewed as a function of variables ω 1 and ω 2 , Γ ∞ ij,kl is a constant, U ij,kl , plus two ridges along the main and the minor diagonal, a structure observed in earlier numerical studies 17 . The cross section of these ridges is given by the local dynamical susceptibility, typically, with a fast decay away from ν = 0. The asymptotic form (10) is used in the numerical treatment of both equations (4,5) as described below.
In this representation the LF block is located in the upper left corner of the matrix and the ω dependent part of Γ ∞ ij,kl (ω 1 , ω 2 ) (10) proportional to χ ph and χ pp has a band diagonal form.
In the following we split the set of Matsubara frequencies into the LF block |ω n | < ω c , denoted with a block index 0, and HF block |ω n | ≥ ω c , denoted with a block index 1. Expansion of a matrix equation of the form (4,5) into the LF and HF blocks is given in Appendix B. The impurity IBS equation (5) is used to compute the LF part of the vertex function, the Γ 00 block. We assume to know the LF (00) part ofχ ij,kl (ω 1 , ω 2 ) (8) from the solution of the impurity problem and to have a complete information aboutχ 0 ij,kl (ω 1 , ω 2 ) (9). In practice we use a numerical representation ofχ 0 ij,kl (ω 1 , ω 2 ) up to a HF cut-off Ω c (typically much larger that the LF cut-off ω c ) and the analytic − 1 ω 2 n tail up to infinity for elements thẽ χ 0 ij,ij elements, while the remaining elements involving off-diagonal G ij (9) are neglected above Ω c . The LF part of the vertex Γ 00 is obtained from (B4), which amounts to inverting equation (5) restricted to the LF block and subtracting a correction. For χ ph (ν) and χ pp (ν) sharply peaked around ν = 0 the correction reduces to a constant shift everywhere except the vicinity of ω c .
Once we have computed Γ 00 we proceed to solve equa-  tion (4) forχ ij,kl (q; ω 1 , ω 2 ) assuming the full knowledge ofχ 0 ij,kl (q; ω 1 , ω 2 ) and Γ ij,kl (ω 1 , ω 2 ). We use the same partitioning into the LF and HF blocks and take advantage of the fact that only ω-summed susceptibility χ ij,kl (q) is of interest. This allows partial ω-summations in the HF block to be preformed at an earlier stage of the calculation and to replace matrix-matrix operations with matrix-vector ones, see Appendix B. Further speed-up comes from treating the ω-dependent (B11) and the constant (B12) parts of Γ 11 in two subsequent steps, which reduces the most demanding part of the HF problem to a repeated band-matrix-vector multiplication.

A. 2D Hubbard model
In the following we present our results for 2D Hubbard model on a square lattice with the nearest neighbor hopping. In particular, we compute the Néel temperature as a function of U and doping. This is an academic example since it is well know that no magnetic order is possible in this model at finite temperature and the magnetic oder studied here is an artefact of the DMFT approximation. Nevertheless, we find it a useful benchmark since i) a comparison to other studied is possible and ii) the model exhibits a crossover from a nesting-driven instability at small U to the Heisenberg model with local moments at large U .
First, we report the results obtained at half filling. After converging the DMFT equations, a finer QMC run of 5×10 6 to 10 7 sweeps is performed to obtain the local twoparticle correlation functionχ ij,kl (ω 1 , ω 2 ) and local susceptibilities χ ph ij,kl (ν), χ pp ij,kl (ν). It is possible to separate the charge and the spin-longitudinal channels from the beginning, but we perform the calculation in the upup, updn, dnup, dndn basis, as would be the case for general orbital indices, and separate the two channels only in the final result. Expression (6) is evaluated on a 20 × 20 grid of q-points (66 irreducible points) using a fine k-mesh of 551 × 551 points, which is more than enough to ensure that the uncertainties associated with k-points sampling as much smaller that other sources of error.
The IBS equations were solved for several LF cut-offs ω c leaving between 5 to 20 lowest Matsubara frequencies for temperatures above 1/20, and 15 to 30 Matsubara frequencies for lower temperatures. In all studies cases, the relative variation of the susceptibility with the LF cut-off was less than 1%. We have also tested the sensitivity to the HF cut-off Ω c . When neglecting the contributions from above Ω c completely, a slow convergence of the susceptibilities close to the AFM instability was observed. The error essentially followed the 1/Ω c dependence which may be expected when truncating a sum over 1/ω 2 n series. The dominant contribution to the error came from the Γ 01 X 11 Γ 10 term in equation (B6). Introduction of an analytic summation of the 1/ω 2 n tails to infinity lead to the converged results already for lowest HF cut-offs corresponding to 80-100 lowest Matsubara frequencies.
To check the consistency of our calculations we have solved equation (4) also for the local bubbleχ 0 ij,kl (ω 1 , ω 2 ), i.e. went fromχ ij,kl (ω 1 , ω 2 ) to Γ ij,kl (ω 1 , ω 2 ) and back, and compared the result with the local susceptibility obtained directly from the QMC simulation. For all reported parameters we have found a relative deviations of less than 2% close to the AFM instability and less than 0.5% deeper in the paramagnetic phase. As a side remark we mention that in the insulating phase for U > 1 it is important thatχ 0 ij,kl (ω 1 , ω 2 ) entering equations (4) and ((5) are exactly the same (e.g. obtained with the same k-point sampling) since in this parameter range the results are extremely sensitive to differences betweenχ 0 's entering equation (4) and (5). Due to the SU(2) symmetry of the model the transverse and the longitudinal spin susceptibilities are identical in principle. This symmetry is not explicitly enforced in the QMC calculation, where the identity is obeyed only approximately within the stochastic error of the simulation. For U < 1 we have found relative deviations between the transverse and longitudinal spin susceptibilities of 0.2-0.5%. For larger U 's (especially at lower temperatures) the transverse susceptibility rapidly becomes useless due to increasingly large errorbars. This problem is specific to the Hirsch-Fye QMC solver and the way the transverse susceptibility is measured in the simulation. We have used the longitudinal susceptibility to determine T N in this parameter range.
In Fig. 1 we show the Néel temperature as a function of U. On the large U side T N approaches the static mean field value of the corresponding Heisenberg model reflecting the insulator with well defined local moments. On the small U side the magnetic order arises from the Fermi surface nesting. For comparison, we show also T N for the RPA condition for the instability, χ 0 (T N ) = 1/U , obtained by replacing the vertex with a constant and taking χ 0 (T ) of the non-interacting electron gas or χ 0 (T ) obtained with the dressed propagators.
In Fig. 2 we show the LF block of the real part of the transverse spin vertex Γ 12,12 for large and small U as a function of the fermionic Matsubara frequencies, which exhibits the two-ridge structure and a rapid convergence towards its asymptotic form. At small U the ridges are only weakly T -dependent and their absolute value is small compared to the constant part of the vertex. At large U s the ridges grow as 1/T and thus become dominant at low enough temperatures. In Fig. 3 we show the same vertices at fixed ω m (still in the LF block) together with the corresponding asymptotic tails Γ ∞ 12,12 . We find a good match between the two dependencies. The remaining mismatch is mainly due to the time discretization error inherent to the Hirsch-Fye impurity solver, which shows up in a spurious curvature of the Γ(ω m , ω n ) for larger ωs (not really visible in the present plots). Using the semi-analytic asymptotic tail Γ ∞ can efficiently remove this spurious curvature. The discretization error arises from performing Greens function convolutions on a discrete time grid and should not be confused with the so called Trotter error, which is another consequence of discretizing the imaginary time and which cannot be removed be restricting oneself to the low frequencies. Impurity solvers using continuous time approaches are not expected to have this problem.
Finally, the spin susceptibility throughout the Brillouin zone is shown in Fig. 4. The plots clearly show that the magnetic instability takes place at (π, π) as expected. The behavior of the vertices reveals that the origin of the instability is quite different for small and large U s. While for large U s the main T -dependence in the problem is in the vertex, related to 1/T behavior of the local susceptibility, for small U x the main T -dependence is in the 'bubble', namely the of thermal smearing of the nesting related peak.

B. Effect of doping
We use doping to break the particle-hole symmetry in the above model. The calculations were performed in the crossover regime between the metallic and insulating phases at U = 1. We have found that magnetic instability survives doping up to about 0.2 electron (or hole) per site. For dopings less or equal to 0.14 the instability appears at [π, π]. For dopings of 0.16 and 0.18 we have noticed that the maximum of χ s (q) moved slightly away from [π, π] as the temperature was lowered. Similar behavior was observed also in a model with the nextnearest-neighbor (diagonal) hopping 18 . The Néel temperature as a function of doping is shown in the inset of Fig. 1. In Fig. 5 we show the detail of the spin susceptibility in the vicinity of [π, π] for the highest studied doping of 0.18. While at high temperatures the susceptibility exhibits a flat maximum at [π, π], at lower temperatures the maximum moves to [1.04π, π] eventually giving rise to an incommensurate order.

C. Bilayer model
As a last example we study a bilayer Hubbard model 19 . We choose this model because: i) it has two bands and the lattice bubbleχ ij,kl (q; ω 1 , ω 2 ) has a more complicated orbital structure than for the single-band model ii) it exhibits a non-trivial temperature dependence of the uniform spin susceptibility characterized by an exponential decrease at low temperatures. In our previous work on a more general version of this model (model of covalent insulator) 20 we have used a different numerical approach to compute the uniform susceptibility without an explicit use of the vertex function and its HF tail, an approach which could not be used for more general models and q-dependent susceptibilities. Obtaining the exponential decay at low temperatures proved then numerically challenging. Therefore we find this model to provide a good test of the performance of the present method. The model is obtained by taking two Hubbard sheets from the previous section (with opposite signs of the in-plane hopping) and adding an inter-plane hopping of 0.2. The different signs of the in-plane hoppings on the two sheets ensures that any non-zero inter-plane hopping opens a hybridization band gap. The inter-plane hopping was chosen such that the interesting temperature range can be easily covered with the QMC simulations. As in the other cases studied so far we were not primarily interested in the physics of the model, but in the performance of the computational method. Therefore we omit that fact that the system actually exhibits an AFM instability and focus on the uniform susceptibility only. While the AFM instability can be removed by introducing an in-plane frustration, the T -dependence of the uniform susceptibility is a generic feature independent of the bandstructure details 20 .
In Fig. 6 we show the local and uniform spin susceptibilities together with the uniform charge susceptibility, which is inversely proportional to the compressibility of the electron liquid. The calculated spin susceptibilities closely resemble those of Ref. 20. The decay of the uniform susceptibilities at low temperatures reflects an evolution towards a band-insulator-like state, characterized by a finite excitation gap for both charge and spin excitations (for more information see Ref. 19). In the present calculations we have not encountered any numerical problems in reproducing the low values of the susceptibility at low T .

IV. CONCLUSIONS
An efficient numerical procedure for calculation of the static two-particle response functions within the DMFT formalism was presented. The main development is the identification of the asymptotic form of the two-particle irreducible vertex function and implementation of a block algorithm for solution of the irreducible Bethe-Salpeter equations. We have studied several test cases including the single-band Hubbard model at and away from half-filling, and a bilayer Hubbard model at half-filling using Hirsch-Fye quantum Monte-Carlo impurity solver. In all cases we have observed a rapid convergence of the vertex function towards its asymptotic form and a very good numerical stability of the computations. This means that the numerical effort can be efficiently reduced by choosing a relatively small low-frequency cutoff, which makes calculations for multi-band systems feasible. Besides the susceptibility calculations the described procedure may be useful in diagrammatic extensions of DMFT, which employ the IBS formalism and two-particle vertices 17,[21][22][23] . δG . The self-energy itself can be obtained as a variation of the Baym-Kadanoff functional Σ = δΦ [G] δG . Diagrammatically this means that Σ is obtained from Φ by cutting a single Green function line and Γ = δ 2 Φ[G] δGδG is obtained by cutting two Green function lines.
We are interested in the high frequency behavior of Σ(ω) and Γ(ω 1 , ω 2 ). The expansion of Σ(ω) to the order 1 ω is well known from the moments of the spectral function 16 . Diagrammatically the zero's order term of Σ is the Hartree-Fock contribution, in which the external frequency does not enter any of the internal lines ( 7a). The first order contribution, 1/ω, comes from diagrams with external vertices connected by single Green function line (Fig. 7b). Their contribution has the form with χ(ν) being the dynamical susceptibility, where the plus sign in front of ω applies for particle-hole and the negative sign for particle-particle susceptibility. Expanding G in (A1) in the powers of 1 (ω−ν) = 1 ω + ν ω(ω−ν) and taking the leading 1 ω term we get the corresponding contribution to the self-energy. The remaining free summation over ν leads to an equal-time correlator of the c † c c † c type.
Our main interest here lies in the high frequency behavior of Γ(ω 1 , ω 2 ). In particular, we are interested in the zero's order contribution, which remains finite in the limit or large ω 1 and/or ω 2 . It is not possible to develop the high frequency expansion by keeping one of the frequencies fixed since the limit ω 1 , ω 2 → ∞ is not uniform, i.e. it depends on the way in which infinity is approached. By inspecting the diagrams for Γ(ω 1 , ω 2 ) we have identified two non-vanishing contributions. The first one, obtained from the variation of the zero's order term in the self-energy (Fig. 7a) is simply the bare vertex Fig. 8a. The other non-vanishing terms are obtained by cutting the Green function line in Fig. 7b, which leads to Fig. 8b and Fig. 8c. The corresponding expression for Γ ∞ reads The expression Γ ∞ reads Γ ∞ ij,kl (ω 1 , ω 2 ) =U ij,kl + U im,kn χ ph mn,pq (ω 1 − ω 2 )U pj,ql The factor 1/4 accounts for the exchanges of the external points of the two vertices. Note, that the second order diagram (obtained by replacing the susceptibility in Fig.  7b with a simple bubble) contributes to both Fig. 8b and Fig. 8c. As a function of ω 1 and ω 2 the high frequency part of the vertex function Γ ∞ has a structure of a constant plus two ridges along the main and minor diagonals, cross section of which is determined by the particle-hole and particle-particle dynamical susceptibilities.
In terms of the Baym-Kadanoff functional, the constant term in Γ ∞ derives from the first order, Hartree-Fock, diagram (Fig. 9a), and the diagrams containing at least one particle-hole or particle-particle bubble (Fig.  9b,c) by cutting the Green function lines corresponding to the same bubble.

Appendix B: Blocked IBS equations
In the following we present the IBS equation written in terms of small-ω and large-ω sectors denoted by indices 0 and 1, respectively. The IBS equation has the form where A, B and Γ are matrices indexed by the discrete Matsubara frequencies. Matrix B is diagonal in the Matsubara frequencies.
In the second step we assume full knowledge of Γ and B and want to compute A. Importantly, we are interested only in the sum over the matrix elements (Matsubara frequencies) of A. To this end we use a set of block equations A 00 = B 00 + B 00 Γ 00 + Γ 01 X 11 Γ 10 A 00 (B6) A 10 = X 11 Γ 10 A 00 (B7) A 01 = A 00 Γ 01 X 11 (B8) A 11 = X 11 + X 11 Γ 10 A 00 Γ 01 X 11 . (B9) The advantage of the blocked scheme is that equation (B5) does not have to be inverted in practice. This is thanks to three factors: i) only sum of elements of X 11 (or a weighted sum in case of X 11 Γ 10 ) is needed, ii) B 11 is small due to 1/ω 2 n decay, iii) Γ 11 is a constant plus a band matrix. We describe in detail the computation of The computation for X 11 Γ 10 is analogous. First, we solve iteratively equation (B5) taking only the band part of Γ 11 where the block indices 11 were dropped for simplicity. AfterX i• is computed the summation over the remaining index is performed to obtainX •• . Next, we should solve equation (B5) taking the constant part of Γ 11 andX in place of B 11 . This is not necessary, since there is a simple, well known, relationship between X •• andX •• Note, that quantities in this equation should be viewed as non-commuting matrices in the orbital indices.