We start with writing down the Anderson impurity Hamiltonian

_{} (1)

describing the
interaction of the impurity level _{} with
bands of conduction electrons _{} via hybridization _{}. _{} is the Coulomb repulsion between different orbitals in the _{}–band. In
order to solve the Anderson
impurity model in general multiorbital case, we use a rational interpolative
formula for the self–energy. This can be encoded into a form

_{} (2)

The coefficients
_{}, _{}, or, alternatively, the poles _{}, zeroes _{} and value
of the self-energy at infinity _{} in this equation are to be determined. The form (2) can be also
viewed as a continuous fraction expansion but continuous fraction
representation will not be necessary for the description of the method.

Our basic
assumption is that only a well distinct set of poles in the rational
representation (2) is necessary to reproduce an overall frequency dependence of
the self–energy. Extensive experience gained from solving Hubbard and periodic
Anderson models within DMFT at various ratios of the on–site Coulomb
interaction _{} to the bandwidth _{} shows the appearance of lower and upper Hubbard bands as well as
renormalized quasiparticle peak in the spectrum of one–electron excitations.

It is clear
that the Hubbard bands are damped atomic excitations and to the lowest order
approximation appear as the positions of the poles of the atomic Green’s
function. In the _{} symmetry case which is described by the Hamiltonian (1), these
energies are numerated by the number of electrons occupying impurity level,
i.e. _{} and the atomic Green’s function takes a simple form

_{} (3)

where _{} are the
probabilities to find an atom in configuration with _{} electrons
while combinatorial factor _{} appear due to equivalence of all states with _{} electrons
in _{}.

We can represent the atomic Green’s function (3) using the rational representation (2), i.e.

_{} (4)

where _{} are all _{} atomic
poles, while _{} denote _{} zeroes with _{} being the total degeneracy. The centers of the Hubbard bands are
thus located at the atomic excitations _{}. Using
standard expression for the atomic Green’s function _{}, we
arrive to a desired representation for the atomic self–energy

_{} (5)

Using this
functional form for finite _{} modifies
the positions of poles and zeroes via recalculating probabilities _{} which is
equivalent to the famous Hubbard I approximation.

We now
concentrate on the description of the quasiparticle peak which is present in
metallic state of the system. For this an extra pole and zero have to be added
in Eq (5). To see this, let us consider the Hubbard model for the _{} case
where the local Green’s function can be written by the following Hilbert
transform _{} If self–energy lifetime effects are ignored, the local spectral
function becomes _{} where _{} is the non–interacting density of states. The peaks of the spectral
functions thus appear as zeroes in Eq. (5) and in order to add the
quasiparticle peak, one needs to add one extra zero (denoted hereafter as _{} to the
numerator in Eq. (5). To make the self–energy finite at _{} one has
to also add one more pole (denoted hereafter as _{}) which
should appear in the denominator or Eq. (5). Furthermore, frequently the
Hartree Fock value for the self–energy can be computed separately and it is
desirable to have a parameter in the functional form (5) which will allow us to
fix _{}. An obvious candidate to be changed is that self–energy pole in (5)
which is closest to _{} equal zero. Let us denote this parameter as _{} and
rewrite the denominator of (5) as _{} where
the product is now extended over all zeroes of the atomic Green’s functions
except the one closest to zero and two extra poles _{} and _{} can
control the width of the quasiparticle peak and _{} Thus, we
arrive to the functional form for the self–energy

_{} (6)

We are now
ready to list all constrains of our interpolative scheme. To fix the unknown
coefficients _{} _{} _{} _{} in Eq. (6) and to write down the linear set of equations for the coefficients_{}, _{} in Eq.
(2) we use the following set of conditions.

*a)* *Hartree
Fock value* _{} In the limit _{} the self–energy takes its Hartree–Fock form

_{} (7)

Mean level
occupancy* *_{}is
defined as a sum over all Matsubara frequencies for the Green’s function, i.e.

_{} (8)

where

_{} (9)

defines the
impurity Green’s function and _{} is the
hybridization function.

*b)* *Zero–frequency
value *_{} The so called Friedel sum rule establishes the relation between the
total density and the real part of the self–energy at zero frequency

_{} (10)

