Renormalization approach to the superconducting Kondo model

Steffen Sykora1$\star$ and Tobias Meng1

1 Institute for Theoretical Physics and Würzburg-Dresden Cluster of Excellence
ct.qmat, Technische Universität Dresden, 01069 Dresden, Germany

⋆<EMAIL_ADDRESS>
## Abstract

An approach to bound states based on unitary transformations of Hamiltonians
is presented. The method is applied to study the interaction between electrons
in a BCS $s$-wave superconductor and a quantum spin. It is shown that known
results from the t-matrix method and numerical studies are reproduced by this
new method. As a main advantage, the method can straightforwardly be extended
to study the topological properties of combined bound states in chains of many
magnetic impurities. It also provides a uniform picture of the interplay
between the Yu-Shiba-Rusinov (YSR) bound states and the Kondo singlet state.

###### Contents

  1. 1 Introduction
  2. 2 Theoretical approach
    1. 2.1 Renormalized Hamiltonian
    2. 2.2 Bound states
    3. 2.3 Expectation values
  3. 3 Numerical results
    1. 3.1 YSR state energy and Kondo resonance
    2. 3.2 Spectral function
    3. 3.3 Influence of temperature and spin size
  4. 4 Conclusion
  5. A Classical spin
  6. B Non-magnetic impurity

##  1 Introduction

The study of individual impurity bound states in electronic systems has become
an important experimental and theoretical tool for the characterization of
correlated quantum materials. A prominent example is a system in which a
magnetic impurity interacts with a superconducting condensate, leading to the
formation of Yu-Shiba-Rusinov (YSR) bound states [1, 2, 3], but also to Kondo
physics [4, 5, 6]. Depending on the strength of the exchange coupling, a
transition from a regime dominated by weakly coupled YSR states to a more
Kondo-like state is found[7, 8, 9]. This physics has at least in parts be
visualized by scanning tunneling spectroscopy experiments [11, 12, 10].

In addition, the controlled extension to more than one local moment has
recently become an important theme. This gives rise to models of molecules
[13, 14, 15, 16] and spin chains [17, 18]. The interest in these systems has
prominently been fueled by predictions of topological superconductivity and
Majorana zero modes [19, 20]. In the discussion of this intriguing physics,
however, the influence of the quantum nature of the impurity spins, and the
resulting Kondo physics has not beet explored extensively. This is in parts
due to the fact that the corresponding theoretical treatment requires true
many-body approaches that properly take into account the quantum nature of an
entire collection of local moments.

In recent years, some powerful theoretical methods like mean field
calculations[7, 21], perturbation theory[22, 23], self-consistent
approaches[24, 26, 25] and the numerical renormalization group (NRG) approach
[27, 28, 29] have been applied to study the formation of bound states in
correlated systems. Many of these methods are restricted to only a very small
number of quantum impurities, or to very small system sizes. In particular,
powerful analytical approaches which are able to tackle a collection of
quantum impurities in an environment of a condensate in the thermodynamic
limit are not available so far. Such an approach would allow for a deeper
understanding of the underlying physical processes and the structure of the
quasiparticle which is responsible for the bound state. In addition, the
quantum nature of the local impurity is often important for the properties of
the spectral function.

Figure 1:  Schematic picture of our theoretical approach to bound states. The
scattering term $\sim J$ in the original Hamiltonian is eliminated by a
unitary transformation $e^{X}$. This technique diagonalizes the Hamiltonian
and generates a bound state with quantum number ${\bf k}_{0}$ (red dot), which
is renormalized in energy. Excitations split into a coherent part and a local
excitation.

Motivated by these open questions, we have developed a new renormalization
scheme for Hamiltonians describing a local moment interacting with a
superconducting condensate. Our diagonalization method shares some basic
concepts with the known flow equation approaches[30, 31, 32]. The method is
technically organized in such a way that the formation of bound states is
highlighted. A schematic picture of the method is shown in Fig. 1. As we
discuss below, bound states lead to singularities in the unitary
transformation on which our method builds. These singularities in turn lead to
additional local contributions to the quasiparticle operators (the bound
state), and to a significant energy renormalization (the bound state energy).
The quantum nature of the impurity spin is taken into account on an
approximate level, but we show that this approximation is sufficiently refined
to allow access to Kondo physics. Due to the readily generalizable unitary
transformations used to diagonalize the Hamiltonian, also more than one
impurity can in principle be implemented (for normal state electronic Kondo
lattices, this has already been worked out in Ref. [33]).

The rest of the paper is organized as follows. In Sec. 2, we introduce the
model Hamiltonian. Sec. 2.1 is devoted to a detailled discussion of our
renormalization method. In particular, we explain the method of integrating
out the electron-impurity spin interaction, and show how the corresponding
renormalization equations are derived. Then, in Sec. 2.2, we explain the
appearance of bound states in our theory, and how the bound states are related
to singularities in the unitary transformation. In Sec. 2.3, we discuss the
access to expectation values within our approach. We apply the developed
concepts to renormalize single-particle operators and present two examples of
expectation values, the occupation number and the single-particle spectral
function. Sec. 3 is devoted to numerical results for the YSR bound state
energy and the single-particle spectral function in the weak and strong
coupling regimes, which we relate to the nature of the excitations to
highlight the coexistence and competition of YSR and Kondo physics. We
furthermore show that our methods compared very well with other methods, in
particular by benchmarking with numerical renormalization group (NRG) studies.
Finally, a conclusion is given in Sec. 4.

##  2 Theoretical approach

Our starting point is a basic BCS-type $s$-wave superconductor coupled to a
single magnetic impurity. Denoting the spin by an index
$\sigma=\uparrow,\downarrow$, this system is modelled by the Hamiltonian

$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{1}\quad\mbox{with}\quad\mathcal{H}_{0}=\sum_{\mathbf{k}\sigma}\varepsilon_{\mathbf{k}}c^{\dagger}_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}+\Delta\sum_{\mathbf{k}}(c^{\dagger}_{\mathbf{k}\uparrow}c^{\dagger}_{-\mathbf{k}\downarrow}+c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow})\quad\mbox{and}\quad\mathcal{H}_{1}=J{\mathbf{S}}\cdot{\mathbf{s}}_{\mathbf{r}_{\text{imp}}}\,,$
(1)

where ${\mathbf{S}}$ and ${\mathbf{s}}_{\mathbf{r}_{\text{imp}}}$ are the
local impurity spin and the electron spin at the impurity site
$\mathbf{r}_{\text{imp}}$, respectively. Here, $\varepsilon_{\mathbf{k}}$ is
the bare electronic dispersion, which is a function of the momentum vector
${\bf k}$, and $\Delta$ is the superconducting gap. The lattice consists of
$N$ sites. We consider the case of a single local magnetic impurity
represented by a quantum spin of arbitrary size that is described by an
angular momentum operator $\mathbf{S}$. This impurity is coupled to the
electrons via a Heisenberg exchange coupling $J$, and the coordinate origin is
chosen to coincide with the impurity position.

Before dealing with the magnetic impurity, our approach requires to
diagonalize the purely electronic Hamiltonian $\mathcal{H}_{0}$. This is
achieved by introducing the usual Bogoliubov quasiparticles,

$\alpha_{\bf k}^{\dagger}=u_{\bf k}c_{\mathbf{k}\uparrow}^{\dagger}-v_{\bf
k}c_{-{\bf k}\downarrow},\quad\beta_{\bf k}^{\dagger}=u_{\bf
k}c_{-\mathbf{k}\downarrow}^{\dagger}+v_{\bf k}c_{{\bf k}\uparrow},$ (2)

with

$u_{\bf
k}^{2}=\frac{1}{2}\left(1+\frac{\varepsilon_{\mathbf{k}}}{\sqrt{\varepsilon_{\mathbf{k}}^{2}+\Delta^{2}}}\right),\quad
v_{\bf
k}^{2}=\frac{1}{2}\left(1-\frac{\varepsilon_{\mathbf{k}}}{\sqrt{\varepsilon_{\mathbf{k}}^{2}+\Delta^{2}}}\right).$
(3)

The electronic Hamiltonian $\mathcal{H}_{0}$ is thus brought to the diagonal
form

$\mathcal{H}_{0}=\sum_{\mathbf{k}}E_{\mathbf{k}}\left(\alpha_{\bf
k}^{\dagger}\alpha_{\bf k}+\beta_{\bf k}^{\dagger}\beta_{\bf
k}\right)+\sum_{\mathbf{k}}\left(\varepsilon_{\mathbf{k}}-E_{\bf k}\right),$
(4)

with the quasiparticle energy
$E_{\mathbf{k}}=\sqrt{\varepsilon_{\mathbf{k}}^{2}+\Delta^{2}}$. Next, we turn
to the electron-impurity coupling. In the original basis, it reads

$\mathcal{H}_{1}=\frac{J}{N}\sum_{{\bf k}{\bf
k}^{\prime}}\sum_{\alpha,\beta}{\mathbf{S}}\cdot\frac{\vec{\sigma}_{\alpha\beta}}{2}c_{{\bf
k}\alpha}^{\dagger}c_{{\bf k}^{\prime}\beta}.$ (5)

Replacing the conduction electron operators with the ones of Bogoliubov
quasiparticles using Eq. (2), and introducing the shorthands

$C_{1,{\bf k}{\bf k}^{\prime}}^{\pm}=u_{\bf k}u_{{\bf k}^{\prime}}\pm v_{\bf
k}v_{{\bf k}^{\prime}},\quad C_{2,{\bf k}{\bf k}^{\prime}}^{\pm}=u_{\bf
k}v_{{\bf k}^{\prime}}\pm v_{\bf k}u_{{\bf k}^{\prime}}$ (6)

for combinations of Bogoliubov coefficients (coherence factors), the electron-
impurity-coupling takes the form

$\begin{split}&\mathcal{H}_{1}=\frac{J}{2N}\sum_{{\bf k}{\bf
k}^{\prime}}\Big{\\{}C_{1,{\bf k}{\bf k}^{\prime}}^{+}\Big{(}\alpha_{\bf
k}^{\dagger}\alpha_{{\bf k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\Big{)}\\\ +&C_{2,{\bf k}{\bf k}^{\prime}}^{-}\left[\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}+\beta_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S_{z}+\frac{1}{2}\left(\beta_{{\bf
k}^{\prime}}\beta_{{\bf k}}-\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\right]\Big{\\}}.\end{split}$ (7)

Here, $S^{+}$ ($S^{-}$) is the spin raising (lowering) operator, while $S_{z}$
denotes the $z$-component of the spin operator $\mathbf{S}$. At this point,
all we demand of these operators is that they satisfy the commutation
relations $[S^{+},S^{-}]=2S_{z}$, $[S_{z},S^{-}]=-S^{-}$, and
$[S_{z},S^{+}]=S^{+}$. In the following, we show that the Hamiltonian
$\mathcal{H}_{1}$ can be diagonalized by a unitary transformation in
combination with a factorization scheme. As we discuss, this approach is able
to describe the presence of the Yu-Shiba-Rusinov (YSR) bound states [1, 2, 3],
but also of Kondo physics.

###  2.1 Renormalized Hamiltonian

Unlike the case of a classical spin in a superconductor, or the even simpler
case of a potential impurity in a spinless electron bath, the case of a
quantum spin immersed into an $s$-wave superconductor constitutes an
interacting quantum problem, and therefore cannot be solved exactly. Our
approach is built to tackle this interaction on a refined approximate level.
The power of our approach can be gauged from the fact that a variant of our
approximation scheme has already used to successfully describe the Kondo
effect of a quantum spin coupled to a normal-state electron system [33]. It
can therefore be expected to also yield excellent results when applied to a
superconducting system. After discussing technical details of the
transformation we use to tackle the problem, we will first compare our results
with recent experimental and numerical studies, and then show original results
for ${\mathbf{k}}$-resolved spectral function. We furthermore discuss that our
approach clearly illustrates how Kondo physics appears. The latter is
associated with an energy scale $k_{\text{B}}T_{\text{K}}$ (with
$k_{\text{B}}$ being the Boltzmann constant and $T_{\text{K}}$ the Kondo
temperature), and competes with the singlet state of the superconductor [34].

The diagonalization method to be used here is reviewed in Ref. [35]. Following
the spirit of this approach, we subject the Hamiltonian (1) to a mapping that
brings the Hamiltonian to a diagonal form $\tilde{\mathcal{H}}$, which in turn
can be used to calculate physical observables. According to Ref. [35], in a
general situation of an interacting many-body problem the diagonalization
procedure of the Hamiltonian needs to be carried out in a stepwise way. In the
present case, we show that this multiple-step procedure can be constructed
such that the coupling term in the Hamiltonian is continously reduced until it
is finally fully eliminated. At each step, the unitary transformation can be
limited to the lowest order with respect to a small variable which accounts
for the difference between the control parameter values of two subsequent
transformation steps. This allows the controlled derivation of renormalization
equations for the parameters of the transformed Hamiltonian, which in practice
are solved numerically. Eventually, we end up with an effectively free system
of Bogoliubov quasiparticles $\tilde{\mathcal{H}}$, whose renormalized
energies account for the effects of the magnetic impurity.

The mapping between the original Hamiltonian ${\mathcal{H}}$ and an arbitrary
intermediate Hamiltonian ${\mathcal{H}}_{\lambda}$ appearing during the
transformation process is implemented by a unitary transformation,

${\mathcal{H}}_{\lambda}=e^{X_{\lambda}}{\mathcal{H}}e^{-X_{\lambda}},$ (8)

where the hermiticity of ${\mathcal{H}}_{\lambda}$ requires
$X_{\lambda}^{\dagger}=-X_{\lambda}$. The dimensionless control parameter
$\lambda$ describes the progress of the diagonalization. It is defined such
that its initial value $\lambda=1$ marks the starting point where
${\mathcal{H}}_{\lambda=1}={\mathcal{H}}$ and its final value $\lambda=0$
describes the fully diagonalized effective Hamiltonian
${\mathcal{H}}_{\lambda=0}=\tilde{\mathcal{H}}$ with

$\tilde{\cal H}=\sum_{\bf k}\tilde{E}_{{\bf k}}\left(\alpha_{\bf
k}^{\dagger}\alpha_{\bf k}+\beta_{\bf k}^{\dagger}\beta_{\bf
k}\right)+\tilde{E},$ (9)

where finding the renormalized dispersion $\tilde{E}_{\bf k}$ and the energy
constant $\tilde{E}$ is the aim of the diagonalization process. Note that
since $\tilde{\cal H}$ is related to ${\cal H}$ by a unitary transformation
(and a factorization scheme, see below), the energy values $\tilde{E}_{\bf k}$
correspond (to a very good approximation) to the eigenvalues of the original
Hamiltonian, including the binding energies of possible bound states. Using
the introduced definition of $\lambda$, we define the transformed Hamiltonian
for an arbitrary $\lambda$ as

$\begin{split}{\mathcal{H}}_{\lambda}=&\sum_{\bf k}E_{{\bf
k},\lambda}\left(\alpha_{\bf k}^{\dagger}\alpha_{\bf k}+\beta_{\bf
k}^{\dagger}\beta_{\bf k}\right)+E_{\lambda}\\\
+&\lambda\frac{J}{2N}\sum_{{\bf k}\neq{\bf k}^{\prime}}\Big{\\{}C_{1,{\bf
k}{\bf k}^{\prime}}^{+}\Big{(}\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\Big{)}\\\ +&C_{2,{\bf k}{\bf k}^{\prime}}^{-}\left[\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}+\beta_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S_{z}+\frac{1}{2}\left(\beta_{{\bf
k}^{\prime}}\beta_{{\bf k}}-\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\right]\Big{\\}}.\end{split}$ (10)

Up to the parameter $\lambda$, the second and third line corresponds to the
original interaction part ${\mathcal{H}}_{1}$.

The renormalization equations for $E_{{\bf k},\lambda}$ and $E_{\lambda}$ are
found from a small transformation step which maps ${\mathcal{H}}_{\lambda}$ to
the corresponding effective Hamiltonian
${\mathcal{H}}_{\lambda-\Delta\lambda}$ referring to a somewhat smaller
parameter $(\lambda-\Delta\lambda)$ with $\Delta\lambda\ll 1$. Thereby we
demand that both ${\mathcal{H}}_{\lambda}$ and
${\mathcal{H}}_{\lambda-\Delta\lambda}$ fulfill Eq. (10) and are related to
each other via a unitary transformation,

${\mathcal{H}}_{\lambda-\Delta\lambda}=e^{X_{\lambda,\Delta\lambda}}{\mathcal{H}}_{\lambda}e^{-X_{\lambda,\Delta\lambda}},$
(11)

where the generator $X_{\lambda,\Delta\lambda}$ has to be constructed such
that ${\mathcal{H}}_{\lambda-\Delta\lambda}$ keeps the structure of Eq. (10)
but has a slightly reduced coupling term. It in that sense is closer to a
diagonalized Hamiltonian than ${\mathcal{H}}_{\lambda}$. For the generator
$X_{\lambda,\Delta\lambda}$ satisfying these requirements, we make the
following ansatz

$\begin{split}&X_{\lambda,\Delta\lambda}=\frac{\Delta\lambda}{N}\sum_{{\bf
k}\neq{\bf k}^{\prime}}\left\\{\left(J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\right)\right.\right.\\\ +&\left.\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}\left(\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}-\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)S_{z}-\frac{1}{2}\left(\beta_{{\bf k}^{\prime}}\beta_{{\bf
k}}+\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}+\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\right)\right]\\\ &-\left.\left(V_{{\bf
k},\lambda}+V_{{\bf k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}+\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}\right)-\frac{C_{2,{\bf k}{\bf k}^{\prime}}^{-}}{E_{{\bf
k},\lambda}+E_{{\bf k}^{\prime},\lambda}}\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\beta_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)\right]\right\\}.\end{split}$ (12)

