# The structure of quasiparticles in a local Fermi liquid

Izak Snyman  Mandelstam Institute for Theoretical Physics, School of Physics,
University of the Witwatersrand, Johannesburg, South Africa

###### Abstract

Conduction electrons interacting with a dynamic impurity can give rise to a
local Fermi liquid. The latter has the same low energy spectrum as an ideal
Fermi gas containing a static impurity. The Fermi liquids’s elementary
excitations are however not bare electrons. In the vicinity of the impurity,
they are dressed by virtual particle-hole pairs. Here we study this dressing.
Among other things, we construct a mode-resolved measure of dressing. To
evaluate it in position representation, we have to circumvent the limitations
of the Numerical Renormalization Group, which discretizes the conduction band
logarithmically. We therefore extend Natural Orbital methods, that
successfully characterize the ground state, to describe excitations. We
demonstrate that the dressing profile shows nontrivial powerlaw decay at large
distances. Our Natural Orbital methodology could lay the foundation for
calculating the properties of local Fermi liquid quasiparticles in nontrivial
geometries such as disordered hosts or mesoscopic devices.

##  I Introduction

Fermi liquids emerge near infrared fixed points in the Renormalization Group
flow of some interacting many-fermion systems [1]. There is something
alchemical about them. Renormalization plays the role of the philosopher’s
stone, so that it often remains a mystery how the independent quasiparticles
(the probverbial gold) are constituted out of the raw ingredients – bare
electrons. Physicists characterize the low energy behavior of Fermi liquids in
terms of a handful of parameters that are determined from experiment. A more
daunting task is to determine how the Fermi liquid parameters depend on the
microscopic parameters of a given physical realization.

Local Fermi liquids are a class of systems in which this challenging problem
has met with success [2]. They occur in dynamic quantum impurity models, where
a small interacting quantum system is coupled to a bulk system of non-
interacting fermions. They are interesting many-body systems in their own
right [3], and also appear as a key ingredient in dynamical mean field theory
[4], an important method in the study of bulk-interacting systems. For these
reasons, their dynamics is an active field of study [5, 6, 7, 8, 9]. Often,
the quantity of interest is the impurity Green function, and nontrivial
methods have been developed, that are geared to calculating it. This exploits
the fact that a very precise description of the impurity can be achieved,
employing a less precise description of the bulk. The lack of precision where
the bulk is concerned, can lead to difficulty, when attempting to study the
structure of correlations that live in the bulk [10, 11, 12, 13, 14, 15].

The fact that a Fermi liquid description applies in certain quantum impurity
models, was established half a century ago, using the Numerical
Renormalization Group (NRG) [16]. This method was able to compute the low
energy many-body spectrum, whose structure was found to be nearly identical to
that of a non-interacting fermion system [17, 18, 19]. The spectrum obtained
for a given set of microscopic parameters was fitted to a non-interacting
fixed point Hamiltonian and leading irrelevant perturbations. In this way, the
microscopic parameters could be mapped to the Fermi liquid parameters they
give rise to [20, 21, 22, 23]. This raises a prospect which seems unfeasible
for other Fermi liquids. Can one explicitly calculate how the quasiparticle
excitations of a local Fermi liquid are constituted out of bare electronic
degrees of freedom? This question is particularly challenging to answer when a
dynamic quantum impurity is imbedded in a host in which electrons experience a
non-trivial potential landscape, such as disorder, or electrostatic gates [24,
25, 26, 27, 28, 29, 30, 31, 32]. Using scanning gate microscopy,
experimentalists have succeeded in obtaining a real space picture of
quasiparticle excitations in such systems [33, 34, 35]. The microscopic
modelling of these experiments may benefit from the development of new
methods. One reason why NRG is inadequate, is that it discretizes the bulk on
what is called a logarithmic energy grid. This throws away short wave-length
information required to achieve good spatial resolution of many-body
correlations in the vicinity of the impurity [11, 12, 36]. A second reason is
that NRG discards all but the very lowest single-quasiparticle excitations
before their long wave length structure is fully resolved [37, 38].

Another single-particle-like picture, that is distinct from Fermi liquid
theory, is important for understanding many-body correlations in the ground
states of quantum impurity problems [39, 40]. Associated with this picture are
a set of bare single-particle orbitals called natural orbitals [41, 42]. They
form the single-particle basis in which a correlated ground state is expressed
as a linear combination of the fewest number of Slater determinants. Natural
orbital methods have an established role in quantum chemistry, where they
successfully account for chemical properties of strongly correlated few-
electron systems [43, 44]. There has recently been increased interest in the
applications of natural orbital methods in condensed matter [45, 46],
including quantum impurity problems [47, 48], where high accuracy in very
large systems have been achieved at modest computational cost [15]. The
present study was inspired by the question as to whether natural orbitals can
shed light on the structure of quasiparticle excitations in a local Fermi
liquid.

Our main results are as follows. We have developed methods for analyzing
single-quasiparticle excitations in a large finite system consisting of a
quantum impurity hybridized with a noninteracting host. It clarifies the
relationship between the single-particle eigenstates of the noninteracting
effective quasiparticle Hamiltonian, and the associated many-body eigenstates
of the original interacting Hamiltonian. This connection allows one to split
the many-body eigenstate into two parts, corresponding respectively to a bare
electron on top of the ground state and to the dressing of the bare electron
by particle-hole pairs. We developed an accurate ansatz for single-
quasiparticle excitations on top of the ground state in terms of natural
orbitals. This allowed us to investigate quasiparticle excitations that are
discarded in NRG, and to investigate the structure of the dressing of bare
electrons in real space. We could study a bulk consisting of thousands of
sites and resolve the wave function with single site precision. We thus
obtained results on the non-trivial power-law decay of dressing with distance,
that NRG is unable to resolve. Our results could lead to advances in the study
of spatial features of low energy excitations in local Fermi liquids,
including the modelling of realistic environments, such as atomic lattices,
static disorder, and mesoscopic electronic devices. It could also lead to the
further refinement of natural orbital methods to study the dynamics of quantum
impurities.

The rest of the Article is structured as follows. In Section II we introduce
the Single Impurity Anderson Model [49], that we will focus on. We review
elementary aspects of the local Fermi liquid theory that applies to it at low
energies, and formulate the questions regarding its quasiparticle excitations
that guided the present study. We explain the limitations of NRG, and review
the natural orbital methods that apply to the study of the model’s ground
state. In Section III we present the theoretical developments that our study
contribute to local Fermi liquid theory. In Section IV we employ the theory
that we developed in Section III to numerically study the structure of single-
quasiparticle excitations. Section V contains a summary of our main results,
and concluding remarks. As evidence that our natural orbital methods are
sufficiently accurate for our purposes, in Appendix A we compare to results
obtained by other means.

##  II Background and aims of the study

###  II.1 Local Fermi liquid theory and the Single Impurity Anderson Model

In this Article, our primary focus will be on the Single Impurity Anderson
Model (SIAM) [49, 2]. It describes $\hat{d}^{\dagger}_{\sigma}$ electrons in a
localized orbital interacting via Coulomb repulsion $U$ and hybridizing with
the non-interacting $\hat{c}_{k\sigma}^{\dagger}$ electrons in a Fermi sea.
The Hamiltonian reads

$\displaystyle\hat{H}=$ $\displaystyle
U\hat{n}_{d\uparrow}\hat{n}_{d\downarrow}+\varepsilon_{d}\left(\hat{n}_{d_{\uparrow}}+\hat{n}_{d\downarrow}\right)+V\sum_{\sigma}\left(\hat{d}^{\dagger}_{\sigma}\hat{\psi}_{0\sigma}+\hat{\psi}_{0\sigma}^{\dagger}\hat{d}_{\sigma}\right)$
$\displaystyle+\sum_{k=1}^{L}\sum_{\sigma}\varepsilon_{k}\hat{c}_{k\sigma}^{\dagger}\hat{c}_{k\sigma},$
(1)

where $\hat{n}_{d\sigma}=\hat{d}_{\sigma}^{\dagger}\hat{d}_{\sigma}$, and
$\hat{\psi}_{0\sigma}^{\dagger}=\frac{1}{\sqrt{L}}\sum_{k=1}^{L}\hat{c}_{k\sigma}^{\dagger}$
creates an electron on the lattice site closest to the impurity. To quantify
the hybridization, it is conventional to cite the spectral density

$\Delta=\pi\nu V^{2}$ (2)

at the Fermi level, of the operator $V\hat{\psi}_{0}$ in the infinite system
uncoupled from the impurity, rather than $V$. In the above expression, $\nu$
is the bulk density of states per unit volume at the Fermi level. We measure
single-particle energies such that the Fermi energy lies at zero. If
$-\varepsilon_{d}>\Delta$ and $\varepsilon_{d}+U>\Delta$, the ground state
occupation probability of the $d$-orbital is approximately one, meaning it
possesses a spin-1/2 magnetic moment. We will focus on the particle-hole
symmetric case, where the energy of the $d$-orbital is tuned such that
$\varepsilon_{d}=-U/2$, and the condition for single occupation of the
$d$-orbital becomes $U/2>\Delta$. Virtual processes in which the singly-
occupied impurity is emptied and then reoccupied, or doubly-occupied and then
returned to single occupancy, lead to an effective antiferromagnetic spin-
exchange interaction between the impurity and the band. Kondo physics results,
associated with an emergent Kondo scale $T_{\text{K}}$, given by the inverse
of the $d$-orbital’s susceptibility to being spin-polarized by a local
magnetic field

$T_{\text{K}}=\frac{1}{4\chi},~{}~{}~{}\chi=\lim_{B\to
0}\left<\left(\hat{n}_{d\uparrow}-\hat{n}_{d\downarrow}\right)\right>_{B}/B.$
(3)

Here $\left<\ldots\right>_{B}$ denotes the ground state expectation value when
the Hamiltonian is perturbed by the local field

$\hat{H}_{B}=\hat{H}-\frac{B}{2}\left(\hat{n}_{d\uparrow}-\hat{n}_{d\downarrow}\right).$
(4)

Krishna-Murthy, Wilkins and Wilson (KWW) [19], showed by means of NRG that the
SIAM flows to strong coupling, and reaches a stable infrared fixed point below
$T_{\text{K}}$. At energies sufficiently lower than the Kondo scale, the many-
body energy spectrum of the model resembles that of a non-interacting Fermi
gas, i.e. it is given by

$E(\\{\nu_{\alpha\sigma}\\})=\text{const.}+\sum_{\alpha}E_{\alpha}n_{\alpha\sigma},$
(5)

where the quantum numbers $n_{\alpha\sigma}\in\\{0,1\\}$ can be interpreted as
the occupation numbers of quasiparticle orbitals. While no explicit mapping
relating the low-energy quasiparticles to the bare fermions appearing in (1)
have ever been derived, a priori one expects the quasiparticle creation and
annihilation operators to be nonlinear functions of the bare fermion operators
$\hat{c}_{k\sigma}$, $\hat{c}_{k\sigma}^{\dagger}$, $\hat{d}_{\sigma}$ and
$\hat{d}_{\sigma}^{\dagger}$. In other words, we expect the quasiparticles to
be fermions dressed with particle-hole pairs.

Hewson [20], building on ideas of Nozières [18], also introduced a
quasiparticle Hamiltonian for the SIAM, as a starting point for a
“renormalized” perturbation theory. Hewson relied on the generic behavior of
the (zero-temperature) retarded self-energy of a class of fermionic models
that includes the SIAM, namely that the self-energy can be expanded in
frequency around $\omega=0$,

$\Sigma(\omega)=\Sigma(0)+\Sigma^{\prime}(0)\omega+\Sigma^{\text{rem}}(\omega),$
(6)

and that $\Sigma(0)$ and $\Sigma^{\prime}(0)$ are both real. The interacting
retarded Green function for the $d$-orbital reads
$G_{\sigma\sigma^{\prime}}(\omega)=\delta_{\sigma\sigma^{\prime}}/(\omega-\varepsilon_{d}+i\Delta-\Sigma(\omega))$.
Hewson rewrote $G_{\sigma\sigma^{\prime}}(\omega)$ as

$G_{\sigma\sigma^{\prime}}(\omega)=\frac{z\delta_{\sigma\sigma^{\prime}}}{\omega-\tilde{\varepsilon}_{d}+i\tilde{\Delta}-\tilde{\Sigma}(\omega)},$
(7)

where $z=1/[1-\Sigma^{\prime}(0)]$,
$\tilde{\varepsilon}_{d}=z[\varepsilon_{d}+\Sigma(0)]$ and
$\tilde{\Sigma}(\omega)=z\Sigma^{\text{rem}}(\omega)\propto\omega^{2}$. (In
the particle-hole symmetric case, $\tilde{\varepsilon}_{d}=0$.) At
sufficiently low frequencies/energies, $\tilde{\Sigma}(\omega)$ makes a very
small contribution to the denominator, and can be neglected. Hewson therefore
defined a quasiparticle Hamiltonian $\hat{H}_{\text{qp}}$ that has the same
form as the non-interacting version of (1), with $U=0$. Furthermore he made
the replacements $V\to\tilde{V}=\sqrt{z}V$ and
$\varepsilon_{d}\to\tilde{\varepsilon}_{d}$. If we denote the retarded
$d$-orbital Green function of Hewson’s quasiparticle Hamiltonian by
$\tilde{G}_{\sigma,\sigma^{\prime}}(\omega)$, and drop
$\tilde{\Sigma}(\omega)$ in the expression (7) for
$G_{\sigma\sigma^{\prime}}(\omega)$, then
$G_{\sigma\sigma^{\prime}}(\omega)=z\tilde{G}_{\sigma\sigma^{\prime}}(\omega)$.
The elastic scattering amplitudes for a bare electrons incident on the
impurity are determined by the transfer matrix

