{'fraction_non_alphanumeric': 0.031396663839411934, 'fraction_numerical': 0.009443030817076618, 'mean_word_length': 4.271704299873314, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational cluster approach to the single impurity Anderson model

We study the single impurity Anderson model by means of cluster perturbation theory and the variational cluster approach (VCA). An expression for the VCA grand potential for a system in a non interacting bath is presented. Results for the single-particle dynamics in different parameter regimes are shown to be in good agreement with established renormalization group results. We aim at a broad and comprehensive overview of the capabilities and shortcomings of the methods. We address the question to what extent the elusive low energy properties of the model are reproducible within the framework of VCA. These are furthermore benchmarked against continuous time quantum Monte Carlo calculations. We also discuss results obtained by an alternative, i.e. self consistent formulation of VCA, which was introduced recently in the context of non equilibrium systems.


I. INTRODUCTION
In recent years both the applications for strongly correlated quantum impurity models and the number of successful approaches to harvest their physical results have grown enormously. Those models were introduced to describe the effects of magnetic transition metal impurities immersed in metallic hosts 1,2 . Originally they were derived to capture remarkable physical properties like the resistance minimum 3,4 at a specific temperature scale T K 5 or the anomalous magnetic susceptibility and specific heat of such materials. Today a whole realm of applications for quantum impurity models has opened. They describe the physics of quantum dots and wires [6][7][8] as well as molecular electronics 9 . Applications range from nanoelectronics all the way to quantum information processing 10 . Their properties are essential for todays technological applications in single electron transistors 11 exhibiting the Coulomb blockade effect 12 or in devices dominated by RKKY interaction [13][14][15] . The behavior of various magnetic phenomena and the fascinating branch of heavy fermion physics is described by strongly correlated quantum impurity models 16,17 . Recent studies have shown that the remarkable material Graphene exhibits Kondo physics 18 which may be investigated theoretically by virtue of quantum impurity models. These models have further been applied to understand the adsorption of atoms onto surfaces [19][20][21] . In addition, they are of theoretical importance as solvable models of quantum field theories 22,23 . A renewed interest in understanding and calculating dynamic quantities of these models was created with the advent of dynamical mean-field theory (DMFT) [24][25][26] . In the foundations of this theory quantum impurity models have to be solved as an auxiliary problem. A wide range of methods and approximations have been suggested for the solution of quantum impurity models. They however prove to be a very delicate subject because standard perturbative approaches diverge 5 . Prominent methods to gain physical conclusions include a self consistent perturbative expansion 27 and Bethe Ansatz techniques 28 for one dimensional problems. The low energy physics is very well described by numerical renormalization group (NRG) 29 and in some limits also by functional renormalization group (FRG) 30,31 , and density matrix renormalization group (DMRG) [32][33][34] . There is a range of slave particle methods 35,36 available as well as methods based on Hubbard's X-operator technique 37,38 and calculations using variational wave functions 21 . Valuable physical insight has been gained by using equation of motion techniques applying different approximation schemes 39 . For moderate system sizes the Hirsch-Fye Quantum Monte Carlo (QMC) 41 algorithm has proven to achieve good results. In the past years different approaches to continuous time QMC 40 have been applied very successfully to solve quantum impurity models especially in application with DMFT. In this context exact diagonalization (ED) methods have been explored to solve small systems 42 . As of today some limits of quantum impurity models are understood with great precision but there appear several gaps to be bridged. The low energy properties of these models are reproduced very well by renormalization group based approaches (i.e. NRG). These approaches in general have trouble to capture the high energy parts of the spectrum. The same may be said about QMC methods which if applicable yield dynamic quantities in imaginary time. The analytic continuation to the real energy axis is ill conditioned. Spectra obtained by, for example, the maximum-entropy method 43,44 have a large uncertainty for higher energies. Exact diagonalization methods, in principle, grant access to low as well as high energy parts of the spectrum at the same time. Due to the prohibitively large Hilbert space however only small systems (about ten to twenty sites) may be treated with this method, whose low-energy behavior is expected to deviate from the one of the infinite lattice significantly. Nevertheless, the advantage consists in the fact that the spectral properties may be determined directly on the real energy axis. Besides the issue of the low energy scale, also the flexibility to adapt to various impurity configurations and geometries is limited in many methods. NRG has been successfully applied only to the one and two impurity case so far. QMC approaches may suffer from the sign problem for more complex multiband models 40 . The region of large interaction strength is naturally difficult to treat in standard perturbative/diagrammatic approaches (i.e. diagrammatic perturbation theory or FRG). In the present work we test cluster perturbation theory (CPT) 45,46 and the variational cluster approach (VCA) [47][48][49] on the single impurity Anderson model 1 . The great flexibility and versatility of CPT/VCA allows for obtaining approximate single-particle dynamic quantities and static expectation values in all parameter regions of any lattice impurity model with local interactions. However these many body cluster methods can not be expected to describe the low energy excitations as accurately as specifically tailored methods do. It is however interesting to see whether the correct low energy behavior may be reproduced at least to some extent. CPT as well as VCA bare several advantages 58 : i) They yield spectra directly on the real axis and ii) also the high energy incoherent part of the dynamics becomes available. iii) They are applicable in all parameter regions and also at high interaction strengths. iv) They have the advantage of comparatively low computational cost for a required resolution. Our main goal in studying the well understood single impurity is to benchmark CPT/VCA for future application to the not so well understood case of multiband impurity models in various spatial geometries. This publication also sets the foundations for a future extension to non equilibrium problems. The text is organized as follows. The single impurity Anderson model is introduced in sec. II. A short review on CPT and VCA in this context is given in sec. III. A self consistent formulation of VCA previously introduced in the context of non equilibrium problems 50 is presented in sec. III 1. Some remarks about the choice of variational parameters are provided in sec. III 2. In sec. IV we discuss the grand potential Ω for infinite fermionic systems in relation with the VCA. Results for the single-particle dynamics of the SIAM are provided in sec. V. In this section also the quality of the low energy Kondo physics is compared to benchmarking results from NRG, DMRG, CT-QMC, Hartree-Fock and Bethe Ansatz calculations. Finally, we summarize and conclude our findings in sec. VI.

II. THE SINGLE IMPURITY ANDERSON MODEL
We consider the single impurity Anderson model (SIAM) 1 in real space, in one dimension H SIAM =Ĥ conduction +Ĥ impurity +Ĥ hybridization . (1) A tight binding band of non interacting s-electrons with nearest neighbor i, j hopping is described bŷ where ǫ s is the on-site energy of the particles, t is the overlap integral between nearest neighbor orbitals and i, j ∈ {1, . . . , N s } where N s is eventually taken to be infinity. The operators c † iσ and c iσ , respectively, create and annihilate electrons in orbital i with spin σ. The impurity Hamiltonian consists of a single f-orbital with local Coulomb repulsion U , with f † σ creating an electron with spin σ and on-site energy ǫ f located at the impurity. The particle number operator is defined asn f σ = f † σ f σ . Finally the coupling between a non interacting s-orbital and the impurity forbital is given bŷ where V is the hybridization matrix element between the s-and the f-orbital of the impurity atom (see fig. 1 for illustration). We have set the chemical potential µ to the center (ǫ s ) of the conduction electron density of states and choose µ = ǫ s = 0. The resonance width ∆ is defined as For the model defined in eq. (2) the local density of states of the conduction electrons ρ s (0) is given by ρ s (0) = 1 π t . In the forthcoming discussion we refer to the particle-hole symmetric case when we furthermore set ǫ f = − U 2 . All calculations are performed with t = 1 and V = 0.3162 which yields ∆ = 0.1. All reported results, except for the CT-QMC data in sec. V F, are for zero temperature. The impurity f-orbital is subjected to a local on-site energy ǫ f and local Coulomb interaction U and is hybridized with one of the s-orbitals (here the one at the beginning of the chain) via a hybridization matrix element V . This maps the impurity f-orbital onto site 0 and the impurity s-orbital onto site 1 in this geometry, the rest of the conduction band s-electron orbitals are mapped onto sites 2 to ∞. The semi-infinite non interacting chain is truncated at some site L. This decomposes the model into two clusters: an interacting cluster of variable size including the interacting impurity f-orbital and a semi-infinite chain of non interacting s-orbitals. In CPT/VCA these decomposed systems are coupled via a hopping element t.
function G σσ ′ ij (ω) which we calculate within cluster perturbation theory (CPT) 45,46 as well as the variational cluster approach (VCA) 47 . CPT and VCA have been previously applied inter alia to the fermionic Hubbard model and VCA also with great success to bosonic systems [51][52][53] with broken symmetry phases. The groundwork of VCA lies within cluster perturbation theory which is a cluster extension of strong coupling perturbation theory, valid to first order in the inter-cluster hopping. The main result of CPT is that the Green's function of the physical system G (which we call full Green's function throughout this text) may be obtained by a Dyson-like equation in matrix form Here g denotes the Green's function of a cluster which comes about by tiling the lattice of the physical system into smaller, numerically exactly solvable patches. This tiling is done by removing the hoppings between sites connecting such clusters. Therefore the matrix T = g −1 0 −G −1 0 contains all single-particle terms connecting clusters (i.e. the inter cluster hopping which will be referred to as T inter below). The subscript 0 denotes the non interacting Green's function. To apply this approach to the SIAM we start by splitting the physical model under consideration (eq. (1)) into appropriate pieces. Here we consider a cluster decomposition consisting of two parts. One part, consisting of a cluster of size L, which contains the interacting impurity f-orbital and a second, infinitely large part, the environment, which contains the rest of the conduction band The original Hamiltonian eq. (1), defined on the semiinfinite lattice, may now be rewritten aŝ Here T inter is the part of T describing the hopping from the interacting cluster to environment "cluster", which is the only term not included in the two clusters. For the SIAM the two bare Green's functions g interacting and g env (, which correspond toĤ interacting (eq. (7)) and H environment (eq. (8)), ) needed for eq. (6) may be evaluated separately. This is a bit different from the usual application of CPT to translationally invariant systems which normally leads to a single cluster having discrete spectra, embedded in a superlattice. Therefore the application of CPT to this problem does obviously not suffer from issues arising due to periodization prescriptions for the Green's function or self-energy 58 . We are dealing with two fundamentally different clusters, where one has a discrete (interacting cluster eq. (7)) and the other a continuous spectrum (environment eq. (8)). Due to the continuous spectrum of the environment a numerically favorable representation of the Green's function of the physical system G in terms of the Lehmann representation (see for example ref. 79) is not possible. For evaluating quantities from the Green's function G one therefore has to use a direct numerical integration, which works best on the Matsubara axis. The cluster Green's function g interacting is determined by exact diagonalization of eq. (7). We apply the Lanczos algorithm 55 to find the ground state and a Band Lanczos method to obtain the Green's function. The Band Lanczos method is initialized with a set of all annihilation and creation operators under consideration applied to the ground state. Thereby we obtain the so-called Qmatrices 56 which are used to calculate the Green's function Essentially this is the Lehmann representation for zero temperature Green's functions. The sum over α denotes a sum over a possibly d-fold degenerate set of ground states. The sum over γ is over a set of orthonormal basis-states having one particle more than the ground state (particle part) and one particle less than the ground state (hole part).
The Green's function of the environment g env is given analytically by the Green's function of a semi-infinite tight binding chain 57 where υ i,j is the retarded / advanced Green's function of the infinite tight binding chain if the infinitesimal imaginary part (0 + ) of ω is positive / negative. VCA, the variational extension of CPT, is based on the self-energy functional approach (SFA) 48,49 . In the SFA one considers the Legendre transformed Luttinger- , which is a universal functional of the self-energy, i.e. it does not depend on G 0 . F generates the Green's function, i.e.
where β denotes the inverse temperature. Introducing the (non-universal) self-energy functional Eq. (13) is an equation for the physical self-energy Σ given the Luttinger-Ward functional F [Σ] and the free Green's function G 0 . The trace Tr is short for Tr ≡ 1 β ωn tr, the sum is over fermionic Matsubara frequencies and the small form trace tr denotes a sum over lattice sites and spin. The idea is that, due to its universality, F [Σ] (and thus Ω[Σ, G 0 ]) can be evaluated exactly by exploiting a different system (so called "reference system") which differs from the physical system by single-particle terms only. This reference systemĤ ′ is simpler, and thus exactly solvable. It is defined on a cluster tiling of the original lattice and has the same interaction as the original systemĤ. The VCA reference system is chosen to be a cluster decomposition of the original lattice, as the one introduced for CPT above. Comparing eq. (12) for the full and the reference system yields where lowercase g denote Green's functions of the reference system. Thus the SFA/VCA approximation consists in solving eq. (13) in a restricted range of self-energies Σ, i.e. the ones produced by the reference system. In this way, the space of allowed Σ is spanned by the set of single-particle parameters of the reference system, x ′ . This means that the functional Ω[Σ, G 0 ] (eq. (14)) becomes a function of those parameters Ω(x ′ ) = Ω ′ (x ′ ) + Tr ln (−G(x ′ )) − Tr ln (−g(x ′ )) , (15) The stationarity condition determining the physical parameters eq. (13) is then given by The Green's function of the physical system is obtained by the CPT relation eq. (6). The matrix T = g −1 0 − G −1 0 (eq. (6)) in VCA contains all single-particle terms not included in the reference system (i.e. T inter ) as well as, in addition, the deviation introduced by VCA, ∆x ≡ x ′ − x of the single-particle parameters of the reference system x ′ with respect to the ones of the original system x.
In the following we fully adopt the zero temperature formalism, in which according expressions for the grand potential and related quantities may be readily evaluated.