Importantly, this generator is proportional to the small parameter
$\Delta\lambda$. It consists of two types of terms. The first two lines
contain terms that occur in the original interaction part ${\mathcal{H}}_{1}$,
and that are composed of products of fermionic and spin operators. The third
line, on the contrary, contains only terms independent of the spin and
quadratic in the fermionic operators. As we will show below, terms of this
form are generated during the diagonalization process, and must therefore also
be included in the generator. The unknown coefficients $J_{{\bf k},\lambda}$
and $V_{{\bf k},\lambda}$, finally, will be determined in such a way that
after evaluating the right-hand side of Eq. (11) using the ansatz (12), the
form of Eq. (10) is exactly reproduced, except that the parameter $\lambda$ is
replaced by $(\lambda-\Delta\lambda)$. The step width $\Delta\lambda\ll 1$ can
be chosen arbitrarily small. We can thus approximate the unitary
transformation in Eq. (11) by

${\mathcal{H}}_{\lambda-\Delta\lambda}\approx{\mathcal{H}}_{\lambda}+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{\lambda}].$
(13)

In order to find the equations for $J_{{\bf k},\lambda}$ and $V_{{\bf
k},\lambda}$, as well as for $E_{{\bf k},\lambda}$ and $E_{\lambda}$, we have
to compute the commutator between $X_{\lambda,\Delta\lambda}$ and
${\mathcal{H}}_{\lambda}$. To facilitate the discussion of the various
resulting terms, let us first formally decompose the renormalized Hamiltonian
into its diagonal and non-diagonal parts,
${\mathcal{H}}_{\lambda}={\mathcal{H}}_{0,\lambda}+{\mathcal{H}}_{1,\lambda}$,
with

${\mathcal{H}}_{0,\lambda}=\sum_{\bf k}E_{{\bf k},\lambda}\left(\alpha_{\bf
k}^{\dagger}\alpha_{\bf k}+\beta_{\bf k}^{\dagger}\beta_{\bf
k}\right)+E_{\lambda}$ (14)

and

$\begin{split}{\mathcal{H}}_{1,\lambda}&=\lambda\frac{J}{2N}\sum_{{\bf
k}\neq{\bf k}^{\prime}}\left\\{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\Big{(}\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\Big{)}\right.\\\ +&\left.C_{2,{\bf k}{\bf
k}^{\prime}}^{-}\left[\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}+\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)S_{z}+\frac{1}{2}\left(\beta_{{\bf k}^{\prime}}\beta_{{\bf
k}}-\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\right]\right\\}.\end{split}$ (15)

The commutator between $X_{\lambda,\Delta\lambda}$ and
${\mathcal{H}}_{\lambda}$ can now be written in the form

$[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{\lambda}]=[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{0,\lambda}]+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}].$
(16)

Since ${\mathcal{H}}_{0,\lambda}$ is diagonal in the fermion operators, the
first part $[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{0,\lambda}]$ is found
directly from Eq. (12) by cancellation of the fermion energies $E_{{\bf
k},\lambda}$ in the denominators. We obtain

$\begin{split}&[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{0,\lambda}]=-\frac{\Delta\lambda}{N}\sum_{{\bf
k}\neq{\bf k}^{\prime}}\Big{\\{}\left(J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}\right)\left[C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\right)\right.\\\ +&C_{2,{\bf k}{\bf
k}^{\prime}}^{-}\left(\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}+\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)S_{z}-\frac{1}{2}\left(-\beta_{{\bf k}^{\prime}}\beta_{{\bf
k}}+\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\right)\Big{]}\\\ &-\left.\left(V_{{\bf
k},\lambda}+V_{{\bf k}^{\prime},\lambda}\right)\left[C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}+\beta_{\bf k}^{\dagger}\beta_{{\bf k}^{\prime}}\right)-C_{2,{\bf
k}{\bf k}^{\prime}}^{-}\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}+\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)\right]\right\\}.\end{split}$ (17)

The second part of the commutator (16) gives rise to an additional internal
summation, and generates new contributions to an effective electron-impurity
scattering. We distinguish three types of terms,

$[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]=[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{1}+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{2}+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{3}.$
(18)

The part $[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{1}$ arises
from commutators involving Bogoliobov quasiparticle operators from the first
and second lines in Eq. (12). Initially, these terms contain products of spin
operators. Due to spin rotation invariance, however, the spin operators
combine to form $(\vec{S}\cdot\vec{S})$, and can be replaced by the number
$S(S+1)$. The term $[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{1}$
is therefore merely a bilinear of fermionic operators. We obtain

$\begin{split}&[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{1}=\\\
&-\frac{J\lambda\Delta\lambda}{2N}S(S+1)\sum_{{\bf k}{\bf
k}^{\prime}}\left(A_{{\bf k},\lambda}^{(1)}+A_{{\bf
k}^{\prime},\lambda}^{(1)}\right)\left[C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}+\beta_{\bf k}^{\dagger}\beta_{{\bf k}^{\prime}}\right)-C_{2,{\bf
k}{\bf k}^{\prime}}^{-}\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}+\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)\right],\end{split}$ (19)

where we introduced the short-hand notation

$A_{{\bf k},\lambda}^{(1)}=\frac{1}{N}\sum_{{\bf q}(\neq{\bf k})}\frac{J_{{\bf
k},\lambda}+J_{{\bf q},\lambda}}{E_{{\bf q},\lambda}^{2}-E_{{\bf
k},\lambda}^{2}}\left[E_{{\bf k},\lambda}+\left(C_{2,{\bf q}{\bf
q}}^{+}-C_{1,{\bf q}{\bf q}}^{-}\right)E_{{\bf q},\lambda}\right].$ (20)

Note that according to Eqs. (3) and (6), the coherence factors with respect to
equal quantum number $\mathbf{q}$ have the solutions $C_{2,{\bf q}{\bf
q}}^{+}=\Delta/E_{\bf q}$ and $C_{1,{\bf q}{\bf q}}^{-}=\varepsilon_{\bf
q}/E_{\bf q}$. It is seen that the operator structure of
$[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{1}$ is different from
${\mathcal{H}}_{1,\lambda}$ since this commutator part only contains
scattering of fermions without a coupling to the spin operator. The natural
appearance of such terms is the reason why we have to take into account such
decoupled fermion scatterings in the ansatz of the generator in Eq. (12)
(third line).

The bilinear terms in the generator give rise to a second type of terms which
we denote as $[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{2}$. It
describes the commutator between the third line of Eq. (12) and all of
${\mathcal{H}}_{1,\lambda}$. We find

$\begin{split}[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{2}&=\frac{J\lambda\Delta\lambda}{2N}\sum_{{\bf
k}{\bf k}^{\prime}}\left(A_{{\bf k},\lambda}^{(2)}+{\cal A}_{{{\bf
k}^{\prime},\lambda}}^{(2)}\right)\Big{\\{}C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\Big{[}\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\Big{]}\\\ &+C_{2,{\bf k}{\bf k}^{\prime}}^{-}\Big{[}\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}+\beta_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S_{z}+\frac{1}{2}\left(\beta_{{\bf
k}^{\prime}}\beta_{{\bf k}}-\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\Big{]}\Big{\\}},\end{split}$ (21)

where we have introduced

$A_{{\bf k},\lambda}^{(2)}=\frac{1}{N}\sum_{{\bf q}(\neq{\bf k})}\frac{V_{{\bf
k},\lambda}+V_{{\bf q},\lambda}}{E_{{\bf q},\lambda}^{2}-E_{{\bf
k},\lambda}^{2}}\left[E_{{\bf k},\lambda}+\left(C_{2,{\bf q}{\bf
q}}^{+}-C_{1,{\bf q}{\bf q}}^{-}\right)E_{{\bf q},\lambda}\right].$ (22)

It is seen that the structure of the original scattering interaction is
maintained.

The remaining contributions, summarised in the term
$[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{3}$, contain all terms
in which commutators between the spin operators are taken. These new terms
would be absent in the case of a classical spin. They describe additional
many-body interactions, and are a crucial new ingredient in our approach.
Concretely, we find the following operator structure

$\begin{split}&[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{3}=\frac{J\lambda\Delta\lambda}{2N^{2}}\sum_{{\bf
k}{\bf k}^{\prime}}\sum_{{\bf k}_{1}{\bf
k}^{\prime}_{1}}\left\\{\left(\frac{J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}+\frac{J_{{\bf k}_{1},\lambda}+J_{{\bf
k}^{\prime}_{1},\lambda}}{E_{{\bf k}_{1},\lambda}-E_{{\bf
k}^{\prime}_{1},\lambda}}\right)C_{1,{\bf k}{\bf k}^{\prime}}^{+}C_{1,{\bf
k}_{1}{\bf k}^{\prime}_{1}}^{+}\right.\\\ &\times\left[\left(\alpha_{\bf
k}^{\dagger}\alpha_{{\bf k}^{\prime}}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}\right)\left(\beta_{{\bf k}^{\prime}_{1}}^{\dagger}\alpha_{{\bf
k}_{1}}S^{+}+\alpha_{{\bf k}_{1}}^{\dagger}\beta_{{\bf
k}^{\prime}_{1}}S^{-}\right)-2\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}\beta_{{\bf k}^{\prime}_{1}}^{\dagger}\alpha_{{\bf
k}_{1}}S_{z}\right]\\\ +&\left(\frac{J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}+\frac{J_{{\bf k}_{1},\lambda}+J_{{\bf
k}^{\prime}_{1},\lambda}}{E_{{\bf k}_{1},\lambda}+E_{{\bf
k}^{\prime}_{1},\lambda}}\right)C_{2,{\bf k}{\bf k}^{\prime}}^{-}C_{2,{\bf
k}_{2}{\bf k}^{\prime}_{1}}^{-}\\\ &\times\left[\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\beta_{{\bf
k}^{\prime}}\alpha_{{\bf k}}\right)\left(\beta_{{\bf
k}^{\prime}_{1}}\beta_{{\bf k}_{1}}S^{+}+\alpha_{{\bf
k}_{1}}^{\dagger}\alpha_{{\bf
k}^{\prime}_{1}}^{\dagger}S^{-}\right)-\frac{1}{2}\left(\beta_{\bf
k}\beta_{{\bf k}^{\prime}}+\alpha_{{\bf k}^{\prime}}^{\dagger}\alpha_{{\bf
k}}^{\dagger}\right)\left(\beta_{{\bf k}^{\prime}_{1}}^{\dagger}\beta_{{\bf
k}_{1}}^{\dagger}-\alpha_{{\bf k}_{1}}\alpha_{{\bf
k}^{\prime}_{1}}\right)S_{z}\right]\\\ +&\left[\left(\frac{J_{{\bf
k},\lambda}+J_{{\bf k}^{\prime},\lambda}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}+\frac{J_{{\bf k}_{1},\lambda}+J_{{\bf
k}^{\prime}_{1},\lambda}}{E_{{\bf k}_{1},\lambda}+E_{{\bf
k}^{\prime}_{1},\lambda}}\right)C_{1,{\bf k}{\bf k}^{\prime}}^{+}C_{1,{\bf
k}_{1}{\bf k}^{\prime}_{1}}^{-}+\left(\frac{J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}+\frac{J_{{\bf k}_{1},\lambda}+J_{{\bf
k}^{\prime}_{1},\lambda}}{E_{{\bf k}_{1},\lambda}-E_{{\bf
k}^{\prime}_{1},\lambda}}\right)C_{1,{\bf k}{\bf k}^{\prime}}^{-}C_{1,{\bf
k}_{1}{\bf k}^{\prime}_{1}}^{+}\right]\\\ &\times\left[\left(\alpha_{\bf
k}^{\dagger}\alpha_{{\bf k}^{\prime}}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}\right)\left(\beta_{{\bf k}^{\prime}_{1}}\beta_{{\bf
k}_{1}}S^{+}+\alpha_{{\bf k}_{1}}^{\dagger}\alpha_{{\bf
k}^{\prime}_{1}}^{\dagger}S^{-}\right)-\frac{1}{2}\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}+\beta_{{\bf
k}^{\prime}}^{\dagger}\alpha_{{\bf k}}\right)\left(\beta_{{\bf
k}^{\prime}_{1}}^{\dagger}\beta_{{\bf k}_{1}}^{\dagger}-\alpha_{{\bf
k}_{1}}\alpha_{{\bf k}^{\prime}_{1}}\right)S_{z}\right].\end{split}$ (23)

It is at this step that our method requires an approximation to be made.
Namely, we apply a factorization approximation [35] that replaces a product of
four fermionic operators by a bilinear times an expectation value. Concretely,
the factorization is carried out as follows. Each time a transformation yields
a term with four Bogoliubov quasiparticle operators of the form
$\alpha^{\dagger}\alpha\beta^{\dagger}\beta$, we make the approximation

$\alpha^{\dagger}\alpha\beta^{\dagger}\beta\approx\alpha^{\dagger}\alpha\langle\beta^{\dagger}\beta\rangle_{\lambda}+\langle\alpha^{\dagger}\alpha\rangle_{\lambda}\beta^{\dagger}\beta-\langle\alpha^{\dagger}\alpha\rangle_{\lambda}\langle\beta^{\dagger}\beta\rangle_{\lambda},$
(24)

where the $\lambda$-dependent expectation value is given by

$\langle{\cal A}\rangle_{\lambda}=\frac{Tr\left({\cal A}e^{-\beta{\cal
H_{\lambda}}}\right)}{Tr\left(e^{-\beta{\cal H_{\lambda}}}\right)},$ (25)