*c) Quasiparticle mass renormalization value *_{} The
slope of the self–energy at zero frequency controls the quasiparticle residue, _{} using
the following relationship

_{} (11)

Formally,
constraints (*b*) and (*c*) hold for zero temperature only but we
expect no significant deviations in many regions of parameters as long as we
stay at low enough temperatures. The behavior may be more elaborated in the
vicinity of Mott transition.

*d)* *Positions
of Hubbard bands.*

As we discussed, in order that the self–energy obeys the atomic limit and places the centers of the Hubbard bands at the positions of the atomic excitations, we demand that

_{} (12)

This condition
fixes almost all self–energy zeroes _{} in Eq.(2)
to the poles _{} However, it alone does not ensure that the weight is correctly
distributed among the Hubbard bands and that the very distant Hubbard bands
disappear. For this to occur, distant poles of Green’s function have to be
canceled out by nearby zeroes of the Green’s function. It is clear that each
pole _{} far from the Fermi level has to be accompanied by a nearby zero _{} in order
the weight of the pole be small. Thus, the self–energy has poles at the
position of Green’s function zeroes which can be encoded into the constrain

_{} (13)

We want to
keep this property of the self–energy for finite _{} and thus
demand that self–energy diverges (when lifetime effects are kept, it only
reaches a local maximum) at the zeros of the functional form (3) of _{} Note
that the relationship (13) holds (approximately) for frequency _{} larger
than the renormalized bandwidth _{} Therefore
the information about one _{} which lies close to _{} is
omitted and replaced by the information about _{} _{} and _{} as it is
done by separating _{} and _{} in the denominator of Eq. (6).

We can now
write down a set of linear equations for all unknown coefficients in the
expression (2). There is total _{} of
parameters _{} and _{} where we can always set _{} The
conditions *a)*,*b),c) *give

_{} (14)

_{} (15)

_{} (16)

_{} (17)

According to
condition *d)* we can use all _{} poles _{} and _{} zeroes _{} The zero
_{} closest to _{} is dropped out. This brings additional _{} equations
for the coefficients and makes _{} as the degree of the rational interpolation which are written below

_{} (18)

_{} (19)

Note that
while _{} may be rather large, the actual number of poles contributing to the
self–energy behavior is indeed very small. We can directly see this from Eq.
(5) which uses all _{} poles _{} fulfilling
Eq. (12) and uses _{} zeroes _{} directly
related to _{} poles _{} Clearly, when the
spectral weight of the atomic excitation becomes small, the corresponding _{} becomes
close to _{} and the cancellation occurs. Therefore in realistic situations when
only the upper and lower Hubbard bands have significant spectral weight along
with the quasiparticle peak, the actual degree of the polynomial expansion is
either two or three. However, it is advantageous numerically and cheap
computationally to keep all poles and zeroes in Eq. (6) because the formula
automatically distributes spectral weight over all existing Hubbard bands.

In the limit
when _{} the self–energy automatically translates to the non–interacting
one. The atomic poles get close to each other but, most importantly, their
spectral weight goes rapidly to zero as it gets accumulated within the
quasiparticle band.

In the Mott
insulating regime, the conditions *b)*
and *c)* drop out while all poles _{} and zeroes _{} can be used to
determine the interpolation. However, in this regime it does not matter whether
one of _{} closest to _{} is dropped out or
kept, since we can always replace this information by information about _{} Therefore
the Mott transition can be studied without changing the constraints.

We thus see
that in the insulating case the self–energy correctly reproduces the well–known
result of the Hubbard I method where the Green’s function is computed after Eq.
(9) with atomic self–energy. If the lifetime effects are computed, the
parameters _{} and _{} become complex and the
Hubbard bands will acquire an additional bandwidth. This effect is evident from
the simulations using various perturbative or QMC impurity solvers and can be
in principle incorporated into the interpolative formulas (2) or (6). However,
in in the current implementation we omit it for simplicity.

We thus see
that the interpolational scheme is defined completely once a prescription for
obtaining parameters such as _{} _{} as well as poles _{},
and zeroes _{} is given. For this purpose we will test two commonly used methods:
SBMF method due to Gutzwiller [1] as described by Kotliar and Ruckenstein [2]
and the well–known Hubbard I approximation [3]. We compare these results
against more accurate but computationally demanding quantum Monte Carlo calculations and establish the
procedure to extract all necessary parameters.