$\displaystyle T_{kk^{\prime}}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{L}\frac{V}{\omega-\varepsilon_{k}}G_{\sigma\sigma}(\omega)\frac{V}{\omega-\varepsilon_{k^{\prime}}}$
(8) $\displaystyle\simeq$
$\displaystyle\frac{1}{L}\frac{\tilde{V}}{\omega-\varepsilon_{k}}\tilde{G}_{\sigma\sigma}(\omega)\frac{\tilde{V}}{\omega-\varepsilon_{k^{\prime}}}.$

The second line of the above equation implies that at low energies, bare
electrons scatter off the impurity as if their dynamics are described by the
non-interacting Hamiltonian $\hat{H}_{\text{qp}}$. Does this mean that our a
priori expectations about KWW’s quasiparticles were wrong? Could the
quasiparticles be bare electrons described by Hewson’s $\hat{H}_{\text{qp}}$?
The answer is “no”. $\hat{H}_{\text{qp}}$ has a $d$-electron Green function
that is off by a factor $z$ compared to the original Hamiltonian (1). Thus at
least the $d$-orbital in $\hat{H}_{\text{qp}}$ is an effective, rather than a
bare degree of freedom. The fact that low energy bare electrons scatter
elastically does however imply the following. If the system is prepared in an
eigenstate corresponding to a single quasiparticle on top of the ground state,
no particle-hole pairs escape to infinity. Instead, the particle-hole pairs
that dress a bare electron to make up a quasiparticle excitation, are confined
to a finite region around the impurity. If we write the creation operator for
a quasiparticle excitation, as

$y_{d}\hat{d}^{\dagger}+\sum_{k}y_{k}\hat{c}_{k\sigma}^{\dagger}+\mbox{higher
order terms
in}\\{\hat{d}^{\dagger}_{\sigma},\hat{d}_{\sigma},\hat{c}^{\dagger}_{k,\sigma},\hat{c}_{k\sigma}\\},$
(9)

then the higher order terms create particle-hole pairs in the vicinity of the
impurity.

###  II.2 Aims of this Article

In this study we want to shed further light on the structure of the
quasiparticles associated with the SIAM. Our aims are (1) to calculate the
quasiparticle wave function $y_{d},\,y_{k},\,k=1,2,\ldots,L$ in (9)
explicitly, (2) to visualize in real space the particle-hole dressed part in
the vicinity of the impurity, described by the higher order terms in (9), and
(3) see the onset of effects associated with the remnant self-energy
$\Sigma^{\text{rem}}(\omega)$ at increased excitation energy. These questions
require the development of new methods, as we explain next.

###  II.3 Limitations of NRG

For the purposes of numerics, one has to study a system with a finite-
dimensional Hilbert space. One possibility is to study a system defined on a
regular finite lattice. However, the many-body Hilbert space becomes too large
for brute-force methods already at $\sim 10$ sites. A more sophisticated
option is to take an infinite system, and rediscretize it. There is
considerable freedom in how to discretize, and particular choices may be
better suited to subsequent approximations than others [16]. The general
rediscretization procedure works as follows. The energy band is partitioned
into intervals $\mathcal{E}_{k}$ called energy shells. One discrete mode per
shell, per channel, is retained, with associated creation operator
$\hat{c}_{k\sigma}^{\dagger}=\int_{\mathcal{E}_{k}}d\varepsilon\,\hat{c}_{\varepsilon\sigma}^{\dagger}/\sqrt{\int_{\mathcal{E}_{k}}d\varepsilon}$.
Here $\hat{c}_{\varepsilon\sigma}^{\dagger}$ creates a particle with energy
$\varepsilon$ and spin $\sigma$ (in a particular channel). Impurity problems,
such as the SIAM, involve local coupling to a point impurity and scattering is
therefore s-wave, so that only a single channel is involved. Discrete mode
$k$’s energy is taken as the average energy of the interval $\mathcal{E}_{k}$.
The NRG method, to which we owe many results on fermionic quantum impurity
problems, relies on the fact that accurate results can be obtained for low-
energy properties, from a cruder description at higher energies. It
discretizes the band “logarithmically” so that the density of discrete states
scales like $1/\left|\varepsilon_{k}-E_{\text{F}}\right|$. This is done by
letting the width of each new interval decrease by a constant factor
$1/\Lambda$ compared to the previous one, from the band edges to the Fermi
energy. NRG proceeds with a sequence of RG transformations [17]. Step $k$
involves an approximate diagonalization that resolves energies up to shell
$k$. Before the modes in lower energy shells are included, high energy states
found in the current iteration are discarded, thus avoiding the dimension of
the many-body Hilbert space becoming too large to handle. It turns out that
the logarithmic discretization ensures accuracy at lower energies being
maintained. This accuracy extends to thermodynamic quantities such as specific
heat, and impurity spectral quantities such as the local density of states of
the $d$-orbital. However, due to the decimation process that keeps the size of
the Hilbert space manageable, in practice, only the lowest few single-
quasiparticle excitations are found in the end. Furthermore, NRG inherently
suffers from poor spatial resolution of correlations in the Fermi sea, because
large shells of high-wavelength modes are crudely lumped together into a few
discrete modes [11, 12]. (Often only two modes are retained to represent the
top and bottom quarters of the band, i.e. $\Lambda=2$.) Our aim is to study
the structure of single-quasiparticle modes in terms of bare electron degrees
of freedom. This requires us to resolve spatial features all the way down to
the Fermi wavelength, which NRG cannot do. Natural orbital methods have
recently been shown capable of the required spatial resolution, where ground
state correlations are concerned [15]. In the next subsection we review the
application of natural orbital methods to impurity ground state problems.

###  II.4 Covariance matrix and natural orbitals

Given a set of fermionic creation and annihilation operators, associated with
an orthonormal single-particle basis, and an arbitrary state $\left|X\right>$
in the associated Fock space, the covariance matrix $C(\left|X\right>)$ is
defined as

$C_{ij}\left(\left|X\right>\right)=\left<X\right|\hat{c}_{i}^{\dagger}\hat{c}_{j}\left|X\right>.$
(10)

We will denote the eigenvectors and eigenvalues of the ground state covariance
matrix as follows.

$\bm{x}_{\alpha}=(x_{1,\alpha},\ldots,x_{L,\alpha})^{\text{T}},~{}~{}~{}C(\left|\text{GS}\right>)\bm{x}_{\alpha}=\lambda_{\alpha}\bm{x}_{\alpha}.$
(11)

We label eigenvectors and eigenvalues such that
$\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{L}$. The eigenvalues of the
covariance matrix are independent of the single-particle basis associated to
$\hat{c}_{i}$. The fermionic creation operators

$\hat{q}_{\alpha}^{\dagger}=\sum_{i=1}^{L}\hat{c}_{i}^{\dagger}x_{i\alpha},$
(12)

associated to eigenstates of the covariance matrix are known as the natural
orbital basis [42]. This basis is independent of the the single-particle basis
associated to $\hat{c}_{i}$. The eigenvalues $\lambda_{\alpha}$ are the ground
state occupation probabilities
$\left<\text{GS}\right|\hat{q}_{\alpha}^{\dagger}\hat{q}_{\alpha}\left|\text{GS}\right>$
of the natural orbitals. When fermions are held at a finite density, the
exclusion principle forces many natural orbitals to have ground state
occupation probabilities near unity. Orbitals that are nearly filled or empty,
are inert and cannot participate in many-body correlations. For a given small
positive number $\epsilon$, this motivates us to define three sets

$\displaystyle\mathcal{O}$ $\displaystyle=$
$\displaystyle\left\\{\alpha|\lambda_{\alpha}>1-\epsilon\right\\},$
$\displaystyle\mathcal{C}$ $\displaystyle=$
$\displaystyle\left\\{\alpha|\epsilon<\lambda_{\alpha}<1-\epsilon\right\\},$
$\displaystyle\mathcal{U}$ $\displaystyle=$
$\displaystyle\left\\{\alpha|\lambda_{\alpha}<\epsilon\right\\},$ (13)

which we call the occupied, correlated, and unoccupied sectors respectively.
(Quantum chemists refer to these sets as the inactive, active and virtual
spaces.)

In a quantum impurity problem, the bulk remains non-interacting, and the
correlated sector only contains a vanishing fraction of the total number of
orbitals in a large system. It was recently realized that for a generic
fermionic impurity problem, the number of orbitals in the correlated sector is
proportional to $-{\rm ln}\,\epsilon$ at sufficiently small $\epsilon$, with a
proportionality constant that remains finite in the thermodynamic limit. As a
consequence, the ground state of a generic fermionic impurity model can be
approximated as follows, if the natural orbitals are known [40]. Consider a
model with $N$ particles. At given $\epsilon$, let $N_{\mathcal{O}}(\epsilon)$
and $N_{\mathcal{C}}(\epsilon)$ be the numbers of orbitals in respectively the
occupied and correlated sectors. For sufficiently small $\epsilon$, we can
approximate the orbitals in the occupied sector as fully occupied, and those
in the unoccupied sector as completely empty, meaning that in the ground state
the correlated sector contains $m=N-N_{\mathcal{C}}(\epsilon)$ particles
distributed among $N_{\mathcal{C}}(\epsilon)$ orbitals. Let

$\left|\text{F}\right>=\hat{q}_{N_{\mathcal{O}}}^{\dagger}\ldots\hat{q}_{1}^{\dagger}\left|0\right>,$
(14)

be the Fermi sea corresponding to the completely filled occupied sector. Given
a set of $n$ orbitals ${\left\\{\alpha_{i}\right\\}_{n}}\subset\mathcal{C}$
belonging to the correlated sector, we define

$\left|{\left\\{\alpha_{i}\right\\}_{n}}\right>=\left({\prod}^{\prime}_{\alpha\in{\left\\{\alpha_{i}\right\\}_{n}}}\hat{q}_{\alpha}^{\dagger}\right)\left|\text{F}\right>.$
(15)

Here the prime denotes a fixed ordering of operators (say $\alpha$ decreasing
from left to right), to remove ambiguity about the phase of the state. We take
as ground state ansatz, an arbitrary linear combination of $N$-particle states
of the form (15):

$\left|\text{GS}\right>\simeq\sum_{\left\\{\alpha_{i}\right\\}_{m}}v_{{\left\\{\alpha_{i}\right\\}_{m}},\text{GS}}\left|{\left\\{\alpha_{i}\right\\}_{m}}\right>.$
(16)

The optimal state is found by minimizing the expectation value of the energy
over the coefficients $v_{{\left\\{\alpha_{i}\right\\}_{m}},\text{GS}}$. Thus
it is found that the optimal expansion coefficients
$v_{{\left\\{\alpha_{i}\right\\}_{m}},\text{GS}}$ correspond to the ground
state eigenvector of the effective few-body Hamiltonian

$\left[H_{\text{few\,body}}\right]_{{\left\\{\alpha_{i}\right\\}_{m}},{\left\\{\beta_{i}\right\\}_{m}}}=\left<{\left\\{\alpha_{i}\right\\}_{m}}\right|\hat{H}\left|{\left\\{\beta_{i}\right\\}_{m}}\right>.$
(17)

This type of approximation has a long history in quantum chemistry, where it
is called the Complete Active Space (CAS) approach [43]. In the context of
quantum impurity problems, stronger results regarding accuracy apply than in
quantum chemistry, thanks to the proven scaling of the size of the correlated
sector with $\epsilon$: The dimension of the effective few-body Hamiltonian is
$\left(\begin{array}[]{c}N_{\mathcal{C}}(\epsilon)\\\ m\end{array}\right)$,
with $m\sim N_{\mathcal{C}}(\epsilon)/2$, which scales exponentially with
$N_{\mathcal{C}}$. However, $N_{\mathcal{C}}$ only scales logarithmically with
$\epsilon$, and the dimension of the Hamiltonian that has to be diagonalized
therefore scales polynomially with $1/\epsilon$. These features of fermionic
quantum impurity problems have been exploited to prove that the computational
complexity of finding the ground state of a fermionic quantum impurity problem
scales quasi-polynomially with the inverse of the required accuracy and
polynomially with the system size [39]. What makes the proof non-trivial is
the fact that the natural orbitals are defined relative to the ground state,
and therefore not known beforehand. In practice, iterative algorithms are
found to work well. In these algorithms, a guess for the natural orbitals is
recursively improved from the previous iteration’s result for the approximate
ground state. This is known as the Recursive Generation of Natural Orbitals
(RGNO) [44]. An early application of the method in a Condensed Matter context
considered a multi-impurity system [48]. Below we find that for the SIAM, the
ground state energy can be determined to an accuracy of 5% of the Kondo
temperature, for realistic Kondo temperatures of 1 % of the band width, and
well-developed correlations (quasiparticle weight $z=0.2$), using very modest
resources, namely a correlated sector consisting of 12 orbitals (six spin up
and six spin down), containing 6 particles (three spin up and three spin
down).

###  II.5 Discretizations employed

Below, we will study quasiparticles using both NRG and natural orbital
methods. As stated above, NRG requires a logarithmic discretization. Natural
orbital methods allow more freedom. Here we provide details regarding the
different discretizations we will employ. For simplicity, we will work with a
continuum model that has a half-bandwidth $D$ throughout, and assume that the
continuum model has a flat density of states per unit length so that

$\Delta=\frac{\pi V^{2}}{2D}.$ (18)

For a logarithmic discretization, the conduction band is divided into
intervals