Alternative: Self consistent VCA
In ref. 50 we explored an alternative version of VCA whereby the variational parameters x ′ were determined by a suitable self consistent criterion, instead of looking for the stationary point of the grand potential Ω eq. (16). This alternative approach was introduced to treat systems out of equilibrium, although it can equally be adopted in equilibrium. The advantage of this approach is that the solution of a self consistent equation is numerically easier than the search for a saddle point. The idea of this self consistent approach is to use a reference system which resembles the full system best. The strategy is to find those values x ′ for the set of parameters of the reference system which let the expectation values of their corresponding single-particle operators Ô cluster,x ′ coincide with those of the full system Ô CPT,x,x ′ . Here, the angle brackets denote expectation values in the cluster and the full system coupled by CPT or VCA respectively. Consider the on-site energies ǫ ′ f and ǫ ′ s as variational parameters. We will look for those cluster parameters ǫ ′ f and ǫ ′ s which fulfill the relations The sum is over all non interacting sites included in the cluster. This amounts to solving a system of non-linear equations in each step of the self consistency cycle. In general it is possible to vary each single-particle parameter individually. For reasons of keeping the numerics tractable we use one ǫ ′ s only, which we take to be the same for each orbital in the chain. Extension to a larger number of ǫ ′ s is straightforward. To fix this parameter we require the average particle density on the non interacting sites to fulfill the condition eq. (17). This corresponds to the condition presented in ref. 50 (see eq. (13) therein). In some situations, (see below,) we will alternatively consider the hybridization matrix element V ′ and the intra-cluster hopping t ′ as variational parameters, and proceed in an analogous way. Specifically, the particle number expectation values in eq. (17) are replaced by hopping expectation values. Again for t ′ we use a single variational parameter for hopping between all uncorrelated sites and fix it by requiring the mean value of hopping in the cluster and the full system to coincide. A discussion of this self consistency condition in connection with (cluster) DMFT is given in ref. 50. We use an improved multidimensional Newton-Raphson algorithm to find the roots of the system eq. (17). In some parameter regions no solution is found. A comparison between results obtained via the usual SFA, i.e. as stationary points of the grand potential Ω, which we will now refer to as VCA Ω , to the ones obtained by the above-mentioned self consistent condition (VCA SC ) will be given in the results section (sec. V).

