A real-frequency solver for the Anderson impurity model based on bath optimization and cluster perturbation theory
Manuel Zingl Martin Nuss Daniel Bauernfeind Markus Aichhorn
cond-mat.str-el





[cor1]Corresponding author
add1]Manuel Zingl
[add1]Institute of Theoretical and Computational Physics,
Graz University of Technology, NAWI Graz, Petersgasse 16, 8010 Graz, Austria
add1]Martin Nuss
add1]Daniel Bauernfeind
add1]Markus Aichhorncor1
aichhorn@tugraz.at




Recently solvers for the Anderson impurity model (AIM) working directly on the real-frequency axis have gained
much interest. A simple and yet frequently used impurity solver is exact diagonalization (ED), which is based on a
discretization of the AIM bath degrees of freedom. Usually, the bath parameters cannot be obtained directly on the
real-frequency axis, but have to be determined by a fit procedure on the Matsubara axis. In this work we present
an approach where the bath degrees of freedom are first
discretized directly on the real-frequency axis using a large number 
of bath sites (≈ 50). Then, the bath is optimized by unitary transformations such that it separates into
two parts that are weakly coupled. One part contains the impurity site and its interacting Green's functions can be
determined with ED. The other (larger) part is a non-interacting system containing all the remaining bath sites. Finally,
the Green's function of the full AIM is calculated via coupling these two parts with cluster perturbation theory.



Anderson impurity model exact diagonalization cluster perturbation theory bath optimization






§ INTRODUCTION

The single-orbital Anderson impurity model (AIM) <cit.> can be
represented exactly by an interacting site coupled to a bath of
infinitely many non-interacting sites. In approaches based on exact
diagonalization (ED), the number of sites in the interacting system is
restricted, and thus the bath needs to be
truncated <cit.>. This is a delicate
step, because no unique procedure exists. Different ways are used,
e.g., fits on the Matsubara axis or continuous fraction
expansions <cit.>. 

Various methods improving on ED have been presented in
recent years, e.g., the variational exact
diagonalization <cit.>, the distributional exact
diagonalization <cit.> and methods based on a restriction
of the basis states <cit.>.
Another way of going beyond ED is the use of cluster perturbation
theory (CPT) <cit.>, i.e. the more advanced
variational cluster approximation (VCA) <cit.>, as a
solver for the AIM <cit.>. 

From now on, we assume a single-orbital AIM coupled to a finite but
large bath of L-1 non-interacting sites. The basic idea of using CPT
as an impurity solver is to separate the L-site AIM into a cluster
of size L_C, which includes the impurity site and L_C-1 bath
sites, and a non-interacting system consisting of the remaining bath
sites. In general, the non-interacting Green's function is specified by
the Hamiltonian H^0, that is a matrix in orbital space of size L
× L. For illustration purposes (see the sketch in
Fig. <ref>), we denote the upper left L_C× L_C block
in H^0 as the interacting cluster, subsequently H^0_C. The
remaining, lower (L-L_C)× (L-L_C) block describes the remainder
of the bath, subsequently H^0_R. Additionally, there are two off-diagonal 
blocks T connecting H^0_C and H^0_R. 
The onsite Hubbard interaction H_U = U n_I,↑ n_I,↓,
where I labels the impurity site, is now added to the cluster
Hamiltonian, H_C=H_C^0+H_U. There are no interactions in the bath
degrees of freedom, hence H_R^0 remains unchanged.

In CPT both Hamiltonians (H_C and H^0_R) are solved exactly
for their single-particle Green's functions G_C(ω) and
G_R(ω). To obtain G_C(ω) we use the Lanczos procedure at zero temperature <cit.>. 
Note that G_R(ω) = G^0_R(ω), as the
remainder of the bath is a non-interacting system. Subsequently, the
two systems are joined to yield the single-particle Green's function
of the full system G(ω) via the CPT relation <cit.> 

    G^-1(ω)    =
    ([ G_C(ω)      0;      0 G_R(ω);        ])^-1-V

where V is a L× L coupling matrix consisting only of the T blocks.
Eq. <ref> is exact in the case of a non-interacting system (U=0). In the
interacting case, the CPT relation is no longer exact, but a result of
perturbation theory in V. CPT approximates the self-energy of the
full system by the self-energy of the interacting cluster.  