Note that once
the constraints such as _{} are computed
from a given approximate method, some of the quantities such as the total
number of particles, _{}, and the value of the
self–energy at zero frequency, _{} can be computed fully
self–consistently. They can be compared with their non–self–consistent values. If
the approximate scheme already provides a good approximation for _{} and
satisfies the Friedel sum rule, the self–consistency can be avoided hence
accelerating the calculation. Indeed we found that inclusion of the
self–consistency improves the results only marginally except when we are in a
close vicinity to the Mott transition but here we do not expect that our simple
interpolative algorithm is very accurate.

[1] M.
Gutzwiller, Phys. Rev. **134**, A923
(1964).

[2] G. Kotliar
and A. E. Ruckenstein, Phys. Rev. Lett. **57**,
1362 (1986).

[3] J.
Hubbard, Proc. Roy. Soc. (London) A**281**, 401 (1964).

· Compile program typing "make" at the source files location. Makefile was tested on Linux OS using PGI compiler. Please adjust it depending on operating system and compiler used.

·
Edit “**input**" file which has the following structure:

1 |
IMOD |

0.0 |
EF |

2.0 |
U |

0.016 |
TEMP |

4 |
NDEG |

500 |
NMSB |

500 |
OMAX |

500 |
NOMG |

50 |
WEND |

F |
COMPUTE_REAL |

where IMOD = 1 is the Hubbard Model, IMOD=3 the Anderson impurity model, EF is the impurity level, U is value of Coulomb repulsion, TEMP is temperature to be used used to create for Matsubara points grid, NDEG is the degeneracy, NMSB is the number of Matsubara point, OMAX is the imaginary frequency cutoff, NOMG is the number is real frequency points, WEND is the real frequency cutoff, COMPUTE_REAL is flag once is “True" tell program to produce the self-consistent solution (last iteration) on real axis.

1.
Program’s input consists from one
more file (provided the Anderson impurity model item is chosen in “**input**"
file):

“**delta.dat**" containing the
hybridization function.

The structure is the following:

Column # |
Value |

1 |
Frequency |

2 |
Re Delta |

3 |
Im Delta |

2. Run the program executable (“main" is the default name).

3. Program’s output consists from the following files:

1) “**gfsig_iw.dat**" containing Green’s
function (GF) and the self-energy (SE) on Matsubara axis

2) “**gfsig_re.dat**" containing Green’s
function (GF) and the self-energy (SE) on real axis provided flag
"COMPUTE_REAL" is “TRUE".

They have the same structure:

Column # |
Value |

1 |
Frequency |

2 |
Re GF |

3 |
Im GF |

4 |
Re SE |

5 |
Im SE |

3) “**grid_re.dat**" containing real frequency
grid for the hybridization function (delta).

4) “**grid_im.dat**" containing Matsubara
frequency grid for the hybridization function (delta).

Both files have the same structure:

Column # |
Value |

1 |
Frequency |

Makefile |
main make file |

dmf_cmpdiag.f |
The solution of the generalized eigenvalue problem |

imp_sunatm.f |
Solves Anderson impurity model, returns GF and SE. |

imp_sunmas.f |
Computes effective mass renormalization |

imp_sunrat.f |
Finds rational interpolation or self-energy |

lib_broy6.f |
Broyden mixing. |

lib_cinv.f |
Finds inverse of square matrix. |

lib_csplines.f |
Splines complex function from one to another grid. |

lib_deriv1.f |
Calculates radial derivative. |

lib_morefun.f |
Calculations of a few auxiliary functions |

lib_pade.f |
Realization of Pбde procedure. |

lib_pcoefs.f |
Finds poles |

lib_pzeros.f |
Finds zeros |

lib_simq.f |
Finds solution of a set of linear equations |

lib_splin3.f |
Computed a natural spline approximation of third order. |

mod_common.f |
File containing common modules used across the program. |

mod_dimart.f |
File containing common modules used across the program. |

mod_init.f |
File containing common modules used across the program. |

sun.f |
Main program. |