Choice of variational parameters
In VCA one can, in principle, optimize all possible single-particle parameters which are present in the original model, as well as additional ones. By adding bath sites not present in the original model, one includes dynamical contributions to the cluster Green's function 48 . The numerical difficulty increases with the number of variational parameters. For the VCA SC case a multidimensional root finding algorithm has to be adopted. For the VCA Ω case, a saddle point in many dimensions has to be located. Since the allowed set of variational parameters limits the search space for the self-energies one will find a solution in this restricted space only. It is therefore desirable to vary as many single-particle quantities as possible. A balance has to be found between a large space of available self energies and numerically feasible multidimensional algorithms. Many works have addressed the question of which parameters are the most important to vary and how the choice of variational parameters will influence or limit the results 58 . As discussed in refs. 59,60, it is important to include an overall chemical potential as a variational parameter in order to preserve thermodynamic consistency. As a compromise, we will take two variational parameters x = {ǫ f , ǫ s }, which cover the overall chemical potential. Note that this amounts to shifting an overall on-site energy in the whole cluster plus an extra independent shift at the correlated site. For the variation of on-site energies we observe the grand potential Ω to be maximal at the stationary point which is in agreement with results for other models. Further parameters in the SIAM are the hopping t and the hybridization V . As discussed for example in ref. 61, the variation of hopping parameters is not straightforward. For the VCA Ω approach, we observe a maximum of Ω at ∆V = −V in the center of two symmetric stationary points. The two symmetrically lying minima are equivalent due to the fact that the self-energy is an even power of V . As one tunes the parameters away from particle-hole symmetry this stationary point is lost in the crossover region from the Kondo plateau to a doubly or unoccupied impurity (see sec. V C). In this parameter region the hopping t and hybridization V are probably not appropriate to be used as variational parameters within VCA Ω . In the following, we always choose the set x = {V } or x = {V, t} for calculations at particle-hole symmetry, which also includes {ǫ f , ǫ s }, since the variation of on-site energies will always yield zero deviations from the physical parameters and thus reproduce the CPT result here. For all other parameter regions it is sufficient to consider x = {ǫ f , ǫ s } as variational parameters.

IV. GRAND POTENTIAL FOR REFERENCE SYSTEMS OF INFINITE SIZE
The reference system consists of two parts, a finite interacting system and a non interacting system of infinite size, the environment. Ω ′ (x ′ ) is given by the sum of the grand potentials of the interacting cluster (Ω ′ interacting ) and of the environment (Ω ′ env ) (, which correspond tô H interacting (eq. (7)) andĤ environment (eq. (8))). Here we outline how to determine the grand potential for such kinds of reference systems. For the Green's function G within the CPT/VCA approximation the Dyson equation is given in eq. (6). The Green's function and T have the block structure Up to this point all matrices involving environment indices have infinite size. As far as the Green's function itself is concerned this is no problem as we are primarily interested in G cc for which the Dyson equation reduces to G cc = g cc + g cc T cc G cc + g cc T ce G ec G ec = g ee T ec G cc , and therefore G cc = g cc + g ccΣcc G cc Σ cc := T cc + T ce g ee T ec .
A bit more tedious is the elimination of the explicit dependence on the environment part of G, as far as the grand potential eq. (15) is concerned. We start out from a form of the grand potential functional given by Sénéchal 58 In app. A it is shown that ∆Ω can be expressed solely in terms of cluster quantities ∆Ω = −Tr ln 11 cc −Σ cc g cc .
Along the lines outlined in ref. 58, the resulting integral can be regularized and expressed as The quantities Ω ′ are the grand potentials of the uncoupled reference system. The constant infinite contribution Ω ′ env is absorbed into the definition of Ω. It plays no further role as it does not depend on the variational parameters. This integral may be evaluated as suggested in ref. 58 by integrating from 0 to Λ 1 , from Λ 1 to Λ 2 and from Λ 2 to ∞. Λ 1 and Λ 2 represent two characteristic scales in the problem (for example the smallest/largest eigenvalue of the Hamiltonian matrix). For the last part of the integral a substitutionω = 1 ω is performed. We use an adaptive Gauss Legendre integrator for the evaluation.

V. RESULTS
We have evaluated several benchmarking dynamic quantities of the SIAM. In the following, results for the impurity density of states will be presented and compared to ED, NRG and DMRG data. We will elaborate on the strengths and weaknesses of the methods as well as the comparison of CPT to VCA. Furthermore, we will discuss the relation between VCA SC , where the variational parameters are determined self consistently via eq. (17) and VCA Ω , where the variational parameters are defined at the stationary point of the grand potential. We will show that the Kondo resonance is reproduced within the framework of CPT/VCA and that the variational results fulfill certain analytic relations like the Friedel sum rule (eq. (23)). The method will be shown to provide reasonably accurate results in a wide range of parameter regimes of the model. Low energy properties related to the Kondo temperature T K will be discussed in context with renormalization group results. The imaginary frequency Green's function and self-energy will be compared to CT-QMC results.
A. Even-Odd Effect -choice of the impurity position CPT/VCA rely on the Green's function of an interacting cluster of size L which is obtained by exact diagonalization. Due to this fact, it is unavoidable that some effects of the finite size cluster influence the solution of the full system. (Except in the case of vanishing interaction strength, i.e. U = 0.) Therefore suitable clusters have to be chosen on a basis of physical results. Some aspects of this are discussed by Balzer et al. 59 in the context of DMFT and VCA and by Hand et al. 62 in the context of DMRG. In this work we consider interacting clusters of even size only. For these systems the ground state does in general not suffer from spin degeneracy. Furthermore, the spatial position of the impurity is important. This can be inferred from the bath's density of states, which vanishes for ω = 0 at every second site. It may also be seen in the structure of the ground state, for which the size of the degenerate sectors alternates with the geometrical size of the cluster. Throughout this work we position the impurity f-orbital at the beginning of the infinite chain, although essentially the same results are achieved by attaching it to an s-orbital at site two, four, etc. inside the chain.

B. Spectral properties
The single-particle spectral function A σ ii is obtained from the retarded Green's function G σ,ret ii (see e.g.