In general, the non-interacting bath can always be transformed to a
tridiagonal representation via a Lanczos tridiagonalization, yielding
a chain representation of the AIM. This representation straight
forwardly allows to define the separation of the interacting cluster
and the remainder of the bath. However, the situation is not so clear in
other representations. Consider for example the case of a star
geometry, where all bath sites couple directly to the impurity site. 
Incorporating just a random set of these star sites into the
interacting cluster will lead to a poor discretization of the bath,
and hence a poor self-energy. 

Any unitary transformation on the non-interacting bath degrees of
freedom leaves the physics of the interacting AIM invariant. However,
such a transformation will influence the self-energy of the
interacting cluster significantly, since it changes the cluster
Hamiltonian H_C. Additionally, such transformations will also alter the
off-diagonal block T, rendering the resulting perturbation in some
cases larger than in others. There exist an infinite number of
representations which all describe the non-interacting bath exactly and
which are related via unitary transformations. However, the CPT method
itself suggests which baths might be the best: Those which
minimize the off-diagonal perturbative elements in T. The key
idea of this work is to use unitary transformations to find those
bath representations with minimal couplings between the cluster and the remainder of the bath.  

In the following, we outline a way to construct CPT-favorable bath
representations in Sec. <ref>, and present results for a L=64 AIM with a semi-circular particle-hole
symmetric bath in Sec. <ref>.



§ METHOD



The general form of the non-interacting Hamiltonian for an L-site AIM is

    H_0    = ϵ_I c^†_I c_I^ + ∑_i=1^L-1(t_iI c^†_i c_I^ +  h.c.) +
     ∑_i, j = 1^L-1 t_ij c^†_i c_j^ ,

where the impurity is denoted by the index I and the L-1 bath
sites by i and j. We omit the spin indices. 
To obtain H_0 for an L-site system one can use a star
representation, where each bath site couples only to the impurity
site. Then, the parameters of H_0 can be determined by a
discretization of the non-interacting bath DOS into equally spaced
intervals. Each interval is represented by a delta peak, where the
energy positions of the delta peaks correspond to the on-site energies
and the hopping parameters are obtained from the spectral weights in the
intervals. Of course, the higher the number of bath sites the better the result
of this discretization.

Under a unitary transformation R, performed in the bath only, with
 c_i^ =  ∑_α R_iα d_α^ and c_i^† = ∑_α d_α^† R^*_α i,
where R R^† =, the transformed Hamiltonian reads

    H'_0    = ϵ_I c^†_I c_I^ + ∑_α=1^L-1(h_α I d_α^† c_I^ +  h.c.) + 
     ∑_α, β = 1^L-1 h_αβ
       d_α^† d_β^ .

The parameters of the Hamiltonian transform like
 h_α I = ∑_i R^*_α i t_i I and
 h_αβ =∑_i,j R^*_α i t_ij R_jβ.
Such a transformation leaves the impurity state I and consequently ϵ_I invariant.

We define an energy of a certain bath representation via the 2-norm of the off-diagonal blocks T

    E = 1/N_T∑_i,j|T_ij|^2

where the number of elements in T is N_T = L_C·
(L-L_C). Transformations on the bath degrees of freedom included in
the interacting cluster do not influence the resulting
self-energy. The same is true for transformations performed only in the remainder of
the bath. This imposes a constraint on the energy E, namely, it has
to be invariant with respect to such transformations, which is indeed
fulfilled by the 2-norm. 

The aim is now to find an optimal bath representation
for CPT by minimizing the energy E. Since the configuration 
space of T_ij is high dimensional, we use 
a Monte Carlo procedure. Initially, we perform
global updates in all dimensions with random
rotation matrices to obtain a randomized starting representation of
H_0. Then, we move through the space of possible H_0 by proposing
random local updates R. In general, any unitary update would be allowed, 
but here we restrict ourselves to two-dimensional rotation matrices for the 
local updates 

    R = ([       1       0       …       0       0       …       0       …       0;       0       1       …       0       0       …       0       …       0;       …       …       …       …       …       …       …       …       0;       0       0       …  cos(ϕ)       0       … -sin(ϕ)       …       0;       0       0       …       0       1       …       0       …       0;       …       …       …       …       …       …       …       …       0;       0       0       …  sin(ϕ)       0       …  cos(ϕ)       …       0;       …       …       …       …       …       …       …       …       …;       0       0       …       0       0       …       0       …       1;         ])