$\left[-D\Lambda^{-k},-D\Lambda^{-k-1}\right)\text{ and
}\left(D\Lambda^{-k-1},D\Lambda^{-k}\right],$ (19)

with $\,k=0,1,2,\ldots$. The energies of the rediscretized conduction band
orbitals are

$\varepsilon_{k\pm}=\pm\frac{D}{2\Lambda^{k}}\left(\frac{1}{\Lambda}+1\right).$
(20)

After the above discretization, the discrete system still possesses an
infinite number of modes. NRG proceeds by turning the diagonal kinetic term
into a tri-diagonal form

$\displaystyle\hat{H}=$ $\displaystyle
U\hat{n}_{d\uparrow}\hat{n}_{d\downarrow}+\varepsilon_{d}\left(\hat{n}_{d_{\uparrow}}+\hat{n}_{d\downarrow}\right)+V\sum_{\sigma}\left(\hat{d}^{\dagger}_{\sigma}\hat{\psi}_{0\sigma}+\hat{\psi}_{0\sigma}^{\dagger}\hat{d}_{\sigma}\right)$
$\displaystyle+\sum_{j=0}^{\infty}\sum_{\sigma}t_{j}\left(\hat{\psi}_{j\sigma}^{\dagger}\hat{\psi}_{j+1\sigma}+\hat{\psi}_{j+1\sigma}^{\dagger}\hat{\psi}_{j\sigma}\right),$
(21)

called the Wilson chain, with the off-diagonal elements $\propto t_{j}$
describing hopping along the chain. Their explict form can be looked up in
[16]. The first chain site corresponds to the lattice site $0$ that is
directly coupled to the impurity. Hopping amplitudes decrease like
$\Lambda^{-j/2}$ along the chain, with site $j$ associated with energy scale
$D\Lambda^{-j/2}$. An infrared cut-off is imposed by truncating the chain
after $L$ sites, which yields a finite system, and a smallest energy scale
$\sim D\Lambda^{-L/2}$. We will perform both NRG and natural orbital
calculations on the Wilson chain.

The Wilson chain is not suitable for resolving the real space structure of
quasiparticle excitations. We will therefore also use a second, more suitable
discretization, in conjunction with natural orbital methods. For this purpose,
we take the discrete energies of the finite system to be

$\varepsilon_{k}=\frac{k-\Omega-1}{\Omega}D,~{}~{}~{}L=2\Omega+1.$ (22)

To study spatial structure, we interpret the $\hat{c}_{k\sigma}^{\dagger}$
operators associated with these energies as the even modes of a one-
dimensional lattice with $4\Omega+1$ sites, the central site ($j=0$) of which
is side-coupled to the $d$-orbital. We introduce symmetrized position
representation operators

$\hat{\gamma}_{j\sigma}=\sqrt{\frac{1}{2\Omega+1}}\sum_{k=1}^{2\Omega+1}\cos\left[\frac{\pi
j(2k-1)}{4\Omega+2}\right]\hat{c}_{k\sigma},$ (23)

for $j=0,\ldots,2\Omega$. (Symmetrization means that $\hat{\gamma}_{j\sigma}$
is one half times the sum of the operators that respectively annihilate an
electron on sites $j$ and $-j$. This way we avoid having to introduce creation
operators for the odd parity modes that do not couple to the $d$-orbital.) The
$\hat{\gamma}_{j\sigma}$ operators obey

$\left\\{\hat{\gamma}_{j,\sigma},\hat{\gamma}_{j^{\prime},\sigma^{\prime}}^{\dagger}\right\\}=\frac{\delta_{j,j^{\prime}}\delta_{\sigma,\sigma^{\prime}}}{(2-\delta_{j,0})}.$
(24)

Since the infrared energy cutoff $D/\Omega$ scales inversely rather than
exponentially with system size, very large systems are required to resolve the
emergent infrared physics of quantum impurity models. Natural orbital methods
have proved capable of this task, where ground state properties are concerned.

##  III Theoretical developments

In this section, we present two theoretical developments, that we made, and
which allows us to learn more about the the structure of single-quasiparticles
in a local Fermi liquid, than was known before. The first is a wave function
picture of local Fermi liquid theory, applicable to finite systems. This
complements Hewson’s Green function picture for infinite systems. It provides
us with tools to analize the structure of single-quasiparticle excitations.
The second development is an ansatz for single-quasiparticle excitations in
terms of natural orbitals. We will find that results are accurate in a regime
where correlations are sufficiently strong to study non-trivial local Fermi
liquids, and amenable to discretization on a large regular energy grid.

###  III.1 Wave function picture of local Fermi liquids

Suppose the ground state $\left|\text{GS}\right>$, as well as excited states
$\left|\text{p},n\sigma\right>$ with a single quasiparticle on top of the
ground state, are known to good accuracy. How do we find the (dominant) linear
part, cf. (9) of a quasiparticle operator $\hat{a}_{n\sigma}$ such that

$\left|\text{p},n\sigma\right>{\simeq}\hat{a}_{n\sigma}^{\dagger}\left|\text{GS}\right>?$
(25)

We can answer this question by setting

$\hat{a}_{n\sigma}=y_{dn}\hat{d}_{\sigma}+\sum_{k=1}^{L}y_{kn}\hat{c}_{k\sigma},~{}~{}~{}|y_{dn}|^{2}+\sum_{k=1}^{L}|y_{kn}|^{2}=1,$
(26)

and maximizing the object function

$p_{n}=\left|\left<\text{GS}\right|\hat{a}_{n\sigma}\left|\text{p},n\right>\right|^{2},$
(27)

over all unit-length vectors

$\bm{y}_{n}=(y_{dn},y_{1n},\ldots,y_{Ln})^{\text{T}}.$ (28)

The optimal solution is found to be

$y_{dn}=\frac{\left<\text{p},n\right|\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>}{\sqrt{p_{n}}},~{}~{}~{}y_{kn}=\frac{\left<\text{p},n\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>}{\sqrt{p_{n}}},$
(29)

The optimal $p_{n}$ that results is

$p_{n}=\left|\left<\text{GS}\right|\hat{d}_{\sigma}\left|\text{p},n\right>\right|^{2}+\sum_{k=1}^{L}\left|\left<\text{GS}\right|\hat{c}_{k\sigma}\left|\text{p},n\right>\right|^{2}.$
(30)

The quantity $p_{n}$ is the probability to measure only a single bare electron
excitation on top of the ground state when the state
$\left|\text{p},n\sigma\right>$ is prepared. The quantity

$\delta\rho(E)=\sum_{n}(1-p_{n})\delta(E-E_{n}),$ (31)

equals the difference in density of single-particle states between a system of
bare non-interacting electrons with the same spectrum as the quasiparticles,
and the actual system. (Here we do not mean “density of states per unit
volume”, but the actual denisty of states, that diverges in the thermodynamic
limit.) The density of states difference remains finite in the thermodynamic
limit, and quantifies the extent to which bare electrons are dressed in order
to form quasiparticles. It is tempting to use Fermi liquid green functions,
with the remnant self-energy $\Sigma^{\text{rem}}(\omega)$ neglected, to
evaluate $\delta\rho$ at low energies in the thermodynamic limit. Since the
quasiparticle Green function and the actual single-electron Green function
only differ on the $d$-orbital, one finds
$\delta\rho(\omega)=\text{Im}\left[\tilde{G}_{\sigma\sigma}(\omega)-G_{\sigma\sigma}(\omega)\right]/\pi$.
At the Fermi level, where $\Sigma^{\text{rem}}(\omega)=0$, this is certainly
valid, and gives

$\delta\rho(0)=\frac{1}{\pi}\left(\frac{1}{\tilde{\Delta}}-\frac{1}{\Delta}\right)=\frac{1}{\pi\Delta}\left(\frac{1}{z}-1\right).$
(32)

In a large but finite system, in which the spacing of single-particle levels
near the Fermi energy is $\delta E$, this gives

$z\simeq\frac{1}{1+p_{n}\pi\Delta/\delta E}.$ (33)

Thus the wave function picture of a local Fermi liquid furnishes us with an
interpretation of the quasiparticle weight in terms of the bare electron
occupation probability.

Above the Fermi level, the reasoning that led to (32) would predict that
$\delta\rho(\omega)$ decreases proportional to
$1/(\omega^{2}+\tilde{\Delta}^{2})$, whereas we expect $p_{n}$ to decrease and
hence $\delta\rho$ to increase as the quality of quasiparticles deteriorates
with increasing excitation energy. Below we perform numerics to resolve the
behaviour of $\delta\rho(\omega)$ at $\omega>0$.

A mode-resolved measure $D_{kn}$ of how a bare electron is dressed to form the
excitation $\left|\text{p},n\sigma\right>$ is obtained as follows. Since
$\left|\text{p},n\sigma\right>=\sqrt{p_{n}}\hat{a}_{n\sigma}^{\dagger}\left|\text{GS}\right>+\left|\mbox{dressing}\right>$,
we set

$\displaystyle D_{kn}$
$\displaystyle\equiv\left<\mbox{dressing}\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
$\displaystyle=\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>-\sqrt{p_{n}}\left<\text{GS}\right|\hat{a}_{n\sigma}\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
$\displaystyle=\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>-\sqrt{p_{n}}\left<\text{GS}\right|\\{\hat{a}_{n\sigma},\hat{c}_{k\sigma}^{\dagger}\\}\left|\text{GS}\right>$
$\displaystyle~{}+\sqrt{p_{n}}\left<\text{GS}\right|\hat{c}_{k\sigma}^{\dagger}\hat{a}_{n\sigma}\left|\text{GS}\right>$
$\displaystyle=\underbrace{\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>-\sqrt{p_{n}}y_{kn}}_{=0}+\sqrt{p_{n}}\left<\text{GS}\right|\hat{c}_{k\sigma}^{\dagger}\hat{a}_{n\sigma}\left|\text{GS}\right>$
$\displaystyle=\left[\sqrt{p_{n}}C(\left|\text{GS}\right>)\bm{y}_{n}\right]_{k},$
(34)

where we used the definition of $y_{kn}$ (29) twice in the last two lines.
This measure equals zero when $\left|\text{p},n\sigma\right>$ is associated
with a bare electron on top of the ground state (i.e. $p_{n}=1$). When the
index $k$ refers to wave number, a spatially resolved measure of dressing is
obtained by Fourier transforming to real space. Given the arguments presented
above, we expect to see a signal that decays to zero as we move away from the
impurity.

Intuitively, we suspect a close link between $\bm{y}_{n}$ and the eigenvectors
of the type of single-quasiparticle Hamiltonian identified by Hewson. To
explore this link, consider the quantities
$\left<\text{p},n\sigma\right|[\hat{H},\hat{d}^{\dagger}_{\sigma}]\left|\text{GS}\right>$
and
$\left<\text{p},n\sigma\right|[\hat{H},\hat{c}_{k\sigma}^{\dagger}]\left|\text{GS}\right>$,
where $\hat{H}$ is the particle-hole symmetric SIAM. This leads to the
equations

$\displaystyle(E_{n}-E_{\text{GS}})$
$\displaystyle\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
$\displaystyle=\frac{V}{\sqrt{L}}\left<\text{p},n\sigma\right|\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>+\varepsilon_{k}\left<\text{p},n\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>,$
(35)

and

$\displaystyle(E_{n}-E_{\text{GS}})$
$\displaystyle\left<\text{p},n\sigma\right|\hat{d}_{\sigma}^{\dagger}\left|\text{GS}\right>$
$\displaystyle=U\left<\text{p},n\sigma\right|(\hat{n}_{d-\sigma}-1/2)\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>$
$\displaystyle~{}~{}+\frac{V}{\sqrt{L}}\sum_{k=1}^{L}\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>.$
(36)

If

$\frac{U\left<\text{p},n\sigma\right|(\hat{n}_{d-\sigma}-1/2)\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>}{\frac{V}{\sqrt{L}}\sum_{k=1}^{L}\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>}\equiv(1-z)$
(37)

does not vary much as a function of $n$, we can employ it as a definition of
the quasiparticle weight $z$. That the above-defined $z$ is approximately
independent of $n$ for low energy states can be made plausible by noting that
$(\hat{n}_{d-\sigma}-1/2)\hat{d}^{\dagger}_{\sigma}$ and
$\sum_{k=1}^{L}\hat{c}_{k\sigma}^{\dagger}/\sqrt{L}$ probe the system locally
at one end, while $n\to n+1$ affects the phase difference between amplitudes
on adjacent sites by an amount $\sim 1/L\ll 1$ only. Note that via (36), (37)
is equivalent to

$z=\frac{(E_{n}-E_{\text{GS}})\left<\text{p},n\sigma\right|\hat{d}_{\sigma}^{\dagger}\left|\text{GS}\right>}{\frac{V}{\sqrt{L}}\sum_{k=1}^{L}\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>}.$
(38)

We substitute (37) into (36) and define

$\psi_{dn}=\frac{\left<\text{p},n\sigma\right|\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>}{\sqrt{z}},~{}~{}~{}\psi_{kn}=\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>,$
(39)

and $\tilde{V}=\sqrt{z}V$, to obtain the effective single-quasiparticle
Schrödinger equation

$\displaystyle(E_{n}-E_{\text{GS}})\psi_{kn}$ $\displaystyle=$
$\displaystyle\frac{\tilde{V}}{\sqrt{L}}\psi_{dn}+\varepsilon_{k}\psi_{kn},$
$\displaystyle(E_{n}-E_{\text{GS}})\psi_{dn}$ $\displaystyle=$
$\displaystyle\frac{\tilde{V}}{\sqrt{L}}\sum_{k=1}^{L}\psi_{kn},$ (40)