ref. 63)
The diagonal element at the impurity f-orbital A σ f f (ω) describes the impurity density of states ρ σ f (ω). A physical property of the SIAM which poses a challenge to numerical methods is the Kondo-Abrikosov-Suhl resonance often referred to as Kondo peak 4 . It arises in the parameter regime where the magnetic moment of the impurity is screened by the conduction electrons to form a singlet state 64 . The particle-hole symmetric model lies in the center of this Kondo region. This quasi particle excitation is for example not captured in mean field approaches. With increasing interaction strength U the numerical solution becomes increasingly challenging. In this section we elaborate on the results for the density of states in the particle-hole symmetric case. A comparison of the local single-particle spectral function at the impurity f-orbital as obtained by exact diagonalization and VCA Ω is shown in fig. 2. The ED result is for a ten site system with open boundary conditions. The VCA Ω result is for an infinite reference system, where the interacting part of the reference system was taken to be of size L = 10. An ED treatment Comparison of the local single-particle spectral function at the impurity f-orbital at particle-hole symmetry as obtained by exact diagonalization (ED) of a ten site chain (black) and VCAΩ (cyan). VCAΩ was used with two variational parameters: the hopping t and the hybridization V considering a length of the interacting part of the reference system of L = 10. All data shown is for an interaction strength of U/∆ = 12. All results have been obtained for a large numerical broadening 0 + = 0.05. The inset shows a zoom to the low energy region.
of a finite-size SIAM can not reproduce the low energy resonance at zero energy in the single-particle spectral function (see also app. B). It consists, in the particle-hole symmetric case, for an even number of sites (and open boundary conditions), of symmetrically lying excitations which shift closer to zero for increasing system size and represent a large energy scale. For an odd number of sites a pole in the local single-particle Green's function of the impurity f-orbtial will be present at ω = 0. CPT as well as VCA are able to reproduce finite spectral weight at ω = 0 (even for 0 + → 0), since these methods are formulated for an infinite system. The finite-size structure in the high energy incoherent part of the spectrum, owing from the excitations of the interacting part of the reference system, is strongly reduced in VCA Ω .
Results for the single-particle spectral function eq. (21) of the impurity f-orbital are shown in fig. 3 for four different interaction strengths U/∆ = 4, 8, 12 and 20. As a reference, the spectra obtained with NRG and DMRG from Peters 34 are plotted. Renormalization group approaches like NRG are especially suited to reproduce the low energy quasi particle excitations of this model and therefore serve as a reference for our data. The spectra of Peters were obtained for a flat conduction electron density of states, which was mapped by linear discretization in energy space onto the corresponding orbitals of a semi-infinite chain. Our model is based on a semi-circular density of states of the conduction electrons. The low energy part of the spectra is comparable because we have chosen the only relevant parameter for the low energy part of the spectrum: ∆ accordingly. This parameter fully determines the influence of the conduction electrons on the impurity f-orbital for low energies and therefore the low energy part of the spectrum. The high energy part of the spectrum may deviate slightly and is not directly comparable but yields a crude reference. In addition, we have chosen here a very large numerical broadening of 0 + = 0.05 for reasons of comparison only. This value was used in the DMRG calculations and is needed there to obtain spectra using the correction vector method. This influences the width and the height of the Kondo resonance, located at ω = 0. The CPT spectral weight at ω = 0 appears too broad in the plot in comparison with the NRG result. This is only partly due to a large numerical broadening. Due to the nature of the CPT method we cannot expect it to reproduce the low energy spectrum as well as RG calculations do. The height of the Kondo resonance appears too small in this figure because of the large 0 + . It converges with 0 + → 10 −6 to the result predicted by scattering theory (see fig. 5 and sec. V C). The high energy incoherent parts of the spectrum located at ω ≈ −ǫ f and ω ≈ −ǫ f + U develop more and more with increasing length of the interacting part of the reference system L. A comparison of the center of gravity of the high energy incoherent part of the spectrum of the L = 14 site CPT result and the fifty site DMRG result is in reasonable agreement. There are spurious structures in the spectral density, originating from the cluster Green's function of the finite system, preventing continuous spectra to form. We would like to note that the accurate determination of the Green's function of the reference system is of prime importance. An inaccuracy in pole-positions or pole-weights for very small but non-vanishing weights will yield spurious artifacts in the spectra in the vicinity of ω = 0. To improve on the result of CPT we considered the hopping matrix element t and the hybridization matrix element V as variational parameters. The parameters used for the evaluation of the reference system were determined with two different methods. VCA Ω results are depicted in the plot for a length of the interacting cluster of L = 10. As shown in the figure this method strongly reduces the finite size peaks in the high energy incoherent part of the spectrum. The width of that part of the spectrum is reproduced correctly for high values of U where the full width at half maximum (FWHM) within VCA is given by ≈ 1.9∆. This comes very close to the expected 2∆ 21,65 of the high energy atomic excitations. VCA improves the spectral properties of the Kondo resonance with respect to CPT, bringing it closer to the fifty site DMRG result. The data obtained using the self consistent VCA approach VCA SC agree very well with the result based on VCA Ω on the position of the spectral features. The respective weight however disagrees for low values of interaction strength U , which is due to the different values predicted for the variational parameters by the two procedures. One should note that the two broad Lorentzian high energy peaks (in VCA Ω as well as VCA SC ) consist of many excitations which will be revealed upon repeating this calculation with smaller 0 + .
A more detailed look on the spectral region of the Kondo resonance is provided in fig. 5. The CPT/VCA Ω data is compared to NRG and FRG data as well as results obtained from a restricted Hartree-Fock calculation from Karrasch et al. 31 . The CPT/VCA results are plotted for lengths of the interacting part of the reference system L = 2, 4, 6, 8 and 10 for two different sets of parameters. The results for higher L are always located towards the center of the figure. The results corresponding to the resonance at ω = 0 were obtained for the particle-hole symmetric model. For this set of parameters we used the hybridization V as a variational parameter. The second peak shown centered around ω/∆ ≈ 0.8 corresponds to a parameter set right at the border of the Kondo region. The variational parameters used away from particle-hole symmetry are x = {ǫ f , ǫ s }.
One can see that the CPT result is not converged for L = 10 site interacting clusters yet. In contrast, the VCA Ω result seems to converge much faster. Although in the plot it looks like the VCA result does not improve much upon a restricted Hartree Fock calculation, we will show in the following that CPT/VCA yields results in all parameter regimes of the SIAM which cannot be reproduced within a mean field treatment.
The variational parameters obtained for the two sets of parameters shown in fig. 5 are presented in fig. 6. In addition to the VCA Ω parameters, which were used for the results above, the variational parameters obtained in VCA SC are also depicted. We plotted the difference of the parameter of the reference system x ′ to the physical parameter x: ∆x. All deviations ∆x appear to converge to zero with increasing length L of the interacting part of the reference system. Notice that the self consistent approach always leads to a ∆x of greater magnitude with respect to VCA Ω . Remarkably, the spectra obtained by VCA Ω and VCA SC for the parameter set x = {ǫ f , ǫ s } are in very good agreement even though the variational parameters are rather different. The most striking difference is that the self consistent approach yields a negative ∆ǫ f while the Ω based VCA yields a positive ∆ǫ f . This is however compensated by the different ∆ǫ s . Using the hybridization V as a variational parameter, the ∆V obtained by VCA Ω and VCA SC agree rather well. Remarkably, the resulting density of states is very different, which shows that the calculation is extremely sensitive to this parameter.
A spatially resolved image of the spectral function, calculated with CPT, for the parameter set used in fig. 3 (c) is shown in fig. 4. The qualitative picture would be the same in VCA, merely the structures are slightly shifted. This view reveals how the perturbation, introduced by the impurity, is fading away slowly in an alternating fashion. At every second site away from the impurity a dip at ω = 0 is present, which is usually referred to as Fano dip. The low computational effort of CPT/VCA proofs advantageous for calculating spectra. The VCA procedure (for a twelve site interacting cluster) usually converges in minutes to hours on a standard workstation PC, while more demanding numerical methods often need days to a week to converge. Furthermore, the spectra are exactly determined from the Lehmann representation and no ill-posed analytical continuation is required in comparison to methods working in imaginary time or imaginary frequency space. To our knowledge the most accurate spectra available for this model so far are published in ref. 66. Overall one may conclude that CPT, VCA Ω and VCA SC reproduce a Kondo resonance, which fulfills the Friedel sum rule (eq. (23)) for 0 + → 10 −6 . The VCA results improve drastically upon the CPT data which may be seen in a much faster convergence in L and a suppression of finite size effects especially in the high energy part of the spectrum which in addition has the expected width within VCA. VCA Ω and VCA SC agree rather well on the position of the spectral features. However they assign very different spectral weight to them at low values of interaction strength U .