A local update matrix R(i,j,ϕ) is drawn by choosing two random
integers i,j ∈ [1,L-1] representing the plane of rotation and one
rotation angle ϕ∈ [0,2π[. A new
representation with energy E' is accepted with probability p =
min(1, e^-γ(E'-E)). We use simulated
annealing to obtain low-energy CPT bath representations by increasing
the parameter γ. 

Although bath rotations leave the particle-hole symmetry invariant on
the L-site H_0, they destroy it on the L_C-site cluster. Therefore, 
as shown in Fig. <ref>, we split the bath sites into an equal 
amount of positive (blue elements)
and negative energy (red elements) sites and one zero mode (green
0). Updates are performed simultaneously on the positive and negative
modes which leaves the whole bath, the bath in the cluster as well as
the remaining bath particle-hole invariant.
To avoid a Kramers-degenerate ground state, clusters with an even
number of sites L_C are chosen. This implies that one bath site (the
zero mode) is exactly located at zero energy. Zero mode updates cannot be achieved by two-dimensional rotations without breaking the
particle-hole symmetry of the cluster, but would rather require
special unitary transformations involving at least three bath
sites. For the proof of concept presented here, we refrain from updating 
the zero mode, i.e. the green elements in Fig. <ref> do not change.
Hence, the zero mode coupling is determined by the initial
discretization of the system. Although this restricts the space of
trial bath representations, we leave the zero mode updates for future
works.  

Next to the energy E (Eq. <ref>), which reflects the
magnitude of the perturbation, we evaluate the influence of the CPT
truncation by comparing the non-interacting single-particle impurity
Green's function of the full system G^I_0(ω) to the one
considering only the sites in the cluster G^I_0,C(ω). The
resulting quantity 

    χ^2_C=∫ dω|G^I_0(ω+iη)-G^I_0,C(ω+iη)|^2

reflects the ability of the L_C cluster sites to represent the bath
degrees of freedom. A-priori a positive correlation of E and
χ^2_C is not ensured but expected. We emphazise that
χ^2_C is not used in the algorithm, but only serves as a
measurement for the quality of the bath optimization. A numerical broadening of η =0.02 is used to evaluate Eq. <ref>.

To asses the quality of the optimization scheme we also perform plain ED calculations for a truncated L=10 system. The Hamiltonian of this 10-sites system (with the Green's function G^I_0,ED(iω_n)) is obtained by fitting G^I_0(iω) on the Matsubara axis with the cost function 

    χ^2_ED=∑_n W_n |G^I_0(iω_n)^-1-G^I_0,ED(iω_n)^-1|^2

We employ the simplex search method by Lagarias et al. <cit.>
and impose particle-hole symmetry to reduce the number of fit
parameters. A Matsubara grid with 1024 points at a fictitious
temperature corresponding to β =100^-1 is
used, but the ED solution itself is obtained at zero temperature. 
The cost function Eq. <ref> is a heuristic choice and can also
take various other forms, e.g. with a different definition of the
distance or a different weight W_n (we set W_n = 1/ω_n) <cit.> . 
This leads to an ad-hoc determination of the bath
parameters and introduces some ambiguity to the solution of the
AIM. For a given number of bath sites also the discretization on the real-frequency axis 
is not uniquely defined, e.g. we could use unequally spaced energy intervals. 
In contrast to ED, where usually only a small number of bath sites is fitted on the Matsubara axis via the
minimization of Eq. <ref>, the discretization on the real-frequency axis can be easily 
performed for a large number of sites. Indubitably, the ambiguities in determining the AIM parameters are less severe for larger system sizes.



§ RESULTS


Here, we discuss a single interacting impurity in a particle-hole symmetric semi-circular bath with a half-bandwidth of 1. For the bath optimization we use a discretized system with a total number of L=64 sites and an interacting cluster of L_C=10. The interacting cluster includes the impurity site, one zero mode, and four additional positive and negative bath sites each. 

Fig. <ref> shows that smaller perturbative elements E correlate positively with smaller χ^2_C, indicating a better representation of the non-interacting bath by the sites contained in the interacting cluster. We compare intermediate bath representations (black crosses) to the random initial representation (orange star), a star representation by choosing ten sites at random to enter the interacting cluster (cyan triangle), and the chain representation cut after 10 sites (green circle). Note that the chain representation hosts only one perturbative matrix element, which is however large. In all cases one finds a higher χ^2_C with respect to the final result of the optimization (blue cross).  Of course, due to the optimization the number of elements in T grows, but their magnitude becomes tiny compared to the chain or the star representation.






In Fig. <ref> the self-energy Σ(iω_n) of the optimized system on the Matsubara axis is shown for different values of U. We compare it to results obtained at an inverse temperature of β=100^-1 with the TRIQS/CTHYB solver <cit.>, which is based
on continuous-time quantum Monte Carlo in the hybridization expansion (CTHYB) <cit.>. Additionally, we show the self-energies calculated with plain ED for a L=10 system.
As the CPT self-energy of the full AIM is given by the self-energy of the 10-site 
cluster, we can only expect it to be compatible with CTHYB on a similar level as the ED self-energies are. For U<=4 the ED 
self-energies are even slightly closer to the CTHYB result (see Fig. <ref>).
Nevertheless, the important message conveyed by Figs. <ref> 
and <ref> is that the minimization of the off-diagonal
block T is a proper procedure to obtain a good representation of the 
full AIM system by the 10-site cluster.

The actual advantage of CPT is revealed on the level of the spectral function A(ω), which results from coupling the interacting cluster 
to the remainder of the non-interacting bath.
In Fig. <ref> we show the CPT spectral function (red) of the
optimized system (L=64) for interaction values U=0, 1, 2,
 and 3. We compare our results to spectral
functions obtained with the density matrix renormalization group (DMRG) <cit.> and real-time evolution (as in 
Ref. <cit.>) using matrix product states (MPS), colored in black. 
This approach is able to provide an excellent spectral resolution on the whole real-frequency axis <cit.>. We use the star representation of the L=64 AIM for the MPS calculation with a truncated weight of 10^-10. Additionally, we show the results of the 10-site ED system (gray).
In contrast to ED, which exhibits strong finite size effects, the CPT optimization scheme provides smoother spectral functions which are in much better agreement with the MPS results. This difference is particularly pronounced for the smaller 
U values, which is a consequence of CPT becoming exact for U→ 0.
Up to the influence of the discretization, CPT
reproduces the non-interacting spectral function for U=0 
(see top grap in Fig. <ref>). Of course, for higher values of U the energy resolution is limited by the size of the exactly solved system (L_C), which can be improved on by including more sites in the interacting cluster.






§ CONCLUSION

In this work we introduced a bath optimization scheme for the Anderson
impurity model. Using unitary transformations in the bath degrees of
freedom, we minimize the coupling between a small cluster
containing the interacting impurity site and the remaining sites of the
bath. These transformations leave the impurity DOS of the
non-interacting bath invariant. In general, the proposed scheme can be
useful for all CPT-based methods when parts of the entire system are
non-interacting, but it does in principle also provide a guideline to
construct finite-size representations of hybridization functions as
needed, e.g., in the framework of dynamical mean-field theory. For a
large enough number of bath sites, the initial AIM can be obtained
directly on the real-frequency axis, and thus a fit on the Matsubara
axis can be avoided. In this work we have presented a proof of
concept, but anticipate to explore the bath optimization scheme for
systems without particle-hole symmetry and multi-orbital impurity models.



§ ACKNOWLEDGMENTS

The authors acknowledge financial support from the Austrian Science Fund
(FWF) through SFB ViCoM F41 (P04 and P03), project P26220, and through 
the START program Y746, as well as from NAWI-Graz. The CTHYB results were calculated using the TRIQS library <cit.> and the TRIQS/CTHYB solver <cit.>. 
The MPS results were obtained using the ITensor library <cit.>.


§ REFERENCES

elsarticle-num