associated with a non-interacting SIAM ($U=0$) and a renormalized
hybridization $\tilde{\Delta}=z\Delta$. This is precisely the form identified
by Hewson, and implies that $\psi_{dn}=\sqrt{p_{n}/z}\,y_{dn}$,
$\psi_{kn}=\sqrt{p_{n}}\,y_{kn}$.

###  III.2 Ansatz

Next, we explain our approach to generalizing natural orbital methods to
excited states. We consider $N\pm 1$ particle excited states that are obtained
by adding or removing a particle from the non-interacting system, and
adiabatically switching on the interaction. If the low-energy physics is that
of a local Fermi liquid, these excitations will be single-quasiparticle- or
-hole-like in nature.

Based on the notion that the natural orbital basis minimizes the Hamiltonian’s
ability to generate particle-hole pairs, we make the following ansatz for
excited states consisting of one particle on top of the ground state.

$\left|\text{p},k\right>_{1}=\sum_{\alpha\in\mathcal{U}}u_{\alpha,k}\,\hat{q}_{\alpha}^{\dagger}\left|\text{GS}\right>+\sum_{\left\\{\alpha_{i}\right\\}_{m+1}}v_{{\left\\{\alpha_{i}\right\\}_{m+1}},k}\left|{\left\\{\alpha_{i}\right\\}_{m+1}}\right>.$
(41)

This amounts to assuming that

  1. 1.

adding a particle leaves the occupied sector undisturbed,

  2. 2.

there is at most one particle in the unoccupied sector,

  3. 3.

and if the added particle ends up in the unoccupied sector, the particles in
the correlated sector are are not disturbed from the way they were configured
in the ground state.

(Of these assumptions, the last one (3) seems the most arbitrary. In Appendix
A we therefore investigate the consequences of not making this assumption. The
conclusion is that assumption (3) does in fact apply in the local Fermi liquid
regime.) We find the coefficients $u_{\alpha,k}$ and
$v_{{\left\\{\alpha_{i}\right\\}_{m+1}},k}$ variationally. Formally, we first
find the lowest energy state of this form, then we vary orthogonal to the
lowest energy state, to get the second lowest energy state and so on. Luckily,
one does not ever explicitly have to parametrize the space orthogonal to the
states found already: it is easy to show that the coefficients of the states
obtained in the procedure are simply the eigenvectors of the effective
Hamiltonian

$H_{\text{p}}=\left(\begin{array}[]{cc}H_{\mathcal{U}\mathcal{U}}&H_{\mathcal{U}\mathcal{C}}\\\
H_{\mathcal{U}\mathcal{C}}^{\dagger}&H_{\mathcal{C}\mathcal{C}}\end{array}\right),$
(42)

consisting of the blocks

$\displaystyle\left[H_{\mathcal{U}\mathcal{U}}\right]_{\alpha,\beta}$
$\displaystyle=$
$\displaystyle\left<\text{GS}\right|\hat{q}_{\alpha}\hat{H}\hat{q}^{\dagger}_{\beta}\left|\text{GS}\right>,$
$\displaystyle\left[H_{\mathcal{U}\mathcal{C}}\right]_{\alpha,{\left\\{\beta_{i}\right\\}_{m+1}}}$
$\displaystyle=$
$\displaystyle\left<\text{GS}\right|\hat{q}_{\alpha}\hat{H}\left|{\left\\{\beta_{i}\right\\}_{m+1}}\right>,$
$\displaystyle\left[H_{\mathcal{C}\mathcal{C}}\right]_{{\left\\{\alpha_{i}\right\\}_{m+1}},{\left\\{\beta_{i}\right\\}_{m+1}}}$
$\displaystyle=$
$\displaystyle\left<{\left\\{\alpha_{i}\right\\}_{m+1}}\right|\hat{H}\left|{\left\\{\beta_{i}\right\\}_{m+1}}\right>.$
(43)

The dimension of this Hamiltonian is
$N_{\mathcal{U}}+\left(\begin{array}[]{c}N_{\mathcal{C}}\\\
m+1\end{array}\right)$, where $N_{\mathcal{U}}$ is the number of orbitals in
the unoccupied sector. Similar to the ansatz for particle-like excitations, we
make an ansatz for hole-like excitations.

$\left|\text{h},k\right>_{1}=\sum_{\alpha\in\mathcal{O}}u_{\alpha,k}\,\hat{q}_{\alpha}\left|\text{GS}\right>+\sum_{\left\\{\alpha_{i}\right\\}_{m-1}}v_{{\left\\{\alpha_{i}\right\\}_{m-1}},k}\left|{\left\\{\alpha_{i}\right\\}_{m-1}}\right>.$
(44)

We note that variational trial states with limited numbers of particles in the
unoccupied sector and holes in the occupied sector have been employed before
in the context of quantum impurity problems [47, 50], but as far as we know,
not to calculate excited quasiparticle states.

To gauge the accuracy of the ansatz, we benchmarked it against NRG results.
Details can be found in Appendix A. This allowed us to identify SIAM
parameters that produce a well-developed Kondo regime, and in which the ansatz
gives accurate results at an affordable computational cost.

##  IV Numerical Results

We now present numerical results that we obtained by applying the tools
developed in Section III.1 to eigenstates of the SIAM that comprise a single
quasiparticle on top of the ground state. We obtained the latter by means of
NRG, or more approximately, using the ansatz (41).

###  IV.1 The bare electron occupation probability $p_{n}$ and the DOS
difference $\delta\rho$.

Figure 1: Left panel: The probaility $p_{n}$ (27) to detect only a bare
electron above the ground state for the symmetric SIAM with $\Delta=0.09D$ and
$U=0.6D$. Right panel: The difference $\delta\rho$ (31) between the density of
states of the single-quasiparticle Hamiltonian (40) and the single-particle
density of states of the actual Hamiltonian, calculated from the same data.
Results are shown for all states with $p_{n}>0.5$. The grey data represent NRG
results for Wilson chains of lengths of $L=2\Omega+1$ with
$\Omega=1,8,\ldots,20$ and $\Lambda=2$. The black solid curve in the left
panel connects the data points corresponding to the lowest energy states with
one particle on top of the ground state, at different $L$. In the right panel,
the states connected by this black line were omitted. The dashed vertical line
indicates the Kondo temperature (3). NRG results at representative $L$ are
highlighted as indicated in the legend. Results obtained with the Ansatz, on
the same Wilson chain as used in NRG, are plotted with the symbols indicated
in the legend. Figure 2: Hybridization $\Delta$ times the density of states
difference $\delta\rho$ (31), for different $U$ and three values of $\Delta$,
versus the corresponding quasiparticle weight $z$. A dashed line indicates the
estimate (32).

We have performed numerical investigations on the particle-hole symmetric SIAM
using NRG and the ansatz (41). In this subsection we show results for
$\Delta=0.09D$ and $U=0.6D$, which correspond to $T_{\text{K}}=0.013D$. In
Figure 1 we present results for the bare electron occupation probability
$p_{n}$ (27) and the difference in single-particle DOS (31) between the
quasiparticle Hamiltonian (40) and the actual system. We present NRG results
for different Wilson chain lengths $L$ (infrared energy cut-offs) for the same
SIAM parameters. For each $L$, NRG could find the lowest $n_{L}$ single-
particle excitations, with $n_{L}$ between 2 and 4. (Higher excited single-
quasiparticle excitations get decimated during the renormalization flow.) In
order to access higher energy single-quasiparticle excitations, we employed
the ansatz (41) on the same Wilson chain Hamiltonian as NRG. We used a
correlated sector with six spin up and six spin down orbitals, that is half-
filled for the ground state. At each $L$ we see that there is a slight
decrease in $p_{n}$ going from the lowest single-quasiparticle excitation to
the second lowest. There is a steeper decrease of $p_{n}$ going from the
second to the third, third to fourth excitation, and so on. The fact that the
lowest single-quasiparticle excitation behaves differently than the rest has
to do with the energy discretization employed by NRG, in which the inverse of
the level spacing between single-particle states of the Wilson chain has a
similar behaviour. At fixed $L$, $1-p_{n}$ is proportional to $E_{n}$, for
$n=2,3,4$, with an $L$-independent proportionality constant. For the lowest
single-quasiparticle state ($n=1$) at different $L$, we also find
$1-p_{1}\propto E_{1}$. We see good agreement between NRG and the ansatz. When
$L$ is varied, the $1-p_{n}\propto E_{n}$ behaviour breaks down when the
infrared cut-off scale is $\sim 0.1T_{\text{K}}$. On the other hand, when $L$
is held fixed such that the corresponding infrared cutoff scale is much
smaller than $T_{\text{K}}$, the $1-p_{n}\propto E_{n}$ behaviour persists up
to nearly $T_{\text{K}}$. We will explore the reason behind these contrasting
behaviours further below in Figure 5.

We obtain the DOS difference $\delta\rho$ by multiplying $1-p_{n}$ by the
density of single particle energies of the Wilson chain
$=1/\text{ln}(\Lambda)E_{n}$, cf. (31). To avoid having to deal with the
irregularity of the level spacing at the lowest single-quasiparticle state, we
exclude it from the presented results in the right panel of Figure 1. At
$E_{n}\ll T_{\text{K}}$, $\delta\rho$ is constant. NRG is sufficient to show
that $\delta\rho$ increases above the low energy plateau, if the system’s
infrared cut-off energy scale is larger than $0.1T_{\text{K}}$. However, the
ansatz is required to study single-quasiparticle excitations close to
$T_{\text{K}}$ when the system’s infrared cut-off scale is much lower than
$T_{\text{K}}$. The ansatz indicates that in this case, $\delta\rho$ remains
constant up to nearly $T_{\text{K}}$, after which it increases. An increase of
$\delta\rho$, rather than a decrease, validates the intuition articulated in
Section III, that $\delta\rho$ reveals deviations from perfect Fermi liquid
behaviour when the system is probed at energies approaching $T_{\text{K}}$.

In Section III, we predicted (32) that the value of $\delta\rho$ on the low
energy plateau should be $(1/z-1)/\pi\Delta$ where $z$ is the quasiparticle
weight. We test this prediction as follows. Using NRG, we calculate the
plateau value of $\delta\rho$ for different $U$ and $\Delta$. For each
calculation, we also determine $z$ by fitting the low energy spectrum to that
of a non-interacting SIAM ($U=0$), with renormalized $\tilde{V}=\sqrt{z}V$
(40). The results of this analysis are shown in Figure 2. We see very good
agreement between NRG results and our prediction. This confirms the validity
of the interpretation given to the quasiparticle weight $z$ within the wave
function picture of local Fermi liquid theory, namely that it is connected to
the bare electron occupation probability via (33).

Figure 3: Left panel: The quasi-particle weight $z$, for the symmetric SIAM
with $\Delta=0.09D$ and $U=0.6D$. Results were obtained with NRG on Wilson
chains of lengths of $L=2\Omega+1$ with $\Omega=7,8,\ldots,20$ and
$\Lambda=2$. Different colors represent different $L$. Each plot marker
represents a single-quasiparticle excitation. Dots represent the estimate (37)
while open disks represent the estimate (38). The dashed line was obtained
from $z=\tilde{V}^{2}/V^{2}$, where $\tilde{V}$ was determined by fitting the
spectrum of the quasiparticle Hamiltonian (40) to the single-particle spectrum
found in NRG, at $L=29$. Rigth panel: Square of the quasiparticle wave
functions $\psi_{j}$, with unit-normalization, for the symmetric SIAM with
$\Delta=0.09D$ and $U=0.6D$ on a Wilson chain of length of $L=29$. (Here
$j=-1$ refers to the $d$-orbital.) The lowest four single-quasiparticle wave
functions are shown. The dots were calculated with NRG. The open circles
represent the wave functions calculated using the single-quasiparticle
Hamiltonian (40). The quasiparticle weight represented by the dashed line in
the left panel was employed. The stars show the single-particle wave functions
of the non-interacting $\Delta=0.09D$, $U=0$ problem, for comparison.

###  IV.2 Quasiparticle weight and wave function

In Section III, we made a connection (39) between the matrix elements
$\left<\text{p},n\sigma\right|\hat{d}_{\sigma}^{\dagger}\left|\text{GS}\right>$
and
$\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
on the one hand, and eigenvectors of a single-quasiparticle Hamiltonian (40)
on the other. This rests on the assumption that a quasiparticle weight $z$,
approximately independent of excitation energy, can be defined through either
(37) or equivalently (38). In Figure 3 we validate this assumption. In the
left panel, we compare several quantities calculated with NRG, for the
particle-hole symmetric SIAM. The dashed line was obtained by fitting the
spectrum of single-quasiparticle Hamiltonian (40) for given quasiparticle
weight $z$, to the low energy spectrum obtained by NRG. (This is the same
procedure as was followed to obtain $z$ in Figure 2.) The symbols represent
the two equivalent expressions (37) and (38) evaluated for every low energy
eigenstate found by NRG, that has $p_{n}>0.5$. Results for Wilson chains of
different lengths $L$ are shown together. There is a clear agreement between
the conventional definition of the quasiparticle weight, represented by the
dashed line, and the symbols representing (37) and (38). However, there also
appear some outliers. Given any exact eigenstate, the quantity resperesnted by
the dots should equal the quantity represented by the disks, regardless of
whether they equal the energy-independent quasiparticle weight. At the lowest
energies (corresponding to the longest Wilson chains) the dots and the open
disks do not lie on top of each other. This indicates that these low energy
deviations are an artefact due to the inherent numerical instability of
calculating matrix elements of irrelevant operators near the infrared fixed
point. Another branch of outliers are indicated by a solid grey curve. Tracing
this curve from right to left, it eventually merges with the correct value of
$z$. We have identified that this branch is associated with hybridization on
the Wilson chain between the second single-quasiparticle eigenstate and the
lowest excited state containing two quasiparticles and one quasihole. During
the renormalization flow, these two many-body levels approach each other.
Initially there is some level repulsion and associated hybridization. However,
eventually (for sufficiently long Wilson chains) the two become degenerate, so
that a well-developed single-quasiparticle excitation can be distinguished.
This is the point where the branch merges with the plateau indicated by the
dashed line. This branch therefore does not represent a failure of the theory
developed in Section III, but rather a genuine interaction between
quasiparticles at intermediate energies. In other words, if the finite Wilson
chain with $L\leq 31$ could be weakly connected to leads, the considered
outliers would represent a scattering resonance, due to an accidental
degeneracy, in which an incoming electron really does scatter into two
electrons and a hole.