C. Impurity density of states and occupation
The occupation of the impurity f-orbital is given at T = 0 by This integral may be evaluated from the imaginary frequency Green's function, which in turn is directly accessible within CPT/VCA. To see whether CPT/VCA are good approximations in all parameter regions of the SIAM, we vary the on-site energy of the impurity ǫ f at fixed interaction strength U . The local impurity density of states at the chemical potential (ω = µ = 0) and the impurity occupation number are plotted for various lengths of the interacting part of the reference system L = 2, 4, 6 and 8 for the same model parameters. The VCA Ω result is shown in fig. 7, a VCA SC calculation in fig. 8 and the CPT data in fig. 9. We start out by discussing the VCA Ω result ( fig. 7). The variational parameters x used within VCA Ω are the onsite energy of the impurity ǫ f and the on-site energies of the uncorrelated cluster sites ǫ s . The density of states ρ f (0) displays a pronounced plateau which is related to the existence of a quasi particle peak (Kondo resonance) pinned at the chemical potential. The parameter regions leading to an empty (−ǫ f < 0) or to a doubly occupied (−ǫ f > U ) impurity do not show a pinning of the Kondo resonance at the Fermi energy, as expected. In the half filled region which lies in between, virtual spin fluctuations lead to a pronounced quasi particle peak at the chemical potential. We observe that the result converges with increasing length of the interacting part of the reference system L to the physically expected result. Due to the variational parameters considered, the deviations of the results as a function of L are rather small as compared to CPT where the results change significantly with increasing size of the reference system (see fig. 9). We expect CPT calculations in the empty or doubly occupied regions to converge rather fast (within a few sites) while calculations in the Kondo regime, and particularly in the crossover region, may fully converge only at very large (i.e. exponentially) sizes of the reference system 62 . This is inferred from the spin-spin correlation function in the cluster which is observed to decay sufficiently fast outside the Kondo plateau (i.e. it is effectively zero at the boundary of the cluster) but shows long range correlations inside the plateau. The VCA SC results are obtained with one variational parameter x = {ǫ f }. The reason for not using x = {ǫ f , ǫ s } again is that the result is almost the same as the one obtained with VCA Ω (see fig. 7). However in some (small) parameter regions the numerical evaluation becomes difficult. The VCA SC data shown in fig. 8 shows a clear improvement as compared to CPT but does not reach the quality of the VCA Ω result in terms of convergence in system size. The Friedel sum rule (FSR) 5,67,68 provides an exact relation between the extra states induced below the Fermi energy by a scattering center and the scattering phase shift. It also holds true for interacting systems. This gives a relation between the f-orbital occupation n f , and the density of states at the Fermi energy: In our case the mean occupation in the Kondo regime is n f ≈ 1. Since both, the occupancy of the f-orbital and the magnitude of the local density of states at the Fermi energy, can be evaluated independently, we can check the validity of the Friedel sum rule in our approximation. Results are shown in fig. 7 applied to the L = 8 site VCA Ω data. The VCA Ω results fulfill the Friedel sum rule almost in the whole Kondo region. At the crossover to an empty or doubly occupied impurity the Friedel sum rule is not fulfilled exactly any more but approximated very well. Further outside the agreement is again excellent. The variational parameters of VCA are crucial to fulfill the Friedel sum rule as can be seen from a CPT calculation ( fig. 9) which violates it in all parameter regions. It appears that VCA Ω with variational parameters x = {ǫ f , ǫ s } naturally drives system L are shown. The CPT results are by far not converged for the interacting cluster sizes considered here. This emphasizes the importance of the variational parameters.
Our results in fig. 7 and fig. 10 agree very well with those of calculations based on X-operator technique exercised by Lobo et al. 38 . In their work a strong coupling perturbation theory is applied starting from the Anderson molecule as a basis and using the Friedel sum rule as a condition to fix the position of an infinitely narrow conduction band. Analytic considerations (see app. B) allow insight into the behavior of the Friedel sum rule in ED, CPT and VCA. There it is shown that ED always has to violate the Friedel sum rule while CPT always fulfills it in the particle-hole symmetric case. This comes about in the first place because the height of the Kondo resonance at ω = 0 does not depend on the self-energy. A pinning of the Kondo resonance however can only be achieved via the improved self-energy contributions obtained within VCA. The results of this section clearly show that VCA is able to capture the basic physics of the SIAM in every parameter region. The improvement obtained by going over from CPT to VCA is crucial to fulfill exact analytic relations. Moreover we have shown that CPT/VCA is clearly superior to ED calculations. VCA is, even at small L, capable of fulfilling the FSR also away from particle-hole symmetry due to a pinning of the Kondo resonance at the Fermi energy. This pinning can be The numerical broadening used is 0 + = 10 −6 . The set of single-particle parameters considered for variation within VCAΩ is x = {ǫ f , ǫs}. Note that here the point ǫ f = − U 2 corresponds to the particle-hole symmetric case. The Friedel sum rule (eq. (23)) was applied to the L = 8 result (dotted-violet). It is fulfilled to a very good approximation in the Kondo region and far outside of it. Small deviations from the Friedel sum rule arise at the crossover region to an empty or doubly occupied impurity. The inset shows a zoom to the Kondo plateau.
attributed to the better approximation of the self-energy of VCA with respect to CPT.