for an arbitrary operator ${\cal A}$. The factor $\beta=1/(k_{B}T)$ contains
the temperature $T$ and the Boltzmann constant $k_{B}$. A similar replacement
is done for terms with four Bogoliubov operators of the same type
($\alpha^{\dagger}$ or $\beta^{\dagger}$) by taking into account normal
ordering. Since the factorization is carried out within the transformation
step from $\lambda$ to $\lambda-\Delta\lambda$, the expectation value is
formed with respect to $\mathcal{H}_{\lambda}$. This factorization restores an
operator structure that is quadratic in fermionic operators, and therefore
allows us to derive renormalization equations for $J_{{\bf k},\lambda}$ and
$A_{{\bf k},\lambda}$. The price to pay is that these renormalization
equations contain expectation values that have to be calculated self-
consistently.

The above factorization approximation looks somewhat similar to a mean field
approximation. However, it differs from such a treatment since the
factorization is applied at each order of the unitary transformation. The
resulting series is summed up to infinite order since all orders have a
similar structure. Our approach is thus different from a simple mean field
approximation in which a single factorization is applied. As discussed in Ref.
[33], our factorization approach has already been employed successfully for a
treatment of the Kondo effect in a normal-state metal, which illustrates its
non-perturbative power. Moreover, this approximation scheme is the standard
technique to treat many-body interactions in common flow equation approaches
[31, 32]. Applying the factorization approximation to
$[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{3}$, we obtain

$\begin{split}&[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{3,\text{factorized}}=\\\
&\frac{J\lambda\Delta\lambda}{2N}\sum_{{\bf k}{\bf k}^{\prime}}\left(A_{{\bf
k},\lambda}^{(3)}+A_{{{\bf
k}^{\prime}},\lambda}^{(3)}\right)\Big{\\{}C_{1,{\bf k}{\bf
k}^{\prime}}^{+}\left[\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}S_{z}-\beta_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S_{z}+\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}S^{-}+\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\right]\\\ &+C_{2,{\bf k}{\bf k}^{\prime}}^{-}\Big{[}\left(\alpha_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}+\beta_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S_{z}+\frac{1}{2}\left(\beta_{{\bf
k}^{\prime}}\beta_{{\bf k}}-\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}\right)S^{-}+\frac{1}{2}\left(\beta_{{\bf
k}}^{\dagger}\beta_{{\bf k}^{\prime}}^{\dagger}-\alpha_{{\bf
k}^{\prime}}\alpha_{\bf k}\right)S^{+}\Big{]}\Big{\\}},\end{split}$ (26)

where we have introduced

$\begin{split}A_{{\bf k},\lambda}^{(3)}=&\frac{1}{N}\sum_{{\bf q}(\neq{\bf
k}),{\bf k}^{\prime}}\frac{J_{{\bf k},\lambda}+J_{{\bf q},\lambda}}{E_{{\bf
q},\lambda}^{2}-E_{{\bf k},\lambda}^{2}}\left[\left(C_{1,{\bf q}{\bf
k}^{\prime}}^{+}+C_{2,{\bf q}{\bf k}^{\prime}}^{-}\right)E_{{\bf
q},\lambda}-\left(C_{2,{\bf q}{\bf k}^{\prime}}^{+}+C_{1,{\bf q}{\bf
k}^{\prime}}^{-}\right)E_{{\bf k},\lambda}\right]\\\ &\times\langle\alpha_{\bf
q}\alpha_{{\bf k}^{\prime}}^{\dagger}-\beta_{\bf q}^{\dagger}\beta_{{\bf
k}^{\prime}}+\beta_{\bf q}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}+\alpha_{{\bf k}^{\prime}}\beta_{\bf
q}\rangle_{\lambda}.\end{split}$ (27)

The effective scattering in Eq. (26), which was obtained after the
factorization, has the same operator structure as the scattering part
proportional to $\lambda$ in Eq. (10). This property allows us to find
equations for the coefficients $J_{{\bf k},\lambda}$ and $A_{{\bf k},\lambda}$
by adding together all the four commutator parts from Eqs. (17),
(LABEL:Commutator_2), (21), and (26). From Eq. (13), we find

${\mathcal{H}}_{\lambda-\Delta\lambda}={\mathcal{H}}_{\lambda}+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{0,\lambda}]+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{1}+[X_{\lambda,\Delta\lambda},{\mathcal{H}}_{1,\lambda}]_{2}+[X_{\lambda,\Delta\lambda},\mathcal{H}_{1,\lambda}]_{3,\text{factorized}}.$
(28)

Using Eq. (10) and the calculated commutators, we compare both sides of the
Eq. (28) and obtain the following set of renormalization equations,

$\displaystyle-\frac{J}{2}$ $\displaystyle=-J_{{\bf
k},\lambda}+\frac{J\lambda}{2}\left(A_{{\bf k},\lambda}^{(2)}+A_{{\bf
k},\lambda}^{(3)}\right),$ (29) $\displaystyle 0$ $\displaystyle=V_{{\bf
k},\lambda}+\frac{J\lambda}{2}S(S+1)A_{{\bf k},\lambda}^{(1)}.$ (30)

This set of linear equations determines the coefficients $J_{{\bf k},\lambda}$
and $V_{{\bf k},\lambda}$ as a function of $\lambda$. In addition, we find the
renormalization equation for $E_{{\bf k},\lambda}$,

$\begin{split}E_{{\bf k},\lambda-\Delta\lambda}=E_{{\bf
k},\lambda}-\frac{J\lambda\Delta\lambda}{N^{2}}S(S+1)\sum_{\bf q}\frac{J_{{\bf
k},\lambda}+J_{{\bf q},\lambda}}{E_{{\bf q},\lambda}^{2}-E_{{\bf
k},\lambda}^{2}}\left[E_{{\bf k},\lambda}+\left(C_{2,{\bf q}{\bf
q}}^{+}-C_{1,{\bf q}{\bf q}}^{-}\right)E_{{\bf q},\lambda}\right].\end{split}$
(31)

A detailled account of how we solve these equations in practice is given
further below. Before turning to implementations, however, let us dwell on the
above equations to show how bound states emerge in our approach.

###  2.2 Bound states

From the renormalization equation (31), one finds that in the thermodynamic
limit a significant change from $E_{{\bf k},\lambda}$ to $E_{{\bf
k},\lambda-\Delta\lambda}$ only occurs if the coefficient $J_{{\bf
k},\lambda}$ is singular at a particular ${\bf k}$, i. e. if $J_{{\bf
k},\lambda}\propto N$. The reason is the factor $1/N^{2}$ which suppresses any
renormalization in the thermodynamic limit if $J_{{\bf k},\lambda}$ is non-
singular. Let us in the following explicitly consider the renormalization
equations at such a singular point ${\bf k}_{0}$ and therefore set ${\bf
k}={\bf k}_{0}$. Specializing to a system with only a single bound state, we
expect that no other energy is renormalized significantly, and therefore set
$E_{{\bf q},\lambda}=E_{{\bf q}}$ in the internal summations in $A_{{\bf
k}_{0},\lambda}^{(1-3)}$. This assumption has been checked to be valid in our
explicit implementations. Eqs. (30) and (20) imply that also $A_{{\bf
k}_{0},\lambda}^{(1)}$ and $V_{{\bf k}_{0},\lambda}$ become singular at ${\bf
k}_{0}$. Moreover, we can neglect $J_{{\bf q},\lambda}$ and $V_{{\bf
q},\lambda}$ with ${\bf q}\neq{\bf k}_{0}$ compared to $J_{{\bf
k}_{0},\lambda}$ and $V_{{\bf k}_{0},\lambda}$. Taking into account all these
simplifications, Eqs. (29) and (30) at ${\bf k}_{0}$ become

$\displaystyle-\frac{J}{2}=$ $\displaystyle-J_{{\bf
k}_{0},\lambda}+\frac{J\lambda}{2}\frac{V_{{\bf k}_{0},\lambda}}{N}\sum_{{\bf
q}(\neq{\bf k}_{0})}\frac{E_{{\bf k}_{0},\lambda}+\Delta-\varepsilon_{\bf
q}}{E_{{\bf q}}^{2}-E_{{\bf k}_{0},\lambda}^{2}}$ (32)
$\displaystyle+\frac{J\lambda}{2}\frac{J_{{\bf k}_{0},\lambda}}{N}\sum_{{\bf
q}(\neq{\bf k}),{\bf k}^{\prime}}\frac{\left(C_{1,{\bf q}{\bf
k}^{\prime}}^{+}+C_{2,{\bf q}{\bf k}^{\prime}}^{-}\right)E_{{\bf
q}}-\left(C_{2,{\bf q}{\bf k}^{\prime}}^{+}+C_{1,{\bf q}{\bf
k}^{\prime}}^{-}\right)E_{{\bf k}_{0},\lambda}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda}^{2}}$ $\displaystyle\times\langle\alpha_{\bf q}\alpha_{{\bf
k}^{\prime}}^{\dagger}-\beta_{\bf q}^{\dagger}\beta_{{\bf
k}^{\prime}}+\beta_{\bf q}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}+\alpha_{{\bf k}^{\prime}}\beta_{\bf
q}\rangle_{\lambda},$ $\displaystyle V_{{\bf k}_{0},\lambda}=$
$\displaystyle-\frac{J\lambda}{2}S(S+1)\frac{J_{{\bf
k}_{0},\lambda}}{N}\sum_{{\bf q}(\neq{\bf k}_{0})}\frac{E_{{\bf
k}_{0},\lambda}+\Delta-\varepsilon_{\bf q}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda}^{2}}.$ (33)

Let us now consider at first the renormalization in the very first
renormalization step from $\lambda=1$ to $\lambda=1-\Delta\lambda$
($\Delta\lambda\ll 1$), where we can set the $\lambda$ factor equal to 1.
Combining the two equations by elimination of $V_{{\bf k}_{0},\lambda}$ and
solving the resulting equation for $J_{{\bf k}_{0},\lambda}$, we obtain

$\begin{split}J_{{\bf
k}_{0},\lambda=1}&=\frac{J}{2}\bigg{[}1-S(S+1)\left(\frac{J}{2}\frac{1}{N}\sum_{{\bf
q}(\neq{\bf k}_{0})}\frac{E_{{\bf k}_{0},\lambda=1}+\Delta-\varepsilon_{\bf
q}}{E_{{\bf q}}^{2}-E_{{\bf k}_{0},\lambda=1}^{2}}\right)^{2}\\\
&+\frac{J}{2}\frac{1}{N}\sum_{{\bf q}(\neq{\bf k}_{0}),{\bf
k}^{\prime}}\frac{\left(C_{1,{\bf q}{\bf k}^{\prime}}^{+}+C_{2,{\bf q}{\bf
k}^{\prime}}^{-}\right)E_{{\bf q}}-\left(C_{2,{\bf q}{\bf
k}^{\prime}}^{+}+C_{1,{\bf q}{\bf k}^{\prime}}^{-}\right)E_{{\bf
k}_{0},\lambda=1}}{E_{{\bf q}}^{2}-E_{{\bf k}_{0},\lambda=1}^{2}}\\\
&\times\langle\alpha_{\bf q}\alpha_{{\bf k}^{\prime}}^{\dagger}-\beta_{\bf
q}^{\dagger}\beta_{{\bf k}^{\prime}}+\beta_{\bf q}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}+\alpha_{{\bf k}^{\prime}}\beta_{\bf
q}\rangle\bigg{]}^{-1}\end{split}$ (34)

Note that the expectation value now also refers to $\lambda=1$, i. e.
according to Eq. (25) it is taken with the full Hamiltonian ${\cal H}$. We
therefore omitted the index $\lambda$ in the expectation value. According to
Eq. (34), a singularity of $J_{{\bf k}_{0},\lambda=1}$, which in turn implies
a finite renormalization of $E_{{\bf k}_{0},\lambda=1}$, is found if the
equation

$\begin{split}1&=S(S+1)\left(\frac{J}{2}\frac{1}{N}\sum_{{\bf q}(\neq{\bf
k}_{0})}\frac{E_{{\bf k}_{0},\lambda=1}+\Delta-\varepsilon_{\bf q}}{E_{{\bf
q}}^{2}-E_{{\bf k}_{0},\lambda=1}^{2}}\right)^{2}\\\
&+\frac{J}{2}\frac{1}{N}\sum_{{\bf q}(\neq{\bf k}_{0}),{\bf
k}^{\prime}}\frac{\left(C_{1,{\bf q}{\bf k}^{\prime}}^{+}+C_{2,{\bf q}{\bf
k}^{\prime}}^{-}\right)E_{{\bf q}}-\left(C_{2,{\bf q}{\bf
k}^{\prime}}^{+}+C_{1,{\bf q}{\bf k}^{\prime}}^{-}\right)E_{{\bf
k}_{0},\lambda=1}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda=1}^{2}}\langle\alpha_{\bf q}\alpha_{{\bf
k}^{\prime}}^{\dagger}-\beta_{\bf q}^{\dagger}\beta_{{\bf
k}^{\prime}}+\beta_{\bf q}^{\dagger}\alpha_{{\bf
k}^{\prime}}^{\dagger}+\alpha_{{\bf k}^{\prime}}\beta_{\bf
q}\rangle\end{split}$ (35)

is fulfilled. The solution $E_{{\bf k}_{0},\lambda=1}$ of this equation is the
renormalized energy value due to the singularity. Note that the summation over
the expectation values in the second line in general prevents an analytical
solution. Therefore, this equation must be solved self-consistently. Only in
the special case of a classical spin, where the second line is absent, is an
analytical solution possible and yields the well-known YSR states (see
Appendix A). In general, one finds two solutions of Eq. (35), where the first
solution is a part of the quasiparticle band, i.e. $E_{{\bf
k}_{0}}\geq\Delta$, and the second solution lies inside the gap, $E_{{\bf
k}_{0},\lambda}<\Delta$. It follows that the renormalization starts at this
particular solution $E_{{\bf k}_{0}}$ and ends up with the second solution
$E_{{\bf k}_{0},\lambda}$ at some finite value of $\lambda$. All other states
at ${\bf q}\neq{\bf k}_{0}$ are unaffected and the corresponding energies
remain unchanged, $E_{{\bf q},\lambda}=E_{\bf q}$.

In summary, at the starting point $\lambda=1$, a particular state with quantum
number ${\bf k}_{0}$ from the quasiparticle band experiences a finite
renormalization as the corresponding transformation coefficients $J_{{\bf
k}_{0},\lambda}$, $V_{{\bf k}_{0},\lambda}$ become singular, i. e. they take
values proportional to $N$. Finally, we emphasise that the formation of a
bound state as described here is a general property of an impurity which is
embedded in a bath of fermions. This is shown in Appendix B for a toy model of
a non-magnetic impurity which is coupled to spinless fermions.

###  2.3 Expectation values

In addition to the energy parameters considered so far also a renormalization
equation for the $\lambda$-dependent expectation value in Eq. (27) has to be
found. In this expectation value, the quantum numbers ${\bf k}$ and ${\bf q}$
may also be equal (the expectation value then describes an occupation number
of the Bogoliubov quasiparticles). Let us begin by considering the expectation
value $\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf q}\rangle_{\lambda}$ for the
$\alpha$-type of the two quasiparticles (the corresponding expectation value
of the $\beta$ type can be found analogously). To derive the $\lambda$
dependence for this expectation value, we use that the trace of any operator
is invariant under a unitary transformation. According to Eq. (25), we can
thus write

$\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf q}\rangle_{\lambda}=\langle
e^{X_{\lambda,\Delta\lambda}}\alpha_{\bf
k}^{\dagger}e^{-X_{\lambda,\Delta\lambda}}e^{X_{\lambda,\Delta\lambda}}\alpha_{\bf
q}e^{-X_{\lambda,\Delta\lambda}}\rangle_{\lambda-\Delta\lambda},$ (36)

where Eq. (11) was used. The transformed single-particle operator up to linear
order in the small parameter $\Delta\lambda\ll 1$ is found from the expression
(12) of $X_{\lambda,\Delta\lambda}$,