In the right panel of Figure 3, we plot the square of
$\psi_{dn}=\left<\text{p},n\sigma\right|\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>/\sqrt{z}$
and
$\psi_{jn}=\left<\text{p},n\sigma\right|\hat{c}^{\dagger}_{j\sigma}\left|\text{GS}\right>$
(39) with the matrix elements calculated by means of NRG, and $z=0.197$ as
indicated by the dashed line in the left panel. Here $j$ is the site-index
along the Wilson chain, and $j=-1$ corresponds to the $d$-orbital. (We
normalized $\psi_{n}$ to unity.) We compare this to the single-particle
eigenvectors of the quasiparticle Hamiltonian (40), at the same quasiparticle
weight. Results are shown for the four lowest quasiparticle states of the
$L=29$ chain that are indicated by arrows in the left panel. For comparison,
we also show the single-particle wave functions of the non-interacting ($U=0$,
$z=1$) chain, with the same bare hybridization $\Delta=0.09D$. We see near
perfect coincidence between the NRG results and the single-quasiparticle wave
functions at $z=0.197$. We also note that amplitudes on low energy shells
(large $j$) are nearly independent of $z$. The quasiparticle weight (or
renormalized hybridization) only affects the very small amplitudes on energy
shells $j=0,2,\ldots$ up to the shell at the Kondo scale. This is a
manifestation of the fact that the hybridization term in the Hamiltonian
represents an irrelevant perturbation around the infrared fixed point (albeit
a leading one).

Figure 4: Left panel: NRG results for the dressing $D_{j}$, (34) normalized so
that $\text{Max}(D_{j})=1$, for the symmetric SIAM with $\Delta=0.09D$ and
$U=0.6D$. Results were obtained with NRG on Wilson chains of lengths as
indicated in the legend and $\Lambda=2$. At each length, the lowest four
single particle states are plotted. They lie on top of each other, except for
$L=23$. The dashed vertical line indicates the position of Kondo temperature
(3) relative to the energy shells $j$ of the Wilson chain. Right panel: Ansatz
results for the dressing $D_{j}$, (34) normalized so that
$\text{Max}(D_{j})=1$, for the same symmetric SIAM as in the left panel.
Results were obtained on Wilson chains of length $L=41$. Different data sets
correspond to different states calculated at the same $L$. Results are shown
for the two highest single quasiparticle states below $T_{\text{K}}$ and the
lowest one above $T_{\text{K}}$. The black data represents the NRG result for
the lowest single quasiparticle state at $L=41$.

###  IV.3 Mode-resolved dressing.

Next we present results for the mode-resolved dressing measure $D_{jn}$, cf.
(34). We work on the Wilson chain, in the site basis $j$, i.e. we resolve
dressing per energy shell. We have already investigated the total amount of
dressing of a single-quasiparticle state, by calculating $p_{n}$. Here we
therefore normalize each $D_{jn}$ such that
$\text{Max}\\{D_{jn}|j=d,0,1,\ldots,L-1\\}=1$. In the left panel of Figure 4,
we present NRG results for $|D_{jn}|^{2}$, for the lowest four single-
quasiparticle states $n=1,2,3,4$, for Wilson chains of different length. For
chains long enough that the infrared cut-off scale is less than $\sim
0.1T_{\text{K}}$, we see that $D_{jn}$ has the same shape for $n=1,2,3,4$, and
for different $L$. This is a non-trivial observation, since the associated
$\bm{y}_{n}$, $n=1,2,3,4$ used to calculate $D_{jn}$ are clearly distinct. It
means that low energy quasiparticles associated with different single-particle
orbitals are dressed in the same form. We find that $D_{jn}$ is peaked at an
energy shell above $T_{\text{K}}$. This reflects the fact that one has to
probe the system at ultraviolet scales to distinguish a quasiparticle from a
bare electron. Regarding $L$-dependence, we see that $D_{jn}$ starts changing
when the Wilson chain is so short that the tail of $D_{jn}$ extends to the
lowest energy shell. This behaviour confirms that the proposed measure really
does give the correct mode-resolved picture: as the length $L$ of the Wison
chain is decreased, the nature of quasiparticles starts changing when the low
energy dressing cannot fit on the Wilson chain any more. In the right panel of
Figure 4, we investigate the behaviour of single-quasiparticle excitations
close to $T_{\text{K}}$ for a chain of length $L=41$, for which the infrared
cut-off scale safely lies well below $T_{\text{K}}$. Results obtained with the
ansatz are presented. The profile of $D_{jn}$ for low-energy excitaions
obtained in the left panel via NRG is shown as a black dashed line. The two
highest single-quasiparticle excitations below $T_{\text{K}}$ still fit this
profile very well. For the lowest excitation above $T_{\text{K}}$ we see the
profile starting to shift, and becoming significant on chain sites
$j=1,3,\ldots$ where previously it was negligible. (For higher excited states,
the local Fermi liquid picture rapidly fails.)

Figure 5: Left panel: Comparison of NRG results for the dressing $D_{j}$, (34)
normalized so that $\text{Max}(D_{j})=1$, for the symmetric SIAM for different
$\Delta$ and $U$. The three data sets are for (a) $\Delta=0.09D$, $U=1.4D$ (b)
$\Delta=0.09D$, $U=2.4D$, and (c) $\Delta=0.0009D$, $U=0.006D$. The
corresponding Kondo temperatures $T_{\text{K}}$ are respectively $4.7\times
10^{-4}D$, $6.6\times 10^{-6}D$, and $1.5\times 10^{-4}D$. The lowest state
with a single quasiparticle on top of the ground state is plotted in each
case. Results were obtained with NRG on Wilson chains with $\Lambda=2$. The
dashed vertical lines indicate the positions of Kondo temperatures (3)
relative to the energy shells $j$ of the Wilson chain. Right panel: Same data
as left panel, but plotted in log-scale on the vertical axis, and shifted
horizontally so the maxima occur at zero. The black line is proportional to
$e^{-j/\pi}$.

We now return to the low-energy behaviour of (normalized) $D_{jn}$. In the
left panel of Figure 4, we show results for the lowest eigenstate with one
quasiparticle on top of the ground state ($n=1$) for different combinations of
$\Delta$ and $U$, all chosen such that $T_{\text{K}}\ll D$. We see that the
general shape is rather similar for vastly different parameters. (Of course,
if we did not normalize, the overall amplitude would depend strongly on
$T_{\text{K}}$.) In the right panel, we replot the same data in log-scale on
the vertical axis, and shifted horizontally to line up the maxima of the
different data sets at $j=0$. At small $j$ (i.e. in the ultraviolet) we see
exponential behaviour that depends on $\Delta$, but not on $U$. In the
infrared (large $j$), we obtain universal behaviour $D_{jn}\sim e^{-j/2\pi}$.
Translating from shell index to energy and then to momentum, this implies
$D_{jn}\sim k^{1/\pi\text{ln}(\Lambda)}$, where $k$ measures momentum relative
to the Fermi momentum. To extrapolate NRG results obtained with logarithmic
discretization, to the original continuum system, one is supposed to take the
$\Lambda\to 1$ limit, but here the limit does not exist. It therefore seems
that the low momentum- , or long-distance real space behaviour of the dressing
profile $D_{kn}$, is beyond the reach of NRG.

###  IV.4 Real-space quasiparticle structure

Figure 6: Left panel: A representative position space quasiparticle wave
function (the 12’th excited single quasiparticle state) for the symmetric SIAM
with $\Delta=0.09D$ and $U=0.2D$ calculated on a regular 1D lattice of 9601
sites, with the $d$-orbital coupled to the central site ($j=0$). The solid
curves represents the result obtained with the Ansatz, while the dashed curves
represent the corresponding eigenvector of the single-quasiparticle
Hamiltonian (40), with the quasiparticle weight estimated as indicated by the
dashed line in the left panel of Figure (3). Each curve was obtained by
connecting amplitudes on every fourth site, starting at sites 0,1,2 or 3.
Right panel: The phase shift. Results extracted from the Ansatz, using (46)
are indicated by dots. The open disks were obtained from the single-
quasiparticle Hamiltonian on the regular 1D lattice (40), with the
renormalized hybridization $\tilde{V}$ inferred from NRG calculations similar
to those in Subsection IV.2.

To access the properties of single-quasiparticle excitations in real space, we
therefore proceed to study the structure of single-quasiparticle states on a
regular energy grid (22) of $L=2\Omega+1$ levels with $\Omega=2400$, using
natural orbital methods. We take $\Delta=0.09$ as above, and work at the
relatively weak interaction strength $U=0.2D$. According to NRG, the
quasiparticle weight $z=0.77$ and the Kondo temperature is $0.064D$. We have
benchmarked the ansatz at this point in parameter space by calculating the
Kondo resonance of the local density of states of the $d$-orbital. We find
that the ansatz is accurate for excitation energies up to $\sim
0.3T_{\text{K}}$. (See Appendix A. For fixed size of the correlated sector,
the ansatz on the large linear grid seems to be somewhat less accurate than on
the logarithmic grid. On the linear grid and with $U=0.6$ for instance, errors
in the $10\%$ range are found for quantities such as the quasiparticle weight.
This is why we chose to work with $U=0.2D$ here.)

Previously, in the right panel of Figure 1, we confirmed that matrix elements
$\left<\text{p},n\sigma\right|\hat{d}_{\sigma}^{\dagger}\left|\text{GS}\right>/\sqrt{z}$
and
$\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
are eigenvectors of the single-quasiparticle Hamiltonian (40) for the Wilson
chain. Now we perform the equivalent analysis for the large regular energy
grid, to establish whether the equality holds at single real-space lattice
site resolution.

In the left panel of Figure 6, we plot
$y_{kn}=\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$,
Fourier transformed to real space, i.e.

$\tilde{y}_{jn}=\frac{1}{\sqrt{2\Omega+1}}\sum_{k=1}^{2\Omega+1}\cos\left[\frac{\pi
j(2k-1)}{4\Omega+2}\right]y_{kn}.$ (45)

We can check this against the predictions of local Fermi liquid theory as
follows. First we extract the renormalized $\tilde{V}$ by fitting the low
energy NRG spectrum to the Wilson chain version of the single-quasiparticle
Hamiltonian (40). Using this $\tilde{V}$, we compute the eigenvectors
$\bm{\psi}_{n}$ of the single-quasiparticle Hamiltonian (40) discretized on
the regular energy lattice. We Fourier transform and suitably normalize the
bulk amplitudes to obtain $\tilde{\psi}_{jn}$. This is compared to
$\tilde{y}_{jn}$. As a representative example, we look at the 12’th excited
single-quasiparticle state (i.e. $n=13$), and see near-perfect agreement.

Except for the first few sites, $\tilde{y}_{jn}$ is a shifted sinusoidal wave

$A_{n}\cos(\left|j\right|k_{n}-\phi_{n}).$ (46)

We extract the phase shift $\phi_{n}$ by fitting (46) to $\tilde{y}_{jn}$, and
remembering that the boundary condition $\tilde{y}_{2\Omega+2,n}=0$ imposes
the quantization

$k_{n}=k_{\text{F}}+\frac{\pi n-\phi_{n}}{2\Omega+1},$ (47)

of the wave number. We can compare this to the predictions of local Fermi
liquid theory by similarly extracting the phase shift of $\tilde{\psi}_{jn}$.
In the right panel of Fig. 6, we see that the two estimates agree very well up
to $E_{n}=0.02D\sim 0.3T_{\text{K}}$. We conclude that indeed the matrix
elements $y_{dn}/\sqrt{z}$ and $\tilde{y}_{jn}$ represent a real-space
eigenvector of the single-quasiparticle Hamiltonian (40).

If a non-interacting quasiparticle picture holds, we should equivalently be
able to extract the phase shift from the quantization of the wave number and
the dispersion relation as follows

$E_{n}=\frac{2\Omega+1}{\pi\Omega}(k_{n}-k_{\text{F}})\implies\phi_{n}=\pi\left(n-\frac{\Omega
E_{n}}{D}\right).$ (48)

We have checked that applying this method to results obtained with the ansatz
gives the same result as the dots in the right panel of Figure 6.