D. Crossover diagram
To delve into the CPT/VCA results for the whole parameter range of the SIAM a "phase diagram" is presented in this section. This should be understood to be a mere scan of the parameters U and ǫ f because the model does not undergo a phase transition. The density of states of the impurity at the chemical potential ρ f (0) is shown in fig. 11 in a density plot. This figure essentially shows the height of the Kondo resonance as a function of interaction strength and on-site energy of the impurity. The different regimes of the SIAM, as obtained by an atomic limit calculation, are indicated as black lines. These lines divide the physics into regions where the impurity is doubly, singly or not occupied. In the singly occupied region ( U 2 > |ǫ f + U 2 |) local moments and their screening are expected to appear. This region, which bestrides the cone enclosed by black lines, is the region where Kondo physics may take place within this approximation. The parameter regions where the impurity is empty or doubly occupied lie above and below this cone. More sophisticated methods will lead to a smearing out of the border of these regions and introduce a crossover area with competing effects. A boundary expected between a single resonance and a spurious local moment behavior where the single resonance is split into two for spin up and spin down respectively is obtained by mean field theory 64 . In the mean field approach the density-density interaction of the impurity Hamiltonian is replaced by a spin dependent density exposed to the mean contribution of the other spins' density. The mean field boundary is then obtained by replacing the mean field parameters n f ↑ and n f ↓ by the particle number n and magnetization m. Setting the magnetization to m = 0 + in the self consistent equations yields the implicit result U c = π∆(1 + cot(π n f 2 ) 2 ) and ǫ f,c = ∆(cot(π ). The plot shows that the Kondo plateau is reproduced very well by VCA Ω . The results appear almost converged for lengths L ≈ 6 of the interacting part of the reference system. Increasing L yields better results in the crossover region. Results obtained by means of CPT do not reproduce the Kondo plateau very well for small L.
The average impurity occupation for the same pa-  rameter region is shown in fig. 12. The result obtained with VCA Ω clearly shows the Kondo plateau where the impurity is singly occupied. The parameter regions of a doubly occupied or empty impurity lead to a density of states in the impurity which is zero at the chemical potential (compare to fig. 11).
The results of this section have been obtained using VCA Ω with variational parameters x = {ǫ f , ǫ s }. It should be noted that using only x = {ǫ f } already yields good results. As mentioned in sec. V C, CPT needs very large values L to yield the same quality of the results as VCA does with much smaller values of L.

E. Low energy properties, Kondo Temperature
In this section we examine the low energy properties of the symmetric SIAM. In the strong coupling limit a single scale, the Kondo temperature T K , governs the low energy physics 5 . The Kondo temperature T K is known from Bethe Ansatz results for the particle-hole symmetric SIAM 69,70 This scale, which is inversely proportional to the spinflip rate of the impurity, divides the physics of the SIAM into two regions: A local moment behavior of the impurity, where the spin is free, and a low temperature region where the local spin and the conduction electrons become entangled and form a singlet state 71 . Quantities which depend inversely on T K are the effective mass m * and the static spin susceptibility χ m . The Kondo temperature may furthermore be extracted from the width or weight of the Kondo resonance in the local density of states of the impurity f-orbital. We investigate and compare the results for the scale T K obtained from the direct determination of T K (from the FWHM and the spectral weight of the Kondo resonance) and the inverse quantities m * and χ m . We find that the results of all four measurements turn out to yield the correct qualitative behavior in VCA Ω . However, in a region where the dependence of T K is exponentially dependent on the interaction strength U the exponential prefactor is not predicted correctly. Therefore we introduce a scaling factor γ (eq. (24)) which turns out to be the same for all four ways of determining T K . In particular this factor is independent of the set of model parameters used. The scaling factor may be calculated semi-analytically for a reference system consisting of a two site interacting cluster and the semi-infinite environment within VCA Ω and VCA SC (x = {V }). The calculation for VCA Ω leads to an integral expression for the stationary point of the grand potential Ω with respect to ∆V from which the optimal ∆V can be obtained numerically (see app. C). The Kondo scale may be determined from the so obtained values of V ′ (U ) = ∆V (U ) + V by This leads to a perfect exponential behavior as defined in eq. (24) with γ = 0.6511 .
The issue of obtaining an exponential scale but not the correct exponent for the functional dependence on U is common to various approximate methods (for example variational wave functions where the issue was cured by introducing an extended Ansatz by Schönhammer 72 , saddle-point approximations of a functional integral approach 73 or FRG 74 ). A faint analogy may be drawn here to Gutzwiller approximation, where an exponential energy scale in U arises by a renormalized hybridization parameter V 75 , which is also the case for VCA Ω . The self consistent calculation for VCA SC also leads to an integral expression for the determination of ∆V . This expression is obtained by requiring the expectation values of the hopping from the impurity f-orbital to the neighboring site in the reference system to be the same as the expectation value in the physical system. This procedure does not yield an exponential scale in U . The optimal cluster parameter V ′ shows spurious behavior as a function of U . We conclude that VCA SC with x = {V } cannot reproduce the low energy properties of the SIAM even qualitatively, while VCA Ω yields the correct behavior apart from an exponential factor.

Effective mass -quasi particle renormalization
The effective mass m * is defined as the quasi particle renormalization 31 where we introduced the dependence on the interaction strength U explicitly. In the Kondo regime, this quantity becomes inversely proportional to the Kondo temperature.
We want to answer the question whether the Kondo scale is approximately captured by CPT/VCA or not. Therefore we compare the functional form and the exponent obtained from the effective mass and the analytic result for T K , eq. (24). The result for the effective mass obtained within VCA Ω is shown in fig. 13. The variational parameter used was x = {V }. The functional form is reproduced well by VCA Ω (i.e. it starts out quadratically and goes over to an exponential behavior in the Kondo region). However the exponent ( π 8∆ ) is not reproduced correctly. VCA Ω yields a lower exponent of ≈ (γ π 8∆ ). The factor γ is defined in app. C, determined from a semi-analytical calculation of T K within VCA Ω . This additional factor is the same for all initial parameters (within the Kondo regime), it is particularly independent of ∆. If plotted over a scaled U axis, U ′ = 1 γ U , the VCA Ω result would lie on top of the NRG data. However using larger sizes of the interacting part of the reference system does not lead to much better results regarding γ. It is to be expected, that a significant improvement can only be obtained using exponentially large L. The CPT result shows a very different convergence behavior in L which is rather slow. An attempt was made to extrapolate the CPT data to L → ∞ by a simple 1 L scaling. It is interesting to observe that this extrapolated curve coincides nicely with the VCA result (L = 6) in the low U region. However, we expect this extrapolation based on small L to be insufficient to capture the exponential scaling of the data in L. Note that since in CPT the self-energy is taken from the cluster, the CPT results for the effective mass coincide with ED results for systems of size L. (area) of the peak are proportional to the Kondo temperature T K . Obtaining the spectral weight or FWHM of the Kondo resonance from the spectrum introduces a large uncertainty. Nevertheless we made an attempt, to get an idea of the behavior of T K . We fixed the spectral weight by the first minimum to the left and to the right of the central peak (see also ref. 62). In general the effective mass and static spin susceptibility will yield more reliable results but it is instructive to compare these four ways of determining T K . Shown in fig. 14 is the evolution of the spectral weight and the FWHM of the Kondo resonance with increasing interaction strength U . The data was acquired using VCA Ω with a variational parameter x = {V } for the particlehole symmetric SIAM. Within the uncertainty, the same exponential behavior for the Kondo temperature T K is obtained as by calculating the effective mass in VCA Ω .