$\begin{split}&e^{X_{\lambda,\Delta\lambda}}\alpha_{\bf
k}^{\dagger}e^{-X_{\lambda,\Delta\lambda}}\approx\alpha_{\bf
k}^{\dagger}+\frac{\Delta\lambda}{N}\sum_{{\bf
k}^{\prime}}\left\\{\left(J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}^{\prime}{\bf
k}}^{+}}{E_{{\bf k}^{\prime},\lambda}-E_{{\bf k},\lambda}}\left(\alpha_{{\bf
k}^{\prime}}^{\dagger}S_{z}-\beta_{{\bf k}^{\prime}}^{\dagger}\alpha_{\bf
k}S^{+}\right)\right.\right.\\\ -&\left.\left.\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}\left(\beta_{{\bf k}^{\prime}}S_{z}-\alpha_{{\bf
k}^{\prime}}S^{+}\right)\right]+\left(V_{{\bf k},\lambda}+V_{{\bf
k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}\alpha_{{\bf k}^{\prime}}^{\dagger}-\frac{C_{2,{\bf
k}{\bf k}^{\prime}}^{-}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}\beta_{{\bf k}^{\prime}}\right]\right\\}.\end{split}$
(37)

Plugging this expression into Eq. (36) we find the following relation between
the expectation values at $\lambda$ and the ones at the reduced parameter
$\lambda-\Delta\lambda$,

$\begin{split}\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf
q}\rangle_{\lambda}&\approx\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf
q}\rangle_{\lambda-\Delta\lambda}\\\ &+\frac{\Delta\lambda}{N}\sum_{{\bf
k}^{\prime}}\left(V_{{\bf q},\lambda}+V_{{\bf
k}^{\prime},\lambda}\right)\left(\frac{C_{1,{\bf q}{\bf
k}^{\prime}}^{+}}{E_{{\bf q},\lambda}-E_{{\bf
k}^{\prime},\lambda}}\langle\alpha_{\bf k}^{\dagger}\alpha_{{\bf
k}^{\prime}}\rangle_{\lambda-\Delta\lambda}-\frac{C_{2,{\bf q}{\bf
k}^{\prime}}^{-}}{E_{{\bf q},\lambda}+E_{{\bf
k}^{\prime},\lambda}}\langle\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}\rangle_{\lambda-\Delta\lambda}\right)\\\
&+\frac{\Delta\lambda}{N}\sum_{{\bf k}^{\prime}}\left(V_{{\bf
k},\lambda}+V_{{\bf k}^{\prime},\lambda}\right)\left(\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
k}^{\prime},\lambda}}\langle\alpha_{{\bf k}^{\prime}}^{\dagger}\alpha_{{\bf
q}}\rangle_{\lambda-\Delta\lambda}-\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{E_{{\bf k},\lambda}+E_{{\bf
k}^{\prime},\lambda}}\langle\beta_{{\bf k}^{\prime}}\alpha_{{\bf
q}}\rangle_{\lambda-\Delta\lambda}\right).\end{split}$ (38)

Similar equations describe the $\lambda$-dependence of the anomalous
expectation values $\langle\alpha_{\bf k}^{\dagger}\beta_{\bf
q}^{\dagger}\rangle_{\lambda}$, $\langle\beta_{\bf q}\alpha_{{\bf
k}}\rangle_{\lambda}$. The above set of linear equations provides a system of
renormalization equations for the $\lambda$-dependent expectation values. The
$\lambda$-dependence becomes particularly significant if one of the
transformation coefficients $J_{{\bf k},\lambda}$ and $V_{{\bf k},\lambda}$
becomes singular ($\propto N$), which is the case for the bound state ${\bf
k}_{0}$ (see previous section). The numerical evaluation of the
renormalization equations is carried out in parallel to the equations
(29)-(31) starting from $\lambda=1$ down to $\lambda=0$. The starting value
$\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf
q}\rangle_{\lambda=1}=\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf q}\rangle$ is
the scattering probability with respect to the original Hamiltonian. This
quantitiy has to be calculated self-consistently as explained in Sec. 3.

In addition to appearing in the renormalization equations, expectation values
are also central for physical observables. We therefore conclude this
technical section by detailing the general procedure of calculating
expectation values in our approach. We are interested in the expectation value
of some operator $\cal A$ formed with respect to the original Hamiltonian.
Formally, this corresponds to Eq. (25) at $\lambda=1$,

$\langle{\cal A}\rangle=\langle{\cal A}\rangle_{\lambda=1}=\frac{Tr\left({\cal
A}e^{-\beta{\cal H}}\right)}{Tr\left(e^{-\beta{\cal H}}\right)}.$ (39)

To evaluate Eq. (39), we use the invariance of the trace with respect to
unitary transformations and rewrite Eq. (39) introducing the stepwise unitary
transformations down to $\lambda=0$, where the effective Hamiltonian ${\cal
H}_{\lambda=0}=\tilde{\cal H}$ is diagonal, see Eq. (9). We find

$\langle{\cal A}\rangle=\frac{Tr\left(\tilde{\cal A}e^{-\beta\tilde{\cal
H}}\right)}{Tr\left(e^{-\beta\tilde{\cal H}}\right)}.$ (40)

The practical calculation of observables therefore requires the transformation
of the associated operator ${\cal A}$ using the same procedure as carried out
to find $\tilde{\cal H}$. To simplify the discussion, let us discuss this
method using the example of the occupation number $\langle\alpha_{\bf
k}^{\dagger}\alpha_{\bf k}\rangle$, which is the quantity entering all
renormalization equations. Eq. (40) then reads

$\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf
k}\rangle=\frac{Tr\left(\tilde{\alpha}_{\bf k}^{\dagger}\tilde{\alpha}_{\bf
k}e^{-\beta\tilde{\cal H}}\right)}{Tr\left(e^{-\beta\tilde{\cal H}}\right)}.$
(41)

Similar expressions can be derived for any other expectation value. The main
task is therefore to find the transformed operators $\tilde{\alpha}_{\bf
k}^{\dagger}$ by applying the same series of unitary transformation that was
used to transform the Hamiltonian. As for the latter, this results in
$\lambda$-dependent renormalizations of $\tilde{\alpha}_{\bf k}^{\dagger}$.
Eq. (37) suggests the ansatz

$\begin{split}\alpha_{{\bf k}\lambda}^{\dagger}&=x_{{\bf k}\lambda}\alpha_{\bf
k}^{\dagger}+\frac{y_{{\bf k}\lambda}}{N}\sum_{{\bf k}^{\prime}(\neq{\bf
k})}\left\\{\left(J_{{\bf k},\lambda}+J_{{\bf
k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}^{\prime}{\bf
k}}^{+}}{E_{{\bf k}^{\prime},\lambda}-E_{{\bf k},\lambda}}\left(\alpha_{{\bf
k}^{\prime}}^{\dagger}S_{z}-\beta_{{\bf
k}^{\prime}}^{\dagger}S^{+}\right)\right.\right.\\\
-&\left.\left.\frac{C_{2,{\bf k}{\bf k}^{\prime}}^{-}}{E_{{\bf
k},\lambda}+E_{{\bf k}^{\prime},\lambda}}\left(\beta_{{\bf
k}^{\prime}}S_{z}-\alpha_{{\bf k}^{\prime}}S^{+}\right)\right]-\left(V_{{\bf
k},\lambda}+V_{{\bf k}^{\prime},\lambda}\right)\left[\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{E_{{\bf k}^{\prime},\lambda}-E_{{\bf
k},\lambda}}\alpha_{{\bf k}^{\prime}}^{\dagger}+\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{E_{{\bf k},\lambda}+E_{{\bf k}^{\prime},\lambda}}\beta_{{\bf
k}^{\prime}}\right]\right\\},\end{split}$ (42)

where the initial conditions at $\lambda=1$ are $x_{{\bf k}\lambda=1}=1$ and
$y_{{\bf k}\lambda=1}=0$. The renormalization equations for the
$\lambda$-dependent coefficients $x_{{\bf k}\lambda}$ and $y_{{\bf k}\lambda}$
in Eq. (42) can be found by evaluating a similar equation as Eq. (11) but for
the renormalized single particle operator,

$\alpha_{{\bf k}\lambda-\Delta\lambda}^{\dagger}\approx
e^{X_{\lambda,\Delta\lambda}}\alpha_{{\bf
k}\lambda}^{\dagger}e^{-X_{\lambda,\Delta\lambda}},$ (43)

where Eq. (12) must be used for $X_{\lambda,\Delta\lambda}$. Evaluating again
the right hand side up to the linear order in $\Delta\lambda\ll 1$ and using
the above ansatz for $\alpha_{{\bf k}\lambda}^{\dagger}$, we find

$\begin{split}x_{{\bf k},\lambda-\Delta\lambda}&=x_{{\bf k},\lambda}-y_{{\bf
k}\lambda}\frac{\Delta\lambda}{N^{2}}\sum_{{\bf q}(\neq{\bf
k})}\left[S(S+1)\left(J_{{\bf k},\lambda}+J_{{\bf
q},\lambda}\right)^{2}+\left(V_{{\bf k},\lambda}+V_{{\bf
q},\lambda}\right)^{2}\right]\\\ &\times\left[\left(\frac{C_{1,{\bf k}{\bf
q}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
q},\lambda}}\right)^{2}+\left(\frac{C_{2,{\bf k}{\bf q}}^{-}}{E_{{\bf
k},\lambda}+E_{{\bf q},\lambda}}\right)^{2}\right].\end{split}$ (44)

The corresponding equation for $y_{{\bf k},\lambda}$ is found from a sum rule
which must be fulfilled at each $\lambda$,

$\begin{split}1&=x_{{\bf k}\lambda}^{2}+y_{{\bf
k}\lambda}^{2}\frac{1}{N^{2}}\sum_{{\bf q}(\neq{\bf
k})}\left[S(S+1)\left(J_{{\bf k},\lambda}+J_{{\bf
q},\lambda}\right)^{2}+\left(V_{{\bf k},\lambda}+V_{{\bf
q},\lambda}\right)^{2}\right]\\\ &\times\left[\left(\frac{C_{1,{\bf k}{\bf
q}}^{+}}{E_{{\bf k},\lambda}-E_{{\bf
q},\lambda}}\right)^{2}+\left(\frac{C_{2,{\bf k}{\bf q}}^{-}}{E_{{\bf
k},\lambda}+E_{{\bf q},\lambda}}\right)^{2}\right].\end{split}$ (45)

This sum rule results from the requirement that $[\alpha_{{\bf
k}\lambda}^{\dagger},\alpha_{{\bf k}\lambda}]_{+}=1$ holds. We find

$\begin{split}y_{{\bf k},\lambda-\Delta\lambda}^{2}&=\left(1-x_{{\bf
k}\lambda-\Delta\lambda}^{2}\right)/\frac{1}{N^{2}}\sum_{{\bf q}(\neq{\bf
k})}\left[S(S+1)\left(J_{{\bf k},\lambda-\Delta\lambda}+J_{{\bf
q},\lambda-\Delta\lambda}\right)^{2}+\left(V_{{\bf
k},\lambda-\Delta\lambda}+V_{{\bf
q},\lambda-\Delta\lambda}\right)^{2}\right]\\\
&\times\left[\left(\frac{C_{1,{\bf k}{\bf q}}^{+}}{E_{{\bf
k},\lambda-\Delta\lambda}-E_{{\bf
q},\lambda-\Delta\lambda}}\right)^{2}+\left(\frac{C_{2,{\bf k}{\bf
q}}^{-}}{E_{{\bf k},\lambda-\Delta\lambda}+E_{{\bf
q},\lambda-\Delta\lambda}}\right)^{2}\right].\end{split}$ (46)

To calculate expectation values, the set of renormalization equations (44) and
(46) has to be solved simultaneous with the corresponding equations for the
Hamiltonian up to $\lambda=0$. This yields the fully renormalized single-
particle operator at $\lambda=0$ as obtained from Eq. (42),

$\begin{split}\tilde{\alpha}_{{\bf k}}^{\dagger}&=\tilde{x}_{{\bf
k}}\alpha_{\bf k}^{\dagger}+\frac{\tilde{y}_{{\bf k}}}{N}\sum_{{\bf
k}^{\prime}(\neq{\bf k})}\left\\{\left(\tilde{J}_{{\bf k}}+\tilde{J}_{{\bf
k}^{\prime}}\right)\left[\frac{C_{1,{\bf k}^{\prime}{\bf
k}}^{+}}{\tilde{E}_{{\bf k}^{\prime}}-\tilde{E}_{{\bf k}}}\left(\alpha_{{\bf
k}^{\prime}}^{\dagger}S_{z}-\beta_{{\bf
k}^{\prime}}^{\dagger}S^{+}\right)-\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{\tilde{E}_{{\bf k}}+\tilde{E}_{{\bf
k}^{\prime}}}\left(\beta_{{\bf k}^{\prime}}S_{z}-\alpha_{{\bf
k}^{\prime}}S^{+}\right)\right]\right.\\\ &-\left.\left(\tilde{V}_{{\bf
k}}+\tilde{V}_{{\bf k}^{\prime}}\right)\left[\frac{C_{1,{\bf k}{\bf
k}^{\prime}}^{+}}{\tilde{E}_{{\bf k}^{\prime}}-\tilde{E}_{{\bf
k}}}\alpha_{{\bf k}^{\prime}}^{\dagger}+\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{\tilde{E}_{{\bf k}}+\tilde{E}_{{\bf k}^{\prime}}}\beta_{{\bf
k}^{\prime}}\right]\right\\}.\end{split}$ (47)

Note that a possible bound state ${\bf k}_{0}$ may be also included in the
summation. If we consider a state with ${\bf k}\neq{\bf k}_{0}$, we find

$\begin{split}\tilde{\alpha}_{{\bf k}}^{\dagger}&=\tilde{x}_{{\bf
k}}\alpha_{\bf k}^{\dagger}+\frac{\tilde{y}_{{\bf k}}}{N}\sum_{{\bf
k}^{\prime}(\neq{\bf k},{\bf k}_{0})}\left\\{\left(\tilde{J}_{{\bf
k}}+\tilde{J}_{{\bf k}^{\prime}}\right)\left[\frac{C_{1,{\bf k}^{\prime}{\bf
k}}^{+}}{E_{{\bf k}^{\prime}}-E_{{\bf k}}}\left(\alpha_{{\bf
k}^{\prime}}^{\dagger}S_{z}-\beta_{{\bf
k}^{\prime}}^{\dagger}S^{+}\right)-\frac{C_{2,{\bf k}{\bf
k}^{\prime}}^{-}}{E_{{\bf k}}+E_{{\bf k}^{\prime}}}\left(\beta_{{\bf
k}^{\prime}}S_{z}-\alpha_{{\bf k}^{\prime}}S^{+}\right)\right]\right.\\\
&-\left.\left(\tilde{V}_{{\bf k}}+\tilde{V}_{{\bf
k}^{\prime}}\right)\left[\frac{C_{1,{\bf k}{\bf k}^{\prime}}^{+}}{E_{{\bf
k}^{\prime}}-E_{{\bf k}}}\alpha_{{\bf k}^{\prime}}^{\dagger}+\frac{C_{2,{\bf
k}{\bf k}^{\prime}}^{-}}{E_{{\bf k}}+E_{{\bf k}^{\prime}}}\beta_{{\bf
k}^{\prime}}\right]\right\\}\\\ &+\frac{C_{1,{\bf k}_{0}{\bf
k}}^{+}}{\tilde{E}_{{\bf k}_{0}}-E_{{\bf k}}}\left(\tilde{J}\alpha_{{\bf
k}_{0}}^{\dagger}S_{z}-\tilde{J}\beta_{{\bf
k}_{0}}^{\dagger}S^{+}-\tilde{V}\alpha_{{\bf
k}_{0}}^{\dagger}\right)-\frac{C_{2,{\bf k}{\bf k}_{0}}^{-}}{E_{{\bf
k}}+\tilde{E}_{{\bf k}_{0}}}\left(\tilde{J}\beta_{{\bf
k}_{0}}S_{z}-\tilde{J}\alpha_{{\bf k}_{0}}S^{+}+\tilde{V}\beta_{{\bf
k}_{0}}\right),\end{split}$ (48)

