{'fraction_non_alphanumeric': 0.03240789612583131, 'fraction_numerical': 0.009975720468700518, 'mean_word_length': 4.622255192878338, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Natural orbital impurity solver for real-frequency properties at finite temperature

We extend the natural orbital impurity solver [PRB 90, 085102 (2014)] to finite temperatures within the dynamical mean field theory and apply it to calculate transport properties of correlated electrons. First, we benchmark our method against the exact diagonalization result for small clusters, finding that the natural orbital scheme works well not only for zero temperature but for low finite temperatures. The method yields smooth and sufficiently accurate spectra, which agree well with the results of the numerical renormalization group. Using the smooth spectra, we calculate the electric conductivity and Seebeck coefficient for the two-dimensional Hubbard model at low temperatures which are in the scope of many experiments and practical applications. These results demonstrate the usefulness of the natural orbital framework for obtaining the real frequency information within the dynamical mean field theory.


I. INTRODUCTION
Understanding correlated electron systems is one of the central issues in condensed matter physics. To address the issue, the dynamical mean-field theory (DMFT) [1][2][3] has been developed and extended in its scope in the past decades. The DMFT properly takes into account the temporal fluctuations which are essential to describe the Mott metal-insulator transition while it ignores the spatial fluctuations which could be important for quasi lowdimensional fluctuating systems such as cuprates. Combined with ab initio calculations, the DMFT has also been used to calculate the electronic structure of real materials, where several orbital degrees of freedom are usually involved 4,5 .
In the DMFT, we need to solve a quantum impurity model which is usually simpler than the original lattice model, but still difficult to solve analytically. Many numerical methods have been proposed for solving that model. The continuous-time quantum Monte Carlo (CT-QMC) method 6 is one of the most commonly used methods since it gives a numerically exact result within a statistical error and is applicable to relatively large degrees of freedom. However, because the formalism is based on the imaginary (Matsubara) time axis, an analytic continuation of the numerical data with statistical errors, which is an ill-posed inverse problem, is required to obtain realfrequency properties. While the numerical analytic continuation method has been improved 7-9 , the dependence on the data precision is unavoidable and it is difficult to correctly describe a high-frequency structure or sharp structure like the Drude peak.
The real-frequency data is, however, indispensable to compare with spectroscopic experiments like optics, photoemission and scanning tunneling spectroscopies, as well as electronic Raman and resonant inelastic X-ray scatterings. Besides these, it is also necessary for calculating transport properties like electric and thermal conductivities, as well as a thermoelectric effect measured by the Seebeck or Nernst coefficient. Furthermore, the real-frequency structure of Green's function and self-energy can give crucial information on underlying the mechanism of physical phenomena, as exemplified by the electron-phonon mechanism of conventional superconductivity established through the analysis of the frequency-dependent superconducting gap function 10 .
Thus, to apply the DMFT to these issues of real materials, we need an impurity solver which allows a direct access to real-frequency properties and is potentially applicable to multiorbital systems. One way in this direction is to improve the numerical renormalization group (NRG) 11,12 or density-matrix renormalization group (DMRG) method [13][14][15] for multiorbital systems, as indeed has been done in Refs. 16-18 up to three orbitals. Recently, the fork tensor-product state was proposed and applied to a three orbital system 19 .
Another way is to improve the exact diagonalization (ED) 20 method. The problem in the ED is a limited number of the bath sites, which results in a discrete spectrum of the impurity Green's function. Even with modern computers, it is difficult to deal with more than 20 degrees of freedom (except for spin) including both impurity and bath sites. Because of this limitation, it has been difficult to apply the ED solver to more than four site or orbital degrees of freedom. Also, to mitigate the finite-size effect, the spectrum is usually calculated with an energy-smearing factor, which defines the width of each discrete peak. While this does not alter the position of the spectral peaks, it changes the broadness of the spectra and can considerably affect the values of transport coefficients so that they cannot be determined unambiguously.
Recently, an interesting idea was proposed in Refs.21 and 22 to remedy this problem of the finite bath. Representing the bath sites in the natural orbital basis, Lu et al. constructed efficient algorithms to take account of low-energy excitations beyond the size tractable with the conventional ED. They calculated a real-frequency Green's function in a restricted Hilbert space, where the significant reduction of its dimension allowed them to consider ∼ 1000 bath sites, to yield a sufficiently smooth spectrum. Considering the tractable size of the Hilbert arXiv:2107.06517v1 [cond-mat.str-el] 14 Jul 2021 space in the ED, the improvement by the natural orbitals may enable us to deal with up to five orbitals in near future, where many correlated materials -in particular transition-metal compounds with five d orbitals -await the clarification by theory.
The idea of the natural-orbital solver is in line with the configuration interaction (CI) method which has recently been applied to the DMFT by several groups [23][24][25][26] . The interesting point is that the CI method seems to work well for the quantum impurity model. The CI method is widely used for molecules in quantum chemistry where it is easier to justify assuming fully filled localized molecular orbitals and only consider the mixture of Slater determinants for the rest degrees of freedom. On the other hand, in metallic states in solids, orbitals are more hybridized with each other and it is questionable to apply the CI method. Actually, the ground state of the Hubbard cluster is quite entangled and difficult to describe by a few Slater determinants. Nevertheless, for the quantum impurity model, it turned out that the ground state is well described by a small number of Slater determinants and low-energy particle-hole excitations from it 24 . The natural orbital framework gives us a general way to construct the new basis set where the CI scheme works efficiently. Recently, it is suggested that the natural orbital configuration is also useful for tensor network calculations 27 .
The application of the natural orbital method has, however, been limited to zero temperature. In view of the comparison with experiments, extension to finite temperatures is indispensable. In this paper, we attempt such an extension, by constructing different natural orbitals for each electron-spin number sector. We find that the minimum energy of each sector and the Green function are well approximated within the scheme for low temperatures at least. We benchmark the accuracy of our result for the small cluster against a direct ED calculation and show the smooth spectrum with the large number of bath sites which is consistent with the previous NRG result. Then, we calculate transport properties for the two-dimensional Hubbard model. We obtain consistent results with the previous CT-QMC+Padé studies, but the current scheme can potentially reach much lower temperatures for general models, which are relevant to many experiments and practical applications.
The paper is organized as follows: In Sec. II, we give an overview of the natural orbital framework and introduce our extension for different sector calculations toward finite temperatures. In Sec. III, we present some benchmark results for small clusters and large systems at zero temperature. Then, we show our results for conductivity and Seebeck coefficient for the two-dimensional Hubbard model on a square lattice.

II. FORMALISM
We apply the idea of natural orbitals 21,22 to the finitetemperature ED solver 20,28,29 of the DMFT. In the ED, we first divide the Hilbert space to independent sectors specified by some quantum number, instead of dealing with the whole Hilbert space at once. In this paper, we assume a spin-conserved system, so that (n ↑ , n ↓ ) specifies each sector, where n σ is the number of spin-σ electrons in the whole system. We then use the Lanczos algorithm for each sector.
We basically follow the previously proposed algorithm of the natural orbital solver with the projection approach for calculating the Green function 22 . The new point here is that we need to consider not only the ground state but also all sectors which contribute to the Green function at finite temperatures. For this purpose, we construct different natural orbitals for each sector to evaluate the energy difference between each sector accurately. Even if the spin-orbit coupling exists, the general framework will not change and we just need to change the way of the sector decomposition based on the conserved quantity of the system. In this paper, we focus on the single band model. The extension to the multi-orbital system would be straightforward, see, e.g., Ref. 27.
The formalism is summarized in Fig 1. We start with the standard Anderson impurity model where imp is the energy level of the impurity site, and a † imp,σ (a imp,σ ) and a † i,σ (a i,σ ), respectively, are the electron creation (annihilation) operator at the impurity site and the i-th bath site with spin σ ∈ {↑, ↓}. U is the onsite Coulomb repulsion, n σ ≡ a † imp,σ a imp,σ , and V i is the hybridization to the i-th bath site.
First, we choose a sector (n ↑ , n ↓ ). Then, following Ref. 22, we make the mean-field Hamiltonian (H MF ) whose energy level is adjusted to keep the density at the impurity site. We calculate the density matrix with the lowest-energy state of H MF within the sector (n ↑ , n ↓ ). Natural orbitals are obtained as the basis set diagonalizing the bath density matrix. Notice that even if we assume the paramagnetic system, we need to construct spin dependent natural orbitals to calculate the contribution from the spin imbalanced sector. After that, we can construct the natural orbital representation of the quantum impurity Hamiltonian (H NO AIM ) for each sector in the same way as the previous study in Refs. 21 and 22.
The natural orbital Hamiltonian contains the impurity site (imp), active bath site (b), and valence (v)/ conduction (c) chains. We then determine the cutoff length n L for the ED calculation. Outside this cutoff, we assume the fully occupied (empty) configuration for the valence (conduction) chain. Thus, the whole Hamiltonian is de- coupled into two terms as where , n c(v) are the energy levels, hoppings of conduction (valence) sites and the number of these sites, respectively, and b are the electron creation (annihilation) operators at the active bath site, i-th conduction-bath site, and i-th valence-bath site with spin σ, respectively (c 0 = v 0 ≡ b). The number of electron within the small system, (n I ↑ , n I ↓ ), is automatically determined as and then we perform the exact diagonalization for H I with (n I ↑ , n I ↓ ) electrons. From this, the energy of the whole system (E AIM ) is calculated as where E HI is the direct result from the ED. The contribution to the Green function from each sector is computed by the projection method proposed in Ref. 22 with p = 2 (up to double particle-hole excitations to the H chain ) and n M = 400 (Krylov dimension). Since the calculation for the Green function of each sector is fully independent of each other, we can just sum-up all the sector contributions with the Boltzmann weight e −β(EAIM−EGS) , where E GS is the ground-state energy of the whole system.
After obtaining the Green function, we complete the DMFT self-consistent loop in a usual way. Since we employ a large number of bath sites to represent a real frequency structure, we can perform the bathfitting directly at real frequencies with equal weight configurations 30 to determine the bath parameters in Eq. (3).

A. Benchmarks for a small cluster
In this section, we start with some benchmarks of finite-temperature results for a small cluster. First, we study the 12-site impurity problem (one impurity site and eleven bath sites) and test the accuracy of the natural orbital solver against the direct ED calculation. For the natural orbital calculation, we diagonalize 6-site H I problem (i.e., n L = 2) with assuming a fixed configuration for the rest 6 bath degrees of freedom. Namely, we construct conduction-and valence-bath chains with these 6 degrees of freedom and connect them to H I . We employ the semicircular density of states ρ( ) = 2 (1 − ω/D) 2 /πD and U = D where D is the half of the bandwidth. The electron density is set to n = 0.8, which is a case hole-doped to the half filling.
In Fig. 2, we show the lowest energy at each sector, (n ↑ , n ↓ ), obtained with the ED and natural orbital calculation. The parameters of the impurity model are determined from the converged 12-site DMFT-ED calculations at two different temperatures, (a) T = 0.02D and (b) T = 0.05D. From ED calculations, (n ↑ , n ↓ ) = (6, 6) is the ground state sector and we show n ↑ dependence of the lowest energy of other sectors for n ↓ = 6 (circles) and n ↓ = n ↑ (squares). The results show a good agreement in particular at low energies which are important to calculate the Green function. At n ↑ = 2 and 10, we find some deviation from the ED result which suggests that the natural orbital representation becomes worse for the sectors that the length of the chain becomes nearly zero, i.e., n c σ (n v σ ) ∼ 0. In general, such states would have a relatively high energy and their contribution to the Green function is small. Furthermore, such contribution becomes more and more negligible as the total number of bath sites (chain length) increases, which we can reach a few hundreds now. Thus, we can conclude that the natural orbital configuration has an accuracy enough to describe the energy not only for the GS-sector but also for all energy configurations of the system.
We also compare the spectrum in Fig. 3 which is obtained from the imaginary part of the Green function as A(ω) = − G(ω)/π. The agreement here is also quite good and especially, they match perfectly at low frequency region. This result demonstrates that the projection method 22 with p = 2 (i.e., up to double particlehole excitations to the H chain ) for calculating the Green function works quite well for the temperatures we study.

B. Spectrum for the large system
Now we move to the result with the large number of bath sites which can not be addressed by the original ED method. Here, we take 301 bath sites and 10 sites for the ED, i.e., n I = 10 and n L = 4.
In Fig. 4, we show the DMFT spectrum at T = 0  for semicircular density of states for several interaction strengths U = 1 − 5D at n = 0.8. The result shows a good agreement with the previous NRG result 31 in a low energy regime and in an overall structure. One difference is in the upper Hubbard band, where it shows an additional peak at the gap edge. Such an peak was reported at half filling by the recent study using matrix product states 32,33 and the D-DMRG study 34 . Our result indicates that the current method works well to obtain the smooth spectrum for calculating real-frequency properties.

C. Real-frequency properties at finite temperatures
Finally, we show the result of real-frequency properties at finite temperatures. Here, we employ the square lattice Hubbard model with only the nearest neighbor hopping t = D/4 and U = 3.5D. Here, we take 101 bath sites, 8 sites for the ED calculation and sum up the Green function contribution only from 6 × 6 sectors around the ground state. As shown in Fig. 5(a,b), the spectrum obtained with the present method is consistent with the previous CT-QMC+Padé results 35 . We obtain the clear upper Hubbard band structure which is sometimes difficult to obtain through a numerical analytical continuation because it is located at a high energy, i.e., away from the Matsubara axis. The smooth spectrum allows us to calculate the resistivity ρ and Seebeck coefficient S from 35 where f is the Fermi distribution function and v x (k) ≡ 0.5D sin k x is the Fermi velocity along the x axis. In Fig. 5(c,d), we show these results at around T = 0.01D. The temperature dependence is also consistent with the CT-QMC+Padé result. This behavior of Seebeck coefficient in doped Mott insulator is commonly observed in many lattice models [36][37][38] . Thus, our calculations demonstrate that the natural orbital configuration provide us an efficient way to approach low temperature physics.

IV. CONCLUSION AND OUTLOOK
In conclusion, we have extended the natural orbital impurity solver to finite temperatures by constructing different natural orbitals for each sector. We have found that the minimum energy of each sector and the Green function are well approximated within this scheme at low but finite temperatures. We have examined the accuracy of our result for the small cluster against direct ED calculation and obtained a good agreement. We have also shown that the smooth spectrum can be obtained with a large number of bath sites. Then, we have calculated real-frequency properties for the two-dimensional Hubbard model. Our results demonstrate the usefulness of the natural orbital framework for obtaining the real frequency information at finite temperatures.
The ED itself has no difficulty for the low temperature calculations, but when used as a DMFT impurity solver, it has severe limitations due to the finite bath site fitting. The natural orbital scheme removes this problem and we can take the advantage of the original strength of the ED for low temperature calculations.
In this paper, we constructed the different natural orbitals at each sector for describing the accurate ground state, while it may not be suitable for excited states in each sector. In this sense, the current treatment may still suffer from the finite size effect which would become important at higher temperatures. Constructing the different natural orbitals even for each excited state may accelerate the convergence against n L (i.e., reducing the finite size effect).
While the present work focuses on the extension to finite temperatures, the extension to multi-orbitals is also important for the application to real materials and this is another direction of future study. There, it is expected that we can easily take off-diagonal elements and general type interactions into account.