Static spin susceptibility
The static spin susceptibility χ m is given by the linear response to an applied magnetic field B in z direction In the Kondo regime this quantity too becomes inversely proportional to the Kondo temperature. For the calcu- lations in this section we introduce an additional spin dependent term in the impurity Hamiltonian (eq. (3)) The static spin susceptibility χ m as obtained with VCA Ω is shown in fig. 15. The variational parameter used was x = {V }. As a reference the results of NRG and FRG 31 are shown. The behavior of the VCA result is good for small interaction strength U . The VCA result shown for L = 6 appears already converged while the CPT result would require much larger systems. We would like to highlight that VCA Ω reproduces an energy scale T K . Results from direct calculation of T K , calculation of the effective mass m * and the static spin susceptibility χ m yield the correct functional form but not the right exponent. In this section we compare CPT/VCA results to QMC data. We obtained the Monte Carlo results using the continuous time Quantum Monte Carlo (CT-QMC) code of the TRIQS 76 toolkit and its implementation of the hybridization expansion (CT-HYB) 77 algorithm using Legendre polynomials 78 . This method enables access to very low temperatures and is especially suited to obtain low energy properties 40 . The CT-QMC data provides statistically exact and reliable results to test our data. All CT-QMC calculations were done for a single impurity orbital at U = 0.8 and ǫ f = −0.4. We used a semicircular hybridization function with half bandwidth D = 2 and V = 0.3162. This setup corresponds to the same model under investigation here. The value for the interaction strength U = 0.8 was chosen because of the relatively low expected Kondo temperature of β K = T −1 K ≈ 100. For all calculations 1.2 × 10 9 MC updates where conducted, with a sweep size of 100 updates, plus a 10% thermalization period. To ensure that the Kondo resonance is correctly reproduced by CT-QMC we evaluated the Matsubara Green's function for various values of inverse temperature β. The height of the Kondo resonance is given by the Friedel sum rule eq. (23) to be ℑm(G f f (iω n = 0)) = −10 for the parameters used here (∆ = 0.1). To obtain ℑm(G f f (iω n = 0)) we extrapolate twice, first in iω n → 0 for each β, then we use these results and extrapolate to T → 0. The extrapolation to iω n → 0 is done linearly using the first two Matsubara frequencies. The imaginary part of G f f (iω n ) and the extrapolated value to iω n → 0 are shown in the inset of fig. 16 for β ∈ [10,1200]. Those extrapolated values are plotted as a function of temperature ( fig. 16). These data points are then extrapolated to T → 0 using a fit by a rational model function. The result clearly shows the onset of the Kondo resonance when the temperature is lowered below the Kondo temperature T K . The extrapolation to T = 0 shows very good agreement (ℑm(G f f (iω n = 0)) ≈ −10.1) with the result expected from the Friedel sum rule within the uncertainty. It is important to note that the CT-QMC results converge very nicely in β. Although for higher β lower Matsubara frequencies become available, the overall shape of the Green's function does not change significantly. Therefore we may compare the T = 0 CPT/VCA results for the Green's function and self-energy to the CT-QMC data. The Matsubara Green's functions of the impurity f-orbital G f f (iω n ) obtained by CT-QMC (β = 400), CPT and VCA are shown in fig. 17. We use β = 400 as a compromise between low temperatures and still reliable CT-QMC results (within manageable computation time). The β = 400 result was obtained using 65 Legendre coefficients. A detailed analysis has shown that this number is sufficient to get high frequency moments of the self-energy Σ accurately. The VCA Ω results were obtained with one variational parameter x = {V } for U = 0.8, ∆ = 0.1 and 0 + = 10 −6 in the particlehole symmetric case. For the CPT calculation we used the same parameters. For both methods we considered lengths of the interacting part of the reference system of L = 2, 4, 6, 8 and 10. The VCA result lies near the CT-QMC data but underestimates the slope of the curve at low iω n . The VCA result provides a huge improvement upon CPT for the lengths of the interacting part of the reference system shown here. The real part of G f f (iω n ) is exactly zero within CPT/VCA as it is supposed to be. Note that the value of G f f (iω n = 0) which is fixed by the Friedel sum rule is exactly reproduced within CPT and VCA for the particle-hole symmetric case. The same is shown for the self-energy of the impurity f-orbital Σ f f (iω n ) in fig. 19. From the imaginary part of Σ f f (iω n ) one can infer the convergence of the CPT/VCA result with larger length of the interacting part of the reference system L. The real part of the selfenergy (ℜe(Σ f f (iω n ) = µ = −e f = U 2 = 0.4) is again exactly reproduced within CPT/VCA. In the following, we discuss the self-energy Σ(iω n ) for the two interesting cases of very low and very high Matsubara frequency. We start out by conducting an expansion of the self-energy Σ(z) for high Matsubara frequencies (z = iω n → ∞) which shall be outlined here briefly. The self-energy matrix is defined by Here T is the one-particle part of the Hamiltonian. In the particle-hole symmetric case considered here it contains all the hoppings as well as the on-site energy of the impurity ǫ f = − U 2 . We conduct a series expansion in powers of z −1 of Σ(z). Apart from the real constant T ii all z-dependent terms of Σ ii (z) are anti-symmetric in z. Therefore even powers in z ±2l , l > 0 vanish. Expanding the Green's function G(z) yields for the self-energy Σ(z) where ∆Ĥ =Ĥ − ω 0 and ω 0 is the ground-state energy ofĤ. Collecting powers of z yields a cumulant-like expansion for the self-energy Σ(z) Here we consider the zeroth and first order in z −1 only and obtain for Σ(iω n ) where the self-energy at the impurity f-orbital Σ f f is the only non-vanishing matrix element of Σ ij . This result is plotted as a reference in fig. 19. Due to the nature of the CPT/VCA approximation these methods always yield the exact self-energy for high Matsubara frequency as shown in fig. 18.
The low energy properties examined in the previous section depend basically on the slope of the Matsubara Green's function at (iω n ) = 0 + . The results shown in fig. 17 and fig. 19 show that this slope is underestimated by CPT/VCA in comparison to CT-QMC, at least at the small lengths of the interacting part of the reference system available. The above results suggest a possible application of VCA as an impurity solver for zero temperature DMFT. The results would not suffer from a bath truncation error as in exact diagonalization based DMFT. A big advantage would be the low demand on computational power of VCA as well as the approximate reproduction of the main features of the local density of states (i.e. Kondo resonance and high energy incoherent part of the spectrum).

G. Introducing a symmetry breaking field
We explore the possibility to improve the VCA results achieved by varying the internal single-particle parameters of the model by introducing a symmetry breaking "spin flip field" at the impurity f-orbital. The term added to the impurity Hamiltonian eq. (3) explicitly breaks the conservation of spin in the cluster solution. We are interested in the model with a physical parameter B x = 0 so this variable may only attain a finite value as a variational parameter B ′ x in the reference system. We investigate the particle-hole symmetric model at V = 0.3162 and t = 1. Our findings indicate that any finite value of B ′ x splits the Kondo resonance and has thus to be discarded on physical grounds for the system under investigation. While this prevents the application of this field to improve the VCA results, it gives very nice insight in the physics of the SIAM as described by CPT/VCA. We find that a critical interaction strength U c depending on the length of the interacting part of the reference system exists which separates solutions which would prefer a finite B ′ x from those which would prefer B ′ x = 0. The critical interaction strength for L = 4 is given by U c /∆ ≈ 4.3. The grand potential Ω−Ω ′ 0,env is plotted for various interaction strengths U in fig. 20. For an analogous calculation for L = 6 site interacting clusters a value of U c /∆ ≈ 4.1 is achieved. The mean field result would yield a critical interaction strength U c /∆ = π for the parameters used here. We interpret this value as a signature of the onset of local moment behavior. The values for U c are of course not to be taken literally, they depend very much on the finite size of the cluster under investigation. The splitting of the Kondo resonance caused by a non-  20)) as a function of the interaction strength U/∆ (see legend). The data was obtained by studying a L = 4 site interacting cluster coupled to a semi-infinite lead. The numerical broadening used was 0 + = 10 −6 . The crosses indicate the respective minimum of the grand potential. There exists a critical Uc/∆ ≈ 4.3 above which a finite B ′ x is preferred by the system.