where the parameters $\tilde{J}$ and $\tilde{V}$ are the renormalized coupling
parameters of the bound state. They fulfil the relations $\tilde{J}_{{\bf
k}_{0}}=N\tilde{J}$ and $\tilde{V}_{{\bf k}_{0}}=N\tilde{V}$. It is seen that
the renormalized quasiparticle operator consists of three contributions. The
first part, proportional to $\tilde{x}_{{\bf k}}$, describes a coherent
excitation reminiscent of the original superconducting quasiparticles. As we
show in Eqs. (48) and (54) below, this part of the renormalized quasiparticle
operator also lives at the energy $E_{\bf{k}}$. This implies that the quantum
number $\mathbf{k}$ can then still be interpreted as a momentum (after all, we
consider a system with a single impurity, such that the breaking of
translational invariance by the localized impurity does not dramatically
impact many of the states in our bulk system). The second part, proportional
to $\tilde{y}_{{\bf k}}$, results from the scattering of an incoming
quasiparticle at ${\bf k}$ into an outgoing state ${\bf k}^{\prime}\neq{\bf
k},{\bf k}_{0}$. This part describes a continuum of excitations which is
governed by the renormalized exchange coupling $\tilde{J}_{{\bf k}}$. Its
operator strucutre, e.g. $\alpha_{{\bf k}^{\prime}}^{\dagger}S_{z}-\beta_{{\bf
k}^{\prime}}^{\dagger}S^{+}$, is well-known in the context of Kondo physics
[39], where operators of this form describe composite fermion excitations. As
we further explain below, these terms indeed correspond to a Kondo resonance.
The remaining terms, finally, are associated with ${\bf k}_{0}$, and therefore
describe the formation of a YSR bound state appearing at the (renormalized)
energy $\tilde{E}_{{\bf k}_{0}}$ inside the superconducting gap.

To illustrate this general scheme, let us consider two specific examples of
expectation values which are needed for our discussions. First, we consider
the static occupation number of the superconducting quasiparticles. Using Eqs.
(41) and (48), the occupation number reads

$\begin{split}\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf
k}\rangle&=\tilde{x}_{{\bf k}}^{2}f(\tilde{E}_{\bf k})+\tilde{y}_{{\bf
k}}^{2}\frac{1}{N^{2}}\sum_{{\bf q}(\neq{\bf k},{\bf
k}_{0})}\left[S(S+1)\left(\tilde{J}_{{\bf k}}+\tilde{J}_{{\bf
q}}\right)^{2}+\left(\tilde{V}_{{\bf k}}+\tilde{V}_{{\bf
q}}\right)^{2}\right]\\\ &\times\left[\left(\frac{C_{1,{\bf k}{\bf
q}}^{+}}{E_{{\bf k}}-E_{{\bf q}}}\right)^{2}f(E_{\bf q})+\left(\frac{C_{2,{\bf
k}{\bf q}}^{-}}{E_{{\bf k}}+E_{{\bf q}}}\right)^{2}\left(1-f(E_{\bf
q})\right)\right]\\\
&+\left[S(S+1)\tilde{J}^{2}+\tilde{V}^{2}\right]\left[\left(\frac{C_{1,{\bf
k}{\bf k}_{0}}^{+}}{E_{{\bf k}}-\tilde{E}_{{\bf
k}_{0}}}\right)^{2}f(\tilde{E}_{{\bf k}_{0}})+\left(\frac{C_{2,{\bf k}{\bf
k}_{0}}^{-}}{E_{{\bf k}}+\tilde{E}_{{\bf
k}_{0}}}\right)^{2}\left(1-f(\tilde{E}_{{\bf
k}_{0}})\right)\right],\end{split}$ (49)

where $f(E)=1/(1+e^{\beta E})$ is the Fermi function with respect to an energy
$E$.

As a second example, we turn to the single-particle spectral function of the
original electrons as a function of ${\mathbf{k}}$. In the absence of a
coupling to the quantum spin, $J=0$, this quantum number corresponds to the
electron’s momentum. The presence of the impurity spin breaks translational
invariance, such that ${\mathbf{k}}$ strictly speaking cannot be interpreted
as a momentum anymore. However, as argued above and illustrated by Eqs. (48)
and (54), most of the states in the considered bulk system are not affected by
the presence of a single impurity. We can therefore in practice understand
${\mathbf{k}}$ as a quantum number that for most states corresponds to
momentum. The single-particle spectral function as a function of
${\mathbf{k}}$ in that sense still allows to visualize the excitation
structure in momentum space. It is defined via the imaginary part of the
single-particle Green’s function,

$\Im G({\bf k},{\bf
k},\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\langle[c_{{\bf
k}\sigma}^{\dagger}(-t),c_{{\bf k}\sigma}]_{+}\rangle e^{i\omega t}dt.$ (50)

The time-dependence and the expectation value is calculated using the original
Hamiltonian ${\cal H}$, i.e. the expectation value is defined as in Eq. (39).
Using the Mori-Zwanzig formalism [37, 38], the anticommutator correlation
function can be at first rewritten to the form

$\Im G({\bf k},{\bf k},\omega)=\langle[c_{{\bf
k}\sigma},\delta(\mathbf{L}-\omega)c_{{\bf k}\sigma}^{\dagger}]_{+}\rangle,$
(51)

where $\mathbf{L}\mathcal{A}=[\mathcal{H},\mathcal{A}]$ is the Liouville
operator. For occupied states ${\bf k}$, this correlation function can be
probed by angle resolved photo emission (ARPES). Our method to calculate this
quantity is by using the invariance of the trace under the unitary
transformation as explained above in the context of Eqs. (40) and (41). This
yields

$\Im G({\bf k},{\bf k},\omega)=\langle[\tilde{c}_{{\bf
k}\sigma},\delta(\tilde{\mathbf{L}}-\omega)\tilde{c}_{{\bf
k}\sigma}^{\dagger}]_{+}\rangle.$ (52)

We evaluate the Liouville operator and the expectation value by replacing the
electron operators with the Bogoliubov quasiparticles using Eq. (2) and obtain

$\begin{split}\Im G({\bf k},{\bf k},\omega)=\langle[u_{\bf
k}\tilde{\alpha}_{{\bf k}}+v_{\bf k}\tilde{\beta}_{{\bf
k}}^{\dagger},\delta(\tilde{\mathbf{L}}-\omega)(u_{\bf k}\tilde{\alpha}_{{\bf
k}}^{\dagger}+v_{\bf k}\tilde{\beta}_{{\bf k}})]_{+}\rangle,\end{split}$ (53)

where for the sake of concreteness, we have here considered the spin-up case.
Finally, it remains to replace the fully renormalized single-particle
operators $\tilde{\alpha}_{{\bf k}}$ with Eq. (48) and the corresponding
equation for $\tilde{\beta}_{{\bf k}}$ (not explicitly given). Since the
effective Hamiltonian $\tilde{\cal H}$ which is used to calculate the
expectation value is diagonal with respect to the $\alpha$ and $\beta$
operators, one can simply replace the Liouville operator with the
corresponding eigenenergy. We obtain

$\begin{split}&\Im G({\bf k},{\bf k},\omega)=\tilde{x}_{\bf k}^{2}\left[u_{\bf
k}^{2}\delta(E_{\bf k}-\omega)+v_{\bf k}^{2}\delta(-E_{\bf
k}-\omega)\right]\\\ &+\frac{\tilde{y}_{\bf k}^{2}}{N^{2}}\sum_{{\bf
q}(\neq{\bf k})}\left[S(S+1)\left(\tilde{J}_{{\bf k}}+\tilde{J}_{{\bf
q}}\right)^{2}+\left(\tilde{V}_{{\bf k}}+\tilde{V}_{{\bf
q}}\right)^{2}\right]\left\\{\left[\left(\frac{u_{\bf k}C_{1,{\bf k}{\bf
q}}^{+}}{E_{{\bf k}}-\tilde{E}_{{\bf q}}}\right)^{2}+\left(\frac{v_{\bf
k}C_{2,{\bf k}{\bf q}}^{-}}{E_{{\bf k}}+\tilde{E}_{{\bf
q}}}\right)^{2}\right]\delta(\tilde{E}_{\bf q}+\omega)\right.\\\
&+\left.\left[\left(\frac{u_{\bf k}C_{2,{\bf k}{\bf q}}^{-}}{E_{{\bf
k}}+\tilde{E}_{{\bf q}}}\right)^{2}+\left(\frac{v_{\bf k}C_{1,{\bf k}{\bf
q}}^{+}}{E_{{\bf k}}-\tilde{E}_{{\bf
q}}}\right)^{2}\right]\delta(-\tilde{E}_{\bf q}+\omega)\right\\}.\end{split}$
(54)

The first line describes the usual coherent dispersive excitation of a
superconducting quasiparticle. The summation over ${\bf q}$ in the second and
third line describes a continuum of excitations associated with scattering off
the impurity, and the Kondo effect. It also contains the YSR bound state at
${\bf k}_{0}$, for which the coefficients $\tilde{J}_{{\bf q}}$ and
$\tilde{V}_{{\bf q}}$ are singular (proportional $N$). In this case a
dispersionless excitation appears at energies $\pm\tilde{E}_{{\bf k}_{0}}$.
Moreover, the renormalized exchange energy $\tilde{J}_{{\bf q}}$ may become
very large around the Fermi level due to the contributions from the
expectation values in the term $A_{{\bf q}\lambda}^{(3)}$ in the
renormalization equation (29). These contributions are assigned to the Kondo
resonance [33].

##  3 Numerical results

We conclude by using the above-described self-consistent scheme to calculate a
number of physical observables. The numerical calculation is carried out in
the thermodynamic limit, $N\rightarrow\infty$. Thus, all summations over
vectors ${\bf q}\neq{\bf k}_{0}$ are calculated by replacing the sum with an
energy integral using the BCS density of states,

$\frac{1}{N}\sum_{{\bf q}(\neq{\bf
k}_{0})}\rightarrow\nu_{0}\int_{\Delta}^{W}\Theta(|E|-\Delta)\frac{|E|}{\sqrt{E^{2}-\Delta^{2}}}dE,$
(55)

where $\nu_{0}$ is the density of states in the normal state near the Fermi
level and $W$ is the bandwidth of the conduction electron band. This
description in the continuum limit is possible because for ${\bf q}\neq{\bf
k}_{0}$ the quasiparticle energy remains unchanged. The introduced density of
states fulfills the sum rule condition

$N=\nu_{0}\int_{\Delta}^{W}\Theta(|E|-\Delta)\frac{|E|}{\sqrt{E^{2}-\Delta^{2}}}dE.$
(56)

The renormalization process starts with a guess for the expectation values
$\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf k^{\prime}}\rangle$. Using these
values, the bound state quantum number ${\bf k}_{0}$ and the corresponding
energy $\tilde{E}_{{\bf k}_{0}}$ are found by evaluating Eq. (35). Next, the
renormalization equations (29) and (30) as well as Eqs. (38), (44), and (46)
are evaluated stepwise starting from $\lambda=1$ down to $\lambda=0$. The step
width is chosen to be $\Delta\lambda=10^{-3}$. The starting values are
$E_{{\bf k},\lambda=1}=E_{\bf k}$, $x_{{\bf k},\lambda=1}=1$, and $y_{{\bf
k},\lambda=1}=0$. In each $\lambda$-step, Eqs. (29) and (30) are solved self-
consistently for $J_{{\bf k},\lambda}$ and $V_{{\bf k},\lambda}$,
respectively. At the end of this procedure, we arrive with the fully
renormalized values $\tilde{E}_{{\bf k}_{0}}$, $\tilde{J}_{\bf k}$,
$\tilde{V}_{\bf k}$, $\tilde{x}_{\bf k}$, and $\tilde{y}_{\bf k}$. With all
these quantities, we calculate the renormalized quasiparticle operator using
Eq. (48), and, correspondingly, we are able to recalculate the expectation
values $\langle\alpha_{\bf k}^{\dagger}\alpha_{\bf k^{\prime}}\rangle$ using
Eq. (41).

Next, we restart the entire renormalization process using these improved
expectation values, and again solve all renormalization equations down to
$\lambda=0$. We conclude the iteration by recalculating the expectation
values. This procedure is repeated until self-consistency is obtained. We
define self-consistency to be fulfilled if the relative difference between the
expectation values in two subsequent procedures is less than $10^{-3}$ for all
${\mathbf{k}}$.

For the evaluation of the spectral function which is given by Eq. (54) we
describe the $\delta$-functions with Lorentzian functions with a manually-
imposed broadening of $\eta=0.03\Delta$. This value is taken in all figures
showing spectral functions.

Figure 2:  (a) Energy of the YSR state as a function of the Kondo temperature.
(b) Calculated density of states at larger values of the Kondo temperature. A
strong enhancement of the spectral intensity around the normal state Fermi
level is found, which marks the crossover into a Kondo resonance regime. The
inset shows the spectral weight of the YSR bound state around the values of
$k_{\text{B}}T_{\text{K}}/\Delta$ where the bound state crossed zero energy.
In agreement with literature results[28], a jump-like behaviour is found.

###  3.1 YSR state energy and Kondo resonance

In Fig. 2(a), we show the result of the self-consistent solution for the
energy $\tilde{E}_{{\bf k}_{0}}$ at $T=0$ as a function of the Kondo
temperature $k_{\text{B}}T_{\text{K}}$. Here, the Kondo temperature is varied
by changing the exchange coupling parameter $J$. The relationship between
$T_{\text{K}}$ and $J$ is

$k_{\text{B}}T_{\text{K}}=W\sqrt{J\nu_{0}}\exp\left(-\frac{1}{J\nu_{0}}\right).$
(57)

The energy of the second YSR excitation at $-\tilde{E}_{{\bf k}_{0}}$ is also
shown (red line). One clearly recognizes the behavior known from former
numerical studies [28, 6] including the crossing at $\tilde{E}_{{\bf
k}_{0}}=0$ (dashed line) which marks the transition from a weak coupling
regime at low $T_{\text{K}}$ to a strong coupling regime at larger
$T_{\text{K}}$ where the Kondo physics becomes significant. This regime is
shown in Fig. 2(b) where the integrated spectral function from Eq. (54) is
shown for three larger values of the Kondo temperature. It is seen that the
spectral intensity is strongly enhanced around the superconducting gap edge.
The width of this resonance peak scales with the Kondo temperature, which is a
typical behavior of the Kondo resonance. The YSR states are pinned to the gap
edges since the exchange coupling is already relatively large. In agreement
with early numerical studies [28], our calculations also show that the
spectral weight of the YSR state within the gap (shown in the inset) also
increases with the Kondo temperature, and exhibits a characteristic jump when
crossing zero energy.

###  3.2 Spectral function

To illustrate the behavior of the YSR state in the spectral function in more
detail, we considered three further values of $J/\Delta$ corresponding to a
rather weak coupling in the YSR regime, $J/\Delta=20$ ($k_{B}T_{K}\approx
0.003$), an intermediate coupling $J/\Delta=55$ ($k_{B}T_{K}\approx 0.120$),
and a large value $J/\Delta=120$ ($k_{B}T_{K}\approx 0.476$) in the Kondo
regime.

Figure 3:  Total density of states for three different values of $J$ in an
energy region around the Fermi level. These results where obtained by
$\mathbf{k}$-summation of both parts ${\bf k}={\bf k}_{0}$ and ${\bf
k}\neq{\bf k}_{0}$ of the spectral function. The superconducting gap is
clearly seen. The bound state inside the gap is shown in a different color.