Figure 7: Main panel: The dressing measure $D_{j}$ Fourier transformed to real
space as in (45), for the symmetric SIAM calculated on a regular 1D lattice of
9601 sites, with the $d$-orbital coupled to the central site ($j=0$) at
$\Delta=0.09D$ and two values of $U$. The curves are envelopes to the actual
data, and were obtained by plotting the result on each fourth site, $j=4l$,
$l=0,1,\ldots$. In each case the lowest single-quasiparticle state as obtained
by the Ansatz, was used to calculate $D_{j}$. The dashed lines represent a fit
to a power law $(a/j)^{b}$. For $U=0.2$ we found a best fit $b=0.97$ and for
$U=0.6$ we found $b=0.93$. Inset: The solid curve represents
$\sqrt{L}\tilde{D}_{j}\times(j/a)^{b}$ for $U=0.2D$, with $a=0.43$ and
$b=0.97$ determined in the main panel by fitting $(a/j)^{j}$ to
$\sqrt{L}\tilde{D}_{j}$. The dashed curve represents
$\sqrt{L}\tilde{D}_{j}\times(j/a)$ for $U=0.2D$, with $a=0.53$ determined by
fitting $a/j$ to $\sqrt{L}\tilde{D}_{j}$. From this we deduce that the data is
sufficient to reveal that the powerlaw is $\tilde{D}_{j}\sim j^{-0.97}$ rather
than $\tilde{D}_{j}\sim j^{-1}$ at $U=0.2D$.

Finally, we study the dressing profile in real space. That is, we compute
$D_{kn}$ for the system discretized on a regular energy grid, and Fourier
transform as $(\ref{eqfourier})$ to obtain the site-resolved dressing
$\tilde{D}_{jn}$. The results presented above, particularly in the left panel
of Figure 5, suggests that $\tilde{D}_{jn}$ is a good measure of how different
from a bare electron the quasiparticle excitation looks on real-space lattice
site $j$. In Figure 7, we show results at $U=0.2D$, which are quantitatively
accurate, and results at $U=0.6D$, that may contain small but non-negligible
errors. We compare the lowest eigenstate with one quasiparticle on top of the
ground state in each case, for which the ansatz should give more accurate
results than for higher excited states. We therefore think that the error in
the $U=0.6D$ results are smaller than the $10\%$ errors observed for
quantities such as the quasiparticle weight, that tests the ansatz at higher
excited states. As expected we see that the dressing decays as we move away
from the impurity, and that increasing $U$, which decreases $T_{\text{K}}$,
leads to more significant dressing. At large distances, the dressing profile
seems to satisfy a power law $\sim j^{-b}$ with $b$ discernibly different from
$1$. For $U=0.2D$, we find $b=0.97$ while for $U=0.6D$ we find $b=0.93$. The
decrease of the exponent at larger $U$ is consistent with the expectation that
the particle-hole pairs that dress the bare electron are less confined, the
lower the Kondo temperature, but it has to be kept in mind that the $U=0.6D$
results are less accurate than the $U=0.2D$ results. We are therefore less
confident about the exact value of the exponent at $U=0.6D$ than at $U=0.2D$.
Nonetheless, evidence of a nontrivial power law is clear, and could not have
been obtained with any method other than the ansatz, as far as we can see.

##  V Summary and conclusion

There are two distinct single-particle-like structures associated with dynamic
fermionic quantum impurity problems. There are the quasiparticles that provide
a local Fermi liquid description of low energy excitations in terms of
independent effective degrees of freedom. Then there is also the natural
orbital basis, associated with the bare electronic degrees of freedom, which
is known to provide an economical description of ground state correlations.
The work presented here was inspired by the question: “Can a synergy of these
two pictures shed light on how bare electrons are dressed to form the
quasiparticle excitations of a local Fermi liquid?”. We succeeded in
establishing the following. The quasiparticle weight $z$, and the effective
quasiparticle Hamiltonian

$\hat{H}_{\text{qp}}=\sqrt{z}V\sum_{\sigma}\left(\hat{d}^{\dagger}_{\sigma}\hat{\psi}_{0\sigma}+\hat{\psi}_{0\sigma}^{\dagger}\hat{d}_{\sigma}\right)+\sum_{k=1}^{L}\sum_{\sigma}\varepsilon_{k}\hat{c}_{k\sigma}^{\dagger}\hat{c}_{k\sigma},$
(49)

of the particle-hole symmetric Single Impurity Anderson Model have meaning not
only in the infinite system, where they determine the scattering matrix of
electrons at low energy, but also in the large finite system. If
$\left|\text{p},n\sigma\right>$ is a low energy single-quasiparticle excited
state of the SIAM, then to a good approximation
$\psi_{dn}=\left<\text{p},n\sigma\right|\hat{d}^{\dagger}_{\sigma}\left|\text{GS}\right>/\sqrt{z}$
and
$\psi_{kn}=\left<\text{p},n\sigma\right|\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$
solve the single-particle Schrödinger equation (40) associated with (49).
Furthermore, if we define $y_{dn}=\sqrt{z/p_{n}}\psi_{dn}$ and
$y_{kn}=\psi_{kn}/\sqrt{p_{n}}$ with $p_{n}$ such that
$(y_{dn},y_{1n},\ldots,y_{Ln})$ is normalized to unity, then
$\hat{a}_{n\sigma}^{\dagger}=y_{dn}\hat{d}_{\sigma}^{\dagger}+\sum_{k=1}^{L}y_{kn}\hat{c}_{k\sigma}^{\dagger}$
represents the best bare electron approximation to the quasiparticle
associated with $\left|\text{p},n\sigma\right>$. We call $p_{n}$ the bare
electron occupation probability because it is the occupation probability of
the $\hat{a}_{n\sigma}$-orbital when the system is in state
$\left|\text{p},n\sigma\right>$. In a system where the single-quasiparticle
level spacing is $\delta E$ in the vicinity of excitation
$\left|\text{p},n\sigma\right>$, at low energies, $p_{n}$ is nearly
independent of energy and related to the quasiparticle weight $z$ by
$z=1/(1+p_{n}\pi\Delta/\delta E)$. We define the dressing of the bare electron
as
$\left|\text{dressing}\right>=\left|\text{p},n\sigma\right>-\sqrt{p_{n}}\hat{a}_{n\sigma}^{\dagger}\left|\text{GS}\right>$,
and obtain a mode-resolved measure $D_{kn}$ of dressing by taking the overlap
with $\hat{c}_{k\sigma}^{\dagger}\left|\text{GS}\right>$. This measure can be
computed from knowledge of $p_{n}$, $(y_{dn},y_{1n},\ldots,y_{Ln})$, and the
ground state covariance matrix $C(\left|\text{GS}\right>)$ whose eigenvectors
define the natural orbital basis. Where possible, we used NRG to study these
quasiparticle-structure-related quantities. However, NRG is inadequate to
study more than the first few single-quasiparticle excitations, or to study
structure in real space. We overcame these limitations by constructing an
ansatz for single-quasiparticle excitations in terms of the natural orbital
basis. The important features of this ansatz are (1) that it provides an
explicit expression in terms of bare electronic degrees of freedom, for a
quasiparticle excitation and (2) that it can be implemented directly on a
microscopic model involving an impurity coupled to electrons on a lattice of
several thousand sites, rather than on a continuum model rediscretized on the
logarithmic energy grid used in NRG. The ansatz is accurate at affordable
computational cost in a regime of well-developed correlations. We further
showed that the difference $\delta\rho$ between the density of single-particle
states of $\hat{H}_{\text{qp}}$ and the actual many-body system is an
increasing function of energy at energies approaching or equal to the Kondo
temperature. This is a “beyond ideal Fermi-liquid theory” effect. We finally
showed that the position resolved dressing $\tilde{D}_{jn}$ decays with a
nontrivial powerlaw at large distances.

We envisage that the work we presented here could provide a foundation for the
explicit calculation of the single particle scattering matrix for a quantum
impurity imbedded in a non-trivial geometry such as a disordered host or a
mesoscopic electronic device. This is a challenging task because the
quasiparticle weight $z$ which appears in the effective Hamiltonian, is
determined by the many-body state of the system, which cannot be calculated
with methods such as NRG, with its crudely resolved discrete representation of
the host. In the present work, we considered single quasiparticle excitations.
The methods we developed may be a good starting point to investigate the weak
interactions between two or more quasiparticles. Finally, we also suspect that
simple modifications to our ansatz could extend its regime of validity. For
instance, in the present work, we determine the natural orbital basis at the
start of the calculation, using the ground state as a reference. In our own
ongoing work, we are investigating whether it is profitable to adjust the
natural orbital basis itself when dealing with excitations.

###### Acknowledgements.

Insigthful comments and suggestions from Serge Florens are gratefully
acknowledged.

##  Appendix A Benchmarking

We benchmark our proposed methods against known results. To gauge generality,
we test the natural orbital methods on the SIAM and on another model. The
choice of this model is dictated by the demand that it should be sufficiently
different from the SIAM to give independent confirmation of the accuracy of
the method, yet its structure must be similar enough that computer codes
developed for the SIAM can easily be adapted.

###  A.1 The Interacting Resonant Level Model

The Interacting Resonant Level Model (IRLM) is ideal for this purpose. It
describes a band of non-interacting, spinless fermions hybridizing with a
localized orbital. Additionally there is a local interaction between a
particle in the local orbital and one in the adjacent lattice site of the
band. We study the particle-hole symmetric version of the Hamiltonian, which
reads

$\displaystyle\hat{H}=$ $\displaystyle
U\left(\hat{n}_{d}-\frac{1}{2}\right)\left(\hat{n}_{0}-\frac{1}{2}\right)+V\left(d^{\dagger}\psi_{0}+\psi_{0}^{\dagger}d\right)$
$\displaystyle+\sum_{k=1}^{L}\varepsilon_{k}\hat{c}_{k}^{\dagger}\hat{c}_{k}.$
(50)

Here $\hat{d}^{\dagger}$ creates a particle in the localized orbital and
$\hat{\psi}_{0}^{\dagger}$ creates a particle on the lattice site adjacent to
the impurity. Furthermore

$\hat{n}_{d}=\hat{d}^{\dagger}\hat{d},~{}~{}~{}\hat{\psi}_{0}=\frac{1}{\sqrt{L}}\sum_{k=1}^{L}\hat{c}_{k},~{}~{}~{}\hat{n}_{0}=\psi_{0}^{\dagger}\psi_{0}.$
(51)

We take the Fermi energy to be at zero, in the middle of the band. As in the
SIAM, hybridization is quantified by the spectral density $\Delta=\pi
V^{2}/2D$ at the Fermi level, of the operator $V\hat{\psi}_{0}$ in the
infinite system uncoupled from the impurity. In the $\Delta\ll D$ regime, and
for $U<0$, the IRLM hosts the same physics as the Kondo model, in which there
is a local exchange interaction between the magnetic moment of a spin-1/2
impurity and nearby conduction band electrons. (There is a mapping between the
spin degrees of freedom in the Kondo model and the IRLM, the charge degrees of
the Kondo model being spectators that do not couple to the impurity [51, 52,
53].) The quantum phase transition seen in the Kondo model when the exchange
interaction switches from ferro- to antiferromagnetic, translates to a
critical point in the IRLM at $U\sim-D$. Increasing $U$ beyond the critical
point corresponds to increasing the Kondo temperature in the antiferromagnetic
regime.

###  A.2 Ansatz 2

Apart from the study of Fermi liquid quasiparticles, we are also interested in
establishing if there are any fundamental limitations to applying natural
orbital methods to excited state problems. Even though the ansatz formulated
in Section III will prove sufficiently accurate for our purposes, we introduce
an improved, but more computationally expensive Ansatz 2 for the purpose of
comapring to the ansatz of Section III, which here we refer to as Ansatz 1.

Ansatz 2 is obtained by relaxing the assumption that when a particle ends up
in the unoccupied sector, the particles in the correlated sector are
undisturbed. Instead, a particle in the unoccupied sector may be associated to
an arbitrary configuration of $m$ particles in the correlated sector. Thus we
consider $N+1$ particle states of the form

$\displaystyle\left|\text{p},k\right>_{2}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\in\mathcal{U}}\sum_{{\left\\{\alpha_{i}\right\\}_{m}}}u_{\alpha,{\left\\{\alpha_{i}\right\\}_{m}},k}\,\hat{q}_{\alpha}^{\dagger}\left|{\left\\{\alpha_{i}\right\\}_{m}}\right>$
(52)
$\displaystyle+\sum_{\left\\{\alpha_{i}\right\\}_{m+1}}v_{{\left\\{\alpha_{i}\right\\}_{m+1}},k}\left|{\left\\{\alpha_{i}\right\\}_{m+1}}\right>.$

The optimal coefficients $u_{\alpha,{\left\\{\alpha_{i}\right\\}_{m}},k}$ and
$v_{{\left\\{\alpha_{i}\right\\}_{m+1}},k}$ are found in the same manner as
before. They are thus seen to be eigenvectors of the effective Hamiltonian

$\tilde{H}_{\text{p}}=\left(\begin{array}[]{cc}\tilde{H}_{\mathcal{U}\mathcal{U}}&\tilde{H}_{\mathcal{U}\mathcal{C}}\\\
\tilde{H}_{\mathcal{U}\mathcal{C}}^{\dagger}&H_{\mathcal{C}\mathcal{C}}\end{array}\right),$
(53)

where the blocks $\tilde{H}_{\mathcal{U}\mathcal{U}}$ and
$\tilde{H}_{\mathcal{U}\mathcal{C}}$ are now given by

$\displaystyle\left[\tilde{H}_{\mathcal{U}\mathcal{U}}\right]_{\left(\alpha,{\left\\{\alpha_{i}\right\\}_{m}}\right),\left(\beta,{\left\\{\beta_{i}\right\\}_{m}}\right)}$
$\displaystyle=$
$\displaystyle\left<{\left\\{\alpha_{i}\right\\}_{m}}\right|\hat{q}_{\alpha}\hat{H}\hat{q}^{\dagger}_{\beta}\left|{\left\\{\beta_{i}\right\\}_{m}}\right>,$
$\displaystyle\left[\tilde{H}_{\mathcal{U}\mathcal{C}}\right]_{\left(\alpha,{\left\\{\alpha_{i}\right\\}_{m}}\right),{\left\\{\beta_{i}\right\\}_{m+1}}}$
$\displaystyle=$
$\displaystyle\left<{\left\\{\alpha_{i}\right\\}_{m}}\right|\hat{q}_{\alpha}\hat{H}\left|{\left\\{\beta_{i}\right\\}_{m+1}}\right>,$