zero variational field B ′
x is shown in fig. 21. The value of U/∆ = 12 used for this calculation lies in the region above U c where the system prefers a nonzero field B ′ x .

VI. CONCLUSIONS
In this work we have applied the variational cluster approach (VCA) to the single impurity Anderson model. We devised a cluster tiling applicable to this nontranslationally invariant model which leads to a cluster with a discrete spectrum and an environment having a continuous spectrum. We have derived an expression for the change of the grand potential originating from the coupling of the impurity to the semi-infinite bath. We have compared results for the single-particle dynamics to data obtained by exact diagonalization and cluster perturbation theory (CPT). We found that the variational extension made by the VCA is vital for a good reproduction of the expected behavior of the SIAM. The CPT/VCA spectra both yield a Kondo resonance in the impurity density of states with the correct height as predicted by the Friedel sum rule. A close look at the Kondo resonance shows that the VCA is able to reproduce the resonance and the functional behavior for the Kondo temperature in a remarkable way. The Kondo temperature is expected to show exponential behavior in interaction strength in the Kondo regime. VCA yielding an exponential behavior however tends to underestimate the exponent. Comparison of dynamic x . The plots were obtained using VCA (i.e. the physical field Bx is always zero). Instead of taking the parameter B ′ x at the stationary point of the grand potential (this value would be B ′ x /∆ ≈ 1.9 for the parameters used) we explicitly plug in a fixed value for B ′ x . The length of the interacting part of the reference system used was L = 6 for the model parameters U/∆ = 12. The numerical broadening used was 0 + = 10 −6 . quantities to continuous time Quantum Monte Carlo solidifies the origins of this behavior. The high energy incoherent part of the spectrum shows strong finite size effects within CPT which are partly removed by virtue of the VCA. VCA furthermore reproduces the expected position and width of the high energy part of the spectrum. For the asymmetric model, the Friedel sum rule is fulfilled in all parameter regions implying that the Kondo resonance is pinned at the chemical potential in the Kondo region. In addition a self consistent formulation of the VCA, previously introduced in the context of non equilibrium problems 50 , was explored. Results obtained by the self consistent approach show agreement with results obtained by VCA based on the grand potential for the density of states of the impurity f-orbital. Thereby the positions of the spectral features agree very well with the traditional VCA result, while the spectral weight distribution may deviate especially for small values of interaction strength. Comparison to results obtained from Bethe Ansatz, renormalization group approaches and data obtained from X-Operator based calculations show reasonable agreement for all quantities investigated. In conclusion, while there are certainly more accurate methods to deal with a single quantum impurity model 23 , especially at low energies, our work shows that VCA is a flexible and versatile method which provides reasonably accurate results with modest computational resources. Here, the VCA self consistency condition proves to be crucial. This allows to obtain the same accuracy that CPT would provide with a much larger, inaccessible, interacting cluster size. One of the advantages is the flexibility of the method, i.e. it is straightforward to extend it to many impurities, non equilibrium problems 50 , etc. In the spirit of NRG 29 one could improve on the present results by carrying out an appropriate unitary transformation on the bath such that the bath hoppings decay with increasing distance from the impurity. In this work, we do not include such an improvement and choose a constant hopping sequence (eq. (2)), since our goal is to benchmark VCA/CPT only. Nevertheless, a hybrid approach combining NRG and VCA would be an interesting extension of the present work.
are no matrices T ce and T ec respectively, left. Hence we can as well introduce in eq. (18) the following replacements g ee →g cc , T ce → 11 cc , T ec → 11 cc , 11 → 11 cc , as it leads in the series expansion to the same expressions. The argument in eq. (18) then assumes the form 11 − Tg = 11 − T cc 11 cc 11 cc 0 cc g cc 0 cc 0 ccgcc = b −g cc −g cc 11 , with the abbreviation b := 11 − T cc g cc . Prompted by the Schur complement, the matrix can be factorized into upper and lower triangular block matrices b −g cc −g cc 11 such that the determinant is easily computed, since the determinant of the second matrix is 1 and the determinant of the first matrix is simply the product of the determinants of the diagonal blocks, resulting in det b −g cc −g cc 11 = det(b) det 11 − b −1g cc g cc = det b −g cc g cc = det 11 − (T cc + T ce g ee T ec )g cc .
and for an odd number of sites ℑm(Γ odd (iν)) = −V 2 ν R L (ν, t) , where R L (ν, t) is a rational function which is well behaved upon taking the limit ν → 0 + (i.e. it approaches a constant lim ν→0 + R L (ν, t) = f L (t)) and the real part is always zero. Upon inserting into eq. (B2) one can easily verify that for an even number of sites this yields zero spectral weight at the Fermi energy, while for an odd number of sites it yields −∞ (i.e. there is a pole exactly at ω = 0). Therefore the ED results alternate with even/odd system size between ρ even f (0) = 0 and ρ odd f (0) = ∞. This result shows that the FSR is always violated in ED because a finite value of the impurity density of states at the Fermi energy may only be obtained using an artificial numerical broadening. It furthermore shows that all terms involving m * go to zero and cannot contribute to the sum rule. Now we turn to the case of CPT/VCA, where the conduction electron hybridization takes the form (see eq. (10)) for the model considered in eq. (1) (i.e. a semi-circular density of states of the conduction electrons). The CPT/VCA Green's function of the physical system is then given upon insertion of this Γ into eq. (B2) ℑm(G ret f f (0)) = lim Which yields upon expansion of the square root up to linear order in ν ℑm(G ret f f (0)) = lim and when the limit ν → 0 is taken which is exactly the value predicted by the FSR (eq. (23)) This result is independent of the size of the reference system. Away from particle-hole symmetriy (occupation of the f-orbital not one) the calculation becomes more tedious. The numerical VCA calculations however show that a pinning of the Kondo resonance at the Fermi energy is obtained in contrast to CPT at small L. and the self-energy Σ ′ (z) is given by From this all elements of the cluster Green's function may be obtained by the equation of motion technique.
To be able to calculate the grand potential Ω, the ground state energy of the interacting part of the reference system, ω ′ 0 needs to be obtained (which can be done for example by diagonalization of the Hamiltonian matrix or by an integral over the Green's function) Using the expression for the CPT/VCA Green's function g and the ground state energy and taking the derivative of the grand potential eq. (20) with respect to ∆V one is able to obtain an integral expression which allows to determine V ′ within VCA Ω : d The resulting V ′ (U ) is shown in fig. 22 (left) and is used to calculate the effective mass eq. (C1) shown right in the figure. The effective mass shows exponential behavior but the exponent does not match the Bethe Ansatz result as discussed in sec. V E. Next we attempt to obtain the VCA SC solution for the two-site problem. The only variational parameter is V and therefore we determine the expectation value of σ < f † σ c 1σ > self consistently. Here 1 denotes the impurity's s-orbital. Since we are considering a spinsymmetric model we sum over both spin directions and denote this expectation value as < f † c > in the following. The hopping expectation value is given by Evaluation of this expectation value in the cluster yields < f † c > cluster = − 2 π ∞ 0 dω 4(∆V + V ) 9(∆V + V ) 2 + w 2 36(∆V + V ) 4 + (U 2 + 40(∆V + V ) 2 ) w 2 + 4w 4 . Upon requiring the two expectation values to coincide the optimal value of ∆V is obtained numerically. The resulting V ′ (U ) is plotted in fig. 22 (left) and is used to calculate the effective mass eq. (C1) shown right in the figure. The effective mass does not show exponential behavior in VCA SC .