First, we turn to the $\mathbf{k}$-integrated spectral function including the
bound state quantum number $\mathbf{k}_{0}$. Fig. 3 shows the spectral density
in a narrow region around the Fermi level. The superconducting gap and the
density of states of the conduction electrons are not found to vary noticeably
with $J$. The excitation peak of the YSR state (shown in red color), on the
other hand, moves in energy as $J$ changes. Furthermore, the spectral
intensity of the YSR state changes, in particular, once the Kondo regime is
entered (compare also the inset in Fig. 2(b)). As a crosscheck, we verified
that the energy integral over the total density of states yields the total
number of states $N$, where in the continuum limit considered here $N$ is
given by the integral over the density of states in Eq. (56).

More interesting is the $\mathbf{k}$-resolved spectral function, which for
$\mathbf{k}\neq\mathbf{k}_{0}$ is shown in Fig. 4. In the weak coupling YSR
regime at $J/\Delta=20$, the spectrum is composed of a dispersive excitation
and a YSR bound state located close to the gap edge. The dispersive excitation
loses spectral weight once the peak approaches $k_{F}$, while the spectral
weight of the YSR state increases.

Figure 4:  Spectral functions for the same parameters as in Fig. 3 but only
for ${\bf k}\neq{\bf k}_{0}$. The bound state at $\tilde{E}_{{\bf k}_{0}}$
becomes extended in $\mathbf{k}$, and the distribution of the spectral weight
gets inverted when the Kondo regime is entered.

A larger value of $J$ is considered in Fig. 4(b). The two YSR peaks move
closer together, and become somewhat more intense. At the same time, the
coherent excitation weakens. The YSR state now also has a somewhat larger
extent in $\mathbf{k}$, but still remains relatively close to $k_{F}$.

If, finally, the exchange coupling is increased even more, the YSR states
cross zero energy. As shown in Fig. 4(c), the bound state spectral weight
coming from momenta $\mathbf{k}\neq\mathbf{k}_{0}$ is now much larger than
that in the weak coupling YSR regime.

###  3.3 Influence of temperature and spin size

To further benchmark our approach against purely numerical methods, we now
turn to the dependence of the bound state on temperature and spin size. Fig. 5
shows results of the YSR bound state energy and weight as a function of the
temperature $T$ and the size $S$ of the local spin. Panel (a) shows the
influence of the temperature $T$ in the two different regimes at $J/\Delta=20$
and $J/\Delta=120$. In agreement with the numerical results from Ref. [6], we
find that the YSR bound state energy is shifted to smaller values at small
exchange couplings. This effect becomes inverted in the Kondo regime, where
the shift is smaller but the energy increases with temperature. This is seen
even more clearly in the panels (b,d) where we have directly calculated the
$\mathbf{k}$-integrated spectral function in a region inside the
superconducting gap. If we compare the zero temperature result (blue line)
with the one at the somewhat higher temperature $k_{B}T=0.35$ below $T_{c}$,
we see that the pair of YSR peaks comes closer together at higher temperature
in the weak coupling case, while the opposite is the case at strong coupling.
This is exactly the behavior found in Ref. [6]. Even the slight change of the
intensity with respect to temperature agrees with the mentioned numerical
study.

Fig. 5(c) finally depicts the dependence of the YSR state energy on the size
of the quantum spin at fixed $T$ (below the critical temperature $T_{c}$). We
observe a strong reduction of $\tilde{E}_{{\bf k}_{0}}$ with $S$ in all
regimes. As expected, the most significant effect of quantum fluctuations is
obtained at small values of $S$ where the energy $\tilde{E}_{{\bf k}_{0}}$
deviates most strongly from the classical spin limit (dashed lines). In the
YSR regime, however, the effect of the quantum nature of the spin is indeed
rather weak, but it extends over a wider range of $S$ than in the Kondo
regime.

Figure 5:  Dependence of the YSR bound state energy of a quantum spin on
temperature and spin size. (a) Variation of the temperature. In agreement with
former numerical studies[6], the YSR state energy decreases with temperature
while in the Kondo regime it slightly increases. (b,d)
${\mathbf{k}}$-integrated spectral function showing the pair of YSR bound
states in a narrow energy region inside the superconducting gap. Shown are the
spectra at two temperatures (blue and red line) and two different Kondo
temperatures (upper/lower panel). Again a very good agreement with Ref. [6] is
found. (c) Variation of the spin for the three different values of $J$
discussed above. The corresponding results for the classical spin are also
shown (dashed lines). $\tilde{E}_{{\bf k}_{0}}$ approaches the value of the
classical spin case in the limit of large $S$ in all cases.

##  4 Conclusion

In this manuscript, we have presented a new approach to solve the paradigmatic
Hamiltonian describing a local quantum spin coupled to a superconducting
substrate. The basic idea of our approach is to eliminate the interaction
between electrons and the impurity by unitary transformations that yields a
diagonal renormalized Hamiltonian. In this process, bound states manifest
themselves as singularities in the unitary transformation. Postponing a
detailed derivation to a subsequent publication, we stress that our method can
straightforwardly be extended to larger numbers of impurities. This has been
worked out already for the Kondo lattice model in the normal state [33].

We first illustrated how the spin-fermion exchange interaction is integrated
out using a stepwise unitary transformation. Performing the transformation in
small steps allowed us consider the lowest order commutators in a well-
controlled way. We have derived renormalization equations for the
quasiparticle energy as well as operators. We have shown that our method fully
captures Yu-Shiba-Rusinov bound states and the Kondo resonance within one
uniform theoretical framework. While the YSR state is a $\delta$ peak inside
the superconducting gap, the Kondo resonance appears outside the gap as a peak
of finite width which is proportional to the Kondo temperature. Our approach
hence draws a clear picture of the physics in the superconducting Kondo model
by capturing on equal footing both the classical scattering problem of a
magnetic impurity, and also the complex many-body effect leading to the Kondo
resonance. In addition, we showed that our approach also yields access to
effectively momentum-resolved observables, which for example are important in
the context of ARPES-experiments.

In contrast to purely numerical approaches, our method also allows for
detailed analytical insights. This power of our approach is rooted in the fact
that bound states manifest themselves as singularities, and that it is
therefore sufficient to extract this singular behavior from the full picture
to investigate bound state physics. This for example allows us to explicitely
construct bound state quasiparticle operators, which in turn reveals how
increasing quantum fluctuations of the local spin pushes the bound state from
a YSR-like form to a Kondo-like expression.

## Acknowledgements

We would like to thank K.W. Becker for helpful discussions.

#### Funding information

This project has received funding from the Deutsche Forschungsgemeinschaft via
the Emmy Noether Programme ME4844/1-1 (project id 327807255), the
Collaborative Research Center SFB 1143 (project id 247310070), and the Cluster
of Excellence on Complexity and Topology in Quantum Matter ct.qmat (EXC 2147,
project id 390858490).

##  Appendix A Classical spin

Here we show that the specific case of a classical spin leads to an exactly
solvable scattering problem which gives rise to the well known YSR bound
states in the absence of the Kondo resonance. To describe a classical spin,
the spin operator ${\mathbf{S}}$ is replaced by a constant vector. We choose
the spin to be oriented along $z$ and thus set
${\mathbf{S}}=S\,\mathbf{e}_{z}$ ($\mathbf{e}_{z}$ being the unit vector along
$z$). The electron-impurity Hamiltonian then takes the form

$\begin{split}\mathcal{H}_{1,cl}&=\frac{JS}{2N}\sum_{{\bf k}{\bf
k}^{\prime}}\left[\left(u_{\bf k}v_{{\bf k}^{\prime}}-u_{{\bf
k}^{\prime}}v_{\bf k}\right)\left(\alpha_{\bf k}^{\dagger}\beta_{{\bf
k}^{\prime}}^{\dagger}+\beta_{{\bf k}^{\prime}}\alpha_{\bf
k}\right)+\left(u_{\bf k}u_{{\bf k}^{\prime}}+v_{{\bf k}^{\prime}}v_{\bf
k}\right)\left(\alpha_{\bf k}^{\dagger}\alpha_{{\bf k}^{\prime}}-\beta_{\bf
k}^{\dagger}\beta_{{\bf k}^{\prime}}\right)\right].\end{split}$ (58)

Because a classical spin has no quantum fluctuations, the Hamiltonian is
entirely quadratic in operators.

Applying our approach described above to the scattering problem of the
classical spin simply means that the last part in Eq. (28) is omitted since it
arises from the commutator between spin operators, which is zero in the
classical spin case. Furthermore, the number $S(S+1)$ is replaced with $S^{2}$
according to the convention in Eq. (58). Thus, due to $A_{{\bf
k},\lambda}^{(3)}=0$ the renormalization equations (29) and (30) simplify to

$\displaystyle-\frac{J}{2}$ $\displaystyle=-J_{{\bf
k},\lambda}+\frac{J\lambda}{2}\frac{1}{N}\sum_{{\bf q}(\neq{\bf
k})}\frac{V_{{\bf k},\lambda}+V_{{\bf q},\lambda}}{E_{{\bf
q},\lambda}^{2}-E_{{\bf k},\lambda}^{2}}\left(E_{{\bf
k},\lambda}+\frac{\Delta-\varepsilon_{\bf q}}{E_{\bf q}}E_{{\bf
q},\lambda}\right),$ (59) $\displaystyle 0$ $\displaystyle=V_{{\bf
k},\lambda}-\frac{J\lambda}{2}S^{2}\frac{1}{N}\sum_{{\bf q}(\neq{\bf
k})}\frac{J_{{\bf k},\lambda}+J_{{\bf q},\lambda}}{E_{{\bf
q},\lambda}^{2}-E_{{\bf k},\lambda}^{2}}\left(E_{{\bf
k},\lambda}+\frac{\Delta-\varepsilon_{\bf q}}{E_{\bf q}}E_{{\bf
q},\lambda}\right).$ (60)

The renormalization equations (32) and (33) for the YSR bound state read for
the classical spin case

$\displaystyle-\frac{J}{2}$ $\displaystyle=-J_{{\bf
k}_{0},\lambda}+\frac{J\lambda}{2}\frac{V_{{\bf k}_{0},\lambda}}{N}\sum_{{\bf
q}(\neq{\bf k}_{0})}\frac{E_{{\bf k}_{0},\lambda}+\Delta-\varepsilon_{\bf
q}}{E_{{\bf q}}^{2}-E_{{\bf k}_{0},\lambda}^{2}},$ (61) $\displaystyle V_{{\bf
k}_{0},\lambda}$ $\displaystyle=-\frac{J\lambda}{2}S^{2}\frac{J_{{\bf
k}_{0},\lambda}}{N}\sum_{{\bf q}(\neq{\bf k}_{0})}\frac{E_{{\bf
k}_{0},\lambda}+\Delta-\varepsilon_{\bf q}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda}^{2}}.$ (62)

Since we consider the renormalization in the very first renormalization step
we can set the $\lambda$ factor equal to 1. Then, combining the two equations
by elimination of $V_{{\bf k}_{0},\lambda}$ and solving the resulting equation
for $J_{{\bf k}_{0},\lambda}$ we obtain

$J_{{\bf k}_{0},\lambda}=\frac{J/2}{1-\left(\frac{J}{2}\frac{S}{N}\sum_{{\bf
q}(\neq{\bf k}_{0})}\frac{E_{{\bf k}_{0},\lambda}+\Delta-\varepsilon_{\bf
q}}{E_{{\bf q}}^{2}-E_{{\bf k}_{0},\lambda}^{2}}\right)^{2}}.$ (63)

It is clearly seen that a singularity of $J_{{\bf k}_{0},\lambda}$, which is
the condition for a renormalization of $E_{{\bf k}_{0},\lambda}$, is found if
the equation

$1=\frac{J}{2}\frac{S}{N}\sum_{{\bf q}(\neq{\bf k}_{0})}\frac{E_{{\bf
k}_{0},\lambda}+\Delta-\varepsilon_{\bf q}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda}^{2}}$ (64)

is fulfilled. The solution $E_{{\bf k}_{0},\lambda}$ of this equation is the
renormalized energy value due to the singularity. To explicitly calculate
$E_{{\bf k}_{0},\lambda}$, we use the BCS density of states,

$\nu(E)=\nu_{0}\Theta(|E|-\Delta)\frac{|E|}{\sqrt{E^{2}-\Delta^{2}}},$ (65)

where $\nu_{0}$ is the density of states at the Fermi level in the normal
state. We also express the summation in terms of an energy integration and
calculate the corresponding integral analytically. The result is

$\frac{1}{N}\sum_{{\bf q}(\neq{\bf k}_{0})}\frac{E_{{\bf
k}_{0},\lambda}+\Delta-\varepsilon_{\bf q}}{E_{{\bf q}}^{2}-E_{{\bf
k}_{0},\lambda}^{2}}=\nu_{0}\pi\frac{\Delta+E_{{\bf
k}_{0},\lambda}}{\sqrt{\Delta^{2}-E_{{\bf k}_{0},\lambda}^{2}}}.$ (66)

Thus, the renormalized energy $E_{{\bf k}_{0},\lambda}$ is given by the
solution of the equation

$1-\frac{JS}{2}\nu_{0}\pi\frac{\Delta+E_{{\bf
k}_{0},\lambda}}{\sqrt{\Delta^{2}-E_{{\bf k}_{0},\lambda}^{2}}}=0.$ (67)

This equation has solutions for $E_{{\bf k}_{0},\lambda}$ only in the
superconducting state, i. e. if $\Delta\neq 0$. As described in Sec. 2.2 there
are two solutions of Eq. (67) where one lies in the conduction band, $E_{{\bf
k}_{0}}>\Delta$, which is the starting value of the bound state
renormalization, and the other one, $E_{{\bf k}_{0},\lambda}<\Delta$, is
inside the gap. The latter reads

$E_{{\bf k}_{0},\lambda}=\tilde{E}_{{\bf
k}_{0}}=\Delta\frac{1-\left(\frac{JS\nu_{0}\pi}{2}\right)^{2}}{1+\left(\frac{JS\nu_{0}\pi}{2}\right)^{2}},$
(68)

which is the final value of the renormalization and therefore it corresponds
to the bound state energy. This solution exactly agrees with the original
result obtained by Shiba [2]. Thus, the particular state ${\bf k}_{0}$ where
the singularity appears has to be identified with the well-known YSR state
appearing in the spectral function as a pair of peaks inside the
superconducting gap.

In summary, the whole renormalization process occurs as follows. At the
starting point $\lambda=1$ of the renormalization process the particular state
${\bf k}_{0}$ where $E_{{\bf k}_{0}}=\Delta$ becomes a singular point for the
transformation coefficient $J_{{\bf k}_{0},\lambda}$ (and also $V_{{\bf
k}_{0},\lambda}$), i. e. it takes a value proportional to $N$. Thus, the
solution for the bound state is described by $\tilde{J}_{{\bf
k}_{0}}=N\tilde{J}$ and $\tilde{V}_{{\bf k}_{0}}=-S\tilde{J}_{{\bf k}_{0}}$
due to Eqs. (62) and (64). The corresponding equation for $\tilde{J}$ is found
from Eq. (31) by taking into account again Eq. (64),

$\tilde{J}=\frac{\Delta-\tilde{E}_{{\bf k}_{0}}}{2S}.$ (69)

Using this expression and Eq. (48) we obtain the explicit form of the
quasiparticle operator with respect to ${\bf k}\neq{\bf k}_{0}$ in the
classical spin case,