while the block $H_{\mathcal{C}\mathcal{C}}$ is the same as before (43). The
block $\tilde{H}_{\mathcal{U}\mathcal{U}}$ has dimension
$N_{\mathcal{U}}\times\left(\begin{array}[]{c}N_{\mathcal{C}}\\\
m\end{array}\right)$, which is considerably larger than the corresponding
block (43) in the effective Hamiltonian associated with Ansatz 1.

###  A.3 Local density of states of the $d$-orbital.

We will explore the accuracy of Ansatz 1 and Ansatz 2 by considering the the
local density of states (LDOS) of the $d$-orbital, a quantity that is
sensitive to single-particle excitations on top of the ground state. The two
ansätze that we study were designed to be accurate in the Fermi liquid regime,
i.e. below the Kondo scale. We intentionally choose the LDOS as a benchmark,
not only to demonstrate the accuracy of the ansätze in the expected regime,
but also to explore how they break down. We will see that the LDOS is
accurately reproduced at frequencies up to $\sim T_{\text{K}}$, and start
breaking down at higher frequencies, where the local Fermi liquid picture does
not apply. We stress that the ansätze studied here, are not being presented to
compete with existing methods to compute the LDOS of the $d$-orbital. Instead,
they come into their own right when quasiparticles living in the bulk, rather
than impurity properties, are considered.

For a finite system, in which the band has $L$ orbitals, the LDOS of the
$d$-orbital is given by

$\displaystyle A_{L}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\text{Re}\,\int_{0}^{\infty}dt\,e^{i\omega
t}\left<\text{GS}\right|\left\\{\hat{d}_{\sigma}^{\dagger}(0),\hat{d}_{\sigma}(t)\right\\}\left|\text{GS}\right>$
$\displaystyle=$
$\displaystyle\sum_{n}\left|\left<\text{GS}\right|\hat{d}^{\dagger}_{\sigma}\left|\text{h},n\sigma\right>\right|^{2}\delta(E_{\text{h},n}-E_{\text{GS}}+\omega)$
$\displaystyle+\sum_{n}\left|\left<\text{GS}\right|\hat{d}_{\sigma}\left|\text{p},n\sigma\right>\right|^{2}\delta(E_{\text{p},n}-E_{\text{GS}}-\omega).$

For the SIAM, we assume zero magnetic field, so that a single spin species may
be considered. For the IRLM, the spin index $\sigma$ is dropped. The
summations in respectively the second and third lines are over all excited
states with one more or fewer particles than the number of particles in the
ground state. For results obtained in the logarithmically discretized model,
we use a trick called $z$-averaging to mitigate discretization errors at
finite frequency [54]. It involves replacing $\Lambda^{-k}$ in (19) with
$\Lambda^{-k+z}$ for $k=1,2,\ldots$, and averaging LDOS results over $n_{z}$
separate calculations, each performed with a different $z=(\zeta-1/2)/n_{z}$,
$\zeta=1,2,\ldots,n_{z}$. Note that with this discretization, the hopping
amplitudes $t_{j}$ along the Wilson chain (21) have to be calculated
recursively [16], using arbitrary precision arithmetic.

In order to arrive at a result independent of the discretization scheme, and
that is near the LDOS in the thermodynamic limit, the delta-peaks of the
finite system LDOS must be broadened. We apply the Gaussian broadening
prescription that is standard within NRG

$\displaystyle A(\omega)=\int_{-\infty}^{\infty}d\varepsilon
A_{L}(\varepsilon)\Gamma(\omega,\varepsilon),$
$\displaystyle\Gamma(\omega,\epsilon)=\frac{1}{\sqrt{2\pi\eta(\varepsilon)}}e^{-(\omega-\epsilon)^{2}/2\eta(\varepsilon)^{2}}.$
(56)

For results obtained in the logarithmically discretized model, we use a width
that decreases linearly as the Fermi energy is approached,

$\eta_{\text{log}}(\varepsilon)=\frac{2|\varepsilon|}{n_{z}}.$ (57)

For results obtained in the system discretized on a regular energy grid, we
use a constant width equal to half the single-particle level spacing

$\eta_{\text{lin}}(\varepsilon)=\frac{D}{L}.$ (58)

If an STM tip is held near the $d$-orbital and a potential $\omega/e$ is
maintained between tip and sample, the tunnelling current between the tip and
the $d$-orbital is proportional to $A(\omega)$. We note that the LDOS of the
IRLM is still a topic of active study [55].

Here and further below, all NRG results for the LDOS of the $d$-orbital, were
obtained with the density matrix NRG (DM-NRG) algorithm [56] which retains
full information of the ground state when the LDOS is evaluated at an
arbitrary frequency. In order to use the maximum available information about
excited states, we employ the Full Hilbert space method [37]. Thus the sum
rule $\int_{-\infty}^{\infty}d\omega\,A_{L}(\omega)=1$ is satisfied
identically. We imposed particle number conservation during the iterative
diagonalization, and retained up to 128 states per particle-number sector,
which corresponds to several thousand kept states in total per iteration.

###  A.4 Results for the IRLM

Figure 8: Top panel: Local Density of States (LDOS) at the $d$-orbital, for
the Interacting Resonant Level Model (IRLM), at interaction strength $U=-0.5D$
and hybridization $\Delta=0.035D$. The main panel shows results in log-log
scale. The inset shows the same data in linear scale. The solid curve
represents the result obtained with Ansatz 1. The band was discretized on a
regular energy grid with 4801 levels. The dashed curve represents an NRG
calculation. The dotted curve represents the non-interacting ($U=0$) result
for comparison. Vertical lines indicate the eigenvalues of
$H_{\mathcal{C}\mathcal{C}}$ (43). Bottom panel: LDOS of the IRLM at the same
parameters $U=-0.5D$ and hybridization $\Delta=0.035D$ as in the top panel.
The NRG data is the same as in the top panel. Results for Ansatz 2 (solid
curve) were obtained using the same logarithmic discretization and
$z$-averaging as the NRG results. For comparison, results for Ansatz 1 were
also recalculated on the logarithmic grid, with $z$-averaging.

First we present results for the IRLM. We calculate $A(\omega)$ at
$\Delta=0.035D$ and $U=-0.5D$. We compare NRG results to results obtained with
Ansatz 1 on a regular energy grid of 4801 orbitals. Natural orbitals were
calculated using the RGNO algorithm, with a correlated sector containing four
particles in $N_{\mathcal{C}}=8$ orbitals. NRG results were obtained on a
Wilson chain of length $52$ sites, including the $d$-orbital, with
$\Lambda=1.5$, leading to an infrared cutoff energy $\sim 4\times 10^{-5}D$.
NRG results were $z$-averaged over 8 $z$-values. Results are shown in the top
panel of Figure 8, in logarithmic scale, and in linear scale (inset). The peak
around $\omega=0$ is the well-known Kondo resonance. The non-interacting case,
$U=0$, corresponds to the Toulouse point of the anisotropic Kondo model, and
therefore still shows a Kondo resonance (of half-width $\Delta=0.035D$). We
see that at $U=-0.5D$, the resonance is about 10 times narrower. This
significant downward renormalization of the hybridization proves that we are
in the strongly correlated regime. We see that Ansatz 1 implemented on a
regular energy grid reproduces the LDOS well, from the Fermi energy up to the
point where it has decayed to about 1% of its value at the Fermi energy. It
therefore captures low-energy excitations well, even beyond the Kondo scale
(half-width of the Kondo resonance). The eventual breakdown occurs as follows.
The effective Hamiltonian (42) describes the hybridization of a (near)
continuum of bare single-particle excitations in the unoccupied sector with a
discrete spectrum of few-body correlated states in the correlated sector. This
discrete spectrum comprises the eigenvalues of the lower-right block
$H_{\mathcal{C}\mathcal{C}}$ of $H_{\text{p}}$ in (42) (together with the
eigenvalues of the particle-hole conjugate, representing one-hole
excitations). In Figure 8, we plot vertical lines representing this discrete
spectrum. We see that Ansatz 1 hybridizes the two few-body states that have
energies closest above and below the Fermi energy, with the continuum of bare
particle and bare hole excitations, to produce a smooth spectral density in
the vicinity of the Kondo resonance. The remaining few-body correlated states
have energies that lie in the tail of the Kondo resonance. Ansatz 1 does not
produce the correct hybridization of these states with the bare particle and
hole continua, as seen from the abrupt features of the resulting spectral
density in the vicinity the vertical lines in the tails of the Kondo
resonance. While this analysis indicates that Ansatz 1 is sufficiently
accurate to allow us to investigate the structure of local Fermi-liquid
quasiparticles, we would nonetheless like to see if we cannot improve accuracy
at higher frequencies.

As we explained above, Ansatz 2 is too computationally expensive to deploy on
a large regular energy grid. In the bottom panel of Figure 8, we compare
Ansatz 1 and Ansatz 2 to NRG results, where now all calculations are performed
for the IRLM discretized on the logarithmic grid, and then $z$-averaged. We
again use the RGNO algorithm, involving a correlated sector of
$N_{\mathcal{C}}=8$ orbitals, containing 4 particles. The interaction strength
$U$ and hybridization $\Delta$ are the same as in the top panel. As might be
expected, $z$-averaging, together with the increased broadening (57) at
intermediate energies, compared to (58), smoothes the spurious abrupt features
in the LDOS obtained by means of Anzatz 1. Nonetheless, in the bottom panel of
Figure 3, Ansatz 1 still shows deviations from the NRG result that stand in
one to one correspondence with those seen in the top panel, for
$|\omega|>0.025D$. Ansatz 2 significantly mitigates the lowest frequency
deviations and remains accurate up to $|\omega|\sim 0.1D$, describing the LDOS
of the IRLM well over almost three decades. The benchmarking exercise thus
proves that natural orbital methods are in principle suitable to study low
energy excitations of fermionic quantum impurity models. We therefore proceed
to our primary system of interest, the Single Impurity Anderson Model (SIAM).

###  A.5 Results for the SIAM

As in Section V, we study the particle-hole symmetric version of the model (1)
where $\varepsilon_{d}=-U/2$. The LDOS of the $d$-orbital shows a resonance
with a width $\tilde{\Delta}=z\Delta$, comparable to $T_{\text{K}}\ll U$,
associated with spin fluctuations. Without the Coulomb repulsion
$U\hat{n}_{d\uparrow}\hat{n}_{d\downarrow}$, this Kondo resonance is absent,
and the LDOS has two peaks at $\omega=\pm\varepsilon_{d}$ instead. In the
presence of Coulomb repulsion, these turn into shoulders at the base of the
Kondo resonance, or even (broad) side peaks, at a scale $\varepsilon_{d}=-U/2$
and $U+\varepsilon_{d}=U/2$, associated with charge fluctuations. We take
$\Delta=0.09D$ and $U=0.6D$, for which we compute $T_{\text{K}}=0.013D$. We
perform all calculations on a logarithmic grid, and employ $z$-averaging over
$n_{z}=8$ values of $z$. We employed a Wilson chain of length $L=51$ sites,
with $\Lambda=2$. Natural orbitals were obtained with the RGNO algorithm. A
correlated sector with $N_{\mathcal{C}}=12$ orbitals (6 spin up and 6 spin
down), containing $3$ spin up and $3$ spin down electrons was used. The RGNO
algorithm could determine the ground state energy to an error of $\sim 5\%$ of
the Kondo temperature, meaning that we may expect Kondo physics to be
described accurately. We note that the accuracy of the ground state
calculation is easily improved by increasing the size of the correlated
sector. However, the enlarged correlated sector significantly slows down the
subsequent calculation of excited states. We have therefore settled for the
minimum accuracy that is sufficient for the excited state calculation.

Figure 9: The local density of states (LDOS) at the $d$-orbital, for the
symmetric Single Impurity Anderson Model (SIAM) with $\Delta=0.09D$ and
$U=0.6D$. The main panel shows the result over a frequency window that spans
the whole band. The inset zooms in on the Kondo resonance. The vertical lines
in the inset represent the eigenvalues of $H_{\mathcal{C}\mathcal{C}}$ (43).
Results for Ansatz 1 and Ansatz 2 were obtained employing the same logarithmic
discretization and $z$-averaging as in NRG.

In Figure 9 we compare NRG results for the LDOS of the $d$-orbital, to results
obtained using Ansatz 1 and Ansatz 2. The NRG results show a clear Kondo
resonance of half-width $\sim 0.02D$ and broad shoulders associated with the
scale $U/2=0.3D$. Both Ansatz 1 and Ansatz 2 reproduce the general shape well
and are quantitatively accurate at low frequencies. The inset to Figure 1
shows a zoom of the central peak of the LDOS, with vertical lines indicating
the eigenenergies of $H_{\mathcal{C},\mathcal{C}}$ (43). As in the case of the
IRLM, Anzatz 1 and Anzatz 2 both handle the hybridization of the two few-body
correlated eigenstates of $H_{\mathcal{C},\mathcal{C}}$ that are closest in
energy to the Fermi level well (the two vertical lines at $\pm 0.006D$).
Again, deviations from NRG results are associated with higher excited few-body
correlated states, for instance those represented by the vertical lines at
$\pm 0.02D$. Ansatz 2 mitigates these deviations compared to Ansatz 1, but
errors at this frequency scale are still on the order of 5% for Ansatz 2.