$\begin{split}\tilde{\alpha}_{{\bf k}}^{\dagger}&=\tilde{x}_{{\bf
k}}\alpha_{\bf k}^{\dagger}+\frac{\tilde{y}_{{\bf k}}}{N}\sum_{{\bf
k}^{\prime}(\neq{\bf k},{\bf k}_{0})}\left[C_{1,{\bf k}^{\prime}{\bf
k}}^{+}\frac{\left(\tilde{J}_{{\bf k}}+\tilde{J}_{{\bf
k}^{\prime}}\right)S-\left(\tilde{V}_{{\bf k}}+\tilde{V}_{{\bf
k}^{\prime}}\right)}{E_{{\bf k}^{\prime}}-E_{{\bf k}}}\alpha_{{\bf
k}^{\prime}}^{\dagger}-C_{2,{\bf k}{\bf
k}^{\prime}}^{-}\frac{\left(\tilde{J}_{{\bf k}}+\tilde{J}_{{\bf
k}^{\prime}}\right)S+\left(\tilde{V}_{{\bf k}}+\tilde{V}_{{\bf
k}^{\prime}}\right)}{E_{{\bf k}}+E_{{\bf k}^{\prime}}}\beta_{{\bf
k}^{\prime}}\right]\\\ &+\left(u_{\bf k}+v_{\bf
k}\right)\frac{\Delta-\tilde{E}_{{\bf k}_{0}}}{E_{{\bf k}}-\tilde{E}_{{\bf
k}_{0}}}\alpha_{{\bf k}_{0}}^{\dagger}.\end{split}$ (70)

It is seen that the quasiparticle in the classical spin case consists of three
basic contributions. The first part proportional to $\tilde{x}_{{\bf k}}$
describes the coherent excitation which follows the original dispersion of the
superconducting quasiparticle. The second part in the first line which is
proportional to $\tilde{y}_{{\bf k}}$ contains the scattering to all other
states except ${\bf k}$ itself and ${\bf k}_{0}$. Since in contrast to the
quantum spin case the renormalization of the exchange coupling,
$\tilde{J}_{{\bf k}}$, is always very small (absence of the Kondo resonance)
this term only produces a very small spectral weight ($\propto 1/N$) in the
single-particle spectral function. Finally, the third term in the second line
is the creation of the YSR bound state ${\bf k}_{0}$ appearing at an energy
$\tilde{E}_{{\bf k}_{0}}$ inside the superconducting gap. Thus, the expression
for the single-particle spectral function for $\mathbf{k}\neq\mathbf{k}_{0}$
given be Eq. (54) simplifies in the classical spin case to

$\begin{split}\Im G({\bf k},{\bf k},\omega)&=\frac{\tilde{x}_{\bf
k}^{2}}{2}\left[\left(1+\frac{\varepsilon_{\mathbf{k}}}{E_{\bf
k}}\right)\delta(E_{\bf
k}-\omega)+\left(1-\frac{\varepsilon_{\mathbf{k}}}{E_{\bf
k}}\right)\delta(-E_{\bf k}-\omega)\right]\\\
&+\frac{1}{2}\left(\frac{\Delta-\tilde{E}_{{\bf k}_{0}}}{E_{{\bf
k}}-\tilde{E}_{{\bf k}_{0}}}\right)^{2}\left[\delta(\tilde{E}_{{\bf
k}_{0}}-\omega)+\delta(-\tilde{E}_{{\bf k}_{0}}-\omega)\right].\end{split}$
(71)

In addition to the usual dispersive states (first line), there appears in the
second line a pair of dispersionless excitations which are attributed to the
YSR bound states. Those have the $\mathbf{k}$-independent energy
$\pm\tilde{E}_{{\bf k}_{0}}$ above and below the Fermi level, but inside the
superconducting gap. Since all excitations appear as $\delta$ peaks without
broadening the Kondo resonance arising from many-particle excitations is
absent here.

Figure 6:  (a) Single-particle Green’s function for a classical spin for ${\bf
k}\neq{\bf k}_{0}$ near the Fermi surface. The chosen parameters are
$J=\Delta=0.004W$, $\tilde{E}_{{\bf k}_{0}}=0.002W$, and $S=1/2$ ($W$: band
width). In this figure the level broadening was fixed manially to the value
$\eta=2\cdot 10^{-4}W$. (b) Energy-integrated spectral weight of the different
excitations in the spectral function of panel (a) as a function of
$\mathbf{k}$.

The behavior of the YSR state in the spectral function is visualized in Fig.
6(a) where the spectral function from Eq. (71) is plotted explicitly for a
specific set of parameters, and for a number of momenta
$\mathbf{k}\neq\mathbf{k}_{0}$. One finds the typical dispersionless double-
peak structure of the YSR state inside the superconducting gap where the
intensity is strongest in the vicinity of the gap edge. Furthermore, one
clearly recognizes the two distinct branches of the superconducting
quasiparticles below and above the Fermi level. The spectral weight of the
corresponding excitations which are shown in Fig. 6(b) reveals an intensity
maximum of the YSR state near the Fermi momentum of the normal state. The
$\mathbf{k}$ dependence of the coherent excitation intensity, however, is
mainly influenced by the Bugoliubov coefficients $|u_{\bf k}|^{2}$ and
$|v_{\bf k}|^{2}$.

##  Appendix B Non-magnetic impurity

Here we explain our approach to bound states by applying it to the perhaps
most simple case of a local potential impurity embedded into a spinless,
normal-state electron system. Concretely, we consider a toy-model Hamiltonian
$\mathcal{H}_{V}$ composed of an electron-part $\mathcal{H}_{0}$ and the
electron-impurity-part $\mathcal{H}_{1,V}$ proportional to the coupling
strength $V$,

$\begin{split}\mathcal{H}_{V}=\mathcal{H}_{0}+\mathcal{H}_{1,V}\quad\mbox{with}\quad\mathcal{H}_{0}=\sum_{k}\varepsilon_{k}c^{\dagger}_{k}c_{k}\quad\mbox{and}\quad\mathcal{H}_{1,V}=\frac{V}{2N}\sum_{k\neq
k^{\prime}}\left(c^{\dagger}_{k}c_{k^{\prime}}+c_{k^{\prime}}^{{\dagger}}c_{k}\right).\end{split}$
(72)

Here, $\varepsilon_{k}$ is the bare electronic dispersion, which is a function
of the one-dimensional momentum $k$, and the lattice consists of $N$ sites.

Following the idea of Sec. 2.1, we subject the Hamiltonian (72) to an exact
mapping that brings the Hamiltonian to a new form $\tilde{\mathcal{H}}_{V}$.
The mapping between the original Hamiltonian ${\mathcal{H}}_{V}$ and
$\tilde{\mathcal{H}}_{V}$ is implemented by a unitary transformation of the
form given in Eq. (8). The transformed Hamiltonian should not only be
quadratic in the original fermionic operators $c^{\dagger}_{k}$, but also have
the diagonal form

$\tilde{\mathcal{H}}_{V}=\sum_{k}\tilde{\varepsilon}_{k}c^{\dagger}_{k}c_{k},$
(73)

with a renormalized dispersion $\tilde{\varepsilon}_{k}$. This treatment is
technically slightly different (but equivalent in the way explained in the
main text) from the usual diagonalization in which the quasiparticle operators
in the Hamiltonian are rotated in the operator subspace. For technical
convenience, our approach keeps the single-particle operators in their
original basis while the Hamiltonian is renormalized. Therefore,
$\tilde{\mathcal{H}}_{V}$ differs in its form from ${\mathcal{H}}_{V}$, but
has the same eigenvalue spectrum as the original model.

As explained above in the context of the spin impurity the unitary
transformation to obtain $\tilde{\mathcal{H}}_{V}$ is not carried out in one
single step but within a process of steps which are described by a parameter
$\lambda$. The $\lambda$-dependent Hamiltonian for the present problem is
defined as follows

$\mathcal{H}_{V,\lambda}=\mathcal{H}_{0V,\lambda}+\mathcal{H}_{1V,\lambda}\quad\mbox{with}\quad\mathcal{H}_{0V,\lambda}=\sum_{k}\varepsilon_{k,\lambda}c^{\dagger}_{k}c_{k},\quad\mathcal{H}_{1V,\lambda}=\lambda\frac{V}{2N}\sum_{k\neq
k^{\prime}}\left(c^{\dagger}_{k}c_{k^{\prime}}+c_{k^{\prime}}^{{\dagger}}c_{k}\right).$
(74)

The unitary transformation combining two subsequent parameter values $\lambda$
and $\lambda-\Delta\lambda$ is constructed such that applying the
transformation as in Eq. (11) the structure of $\mathcal{H}_{V,\lambda}$ is
kept for any value of $\lambda$. A convenient ansatz for
$X_{V,\lambda,\Delta\lambda}$ is

$X_{V,\lambda,\Delta\lambda}=\frac{\Delta\lambda}{N}\sum_{k\neq
k^{\prime}}\frac{A_{k,\lambda}}{\varepsilon_{k,\lambda}-\varepsilon_{k^{\prime},\lambda}}\left(c^{\dagger}_{k}c_{k^{\prime}}-c_{k^{\prime}}^{{\dagger}}c_{k}\right).$
(75)

The real coefficients $A_{k,\lambda}$ will again turn out to be key in the
description of bound states. To determine their renormalization equations, we
plug the ansatz (75) into Eq. (11), which allows to determine equations for
$A_{k,\lambda}$ from an expansion as in Eq. (13) in terms of commutators up to
the lowest order in $\Delta\lambda$. Up to first order in the generator
$X_{V}$, the commutators of $X_{V}$ with ${\cal H}_{0}$ and ${\cal H}_{1,V}$
as defined in Eq. (72) read,

$\begin{split}[X_{V,\lambda,\Delta\lambda},{\cal
H}_{0V,\lambda}]&=-\frac{\Delta\lambda}{N}\sum_{k\neq
k^{\prime}}A_{k,\lambda}\left(c^{\dagger}_{k}c_{k^{\prime}}+c_{k^{\prime}}^{{\dagger}}c_{k}\right),\\\
[X_{V,\lambda,\Delta\lambda},{\cal
H}_{1V,\lambda}]&=\frac{V}{N^{2}}\lambda\Delta\lambda\sum_{kk^{\prime}}\sum_{q(\neq
k)}\frac{A_{k,\lambda}+A_{q,\lambda}}{\varepsilon_{k,\lambda}-\varepsilon_{q,\lambda}}\left(c^{\dagger}_{k}c_{k^{\prime}}+c_{k^{\prime}}^{{\dagger}}c_{k}\right).\end{split}$
(76)

Comparing again the both sides of Eq. (13) after plugging in the two above
commutators we find the following renormalization equations,

$\displaystyle-\frac{V}{2}$
$\displaystyle=-A_{k,\lambda}+\frac{V}{N}\lambda\sum_{q(\neq
k)}\frac{A_{k,\lambda}+A_{q,\lambda}}{\varepsilon_{k,\lambda}-\varepsilon_{q,\lambda}},$
(77) $\displaystyle\varepsilon_{k,\lambda-\Delta\lambda}$
$\displaystyle=\varepsilon_{k,\lambda}+\frac{2V}{N^{2}}\lambda\Delta\lambda\sum_{q(\neq
k)}\frac{A_{k,\lambda}+A_{q,\lambda}}{\varepsilon_{k,\lambda}-\varepsilon_{q,\lambda}}.$
(78)

Again $A_{k,\lambda}$ can become singular at a specific value of $k=k_{0}$,
which indicates the formation of a bound state. Furthermore, since
$A_{k_{0},\lambda}\propto N$ such a singularity entails a renormalization of
the single-particle energy from the bare value $\varepsilon_{k_{0}}$ to a new
value – but only at the $k_{0}$ for which $A_{k_{0},\lambda}$ is singular, i.
e. for all $q\neq k_{0}$ we can set $\varepsilon_{q,\lambda}=\varepsilon_{q}$.
The renormalization at $k_{0}$ will again be related to the binding energy of
the bound state. Thus, for the singular state $k_{0}$ we can solve Eq. (77)
for $A_{k_{0},\lambda}$,

$A_{k_{0},\lambda}=\frac{V/2}{1-\lambda\frac{V}{N}\sum_{q(\neq
k_{0})}\frac{1}{\varepsilon_{k_{0},\lambda}-\varepsilon_{q}}}.$ (79)

In the beginning of the renormalization process we can set $\lambda=1$. Thus,
a singularity is found if the equation

$1=\frac{V}{N}\sum_{q(\neq
k_{0})}\frac{1}{\varepsilon_{k_{0},\lambda=1}-\varepsilon_{q}}$ (80)

is fulfilled. The solution $\varepsilon_{k_{0},\lambda=1}$ of this equation is
the renormalized energy value due to the singularity. In general one finds two
solutions, where the first solution $\varepsilon_{k_{0}}$ in general is a part
of the quasiparticle band and the second solution lies below or above this
band depending on the sign of $V$. It follows that the first solution
$\varepsilon_{k_{0}}$ changes to the renormalized value
$\varepsilon_{k_{0},\lambda=1}$ and remains at this value down to $\lambda=0$,
i. e. $\varepsilon_{k_{0},\lambda=1}=\tilde{\varepsilon}_{k_{0}}$. All other
$q\neq k_{0}$ are unaffected and the corresponding energies remain unchanged,
$\varepsilon_{q,\lambda}=\varepsilon_{q}$.

The fully renormalized energy value $\tilde{\varepsilon}_{k_{0}}$ which is a
solution of Eq. (80) exactly agrees with the result of the well known t-matrix
approach to bound states [36]. The bound state energy is calculated by
evaluating the sum using the density of states of the free system and solving
the resulting equation for $\tilde{\varepsilon}_{k_{0}}$. Usually the energy
value $\tilde{\varepsilon}_{k_{0}}$ lies below the bottom of the conduction
electron band and it only exists if $V<0$, i. e. in the case of an attractive
scattering potential.

The whole renormalization process occurs similar to the spin impurity problem
as follows. At the starting point $\lambda=1$ of the renormalization process a
particular state $k_{0}$ from the quasiparticle band becomes a singular point
for the corresponding transformation coefficient $A_{k_{0},\lambda}$, i. e. it
takes a value proportional to $N$. Thus, the solution for the coefficient at
the bound state can be written in the form $\tilde{A}_{k_{0}}=N\tilde{A}$. The
corresponding equation for $\tilde{A}$ is found from Eq. (78),

$\tilde{A}=\frac{\tilde{\varepsilon}_{k_{0}}-\varepsilon_{k_{0}}}{2}.$ (81)

For all other $k\neq k_{0}$, the coefficients $A_{k,\lambda}$ are found by
self consistent solution of Eq. (77).

Let us again visualize the bound state by studying the spectral function,
which is defined via the imaginary part of the single-particle Green’s
function,

$\Im
G(k,k,\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\langle[c_{k}^{\dagger}(-t),c_{k}]_{+}\rangle
e^{i\omega t}dt.$ (82)

The time dependence and the expectation value are now formed with the
Hamiltonian ${\cal H}_{V}$. Correspondingly, the expectation value of an
operator $\mathcal{A}$ is defined as in Eq. (39) but now formed with
$\mathcal{H}_{V}$. Using the Mori-Zwanzig formalism [37, 38], the
anticommutator correlation function can be rewritten to the form

$\Im
G(k,k,\omega)=\langle[c_{k},\delta(\mathbf{L}_{V}-\omega)c_{k}^{\dagger}]_{+}\rangle,$
(83)

where $\mathbf{L}_{V}\mathcal{A}=[\mathcal{H}_{V},\mathcal{A}]$ is the
Liouville operator. Our method to calculate this quantity is again by using
the invariance of the trace under a unitary transformation. This yields

$\Im
G(k,k,\omega)=\langle[\tilde{c}_{k},\delta(\tilde{\mathbf{L}}_{V}-\omega)\tilde{c}_{k}^{\dagger}]_{+}\rangle.$
(84)

To find the transformed operators $\tilde{c}_{k}$, we use the same procedure
that we used above to determine the fully renormalized Hamiltonian
$\tilde{\cal H}_{V}$ from the original Hamiltonian. We start with an ansatz
for the renormalized single-particle operator where the transformation is
already carried out until an arbitrary $\lambda$. This ansatz can be suggested
by forming the first order commutator of the single-particle operator and
$X_{V,\lambda,\Delta\lambda}$ where we use the same operator form but with
$\lambda$-dependent prefactors,

$c_{k,\lambda}^{\dagger}=x_{k,\lambda}c_{k}^{\dagger}+y_{k,\lambda}\frac{1}{N}\sum_{q(\neq
k)}\frac{A_{k,\lambda}+A_{q,\lambda}}{\varepsilon_{q,\lambda}-\varepsilon_{k,\lambda}}c_{q}^{\dagger},$
(85)

where the initial conditions at $\lambda=1$ are $x_{k\lambda=1}=1$ and
$y_{k\lambda=1}=0$. Now we find the renormalization equations for the
$\lambda$-dependent coefficients $x_{k\lambda}$ and $y_{k\lambda}$ in Eq. (85)
by evaluating a similar equation as Eq. (11) but for the renormalized single
particle operator,

$c_{k,\lambda-\Delta\lambda}^{\dagger}=e^{X_{V,\lambda,\Delta\lambda}}c_{k,\lambda}^{\dagger}e^{-X_{V,\lambda,\Delta\lambda}},$
(86)

where the operator expression as given by Eq. (75) must be used for
$X_{V,\lambda,\Delta\lambda}$. Evaluating again the right hand side up to the
linear order in $\Delta\lambda$ and using the above ansatz for
$c_{k,\lambda}^{\dagger}$ we find the following set of renormalization
equations for the $\lambda$-dependent coefficients,

$\begin{split}x_{k,\lambda-\Delta\lambda}&=x_{k,\lambda}-y_{k,\lambda}\frac{\Delta\lambda}{N^{2}}\sum_{q(\neq
k)}\left(\frac{A_{k,\lambda}+A_{q,\lambda}}{\varepsilon_{q,\lambda}-\varepsilon_{k,\lambda}}\right)^{2},\\\
1&=x_{k,\lambda-\Delta\lambda}^{2}+y_{k,\lambda-\Delta\lambda}^{2}\frac{1}{N^{2}}\sum_{q(\neq
k)}\left(\frac{A_{k,\lambda-\Delta\lambda}+A_{q,\lambda-\Delta\lambda}}{\varepsilon_{q,\lambda-\Delta\lambda}-\varepsilon_{k,\lambda-\Delta\lambda}}\right)^{2}.\end{split}$
(87)

This set of renormalization equations is solved simultaneously with the
corresponding equations for the Hamiltonian until $\lambda=0$ is obtained.
Then the Hamiltonian is diagonal and the fully renormalized single-particle
operator at $\lambda=0$ is obtained from Eq. (85),

$\tilde{c}_{k}^{\dagger}=\tilde{x}_{k}c_{k}^{\dagger}+\tilde{y}_{k}\frac{1}{N}\sum_{q(\neq
k)}\frac{\tilde{A}_{k}+\tilde{A}_{q}}{\tilde{\varepsilon}_{q}-\tilde{\varepsilon}_{k}}c_{q}^{\dagger}.$
(88)

Note that a possible bound state $k_{0}$ may be also included in the
summation. If we consider a state $k\neq k_{0}$ we find

$\tilde{c}_{k}^{\dagger}=\tilde{x}_{k}c_{k}^{\dagger}+\tilde{y}_{k}\frac{1}{N}\sum_{q(\neq
k,k_{0})}\frac{\tilde{A}_{k}+\tilde{A}_{q}}{\varepsilon_{q}-\varepsilon_{k}}c_{q}^{\dagger}+\frac{1}{2}\frac{\tilde{\varepsilon}_{k_{0}}-\varepsilon_{k_{0}}}{\tilde{\varepsilon}_{k_{0}}-\varepsilon_{k}}c_{k_{0}}^{\dagger}.$
(89)

Figure 7:  (a) single-particle Green’s function for $k\neq k_{0}$ of a non-
magnetic impurity. The level broadening is introduced by hand and the
broadening value is fixed to $\eta=0.1\tilde{\varepsilon}_{k_{0}}$. (b)
Energy-integrated spectral weights of the spectral function according to Eq.
(90). The black line shows the coherent excitation and the red line shows the
bound state excitation as a function $k$.

Inserting Eq. (89) into the correlation function (84), we can directly
evaluate the operator $\delta(\tilde{\mathbf{L}}_{V}-\omega)$ by simply
replacing the transformed Liouville operator $\tilde{\mathbf{L}}_{V}$ with an
eigenvalue of the corresponding single-particle operator it acts on. This is
possible because it refers to the diagonalized Hamiltonian
$\tilde{\mathcal{H}}_{V}$. Thus, we obtain for $k\neq k_{0}$

$\Im
G(k,k,\omega)=\tilde{x}_{k}^{2}\delta(\varepsilon_{k}-\omega)+\frac{\tilde{y}_{k}^{2}}{N^{2}}\sum_{q(\neq
k,k_{0})}\left(\frac{\tilde{A}_{k}+\tilde{A}_{q}}{\varepsilon_{q}-\varepsilon_{k}}\right)^{2}\delta(\varepsilon_{q}-\omega)+\left(\frac{1}{2}\frac{\tilde{\varepsilon}_{k_{0}}-\varepsilon_{k_{0}}}{\tilde{\varepsilon}_{k_{0}}-\varepsilon_{k}}\right)^{2}\delta(\tilde{\varepsilon}_{k_{0}}-\omega).$
(90)

As in the case of a spin impurity we find the spectral function to exhibit a
state at energy $\tilde{\varepsilon}_{k_{0}}$ in addition to the bare energies
$\varepsilon_{k}$ of the free system. The energy $\tilde{\varepsilon}_{k_{0}}$
of this state is $k$-independent, indicating the state to be local in space –
the hallmark of a bound state. Note that the intensity of the bound state (and
also the other excitations) depends on $k$. Note that the continuum of states
in the second term in Eq. (90) vanishes in the thermodynamic limit due to the
factor $1/N^{2}$. Thus, all excitations exhibit $\delta$-like behavior in the
spectral function.

These properties are shown in Fig. 7(a) where the spectral function was
calculated using Eq. (90) for a quadratic dispersion $\varepsilon_{k}$.
Furthermore, the spectral weights of the bound state and the dispersive states
as calculated from the prefactors in Eq. (90) is shown in Fig. 7(b). We find
that the spectral weight of the bound state increases (at the cost of the
dispersive states) if it comes closer to the conduction band.

## References

  * [1] L. Yu, Acta Physica Sinica 21, 75 (1965).
  * [2] H. Shiba, Classical Spins in Superconductors, Prog. Theoret. Phys. 40, 435 (1968), https://doi.org/10.1143/PTP.40.435.
  * [3] A. I. Rusinov, Superconductivity near a Paramagnetic Impurity JETP Lett. 9, 85 (1969).
  * [4] J. Kondo, Resistance Minimum in Dilute Magnetic Alloys Prog. Theor. Phys. 32, 37 (1964), https://doi.org/10.1143/PTP.32.37.
  * [5] R. Wang, W. Su, J.-X. Zhu, C. S. Ting, H. Li, C. Chen, B. Wang, and X. Wang, Kondo Signatures of a Quantum Magnetic Impurity in Topological Superconductors Phys. Rev. Lett. 122, 087001 (2019), https://doi.org/10.1103/PhysRevLett.122.087001.
  * [6] C. Liu, Y. Huang, Y. Chen, C. S. Ting, Temperature-dependent spectral function of a Kondo impurity in an s-wave superconductor, Phys, Rev. B 99, 245139 (2018), 10.1103/PhysRevB.99.174502.
  * [7] M. I. Salkola, A. V. Balatsky, and J. R. Schrieffer, Spectral properties of quasiparticle excitations induced by magnetic moments in superconductors Phys. Rev. B 55, 12648 (1997), https://doi.org/10.1103/PhysRevB.55.12648.
  * [8] R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, and S. Tarucha, Tunneling Spectroscopy of Andreev Energy Levels in a Quantum Dot Coupled to a Superconductor, Phys. Rev. Lett. 104, 076805 (2010), https://doi.org/10.1103/PhysRevLett.104.076805.
  * [9] K. J. Franke, G. Schulze, J. I. Pasual, Competition of Superconducting Phenomena and Kondo Screening at the Nanoscale, Science 332, 6032 (2011), 10.1126/science.1202204.
  * [10] B. W. Heinrich, J. I. Pascual, K. J. Franke, Single magnetic adsorbates on s-wave superconductors, Progress in Surface Science 93, 1 (2018), https://doi.org/10.1016/j.progsurf.2018.01.001.
  * [11] A. Yazdani, B. A. Jones, C. P. Lutz, M. F. Crommie, and D. M. Eigler, Probing the Local Effects of Magnetic Impurities on Superconductivity, Science 275, 1767 (1997), 10.1126/science.275.5307.1767.
  * [12] S.-H. Ji, et.al., High-Resolution Scanning Tunneling Spectroscopy of Magnetic Impurity Induced Bound States in the Superconducting Gap of Pb Thin Films, Phys. Rev. Lett. 100, 226801 (2008), https://doi.org/10.1103/PhysRevLett.100.226801.
  * [13] M. Flatte and D. Reynolds, Local spectrum of a superconductor as a probe of interactions between magnetic impurities, Phys. Rev. B 61, 14810 (2000), https://doi.org/10.1103/PhysRevB.61.14810.
  * [14] D. Morr and N. Stavropoulos, Quantum interference between impurities: Creating novel many-body states in s-wave superconductors, Phys. Rev. B 67, 020502(R) (2003), https://doi.org/10.1103/PhysRevB.67.020502.
  * [15] T. Meng, J. Klinovaja, S. Hoffman, P. Simon, and D. Loss, Superconducting gap renormalization around two magnetic impurities: From Shiba to Andreev bound states, Phys. Rev. B 92, 125422 (2015), https://doi.org/10.1103/PhysRevB.92.064503.
  * [16] S. Hoffman, J. Klinovaja, T. Meng, and D. Loss, Impurity-induced quantum phase transitions and magnetic order in conventional superconductors: Competition between bound and quasiparticle states, Phys. Rev. B 92, 064503 (2015), https://doi.org/10.1103/PhysRevB.92.125422.
  * [17] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, and A. Yazdani, Proposal for realizing Majorana fermions in chains of magnetic atoms on a superconductor, Phys. Rev. B 88, 020407(R) (2013), https://doi.org/10.1103/PhysRevB.88.020407.
  * [18] F. Pientka, L. I. Glazman, and F. V. Oppen, Topological superconducting phase in helical Shiba chains Phys. Rev. B 88, 155420 (2013), https://doi.org/10.1103/PhysRevB.88.155420.
  * [19] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor, Science 346, 602 (2014), 10.1126/science.1259327.
  * [20] J. Li, H. Chen, I. K. Drozdov, A. Yazdani, B. A. Bernevig, and A. H. MacDonald, Topological superconductivity induced by ferromagnetic metal chains, Phys. Rev. B 90, 235433 (2014), https://doi.org/10.1103/PhysRevB.90.235433.
  * [21] A. Martin-Rodero and A. L. Yeyati, The Andreev states of a superconducting quantum dot: mean field vs exact numerical results, Journal of Physics: Condensed Matter 24, 385303 (2012), 10.1088/0953-8984/24/38/385303.
  * [22] M. Zonda, V. Pokorny, V. Janis, and T. Novotny, Perturbation theory of a superconducting $0-\pi$ impurity quantum phase transition, Scientific Reports 5, 8821 (2015), https://www.nature.com/articles/srep08821.pdf?origin=ppub.
  * [23] M. Zonda, V. Pokorny, V. Janis, and T. Novotny, Perturbation theory for an Anderson quantum dot asymmetrically attached to two superconducting leads, Phys. Rev. B 93, 024523 (2016), https://doi.org/10.1103/PhysRevB.93.024523.
  * [24] T. Meng, S. Florens, and P. Simon, Self-consistent description of Andreev bound states in Josephson quantum dot devices, Phys. Rev. B 79, 224521 (2009), https://doi.org/10.1103/PhysRevB.79.224521.
  * [25] N. Wentzell, S. Florenz, T. Meng, V. Meden, and S. Andergassen, Magnetoelectric spectroscopy of Andreev bound states in Josephson quantum dots, Phys. Rev. B 94, 085151 (2016), https://doi.org/10.1103/PhysRevB.94.085151.
  * [26] R. Maurand, T. Meng, E. Bonet, S. Florens, L. Marty, and W. Wernsdorfer, First-Order $0-\pi$ Quantum Phase Transition in the Kondo Regime of a Superconducting Carbon-Nanotube Quantum Dot, Phys. Rev. X 2 , 011009 (2012).
  * [27] H. Shiba, K. Satori, O. Sakai, and Y. Shimizu, Numerical renormalization group study of the Kondo effect in superconductors Physica B: Condensed Matter 186-188, 239 (1993), https://doi.org/10.1016/0921-4526(93)90540-M.
  * [28] O. Sakai, Y. Shimizu, H. Shiba, K. Satori, Numerical Renormalization Group Study of Magnetic Impurities in Superconductors. II. Dynamical Excitation Spectra and Spatial Variation of the Order Parameter, J. Phys. Soc. Jpn. 62, 3181 (1993), https://doi.org/10.1143/JPSJ.62.3181.
  * [29] J. Bauer, A. Oguri, and A. C. Hewson, Spectral properties of locally correlated electrons in a Bardeen–Cooper–Schrieffer superconductor Journal of Physics: Condensed Matter 19, 486211 (2007), https://iopscience.iop.org/article/10.1088/0953-8984/19/48/486211/meta.
  * [30] S. D. Głazek and K. G. Wilson, Renormalization of Hamiltonians Phys. Rev. D 48, 5863 (1993), https://doi.org/10.1103/PhysRevD.48.5863.
  * [31] F. J. Wegner, Flow-equations for Hamiltonians Ann. Phys. (Leipzig) 506, 77 (1994), https://doi.org/10.1002/andp.19945060203.
  * [32] S. Kehrein, The Flow Equation Approach to Many-Particle Systems, Springer Tracts in Modern Physics (Springer Verlag GmbH, Berlin, Heidelberg, 2006).
  * [33] S. Sykora, K. W. Becker, Many-body approach to Luttinger’s theorem for the Kondo lattice, Phys, Rev. B 98, 174502 (2019), https://doi.org/10.1103/PhysRevB.98.245139.
  * [34] N. Hatter, B. W. Heinrich, D. Rolf, K. J. Franke, Scaling of Yu-Shiba-Rusinov energies in the weak-coupling Kondo regime, Nature Communications 8, 2016 (2017), https://doi.org/10.1038/s41467-017-02277-7.
  * [35] S. Sykora, A. Hübsch, K. W. Becker, Generalized diagonalization scheme for many-particle systems, Phys. Rev. B 102, 165122 (2020), https://doi.org/10.1103/PhysRevB.102.165122.
  * [36] A. V. Balatsky, I. Vekhter, and Jian-Xin Zhu, Impurity-induced states in conventional and unconventional superconductors, Rev. Mod. Phys. 78, 373 (2006), https://doi.org/10.1103/RevModPhys.78.373.
  * [37] R. Zwanzig, Nonequilibrium Statistical Mechanics 3rd ed., Oxford University Press, New York (2001).
  * [38] H. Grabert, Projection operator techniques in nonequilibrium statistical mechanics, Springer Tracts in Modern Physics 95 (1982).
  * [39] M. Maltseva, M. Dzero, P. Coleman, Electron Cotunneling into a Kondo Lattice, Phys. Rev. Lett. 103, 206402 (2009), https://doi.org/10.1103/PhysRevLett.103.206402