Figure 10: The LDOS of the $d$-orbital (LABEL:eqdos) as calculated with the
Ansatz (dots), for the symmetric SIAM with $\Delta=0.09D$ and $U=0.2D$,
calculated on a regular energy grid with $L=4801$ orbitals. The solid line
represents the Lorentzian LDOS $z^{2}\Delta/\pi(\omega^{2}+z^{2}\Delta)$, with
$z$ estimated from the low-energy spectrum obtained with NRG. The Dashed line
represents the unrenormalized $z=1$ Lorentizain LDOS at the same $\Delta$.

At this point we conclude our study of Ansatz 2. Benchmarking shows that
Ansatz 1 is accurate in the Fermi liquid regime below $T_{\text{K}}$. It
remains to check Ansatz 1 for the SIAM discretized on the large regular energy
grid necessary for obtaining real space resolution at the scale of the Fermi
wavelength. Here we find that results are less accurate than on the Wilson
chain. To obtain a quantitatively accurate LDOS on a regular energy grid with
$4801$ levels at $\Delta=0.09$, we had to reduce the interaction strength to
$U=0.2D$. This gives a quasiparticle weight $z=0.77$, which still represents a
finite renormalization. In Figure 10 we show the low-energy LDOS computed with
the anzatz, compared to what it should be for $z=0.77$. We see good agreement
up to $\omega=0.02D$, and a breakdown at $\omega=0.03D$. Here
$T_{\text{K}}=0.064D$. From various trial calculations we performed with the
ansatz, we believe that the energy range in which the ansatz is accurate would
be extended, if we could use a larger correlated sector, but that would make
the calculation too expensive computationally, to perform in a reasonable
time.

## References

  * Shankar [1994] R. Shankar, Renormalization-group approach to interacting fermions, Rev. Mod. Phys. 66, 129 (1994).
  * Hewson [1993a] A. C. Hewson,  _The Kondo Problem to Heavy Fermions_ (Cambridge University Press, 1993).
  * Pustilnik and Glazman [2004] M. Pustilnik and L. Glazman, Kondo effect in quantum dots, Journal of Physics: Condensed Matter 16, R513 (2004).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Eidelstein  _et al._ [2020] E. Eidelstein, E. Gull, and G. Cohen, Multiorbital quantum impurity solver for general interactions and hybridizations, Phys. Rev. Lett. 124, 206405 (2020).
  * Linden  _et al._ [2020] N.-O. Linden, M. Zingl, C. Hubig, O. Parcollet, and U. Schollwöck, Imaginary-time matrix product state impurity solver in a real material calculation: Spin-orbit coupling in $\mathrm{Sr}{}_{2}\mathrm{RuO}{}_{4}$, Phys. Rev. B 101, 041101 (2020).
  * Werner  _et al._ [2023] D. Werner, J. Lotze, and E. Arrigoni, Configuration interaction based nonequilibrium steady state impurity solver, Phys. Rev. B 107, 075119 (2023).
  * Erpenbeck  _et al._ [2023] A. Erpenbeck, W.-T. Lin, T. Blommel, L. Zhang, S. Iskakov, L. Bernheimer, Y. Núñez Fernández, G. Cohen, O. Parcollet, X. Waintal, and E. Gull, Tensor train continuous time solver for quantum impurity models, Phys. Rev. B 107, 245135 (2023).
  * Kloss  _et al._ [2023] B. Kloss, J. Thoenniss, M. Sonner, A. Lerose, M. T. Fishman, E. M. Stoudenmire, O. Parcollet, A. Georges, and D. A. Abanin, Equilibrium quantum impurity problems via matrix product state encoding of the retarded action (2023), arXiv:2306.17216 [cond-mat.str-el] .
  * Gubernatis  _et al._ [1987] J. E. Gubernatis, J. E. Hirsch, and D. J. Scalapino, Spin and charge correlations around an anderson magnetic impurity, Phys. Rev. B 35, 8478 (1987).
  * Barzykin and Affleck [1996] V. Barzykin and I. Affleck, The kondo screening cloud: What can we learn from perturbation theory?, Phys. Rev. Lett. 76, 4959 (1996).
  * Borda [2007] L. Borda, Kondo screening cloud in a one-dimensional wire: Numerical renormalization group study, Physical Review B 75, 041307 (2007).
  * Lechtenberg and Anders [2014] B. Lechtenberg and F. B. Anders, Spatial and temporal propagation of kondo correlations, Phys. Rev. B 90, 045117 (2014).
  * Florens and Snyman [2015] S. Florens and I. Snyman, Universal spatial correlations in the anisotropic kondo screening cloud: Analytical insights and numerically exact results from a coherent state expansion, Phys. Rev. B 92, 195106 (2015).
  * Debertolis  _et al._ [2022] M. Debertolis, I. Snyman, and S. Florens, Simulating realistic screening clouds around quantum impurities: Role of spatial anisotropy and disorder, Phys. Rev. B 106, 125115 (2022).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Numerical renormalization group method for quantum impurity systems, Reviews of Modern Physics 80, 395 (2008).
  * Wilson [1975] K. G. Wilson, The renormalization group: Critical phenomena and the kondo problem, Rev. Mod. Phys. 47, 773 (1975).
  * Nozières [1974] P. Nozières, A “fermi-liquid”description of the kondo problem at low temperatures, Journal of Low Temperature Physics 17, 31 (1974).
  * Krishna-Murthy  _et al._ [1980] H. Krishna-Murthy, J. Wilkins, and K. Wilson, Renormalization-group approach to the anderson model of dilute magnetic alloys. i. static properties for the symmetric case, Physical Review B 21, 1003 (1980).
  * Hewson [1993b] A. C. Hewson, Renormalized perturbation expansions and fermi liquid theory, Phys. Rev. Lett. 70, 4007 (1993b).
  * Hewson  _et al._ [2004] A. C. Hewson, A. Oguri, and D. Meyer, Renormalized parameters for impurity models, The European Physical Journal B - Condensed Matter and Complex Systems 40, 177 (2004).
  * Mora  _et al._ [2015] C. Mora, C. Moca, J. von Delft, and G. Zaránd, Fermi-liquid theory for the single-impurity anderson model, Phys. Rev. B 92, 075120 (2015).
  * Filippone  _et al._ [2018] M. Filippone, C. Moca, A. Weichselbaum, J. von Delft, and C. Mora, At which magnetic field, exactly, does the kondo resonance begin to split? a fermi liquid description of the low-energy properties of the anderson model, Phys. Rev. B 98, 075404 (2018).
  * Dobrosavljević  _et al._ [1992] V. Dobrosavljević, T. R. Kirkpatrick, and B. G. Kotliar, Kondo effect in disordered systems, Physical Review Letters 69, 1113 (1992).
  * Zaránd and Udvardi [1996] G. Zaránd and L. Udvardi, Enhancement of the kondo temperature of magnetic impurities in metallic point contacts due to the fluctuations of the local density of states, Phys. Rev. B 54, 7606 (1996).
  * Aleiner  _et al._ [2002] I. Aleiner, P. Brouwer, and L. Glazman, Quantum effects in coulomb blockade, Physics Reports 358, 309 (2002).
  * Kaul  _et al._ [2006] R. K. Kaul, G. Zaránd, S. Chandrasekharan, D. Ullmo, and H. U. Baranger, Spectroscopy of the kondo problem in a box, Phys. Rev. Lett. 96, 176802 (2006).
  * Ullmo [2008] D. Ullmo, Many-body physics and quantum chaos, Reports on Progress in Physics 71, 026001 (2008).
  * Liu  _et al._ [2012] D. E. Liu, S. Burdin, H. U. Baranger, and D. Ullmo, Mesoscopic anderson box: Connecting weak to strong coupling, Phys. Rev. B 85, 155455 (2012).
  * Ullmo  _et al._ [2013] D. Ullmo, D. E. Liu, S. Burdin, and H. U. Baranger, Mesoscopic fluctuations in the fermi-liquid regime of the kondo problem, The European Physical Journal B 86, 353 (2013).
  * Miranda  _et al._ [2014] V. G. Miranda, L. G. G. V. Dias da Silva, and C. H. Lewenkopf, Disorder-mediated kondo effect in graphene, Physical Review B 90, 201101 (2014).
  * Slevin  _et al._ [2019] K. Slevin, S. Kettemann, and T. Ohtsuki, Multifractality and the distribution of the kondo temperature at the anderson transition, The European Physical Journal B 92, 281 (2019).
  * Brun  _et al._ [2014] B. Brun, F. Martins, S. Faniel, B. Hackens, G. Bachelier, A. Cavanna, C. Ulysse, A. Ouerghi, U. Gennser, D. Mailly, S. Huant, V. Bayot, M. Sanquer, and H. Sellier, Wigner and kondo physics in quantum point contacts revealed by scanning gate microscopy, Nature Communications 5, 4290 (2014).
  * Kolasiński  _et al._ [2016] K. Kolasiński, B. Szafran, B. Brun, and H. Sellier, Interference features in scanning gate conductance maps of quantum point contacts with disorder, Phys. Rev. B 94, 075301 (2016).
  * Brun  _et al._ [2016] B. Brun, F. Martins, S. Faniel, B. Hackens, A. Cavanna, C. Ulysse, A. Ouerghi, U. Gennser, D. Mailly, P. Simon, S. Huant, V. Bayot, M. Sanquer, and H. Sellier, Electron phase shift at the zero-bias anomaly of quantum point contacts, Phys. Rev. Lett. 116, 136801 (2016).
  * Affleck [2010] I. Affleck, The kondo screening cloud: What it is and how to observe it, in  _Perspectives of Mesoscopic Physics_ , edited by A. Aharony and O. Entin-Wohlman (WORLD SCIENTIFIC, 2010) pp. 1–44.
  * Peters  _et al._ [2006] R. Peters, T. Pruschke, and F. B. Anders, Numerical renormalization group approach to green’s functions for quantum impurity models, Phys. Rev. B 74, 245114 (2006).
  * Weichselbaum and von Delft [2007] A. Weichselbaum and J. von Delft, Sum-rule conserving spectral functions from the numerical renormalization group, Phys. Rev. Lett. 99, 076402 (2007).
  * Bravyi and Gosset [2017] S. Bravyi and D. Gosset, Complexity of quantum impurity problems, Communications in Mathematical Physics 356, 451 (2017).
  * Debertolis  _et al._ [2021] M. Debertolis, S. Florens, and I. Snyman, Few-body nature of kondo correlated ground states, Physical Review B 103, 235166 (2021).
  * Löwdin [1955] P.-O. Löwdin, Quantum theory of many-particle systems. i. physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction, Physical Review 97, 1474 (1955).
  * Davidson [1972] E. R. Davidson, Properties and uses of natural orbitals, Rev. Mod. Phys. 44, 451 (1972).
  * Olsen [2011] J. Olsen, The casscf method: A perspective and commentary, International Journal of Quantum Chemistry 111, 3267 (2011).
  * Li and Paldus [2005] X. Li and J. Paldus, Recursive generation of natural orbitals in a truncated orbital space, International Journal of Quantum Chemistry 105, 672 (2005).
  * Aikebaier  _et al._ [2023] F. Aikebaier, T. Ojanen, and J. L. Lado, Extracting electronic many-body correlations from local measurements with artificial neural networks, SciPost Phys. Core 6, 030 (2023).
  * Vanhala and Ojanen [2023] T. I. Vanhala and T. Ojanen, Complexity of fermionic states (2023), arXiv:2306.07584 [quant-ph] .
  * Zgid  _et al._ [2012] D. Zgid, E. Gull, and G. K.-L. Chan, Truncated configuration interaction expansions as solvers for correlated quantum impurity models and dynamical mean-field theory, Physical Review B 86, 165128 (2012).
  * He and Lu [2014] R.-Q. He and Z.-Y. Lu, Quantum renormalization groups based on natural orbitals, Physical Review B 89, 085108 (2014).
  * Anderson [1961] P. W. Anderson, Localized magnetic states in metals, Phys. Rev. 124, 41 (1961).
  * Lin and Demkov [2013] C. Lin and A. A. Demkov, Efficient variational approach to the impurity problem and its application to the dynamical mean-field theory, Physical Review B 88, 035123 (2013).
  * Guinea  _et al._ [1985] F. Guinea, V. Hakim, and A. Muramatsu, Bosonization of a two-level system with dissipation, Phys. Rev. B 32, 4410 (1985).
  * Kotliar and Si [1996] G. Kotliar and Q. Si, Toulouse points and non-fermi-liquid states in the mixed-valence regime of the generalized anderson model, Phys. Rev. B 53, 12373 (1996).
  * Costi and Zaránd [1999] T. A. Costi and G. Zaránd, Thermodynamics of the dissipative two-state system: A bethe-ansatz study, Phys. Rev. B 59, 12398 (1999).
  * Oliveira and Oliveira [1994] W. C. Oliveira and L. N. Oliveira, Generalized numerical renormalization-group method to calculate the thermodynamical properties of impurities in metals, Phys. Rev. B 49, 11986 (1994).
  * Camacho  _et al._ [2022] G. Camacho, P. Schmitteckert, and S. T. Carr, Local density of states of the interacting resonant level model at zero temperature, Phys. Rev. B 105, 075116 (2022).
  * Hofstetter [2000] W. Hofstetter, Generalized numerical renormalization group for dynamical quantities, Phys. Rev. Lett. 85, 1508 (2000).