# Yu-Shiba-Rusinov states, the BCS-BEC crossover, and the exact solution in
the flat-band limit

R. Žitko  Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19,
SI-1000 Ljubljana, Slovenia  L. Pavešić  Jožef Stefan Institute, Jamova 39,
SI-1000 Ljubljana, Slovenia  Faculty of Mathematics and Physics, University of
Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia

###### Abstract

We study the subgap Yu-Shiba-Rusinov (YSR) states in the Richardson’s model of
a superconductor with a magnetic impurity for different electron pairing
strengths from the weak-coupling Bardeen-Cooper-Schrieffer (BCS) regime to the
strong-coupling Bose-Einstein-condensate (BEC) regime. We observe that the
effect of the increasing pairing strength on the YSR excitation spectrum as a
function of hybridisation strength and impurity on-site potential is only
quantitative and, in fact, rather weak when the results are appropriately
rescaled. We furthermore show that the problem is analytically solvable in the
deep BEC limit which is equivalent to flat-band superconductivity. This exact
solution can be related to a zero-bandwidth (ZBW) effective BCS mean-field
Hamiltonian where the superconductor is described by a single electron level
with onsite pairing. The small difference between the BCS and BEC regimes of
the Richardson’s model explains the success of the simple ZBW calculations for
BCS mean-field Hamiltonians. A ZBW model requires only a suitable parameter
rescaling to become useful as a quantitative predictive tool for the full
problem.

##  I Introduction

The pairing interaction is a key concept in nuclear and solid-state physics,
as well as in quantum many-body theory in general. Pairing between fermions
has been studied in bulk materials [1, 2, 3], thin films [4], layered
materials [5, 6, 7], atomic nuclei [8, 9, 10, 11, 12], nuclear matter in
neutron stars [13], cold atom systems [14, 15, 16, 17, 18, 19], nanoscopic
metal grains [20, 21, 22, 23, 24], and ultrasmall superconducting islands [25,
26]. When the pairing is weak, the BCS mean-field approach works well [1].
When the pairing is strong, all fermionic particles pair up into bosons, then
the bosons condense: this produces a BEC. Both regimes are smoothly connected,
despite the fact that the two limits are physically quite different [14, 19,
27].

A magnetic impurity, such as a magnetic dopant or a semiconductor quantum dot
(QD), is well known to induce subgap states in the superconducting gap of a
BCS superconductor [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39], which are
known as the Yu-Shiba-Rusinov (YSR) states in the limit of well-defined local
moment. This raises the questions about the persistence of the YSR states all
the way to the BEC limit and about the evolution of their properties in this
cross-over. These questions are pertinant because the ground state and the
elementary excitations in the two limits have different properties [19].

In this work, we address this subject by coupling a magnetic impurity to a
superconductor (SC) described by the Richardson’s ”picket-fence” pairing
model, and tracking the subgap states through the crossover from weak to
strong pairing. We show that the YSR states not only persist in the cross-over
from the BCS to the BEC regime, but do not even change significantly at the
quantitative level after an appropriate rescaling of the input parameters and
the resulting energies. We also establish that a simple model based on the
flat-band limit of the Richardson’s model (RM) reproduces the full
phenomenology of an interacting QD coupled to a SC island.

The RM first appeared as a model of pairing forces in nuclear physics [40, 11,
41] and was later reintroduced as a description of nanoscale SC grains [42,
24, 43, 44]. The SC is modelled as a set of equidistant energy levels
representing the time-reversal conjugate single-particle states with all-to-
all pairing interaction. While the RM reduces to a BCS system in the
thermodynamic limit, it has important advantages for smaller systems. Namely,
it does not use the mean-field approximation and thus conserves particle
number. This makes it a more suitable choice for modelling of mesoscopic SC
systems. Unlike BCS, it is also applicable for all coupling strengths [45]
which is critical in the present work.

The RM is integrable and analytically solvable in terms of hard-core
bosons[11, 41], but coupling it with an impurity level introduces pair-
breaking processes and breaks integrability. We have previously found a
representation of the QD-SC model as a matrix product operator (MPO), allowing
us to use density matrix renormalization group (DMRG) [46, 47, 48] to obtain
exact solutions for a QD coupled to one [49] or two [50] SC channels. With
this method we obtain exact properties of a magnetic impurity coupled to an
interacting SC with no parameter restrictions, from the BCS to the BEC limits,
for all system sizes including ultrasmall SC islands. Good agreement to
experiment confirms that this is a suitable description of real systems [51].

This article is organized as follows. In Sec. II we present the Hamiltonian
and the methods that permit its numerical solution. In Sec. III we present the
results for the evolution from the BCS to the BEC limit. In Sec. IV we present
an exact solution of the model in the limit of a completely flat band
(equivalent to the deep BEC regime). In Sec. V we discuss how this solution
can also be considered as arising from a ZBW BCS model which leads to exactly
the same matrix representation. We conclude in Sec. VI with a perspective on
the possible extensions of the flat-band model, e.g. to a two-channel
(Josephson junction) situation. Some lengthy derivations are presented in
Appendices A-C. In the Supplemental material we provide Mathematica notebooks
with a computer-algebra verification of all mathematical statements presented
in this work [52]. The full source code of the DMRG solver is available on a
public repository, including a large set of examples [53]. We also provide the
input files for the solver and a corresponding Mathematica notebook with an
exact calculation of the ground state energies in the flat-band limit [52];
this reference calculation is discussed in Appendix D.

##  II Model and method

The Richardson’s model is [24, 49]

$H=\sum_{i\sigma}\epsilon_{i}c^{\dagger}_{i\sigma}c_{i\sigma}-G\sum_{i,j}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{j\downarrow}c_{j\uparrow},$
(1)

where $c^{\dagger}_{i\sigma}$ creates a particle in the level $i$ with spin
$\sigma\in\\{\uparrow,\downarrow\\}$ and $\epsilon_{i}$ are the energy levels
spanning the interval $[-D:D]$ spaced by $d=2D/\mathcal{L}$, where the half-
bandwidth $D\equiv 1$ sets the energy scale, and $\mathcal{L}$ is the total
number of levels. The indexes $i$ and $j$ range from 1 to $\mathcal{L}$. The
coupling $G$ can also be written as

$G=\alpha d=D\frac{2\alpha}{\mathcal{L}}=\frac{g}{\mathcal{L}},$ (2)

where we have introduced $\alpha$ as the dimensionless strength of the pairing
interaction, and $g$ as the corresponding dimensionfull strength; $\alpha$ and
$g$ are intensive quantities.

The BCS to BEC limiting process can be implemented in two equivalent ways: by
increasing the coupling $G$ at constant bandwidth $D$, or by decreasing the
bandwidth $D$ at constant coupling $G$. We opt for the second possibility in
this work. We implement this by appropriately rescaling the coefficients
$\epsilon_{i}$ in the kinetic-energy term of the SC Hamiltonian. Going to the
BEC limit thus corresponds to taking the flat-band limit, i.e., omitting the
kinetic-energy terms from the Hamiltonian.

We describe the coupling of the interacting QD to the SC using an Anderson
impurity Hamiltonian [54, 55, 56, 49]:

$\begin{split}H_{\mathrm{QD}}&=\epsilon\hat{n}+U\hat{n}_{\uparrow}\hat{n}_{\downarrow}+v\frac{1}{\sqrt{\mathcal{L}}}\sum_{i,\sigma}\left(c^{\dagger}_{i,\sigma}d_{\sigma}+d^{\dagger}_{\sigma}c_{i,\sigma}\right)\\\
&=\epsilon\hat{n}+U\hat{n}_{\uparrow}\hat{n}_{\downarrow}+v\sum_{i,\sigma}\left(f^{\dagger}_{\sigma}d_{\sigma}+d^{\dagger}_{\sigma}f_{\sigma}\right),\end{split}$
(3)

where we have defined the bath orbital $f$ through

$f_{\sigma}=\frac{1}{\sqrt{\mathcal{L}}}\sum_{i=1}^{\mathcal{L}}c_{i,\sigma}.$
(4)

The QD filling is controlled through the occupancy parameter $\nu$ defined via

$\epsilon=U(1/2-\nu).$

The particle-hole symmetric point of the QD corresponds to $\epsilon=-U/2$,
i.e., $\nu=1$. The hybridisation is quantified using $\Gamma=\pi\rho v^{2}$
with the density of states

$\rho=1/2D=1/d\mathcal{L}.$

We use the DMRG to find the ground state and a few lowest excited states in
symmetry sectors with conserved particle number $n$ and z-component of spin
$s_{z}$. We choose $n$ close to half filling with $s_{z}=0$ for even $n$ and
$s_{z}=1/2$ for odd $n$.

The Hamiltonian is written in MPO form with $9\times 9$ matrices [49]. The
DMRG results are obtained with system size $\mathcal{L}=100$, with maximal MPS
matrix dimension of $2000$ and the energy convergence criterion set at
$\delta=10^{-8}$. System sizes up to $\mathcal{L}=2000$ with $\delta=10^{-12}$
are realistically achievable, but increasing the system size produces marginal
gains. Using $\mathcal{L}\gtrsim 30$ already gives good results at much
smaller computational expense. We set the pairing strength $\alpha=0.4$ to
ensure that an appropriate amount of energy levels participate in pairing even
when they span the entire energy range $[-1:1]$.

The SC-QD models can also be reliably solved by the numerical renormalization
group [57, 58, 59, 60, 61, 62, 38, 63, 64, 65] as well as with continuous-time
Quantum Monte Carlo (QMC) using $U$ [66] and hybridization [67, 68, 69, 70]
expansion. These methods require a weakly interacting bath and are thus
limited to the mean-field BCS description of superconductivity. The RM without
the impurity has been treated by QMC in the canonical ensemble [71], and its
strong coupling limit has been investigated by expanding the Richardson’s
solution in powers of $1/\alpha$ using pseudospin operators [72]. Our DMRG
approach is currently the only reliable technique to address the full
interacting problem in all parameter regimes.

##  III Numerical results for the BCS-BEC cross-over

We introduce $\gamma$ as a parameter that multiplies the kinetic energy term
of the Hamiltonian. The flattening of the band (and the BCS-BEC crossover)
thus correspond to a variation from $\gamma=1$ to $\gamma=0$. In this section
we present the evolution of the subgap spectra in this crossover.

###  III.1 Excitation energies

Fig. 1 shows the excitation energy of the subgap state at the particle-hole
symmetric point ($\epsilon=-U/2$) as a function of the hybridisation strength
$\Gamma$. We consider two cases, $U/\Delta_{0}=40$ and $U/\Delta_{0}=5$, where
$\Delta_{0}$ is the SC gap in the $\gamma=1$ limit. The first is
representative of the large-$U$ limit where the impurity is a pure exchange
scatterer (Kondo limit) and the subgap states are well-defined YSR states,
while the second is a generic situation where the singlet subgap states have a
mixed character due to a stronger proximity effect which admixes states of
zero and double impurity occupancy. The first row shows the results without
any rescaling. Withing $\Gamma$ increasing from zero, we see a singlet state
detaching from the continuum of Bogoliubov states: this is the YSR singlet. It
crosses the zero energy-difference line at $\Gamma=\Gamma_{c}$ at which point
it becomes the ground state of the system.

We now rescale the results in two steps: 1) we rescale the energies in terms
of the SC gap at given $\gamma$ (second row), 2) we furthermore rescale the
$\Gamma/U$ axis so that the singlet-doublet transition points for all $\gamma$
coincide (third row). These fully rescaled results demonstrate a reasonably
good overlap, showing that during the BCS to BEC crossover the changes in the
YSR excitation energy are only quantitative and actually rather small. This is
one of the key results of this work.

Figure 1: Singlet-doublet energy difference $E_{\mathrm{YSR}}=E_{S}-E_{D}$ as
a function of the hybridisation $\Gamma$ for a range of bandwidths $\gamma$.
Energy difference $E_{S}-E_{D}$ is expressed either in (a) absolute units or
(b,c) rescaled by the SC gap $\Delta$. In (c), the horizontal axis is
furthermore rescaled so that the singlet-doublet transition points coincide.
The parameters are $\alpha=0.4$, strong electron repulsion $U/\Delta_{0}=40$
(left column) and moderate repulsion $U/\Delta_{0}=5$ (right column), where
$\Delta_{0}\approx 0.16$ is the gap at $\gamma=1$, and particle-hole symmetric
tuning $\nu=1$. The results have been corrected for the finite-size effects as
described in Ref. 49.

Fig. 2(a) presents the transition point $(\Gamma/U)_{c}$ extracted from Fig.
1. The hybridization rescaling factor used in the bottom row of Fig. 1 is
shown in the inset. By decreasing the bandwith $\gamma$ we allow more levels
to participate in hopping processes resulting in larger effect of
hybridization $\Gamma$. The transition point thus moves to smaller value of
$\Gamma$ with decreasing bandwidth $\gamma$.

Figure 2:  (a) Bandwidth $\gamma$ dependence of the transition point
$(\Gamma/U)_{c}$. The inset shows the rescaling factor for the hybridisation
strength that would be required to obtain equivalent results. (b) Pairing
strength $\bar{\Delta}$ dependence on $\gamma$. It measures the proportion of
levels which contribute in pairing, with all levels contributing equally when
$\bar{\Delta}=1/2$.

In Fig. 2(b) we plot the $\gamma$ dependence of the pairing correlations
$\sum_{i}u_{i}v_{i}$, where $u_{i}=\langle
c_{i\downarrow}c_{i\uparrow}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}\rangle^{1/2}$
and $v_{i}=\langle
c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{i\downarrow}c_{i\uparrow}\rangle^{1/2}$,
for a range of pairing interaction strengths $\alpha$. The sum measures the
proportion of levels that contribute to pairing, with all levels contributing
equally when $\sum_{i}u_{i}v_{i}=1/2$. Similarly as in Fig. 2(a), the decrease
in bandwidth allows more levels to contribute in pairing processes,
strengthening the correlations. Increasing the pairing strength $\alpha$ at
given bandwidth results in qualitatively the same crossover as decreasing the
bandwidth $\gamma$, illustrating the two possible ways of probing the BEC -
BCS crossover. The pairing correlations in the zero bandwidth limit are
further discussed in section IV.1.

In Fig. 3 we show the gate-voltage dispersion for two fixed values of
$\Gamma/U$, on either side of the singlet-doublet transition value
$\Gamma_{c}$. The energies are scaled by the gap. We find a good overlap,
which could be further improved by adjusting the $\Gamma$ values as described
above.

Figure 3: Gate-voltage dispersion of the subgap states for two values of
hybridisation: (a) $\Gamma/U=0.05$, (b) $\Gamma/U=0.2$. Other parameters are
$\alpha=0.4$, $U/\Delta_{0}=40$. The frequencies are rescaled by $\Delta$.

###  III.2 Subgap state properties

Figure 4:  Impurity occupation probabilities $P_{1}$ and $P_{02}$ versus
binding energy $E_{B}=E_{S}-(E_{D}+\Delta)$ for different values of bandwidth
$\gamma$ for the singlet (top) and doublet (bottom) ground state. $E_{B}$ is
the energy gain of the subgap singlet state compared to the first excited
(decoupled) singlet state and is rescaled by the optical gap $\Delta$ for each
$\gamma$. The pairing interaction strength is $\alpha=0.4$.

In addition to the approximate overlap of the YSR state energies, we also find
a high degree of similarity in the state properties across the crossover. As
an example, in Fig. 4 we plot the probabilities that the impurity level is
singly occupied ($P_{1}$) or doubly occupied/empty ($P_{02}$) versus binding
energy $E_{B}=E_{S}-(E_{D}+\Delta)$. These are obtained from the reduced
density matrix for the impurity site and contain important information about
the nature of the states.

In the top row we plot the singlet ground state. For the YSR state in a
superconductor we expect $P_{\uparrow}=P_{\downarrow}=1/2$ and thus $P_{1}=1$.
With increasing $E_{B}$ (by increasing the hybridization $\Gamma$) the local
magnetic moment of the impurity level decreases due to larger charge
fluctuations. $P_{02}$ measures the contribution of the proximitized ABS
state, which becomes dominant in the limit of large $\Gamma/U$. Fig. 4(a)
shows that this behaviour persists throughout the BCS-BEC transition. In this
case with $U\gg\Delta$ we find good overlap of the curves. The situation away
from this limit, where $U\approx\Delta$, is shown in Fig. 4(b). Here the
magnitude of charge fluctuations depend not only on $\Gamma/U$ but also on
$\Gamma/\Delta$. The fact that two different scales control the charge
fluctuations decreases the quality of the collapse, but the qualitative
agreement for different $\gamma$ remains.

The bottom row shows the doublet ground state, where the mechanism of changing
state character is a bit different. Increasing $E_{B}$ increases the admixture
of the excited doublet state, which is identical to the YSR singlet with an
additional free quasiparticle in the bath. This excited state transforms from
a YSR-like singlet to an ABS-like proximitized state (as discussed in the
previous paragraph), resulting in the decrease of $P_{1}$, which is, however,
more gradual compared to the singlets. This is most obvious when comparing the
two cases with smaller $U$ in Fig. 4(b) and Fig. 4(d).

We note that we find a similar qualitative agreement for different $\gamma$
for many state properties. This is an important result in its own, as it
implies that it is possible to use the flat-band solutions to investigate the
nature of the low energy states in realistic QD–SC systems.

##  IV Exact solution for the flat-band limit

In the flat-band limit ($\gamma=0$, i.e., $D=0$ or $\epsilon_{i}\equiv 0$),
two renormalisation processes associated with the kinetic energy of band
electrons no longer occur. These are the renormalization of the exchange
interaction from the bare Kondo exchange coupling $J_{K}=8\Gamma/\pi\rho U$ to
the scale of the exponentially lower Kondo temperature $T_{K}=\exp(-1/\rho
J_{K})$, and the renormalization of the pairing interaction from the bare
scale of attractive coupling $g$ to the scale of the exponentially lower SC
gap $\Delta=\exp(-1/\rho g)$. Both renormalization processes have exactly the
same origin: the kinetic energy (dispersion) of the conduction band electrons.
In the following we show that the flat-band limit is analytically solvable, in
the sense of being reducible to numerical diagonalisations of very small
matrices. This model provides a conceptually clear picture of all phenomena
occurring in ultrasmall SCs and provides exact benchmark results for testing
other approaches.

###  IV.1 Eigenstates of the superconductor

We take the $\epsilon_{i}\equiv 0$ limit of the RM and write the SC part of
the Hamiltonian as

$H=-\frac{g}{\mathcal{L}}\sum_{i,j}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{j\downarrow}c_{j\uparrow}=-\frac{g}{\mathcal{L}}\sum_{i,j}A^{\dagger}_{i}A_{j},$
(5)

where $A_{j}=c_{j\downarrow}c_{j\uparrow}$ are the hard-core boson operators
with the commutation rule:

$[A_{i},A_{j}^{\dagger}]=\delta_{ij}(1-\hat{n}_{i}),$ (6)

where $\hat{n}_{i}=\sum_{\sigma}c^{\dagger}_{i\sigma}c_{i\sigma}$ is the level
occupancy operator. The following relations also hold:

$\begin{split}A_{i}^{2}&\equiv 0,\\\ (A_{i}^{\dagger})^{2}&\equiv 0,\\\
A_{i}^{\dagger}A_{i}&=\hat{P}_{2,i}=\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},\\\
A_{i}A_{i}^{\dagger}&=\hat{P}_{0,i}=1-\hat{n}_{i}+\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}.\end{split}$
(7)

The operators $\hat{P}_{0,i}$ and $\hat{P}_{2,i}$ are the projectors on the
subspace with occupancy 0 or 2 on level $i$.

The Hamiltonian $H$ is not particle-hole (p-h) symmetric. The symmetry can be
restored by adding a potential energy term:

$\tilde{H}=H+\tilde{\epsilon}\sum_{i}\hat{n}_{i},$ (8)

with

$\tilde{\epsilon}=\frac{g}{2\mathcal{L}}=\frac{\alpha
d}{2}=\frac{\alpha}{\mathcal{L}},$

using $D\equiv 1$ as energy unit.

Let $\mathcal{U}$ be the set of “unblocked levels” (those occupied by either 0
or 2 electrons). The blocked levels $\mathcal{B}$ (those with occupancy 1) do
not participate in pairing and they fully decouple from the problem in the
absence of an impurity, thus in this subsection all the sums that follow are
constrained to the set of unblocked levels. We will also use the same symbol
$\mathcal{U}$ to denote the number of unblocked levels, since this leads to no
ambiguity. In the absence of quasiparticles no level is blocked,
$\mathcal{B}=\emptyset$, and we have $\mathcal{U}=\mathcal{L}$.

We introduce the equal-weight linear combination of hard-core operators
corresponding to all unblocked levels,

$B=\frac{1}{\sqrt{\mathcal{U}}}\sum_{i}^{\mathcal{U}}A_{i},$ (9)

so that

$H=-g\frac{\mathcal{U}}{\mathcal{L}}B^{\dagger}B+H_{\mathrm{blocked}}.$ (10)

Dropping the blocked subspace from this point on, the full Hamiltonian
including the p-h symmetry restoring energy shift is thus

$\tilde{H}=-g\frac{\mathcal{U}}{\mathcal{L}}B^{\dagger}B+\frac{\alpha}{\mathcal{L}}\hat{N},$
(11)

where

$\hat{N}=\sum_{i}^{\mathcal{U}}\hat{n}_{i}$ (12)

is the total number of electrons in the unblocked subset of levels (i.e., the
energy of unpaired electrons in the blocked levels is not included).

The eigenstates for $M$ Cooper pairs in the $\mathcal{U}$ levels are

$\Psi^{\mathcal{L},\mathcal{U}}_{M}=\mathcal{N}_{M}|M\rangle=\mathcal{N}_{M}(B^{\dagger})^{M}|0\rangle,$
(13)

where $|M\rangle=(B^{\dagger})^{M}|0\rangle$ and $\mathcal{N}_{M}$ is a
normalization factor (see Appendix B). In the following, the superscripts
$\mathcal{L},\mathcal{U}$ in $\Psi^{\mathcal{L},\mathcal{U}}_{M}$ will be
dropped when no ambiguity may result; when a single superscript is indicated,
it corresponds to $\mathcal{U}$.

We have $H\Psi_{0}=0$ and $H\Psi_{1}=-g(\mathcal{U}/\mathcal{L})\Psi_{1}$.
General eigenvalues can be computed by recursion (see Appendix A), finding

$H|M\rangle=-g(\mathcal{U}/\mathcal{L})c_{M}|M\rangle$ (14)

with

$c_{M}=\frac{(1+\mathcal{U}-M)M}{\mathcal{U}}.$ (15)

The ground state energy of a system of size $\mathcal{L}$ with $\mathcal{U}$
unblocked levels containing $M$ Cooper pairs is thus

$E^{\mathcal{L},\mathcal{U}}_{M}=-2\alpha\frac{(1+\mathcal{U}-M)M}{\mathcal{L}},$
(16)

while the shifted eigenvalues (corresponding to $\tilde{H}$) are

$\tilde{E}^{\mathcal{L},\mathcal{U}}_{M}=-2\alpha\frac{(\mathcal{U}-M)M}{\mathcal{L}}.$
(17)

The interpretation of this expression is very simple. Each of the $M$ Cooper
pairs resonates with the $\mathcal{U}-M$ empty levels in the system, each
combination contributing $2\alpha/\mathcal{L}=G$ to the total energy. We note
that $\tilde{E}$ does not include the energy of single electrons in the
blocked levels, which is not included in Eq. (11). We therefore define the
total energy $\hat{E}$ as

$\begin{split}\hat{E}^{\mathcal{L},\mathcal{U}}_{M}&=\tilde{E}^{\mathcal{L},\mathcal{U}}_{M}+\frac{\alpha}{\mathcal{L}}(\mathcal{L}-\mathcal{U})\\\
&=\alpha\frac{-2(\mathcal{U}-M)M+\mathcal{L}-\mathcal{U}}{\mathcal{L}}.\end{split}$

For half filling, $n=\mathcal{L}$ and $M=\mathcal{L}/2$, and one has for
$\mathcal{U}=\mathcal{L}$

$\tilde{E}^{\mathcal{L},\mathcal{L}}_{\mathcal{L}/2}=-\alpha\frac{\mathcal{L}}{2}\equiv\mathcal{E}.$
(18)

Also, we find

$\tilde{E}^{\mathcal{L},\mathcal{L}}_{\mathcal{L}/2+1}=\mathcal{E}+\frac{2\alpha}{\mathcal{L}},\,\,\tilde{E}^{\mathcal{L},\mathcal{L}}_{\mathcal{L}/2+2}=\mathcal{E}+\frac{8\alpha}{\mathcal{L}},\ldots$
(19)

thus additional Cooper pairs have an energy cost of order $1/\mathcal{L}$. In
the thermodynamic limit ($\mathcal{L}\to\infty$), the ground state would be
macroscopically degenerate, as is the case in bulk SCs.

We now consider a single quasiparticle in the system, i.e., we make one of the
$\mathcal{L}$ levels singly occupied, so that $\mathcal{U}=\mathcal{L}-1$. The
energy cost in the $n=\mathcal{L}+1$ charge sector is

$E^{qp,+}=\left(\tilde{E}^{\mathcal{L},\mathcal{L}-1}_{\mathcal{L}/2}+\tilde{\epsilon}\right)-\mathcal{E}=\left(1+\frac{1}{\mathcal{L}}\right)\alpha.$
(20)

and in the $n=\mathcal{L}-1$ charge sector it is

$E^{qp,-}=\left(\tilde{E}^{\mathcal{L},\mathcal{L}-1}_{\mathcal{L}/2-1}+\tilde{\epsilon}\right)-\mathcal{E}=\left(1+\frac{1}{\mathcal{L}}\right)\alpha.$
(21)

In the equal-charge ($n=\mathcal{L}$) sector, the lowest excitation contains
two quasiparticles. Its energy is

$\left(\tilde{E}^{\mathcal{L},\mathcal{L}-2}_{\mathcal{L}/2-1}+2\tilde{\epsilon}\right)-\mathcal{E}=2\alpha.$
(22)

The cost of a quasiparticle excitation in flat-band superconductors is equal
to the bare pairing interaction [10, 73, 74].

#### Pairing strength

The strength of pairing correlations may be quantified by evaluating the
following correlators [75]:

$\begin{split}\bar{v}_{i}&=\langle
c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{i\downarrow}c_{i\uparrow}\rangle^{1/2}=\langle
A_{i}^{\dagger}A_{i}\rangle^{1/2},\\\ \bar{u}_{i}&=\langle
c_{i\downarrow}c_{i\uparrow}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}\rangle^{1/2}=\langle
A_{i}A_{i}^{\dagger}\rangle^{1/2},\\\
\bar{\Delta}&=\frac{g}{\mathcal{L}}\sum_{i}\bar{v}_{i}\bar{u}_{i},\\\
\bar{\Delta}^{\prime}&=\frac{g}{\mathcal{L}}\sum_{i}\Bigl{(}\langle
n_{i\uparrow}n_{i\downarrow}\rangle-\langle n_{i\uparrow}\rangle\langle
n_{i\downarrow}\rangle\Bigr{)}^{1/2}.\end{split}$ (23)

$\bar{v}_{i}$ is the amplitude to find a pair of states occupied, while
$\bar{u}_{i}$ is the amplitude to find a pair of states empty. If the system
is homogeneous, all amplitudes are the same, $\bar{v}\equiv\bar{v}_{i}$ and
$\bar{u}\equiv\bar{u}_{i}$, so that

$\bar{\Delta}=g\bar{v}\bar{u}.$ (24)

Both $\bar{\Delta}$ and $\bar{\Delta}^{\prime}$ reduce to
$\Delta_{\mathrm{BCS}}$ in the thermodynamic limit [42].

At half filling, the results is obviously $\bar{v}^{2}=\bar{u}^{2}=1/2$,
because there is probability $1/2$ to find any given level unoccupied and the
same probabilty to find it doubly occupied, and by construction it cannot be
singly occupied, since no level is blocked. Thus

$\bar{\Delta}=\frac{g}{2}=\alpha.$ (25)

In general, the amplitudes $\hat{v}_{i}$ and $\hat{u}_{i}$ are related, since
$A_{i}^{\dagger}A_{i}=A_{i}^{\dagger}A_{i}-(1-\hat{n}_{i})$, thus

$\bar{v}^{2}=\bar{u}^{2}-\langle\Psi_{M}|(1-\hat{n}_{i})|\Psi_{M}\rangle=\bar{u}^{2}-1+\frac{2M}{\mathcal{U}}.$
(26)

In the subspace where all levels are either unoccupied or doubly occupied, we
also have $\bar{v}^{2}+\bar{u}^{2}=1$. Then it follows immediately that

$\bar{v}^{2}=\frac{M}{\mathcal{U}},\quad\bar{u}^{2}=\frac{\mathcal{U}-M}{\mathcal{U}},\quad\bar{\Delta}=2\alpha\frac{\sqrt{M(\mathcal{U}-M)}}{\mathcal{U}}.$
(27)

The other estimate, $\bar{\Delta}^{\prime}$, can also be easily computed:

$\bar{\Delta}^{\prime}=\frac{g}{\mathcal{L}}\mathcal{U}\left(\frac{M}{\mathcal{U}}-\frac{M}{\mathcal{U}}\frac{M}{\mathcal{U}}\right)^{1/2},$
(28)

because in a homogeneous system
$\langle\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\rangle=\langle\hat{P}_{2i}\rangle=M/\mathcal{U}$
and $\langle n_{i\sigma}\rangle=M/\mathcal{U}$. At half filling, this yields

$\bar{\Delta}^{\prime}=\alpha.$

In this flat-band model, the pairing correlations are thus directly
proportional to the coupling-strength $\alpha$. This is a known feature of
strong-coupling SC systems [10, 73, 74]. The exponential reduction in the
standard weak-coupling BCS theory is a consequence of the finite width of the
electron band.

###  IV.2 YSR states for a Kondo impurity

We now include a magnetic impurity into consideration. First we consider the
case of a pure exchange scatterer and take the large-$U$ limit of the impurity
Hamiltonian at $\nu=1$ to obtain an effective Kondo Hamiltonian [76, 77]:

$H_{\mathrm{QD}}=\frac{J}{\mathcal{L}}\sum_{i,j}\mathbf{S}\cdot\mathbf{s}_{i,j},$
(29)

where the Kondo exchange coupling is

$J=8v^{2}/U=\frac{16}{\pi}(\Gamma/U),$ (30)

and the inter-level spin operators are defined as

$\mathbf{s}_{i,j}=\frac{1}{2}c^{\dagger}_{i,\alpha}\boldsymbol{\sigma}_{\alpha,\beta}c_{j,\beta}.$
(31)

Here $\boldsymbol{\sigma}=\\{\sigma^{x},\sigma^{y},\sigma^{z}\\}$ is a vector
of Pauli matrices and the expression is summed over repeated spin indexes
$\alpha$ and $\beta$. Note that in this discrete model $\Gamma$ does not have
the significance of a level width, but it nevertheless quantifies the coupling
strength. Equivalently, we may write

$H_{\mathrm{QD}}=J\mathbf{S}\cdot\mathbf{s},$ (32)

with $\mathbf{s}$ the spin of an electron in the orbital $f$, see Eq. (4),

$\mathbf{s}=\frac{1}{2}f^{\dagger}_{\alpha}\boldsymbol{\sigma}_{\alpha,\beta}f_{\beta}.$
(33)

For half filling with $M=\mathcal{L}/2$ Cooper pairs (assuming even
$\mathcal{L}$) and one electron in the impurity level (odd total number of
electrons $n=\mathcal{L}+1$, i.e., spin-doublet sector), the impurity remains
entirely decoupled from the SC in the ground state. The reason is that the
product state $\psi=|M\rangle\otimes|\sigma\rangle$ is an eigenstate of
$H_{\mathrm{QD}}$. First we note that the state $|M\rangle$ is an eigenstate
of $\hat{s}_{z}$ with $\hat{s}_{z}|M\rangle=0$. We thus only need to consider
the transverse $S^{-}s^{+}$ terms. This operator is composed of pairs

$c^{\dagger}_{i\uparrow}c_{j\downarrow}+c^{\dagger}_{j\uparrow}c_{i\downarrow}.$
(34)

We have

$[c^{\dagger}_{i\uparrow}c_{j\downarrow},A_{j}^{\dagger}]=-c^{\dagger}_{i\uparrow}c^{\dagger}_{j\uparrow},\quad[c^{\dagger}_{j\uparrow}c_{i\downarrow},A_{i}^{\dagger}]=c^{\dagger}_{i\uparrow}c^{\dagger}_{j\uparrow}.$
(35)

This leads to a cancellation for any $i\neq j$ pair. Furthermore, the diagonal
$i=j$ terms are trivially zero due to Pauli exclusion principle. Thus we have
$s^{+}|M\rangle=0$, showing that $|M\rangle$ is indeed an eigenstate.

The same commutation relations also establish that the pure exchange
interaction preserves the decoupling between the blocked (singly occupied)
levels $\mathcal{B}$ and the unblocked (zero or double occupancy) levels
$\mathcal{U}$. Thus the Hamiltonian in the subspace of $\mathcal{B}$ blocked
levels is

$H_{\mathcal{B}}=J\mathbf{S}\cdot\mathbf{s},$ (36)

where $\mathbf{s}$ is now constrained only over the blocked levels. This is an
exact statement. The ground state in the even occupancy $n=\mathcal{L}+2$
sector results from the interaction between the single quasiparticle and the
impurity spin (i.e., the YSR singlet state). This is clearly a simple two-body
problem of two spins coupled by an isotropic exchange interaction, with one
singlet eigenstate at energy $-3/4J$ and one triplet eigenstate at energy
$+1/4J$. The creation of the quasiparticle has an energy cost of
$(1+1/\mathcal{L})\alpha$. The excitation energy is thus

$E_{\mathrm{YSR}}=\left(1+\frac{1}{\mathcal{L}}\right)\alpha-\frac{3J}{4}.$
(37)

The YSR singlet binding energy is

$\Delta E=E^{qp,+}-E_{\mathrm{YSR}}=\frac{3}{4}J=\frac{12}{\pi}(\Gamma/U).$
(38)

There is a YSR triplet state in the “continuum” with an energy surplus of
$1/4J$.

###  IV.3 YSR states for an Anderson impurity

At finite $U$, an electron from a doubly occupied bath level can hop on the
impurity site, and an electron can hop from the impurity site to an unoccupied
bath level. The decoupling between the blocked levels $\mathcal{B}$ and the
unblocked levels $\mathcal{U}$ is lost: the subsectors mix. We need to
separately discus the cases with different fermionic parity (and consequently
total spin).

####  IV.3.1 Odd total number of electrons (spin-doublet)

We first consider the case with an odd total number of electrons. The $z$
component of total spin is a conserved quantity; we will work in the
$s_{z}=1/2$ subspace. The following normalized states are coupled by the
hopping term:

$\begin{split}\psi_{0,b}&=|0\rangle\otimes
c^{\dagger}_{b\uparrow}\Psi^{\mathcal{L}\backslash b}_{M},\\\
\Psi_{1}&=|1\rangle\otimes\Psi^{\mathcal{L}}_{M},\\\
\psi_{1,bb^{\prime}}&=|1\rangle\otimes
c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M-1},\\\ \psi_{2,b}&=|2\rangle\otimes
c^{\dagger}_{b\uparrow}\Psi^{\mathcal{L}\backslash b}_{M-1}.\end{split}$ (39)

The ket $|i\rangle$ represents the impurity state with $i$ electrons, defined
as $|1\rangle=d^{\dagger}_{\uparrow}|0\rangle$ and
$|2\rangle=d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}|1\rangle$. In
$\psi_{0,b}$ and $\psi_{2,b}$ the level $b$ is blocked; here $1\leq
b\leq\mathcal{L}$. The notation $\mathcal{L}\backslash b$ in the superscript
denotes that the set of unblocked levels $\mathcal{U}$ encompasses all levels
except $b$, while the subscript represent the number of Cooper pairs. In a
similar vein, in $\psi_{1,bb^{\prime}}$ two levels are blocked. We note that
states with two quasiparticles with equal spin in the superconductor (and
opposite spin in the QD, such that the total spin is $S=1/2$ and $S_{z}=1/2$)
do not couple with the subset of states in Eq. (39) in the flatband limit.

We now form normalized equal-weight superpositions of $\psi_{0,b}$ and
$\psi_{2,b}$,

$\Psi_{0}=\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}\psi_{0,b},\quad\Psi_{2}=\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}\psi_{2,b},$

as well as

$\begin{split}\Psi^{\prime}_{1}&=\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\psi_{1,bb^{\prime}}.\end{split}$

Note that the double sum is over all $\mathcal{L}$ so the two quasiparticles
form a singlet state. For large $\mathcal{L}$ this simplifies to

$\Psi^{\prime}_{1}\approx\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}\psi_{1,bb^{\prime}}.$

The energies of the doublet states $\Psi$ are

$\begin{split}E^{D}_{0}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-1}_{M}+\tilde{\epsilon}-\epsilon,\\\
E^{D}_{1}&={\tilde{E}}^{\mathcal{L},\mathcal{L}}_{M},\\\
E^{D^{\prime}}_{1}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-2}_{M-1}+2\tilde{\epsilon},\\\
E^{D}_{2}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-1}_{M-1}+\tilde{\epsilon}+\epsilon+U.\end{split}$
(40)

Here we included the energy shifts $\tilde{\epsilon}$ (hence we use energies
$\tilde{E}$), which are important for this consideration, and we have
subtracted $\epsilon$. The superscript $D$ denotes that these are energies in
the spin-doublet subspace.

To compute the matrix elements, we make use of an expression derived in
Appendix B after setting $M=\mathcal{L}/2$:

$f^{\dagger}_{\sigma}\Psi^{\mathcal{L}}_{\mathcal{L}/2}=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{\mathcal{L}}}\sum_{b}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L}\backslash
b}_{\mathcal{L}/2},$ (41)

thus $f^{\dagger}_{\sigma}d_{\sigma}\Psi_{1}=(1/\sqrt{2})\Psi_{0}$. An
analogous expression can be obtained for
$f_{\sigma}\Psi^{\mathcal{L}}_{\mathcal{L}/2}$, giving
$d^{\dagger}_{\sigma}f_{\sigma}\Psi_{1}=(1/\sqrt{2})\Psi_{2}$. In a similar
way we obtain all other matrix elements.

Figure 5: Schematic representation of the states coupled in the doublet
subspace. In the flat-band limit, no other states couple with this subset.
Left line: QD level. Right lines: SC levels. Quasiparticles (singly occupied
SC levels) represent superpositions of states where the single electron
occupies any of the superconductor $\mathcal{L}$ levels.

We obtain a four-level problem with the Hamiltonian in the orthonormal basis
$\\{\Psi_{1},-\Psi_{1}^{\prime},\Psi_{0},\Psi_{2}\\}$ that can be expressed in
the large-$\mathcal{L}$ limit and at half-filling as

$H^{D}_{\mathrm{eff}}=\begin{pmatrix}E^{D}_{1}&0&v/\sqrt{2}&-v/\sqrt{2}\\\
0&E^{D^{\prime}}_{1}&v/\sqrt{2}&v/\sqrt{2}\\\
v/\sqrt{2}&v/\sqrt{2}&E^{D}_{0}&0\\\
-v/\sqrt{2}&v/\sqrt{2}&0&E^{D}_{2}\end{pmatrix}.$ (42)

The schematic representation of this Hamiltonian is shown in Fig. 5. The full
expression for $H^{D}_{\mathrm{eff}}$ for arbitrary $\mathcal{L}$ and $M$ with
exact matrix elements is given in Appendix D. Close to half filling and for
small to moderate $v$, the state $\Psi_{1}^{\prime}$ is a highly excited
states and plays little role (it admixes in the ground state wavefunction only
as a $v^{4}$ correction). We may thus drop the second row and second column in
the matrix in Eq. (42) and write

$H^{D}_{\mathrm{eff}}=\begin{pmatrix}E^{D}_{1}&v/\sqrt{2}&-v/\sqrt{2}\\\
v/\sqrt{2}&E^{D}_{0}&0\\\ -v/\sqrt{2}&0&E^{D}_{2}\end{pmatrix}.$ (43)

The eigenenergies can be obtained from the cubic roots. At half filling in the
particle-hole symmetric point ($\epsilon=-U/2$) the problem simplifies even
further. We then have

$\begin{split}E^{D}_{0}=E^{D}_{2}&=-\alpha\frac{\mathcal{L}}{2}+\alpha+U/2\equiv
E^{D}_{02},\\\ E^{D}_{1}&=-\alpha\frac{\mathcal{L}}{2}.\end{split}$ (44)

The eigenvector $(0,1,1)$ with eigenvalue $E^{D}_{02}$ decouples. The ground
state can then be obtained by diagonalising the $2\times 2$ matrix

$H^{D}_{\mathrm{eff}}=\begin{pmatrix}E^{D}_{1}&v\\\ v&E^{D}_{02}\\\
\end{pmatrix}.$ (45)

Reference calculations demonstrating that the method is exact are presented in
Appendix D. Even using the approximate expressions from this section which are
strictly correct only in the $\mathcal{L}\to\infty$ limit, we find excellent
agreement with the DMRG results even at moderate $\mathcal{L}=50$, see Fig. 6
(top panel) where we plot the lowest eigenvalue of $H^{D}_{\mathrm{eff}}$
using the dashed line and the reference results as full line.

Figure 6: Ground state eigenenergies in the flat-band limit for $U=10\alpha$,
specifically $\alpha=0.1$ and $U=1$ in units of $D$, $\mathcal{L}=50$. Full
line in both panels: DMRG results. (top) Odd-parity (spin-doublet) sector with
$n=\mathcal{L}+1$, $s_{z}=1/2$. Dashed line: lowest eigenvalue of the
Hamiltonian in Eq. (45). (bottom) Even-parity (spin-singlet) sector with
$n=\mathcal{L}+2$, $s_{z}=0$. Dashed line: lowest eigenvalue of the
Hamiltonian in Eq. (55).

####  IV.3.2 Even total number of electrons (spin-singlet)

We now consider the case with an even total number of electrons, specifically
$n=2M+2$, in the spin singlet state. We first enumerate all relevant states in
the subspace with $s_{z}=0$. There are two states with single occupancy of the
impurity, while in the sectors with $n_{\mathrm{imp}}=0$ and
$n_{\mathrm{imp}}=2$ we need to include two types of states, those with no
quasiparticles and those with two quasiparticles:

$\begin{split}\Psi_{0}&=|0\rangle\otimes\Psi^{\mathcal{L}}_{M+1},\\\
\psi^{\prime}_{0,bb^{\prime}}&=|0\rangle\otimes
c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M},\\\ \psi_{1,b,\sigma}&=|\sigma\rangle\otimes
c^{\dagger}_{b\bar{\sigma}}\Psi^{\mathcal{L}\backslash b}_{M},\\\
\Psi_{2}&=|2\rangle\otimes\Psi^{\mathcal{L}}_{M},\\\
\psi^{\prime}_{2,bb^{\prime}}&=|2\rangle\otimes
c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M-1},\end{split}$ (46)

with energies

$\begin{split}E^{S}_{0}&={\tilde{E}}^{\mathcal{L},\mathcal{L}}_{M+1}-\epsilon,\\\
{E^{S}_{0}}^{\prime}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-2}_{M}+2\tilde{\epsilon}-\epsilon,\\\
E^{S}_{1}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-1}_{M}+\tilde{\epsilon},\\\
E^{S}_{2}&={\tilde{E}}^{\mathcal{L},\mathcal{L}}_{M}+\epsilon+U,\\\
{E^{S}_{2}}^{\prime}&={\tilde{E}}^{\mathcal{L},\mathcal{L}-2}_{M-1}+2\tilde{\epsilon}+\epsilon+U.\end{split}$
(47)

For one electron beyond half filling (i.e., for $M=\mathcal{L}/2$) and
$\epsilon=-U/2$, this yields

$\begin{split}E^{S}_{0}&=-\alpha\frac{\mathcal{L}}{2}+U/2+\frac{2\alpha}{\mathcal{L}},\\\
{E^{S}_{0}}^{\prime}&=-\alpha\frac{\mathcal{L}}{2}+2\alpha+U/2+\frac{2\alpha}{\mathcal{L}}=E^{S}_{0}+2\alpha,\\\
E^{S}_{1}&=-\alpha\frac{\mathcal{L}}{2}+\alpha+\frac{\alpha}{\mathcal{L}},\\\
E^{S}_{2}&=-\alpha\frac{\mathcal{L}}{2}+U/2,\\\
{E^{S}_{2}}^{\prime}&=-\alpha\frac{\mathcal{L}}{2}+2\alpha+U/2=E^{S}_{2}+2\alpha.\end{split}$
(48)

The difference between $E^{S}_{0}$ and $E^{S}_{2}$ is a
$\mathcal{O}(1/\mathcal{L})$ correction,
$E^{S}_{0}-E^{S}_{2}=2\alpha/\mathcal{L}$ and vanishes in the thermodynamic
limit.

We now form normalized equal superpositions of $\psi^{\prime}_{0,bb^{\prime}}$
and $\psi^{\prime}_{2,bb^{\prime}}$ states:

$\begin{split}\Psi^{\prime}_{0}&=\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\psi^{\prime}_{0,bb^{\prime}},\\\
\Psi^{\prime}_{2}&=\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\psi^{\prime}_{2,bb^{\prime}},\end{split}$ (49)

we can be simplified in the $\mathcal{L}\to\infty$ limit to

$\begin{split}\Psi^{\prime}_{0}&\approx\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M},\\\
\Psi^{\prime}_{2}&\approx\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M-1},\\\ \end{split}$

These are spin-singlet states. Furthermore, we make a normalized spin-singlet
combination of $\psi_{1,b,\sigma}$:

$\begin{split}\Psi_{1}&=\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}\frac{1}{\sqrt{2}}\left(\psi_{1,b,\uparrow}-\psi_{1,b,\downarrow}\right)\\\
&=\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}\frac{1}{\sqrt{2}}\left(|\uparrow\rangle\otimes
c^{\dagger}_{b,\downarrow}-|\downarrow\rangle\otimes
c^{\dagger}_{b,\uparrow}\right)\Psi^{\mathcal{L}\backslash b}_{M}.\end{split}$
(50)

Making use of Eq. (41) at half-filling ($M=\mathcal{L}/2$), this can also be
expressed as

$\begin{split}\Psi_{1}&\approx\left(|\uparrow\rangle\otimes
f^{\dagger}_{\downarrow}-|\downarrow\rangle\otimes
f^{\dagger}_{\uparrow}\right)\Psi^{\mathcal{L}}_{M}.\end{split}$ (51)

We note that no spin-triplet states are included here.

To find the matrix elements, we make use of the result from Appendix B in the
large $\mathcal{L}$ limit:

$f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\Psi^{\mathcal{L}}_{M}\approx\frac{1}{2}\Psi^{\mathcal{L}}_{M+1}+\frac{1}{2}\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}}_{M}.$
(52)

We find, for example,

$\begin{split}\langle\Psi_{0}|f^{\dagger}_{\uparrow}d_{\uparrow}|\Psi_{1}\rangle&=\langle
0|(\Psi^{\mathcal{L}}_{M+1})^{*}(-d_{\uparrow}f^{\dagger}_{\uparrow})|\uparrow\rangle\otimes
f_{\downarrow}^{\dagger}\Psi^{\mathcal{L}}_{M}\\\ &=\langle
0|d_{\uparrow}|\uparrow\rangle(\Psi^{\mathcal{L}}_{M+1})^{*}f^{\dagger}_{\uparrow}f_{\downarrow}^{\dagger}\Psi^{\mathcal{L}}_{M}\\\
&\approx 1/2,\end{split}$ (53)

and similarly

$\begin{split}\langle\Psi_{0}|f^{\dagger}_{\downarrow}d_{\downarrow}|\Psi_{1}\rangle\approx
1/2,\end{split}$ (54)

giving a total electron hopping matrix element of $1$. Other matrix elements
can be computed in an analogous way. The effective Hamiltonian in the basis
$\\{\Psi_{1},\Psi_{0},\Psi^{\prime}_{0},\Psi_{2},-\Psi^{\prime}_{2}\\}$ in the
large-$\mathcal{L}$ limit is then

$H^{S}_{\mathrm{eff}}=\begin{pmatrix}E^{S}_{1}&v&v&v&v\\\ v&E^{S}_{0}&0&0&0\\\
v&0&{E^{S}_{0}}^{\prime}&0&0\\\ v&0&0&E^{S}_{2}&0\\\
v&0&0&0&{E^{S}_{2}}^{\prime}\\\ \end{pmatrix}.$ (55)

The signs of basis states have been chosen so that the out-of-diagonal matrix
elements are all positive (in the sense of having the same sign as $v$). The
schematic representation of this Hamiltonian is shown in Fig. 7. If the
$1/\mathcal{L}$ corrections in diagonal elements are neglected (as is valid in
the large $\mathcal{L}$ limit), so that

$\begin{split}E^{S}_{02}&=-\alpha\frac{\mathcal{L}}{2}+U/2,\\\
{E^{S}_{02}}^{\prime}&=-\alpha\frac{\mathcal{L}}{2}+2\alpha+U/2=E^{S}_{02}+2\alpha,\\\
E^{S}_{1}&=-\alpha\frac{\mathcal{L}}{2}+\alpha,\\\ \end{split}$ (56)

the problem simplifies further to a $3\times 3$ matrix problem

$H^{S}_{\mathrm{eff}}=\begin{pmatrix}E^{S}_{1}&\sqrt{2}v&\sqrt{2}v\\\
\sqrt{2}v&E^{S}_{02}&0\\\ \sqrt{2}v&0&{E^{S}_{02}}^{\prime}\\\ \end{pmatrix}.$
(57)

Figure 7: Schematic representation of the states coupled in the singlet
subspace. In the flat-band limit, no other states couple with this subset.

The numeric diagonalisation of this Hamiltonian produces a result which agrees
very well with the DMRG results even for $\mathcal{L}=50$, see Fig. 6 (bottom
panel). Exact matrix elements for any values of $\mathcal{L}$ and $M$ are
provided in Appendix D.

###  IV.4 Dependence on hybridisation strength

We now consider the $\Gamma$-dependence of the singlet-doublet excitation
energy for different values of $U$.

####  IV.4.1 Large-$U$ regime

At the particle-hole symmetric point, the lowest eigenvalue of the Hamiltonian
$H^{D}_{\mathrm{eff}}$ in Eq. (45) is

$E_{D}=-\alpha\frac{\mathcal{L}-1}{2}+\frac{U-\sqrt{16v^{2}+(U+2\alpha)^{2}}}{4}.$
(58)

The expansion around the $U\to\infty$ limit is

$\begin{split}E_{D}&\approx-\alpha\frac{\mathcal{L}}{2}-2\frac{v^{2}}{U}+\frac{4v^{2}\alpha}{U^{2}}+\mathcal{O}(1/U^{3})\\\
&=-\alpha\frac{\mathcal{L}}{2}-\frac{J}{4}+\frac{J\alpha}{2U}+\mathcal{O}(1/U^{3}),\end{split}$
(59)

where $J=8v^{2}/U$ as defined in Eq. (30). The first term is the energy of the
half-filled decoupled bath, Eq. (18). The second term is the trivial linear
shift of energies in the Schrieffer-Wolff transformation [76], $-J/4$. The
third term is the leading order $J/U$ correction. As expected, it is
proportional to $\alpha$ and $J$, and inversely proportional to $U$ (the
energy cost of valence fluctuations), since it is due to the perturbation of
the SC state by the hybridisation with the impurity.

The lowest eigenvalue of the Hamiltonian $H^{S}_{\mathrm{eff}}$ in Eq. (57) in
the large-$U$ limit is

$\begin{split}E_{S}&\approx-\alpha\frac{\mathcal{L}-2}{2}-\frac{8v^{2}}{U}+\mathcal{O}(1/U^{3})\\\
&=-\alpha\frac{\mathcal{L}-2}{2}-J+\mathcal{O}(1/U^{3}).\end{split}$ (60)

Now we consider the difference, i.e., the Yu-Shiba-Rusinov excitation energy
$E_{\mathrm{YSR}}=E_{S}-E_{D}$. At large $U$ we then recover the Kondo-limit
results from Sec. IV.2,

$E_{\mathrm{YSR}}=\alpha-\frac{3J}{4},$ (61)

up to $1/\mathcal{L}$ corrections.

For moderate $U$ one finds deviations from the simple linear-in-$\Gamma$
behavior, see Fig. 8. The YSR energy at the same $\Gamma/U$ ratio is found to
be higher, the lower $U$ is. The charge fluctuations thus reduce the effects
of the exchange coupling (at the same value of $J$). Note that at large
$\Gamma$ the eigenvalue reaches the limiting value $-(1-1/\mathcal{L})\alpha$
(“gap edge”) with a finite slope, indicating a level crossing. This is the
main differenece compared with the problem with a finite bandwidth, where the
subgap state experiences a level repulsion and at large $\Gamma$
asymptotically approaches the gap edge which is accompanied by a continuous
transfer of the spectral weight from the subgap states to the continuum of
Bogoliubov excitations [78]. In the flat-band case there is no continuum, only
a set of degenerate excitations which do not hybridize with the subgap state,
thus the hybridized state can cross this set unperturbed. The nonlinearity in
$\Gamma$ can be obtained by expanding the lowest eigenvalue of
$H^{S}_{\mathrm{eff}}$ to the next order in $1/U$:

$\begin{split}E_{S}&\approx-\alpha\frac{\mathcal{L}-2}{2}-\frac{8v^{2}}{U}+\frac{128v^{4}-32v^{2}\alpha^{2}}{U^{3}}\\\
&=-\alpha\frac{\mathcal{L}-2}{2}-\frac{8\Gamma}{\pi\rho
U^{2}}-\frac{32\frac{\Gamma}{U}\left[\left(\pi\rho\alpha\right)^{2}-4\pi\rho\frac{\Gamma}{U}\right]}{(\pi\rho
U)^{3}}.\end{split}$ (62)

####  IV.4.2 Small-$U$ limit

We now consider the expansion for small values of $U$. For Eq. (58) we find

$E_{D}\approx-\alpha\frac{\mathcal{L}-1}{2}+\frac{U}{4}\left(1-\frac{\alpha}{\sqrt{\alpha^{2}+4v^{2}}}\right)-\frac{\sqrt{\alpha^{2}+4v^{2}}}{2}.$
(63)

Taking the small-$\Gamma$ limit gives $E_{D}\approx-\alpha\mathcal{L}/2$,
thereby fully cancelling out the $U$ dependence.

The lowest eigenvalue of the Hamiltonian $H^{S}_{\mathrm{eff}}$ in Eq. (57) in
the small-$U$ limit is

$E_{S}\approx-\alpha\frac{\mathcal{L}-2}{2}+\frac{U}{2}\frac{\alpha^{2}+2v^{2}}{\alpha^{2}+4v^{2}}-\sqrt{\alpha^{2}+4v^{2}}.$
(64)

Taking the $\Gamma\to 0$ limit, we find
$E_{S}\approx-\alpha\mathcal{L}/2+U/2$. At zero $\Gamma$, the excitation
energy is

$\Delta E(\Gamma=0)=\lim_{\Gamma\to 0}\left(E_{S}-E_{D}\right)=U/2,$ (65)

as expected and as seen in Fig. 8. The lowest excitation is not a YSR singlet,
but a ABS-like state with a QD in a proximitized state ($|0\rangle+|2\rangle$
superposition). We note in passing that the order of taking the limits is
important; here we took the limit of $U\to 0$ followed by $\Gamma\to 0$.

The $\Gamma$ dependence to lowest order is

$\Delta
E=E_{S}-E_{D}\approx\frac{U}{2}-\frac{\Gamma}{\pi\rho\alpha}\left(1+\frac{3U}{2\alpha}\right).$
(66)

The energy reduction in the Andreev-bound-state regime is controlled by the
$\Gamma/\alpha$ ratio, a non-zero repulsion $U$ is merely an order $U/\alpha$
correction to the prefactor.

Figure 8: Excitation energy $E_{\mathrm{YSR}}=E_{S}-E_{D}$ in the flat-band
limit for the $\Delta N=1$ transitions. Parameters are $\mathcal{L}=50$,
$\alpha=0.1$. Full lines: DMRG results. Dashed lines: analytical calculation.
The small deviations are due to the $1/\mathcal{L}$ correction that for
simplicity have not been included here.

####  IV.4.3 Special case of $U=2\alpha$

In the limit of $U=2\alpha$, some states in the singlet subspace become
degenerate: $E^{S}_{02}=E^{S}_{1}$. In this regime, we find the following
expansions in the weak hybridisation limit:

$\begin{split}E_{D}&=-\alpha\frac{\mathcal{L}}{2}-\frac{v^{2}}{2\alpha},\\\
E_{S}&=-\alpha\frac{\mathcal{L}}{2}-\frac{v^{2}}{2\alpha}+\alpha-\sqrt{2}v.\end{split}$
(67)

This leads to the interesting result that

$E_{S}-E_{D}=\alpha-\sqrt{\frac{2\Gamma}{\pi\rho}}.$ (68)

This square root singularity is indeed observed in the results shown in Fig.
8.

The singular behaviour observed at $U=2\alpha$ marks the boundary between the
regimes of proximitized (Andreev) subgap states for $U<2\alpha$ and the Yu-
Shiba-Rusinov states for $U>2\alpha$. This is a nontrivial observation because
it is generally believed that the cross-over between these two limits is
perfectly smooth. In fact, there exists a clear anomaly at $U=2\alpha$ in the
low-$\Gamma$ limit.

###  IV.5 Dependence on gate voltage

The formalism also correctly describes the transitions to different impurity
occupancies when $\epsilon=U(1/2-\nu)$ with $\nu\neq 1$. We find good
agreement with the reference DMRG results, see Fig. 9. This is the case even
taking the truncated 3-dimensional basis from Eq. (45). For $\nu$ sufficiently
far away from $1$, the doublet solution enters the quasicontinuum: the black
line in Fig. 9 shoots past the gap edge $\alpha(1+1/\mathcal{L})$. The true
lowest excited state is then a decoupled Bogoliubov quasiparticle sitting at
the bottom of the quasicontinuum.

In the large-$U$ limit, the lowest lying eigenstates are

$\begin{split}E_{S}&\approx-\alpha\frac{\mathcal{L}-2}{2}+\frac{8v^{2}}{u(3+4\nu(\nu-2)},\\\
E_{D}&\approx-\alpha\frac{\mathcal{L}}{2}+\frac{2v^{2}}{u(3+4\nu(\nu-2)},\end{split}$
(69)

so that

$\begin{split}\Delta
E=E_{S}-E_{D}&=\frac{(\mathcal{L}+1)}{\mathcal{L}}\alpha+\frac{6v^{2}}{3+4\nu(\nu-2)}\\\
&\approx\alpha-\frac{3J}{4},\end{split}$ (70)

where we noted that away from half filling, the Kondo exchange coupling is

$\begin{split}J&=\frac{2v^{2}}{-\epsilon}+\frac{2v^{2}}{U+\epsilon}\\\
&=\frac{1}{4\nu(2-\nu)-3}\frac{8v^{2}}{U}.\end{split}$ (71)

As expected, to first order in $1/U$, the only effect of departure from the
particle-hole asymmetry is the modified value of the Kondo exchange. The
leading correction is of order $\alpha J/U$, as found in the previous section.

Figure 9: YSR subgap state dispersion (positive frequency branch) in the flat-
band limit. The lineshape corresponds to the conventional YSR loops.
Parameters are $U=100$, $\mathcal{L}=50$, $\Gamma/U=0.01$,
$\epsilon=U(1/2-\nu)$, where $\nu$ controls the impurity filling.

##  V Comparison to the BCS zero-bandwidth approximation

A commonly used approximation in QD-SC models is the ZBW limit of BCS, where
the SC is modelled by a single level with mean-field pairing interaction [79,
80, 81, 82]. Here we show that this description is equivalent to the RM in the
$\gamma\to 0$ limit. Consider the ZBW Hamiltonian

$H=\epsilon(\hat{n}-1)+U\hat{n}_{\uparrow}\hat{n}_{\downarrow}+v\sum_{i,\sigma}\left(f^{\dagger}_{\sigma}d_{\sigma}+d^{\dagger}_{\sigma}f_{\sigma}\right)-\Delta(f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}+\text{H.c.}).$
(72)

In the doublet subspace with the basis states
$\\{d^{\dagger}_{\uparrow},d^{\dagger}_{\uparrow}f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow},f^{\dagger}_{\uparrow},d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}f^{\dagger}_{\uparrow}\\}$
(applied to the zero occupancy vacuum state) one finds the following matrix
representation:

$H^{D}_{\mathrm{ZBW}}=\begin{pmatrix}0&-\Delta&v&0\\\ -\Delta&0&0&-v\\\
v&0&-\epsilon&0\\\ 0&-v&0&U+\epsilon\end{pmatrix}.$ (73)

Alternatively, one may work in the basis

$\begin{split}&\left(d^{\dagger}_{\uparrow}+d^{\dagger}_{\uparrow}f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2}=d^{\dagger}_{\uparrow}\left(1+f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\\\
&\left(d^{\dagger}_{\uparrow}-d^{\dagger}_{\uparrow}f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2}=d^{\dagger}_{\uparrow}\left(1-f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\\\
&f^{\dagger}_{\uparrow},\\\
&d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}f^{\dagger}_{\uparrow},\end{split}$

where the matrix representation is:

$H^{D}_{\mathrm{ZBW}}=\begin{pmatrix}-\Delta&0&v/\sqrt{2}&-v/\sqrt{2}\\\
0&+\Delta&v/\sqrt{2}&v/\sqrt{2}\\\ v/\sqrt{2}&v/\sqrt{2}&-\epsilon&0\\\
-v/\sqrt{2}&v/\sqrt{2}&0&U+\epsilon\end{pmatrix},$ (74)

which is equivalent to $H^{D}_{\mathrm{eff}}$ in Eq. (42) up to a trivial
shift of all diagonal elements. The combination
$(1+f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow})/\sqrt{2}$ represents a
Cooper-pair state, while the combination
$(1-f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow})/\sqrt{2}$ represents a
two-quasiparticles state with the energy cost of $2\Delta$.

In the singlet subspace, we work in the basis (the signs are chosen so as to
produce positive signs for the out-of-diagonal matrix elements)

$\begin{split}&(d^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}-d^{\dagger}_{\downarrow}f^{\dagger}_{\uparrow})/\sqrt{2},\\\
&\left(1+f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\\\
&-\left(1-f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\\\
&d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}\left(1+f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\\\
&d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}\left(1-f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\right)/\sqrt{2},\end{split}$

where the matrix takes the following form

$H^{S}_{\mathrm{ZBW}}=\begin{pmatrix}0&v&v&v&v\\\ v&-\Delta-\epsilon&0&0&0\\\
v&0&\Delta-\epsilon&0&0\\\ v&0&0&-\Delta+U+\epsilon&0\\\
v&0&0&&\Delta+U+\epsilon\end{pmatrix}.$ (75)

Up to a trivial energy shift, this is exactly the same as
$H^{S}_{\mathrm{eff}}$ from Eq. (55).

The equivalence between the models is exact in the ZBW limit $\gamma\to 0$,
when $\alpha$ plays the role of excitation energy gap, akin to $\Delta$ in
BCS. The demonstration of this equivalence is one of the key results of this
work.

##  VI Conclusion

Even though superconductivity and Bose-Einstein condensation are apparently
quite different states of matter, they emerge as a result of the same
attractive pairing interaction in the opposite limits of pairing strength.
When a SC is coupled to a magnetic impurity, it induces discrete singlet
states inside the SC gap. We find that the subgap states persist throughout
the weak to strong pairing crossover, and we furthermore show that their
nature remains qualitatively the same. The differences are limited to the
scaling of hybridization strength and excitation gap, both due to the
decreasing importance of bath kinetic energy with increased pairing.

In the limit of strong pairing interaction the kinetic energy term can be
omitted. The single-particle levels in the bath are then degenerate, which
lets us obtain a flat-band Hamiltonian in the reduced Hilbert space spanned by
just a few many-particle states connected by the Hamiltonian. The BEC limit of
the problem is in this sense analytically solvable, and we show that the
results agree very well with exact calculations. We also show that the flat-
band Hamiltonian in the limit of large system is formally exactly equivalent
to the BCS zero bandwidth approximation.

ZBW approximations and similar toy models are very commonly applied when
modeling QDs coupled to SCs [83, 82, 84, 85, 86]. Our result explains their
surprisingly strong predictive power and thus provides a solid foundation for
their application in this context. Typically the argument for using the ZBW
approximation is that a single level at the gap edge is sufficient to
encompass the important phenomena, because the SC has a diverging density of
states at the gap edge and because these near-edge are the most important for
the subgap state formation: in fact, the subgap state weight originates
precisely from the depletion of the spectral weight from the superconducting
coherence peaks. The higher excited states largely remain unperturbed by the
impurity, and their small contributions can be discarded. We find that this is
a very reasonable approximation in a variety of situations.

Due to its simplicity, the flat-band model can be straightforwardly extended
in various ways. To accurately describe mesoscopic SC islands, it is important
to take into account the Coulomb repulsion between the island electrons [51].
This can be done by augmenting the model with a charging energy term
$E_{C}\hat{n}_{\mathrm{SC}}^{2}$. Because the RM conserves particle number we
are operating in the canonical ensemble, thus implementing such operator is
very simple and it is straightforward to extend the effective flat-band model
by adding $E_{C}$ to appropriate diagonal matrix elements. For a problem with
a single SC bath, this merely reduces the energy scale for charge fluctuations
from $U/2$ to $U/2+E_{c}$ [49, 51, 50]. The Zeeman splitting can also be
trivially included in the flat-band model. Additional triplet (and quadruplet)
states that become relevant at large magnetic fields can also be included in
consideration.

The analytical approach can also be extended to multichannel impurity
problems. These have additional degrees of freedom, namely the occupation
number differences between channels (i.e. the distribution of Cooper pairs
between the two superconductors). An appropriate linear combination of such
states corresponds to well defined phase differences which would appear in the
equivalent BCS formalism, while a charging term would induce phase
fluctuations. We will pursue this direction in our future work.

###### Acknowledgements.

We acknowledge discussions with Gorm Steffensson, Jens Paaske, and Martin
Žonda. We acknowledge the support of the Slovenian Research Agency (ARRS)
under P1-0044, P1-0416, and J1-3008.

##  Appendix A Eigenenergy calculation

Summing Eq. (6) over $i,j$ in the set $\mathcal{U}$ we obtain

$[B,B^{\dagger}]=1-\frac{1}{\mathcal{U}}\hat{N},$ (76)

where $\hat{N}=\sum_{i}^{\mathcal{U}}\hat{n}_{i}$ is the total number of
electrons in the unblocked levels.

In this appendix we evaluate $B^{\dagger}B$ on state
$|M\rangle=(B^{\dagger})^{M}|0\rangle$. By repeatedly commuting $B^{\dagger}$
and $B$ we find

$\begin{split}&B^{\dagger}B|M\rangle=\left(B^{\dagger}B\right)(B^{\dagger})^{M}|0\rangle\\\
&=\left[B^{\dagger}\left(B^{\dagger}B+1-\frac{1}{\mathcal{U}}\hat{N}\right)(B^{\dagger})^{M-1}\right]|0\rangle\\\
&=\left[(B^{\dagger})^{2}B(B^{\dagger})^{M-1}+\left(1-\frac{2(M-1)}{\mathcal{U}}\right)(B^{\dagger})^{M}\right]|0\rangle\\\
&=\Big{[}(B^{\dagger})^{3}B(B^{\dagger})^{M-2}+\left(1-\frac{2(M-2)}{\mathcal{U}}\right)(B^{\dagger})^{M}\\\
&\quad+\left(1-\frac{2(M-1)}{\mathcal{U}}\right)(B^{\dagger})^{M}\Bigr{]}|0\rangle\\\
&=\left[\sum_{m=1}^{M}\left(1-\frac{2(m-1)}{\mathcal{U}}\right)\right]|M\rangle=c_{M}|M\rangle,\end{split}$
(77)

with

$c_{M}=\frac{(1+\mathcal{U}-M)M}{\mathcal{U}}.$ (78)

##  Appendix B Action of creation operators on $\Psi^{\mathcal{L}}_{M}$

In this Appendix we discuss the action of fermionic operators on the
eigenstates $\Psi^{\mathcal{L},\mathcal{U}}_{M}$ of the Richardson’s model.
This defines the relations between the states from subspaces distinguished by
the different number $M$ of Cooper pairs. The subsectors differing by $\Delta
M=\pm 1$ are coupled due to the Cooper-pair-breaking property of the exchange
interaction induced by the presence of the magnetic impurity. The relations
derived in the following provide the full set of algebraic relations required
to determine the coupling between the accessible states and all matrix
elements in the flat band limit.

###  B.1 Single-quasiparticle state

In this section the number of levels is fixed to $\mathcal{L}$. When a single
superscript is indicated for brevity, the number indicates the number of
unblocked levels, i.e.,
$\Psi^{\mathcal{U}}_{M}\equiv\Psi^{\mathcal{L},\mathcal{U}}_{M}$. When no
superscript is indicated, all levels are unblocked, i.e.,
$\Psi_{M}\equiv\Psi^{\mathcal{L},\mathcal{L}}_{M}$.

Let us consider the spin-doublet state

$\psi_{\sigma}=f^{\dagger}_{\sigma}\Psi^{\mathcal{L},\mathcal{L}}_{M}=\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L},\mathcal{L}}_{M}$
(79)

with $\Psi^{\mathcal{L},\mathcal{L}}_{M}$ being the normalized eigenstate for
$M$ Cooper pairs in a system with $\mathcal{L}$ levels, all of them unblocked.
The application of the creation operator creates a superposition of states
with one blocked level indexed by $b$. This state is thus spanned by the
following orthonormal set of $\mathcal{L}$ eigenstates of the SC bath:

$\phi_{b\sigma}=c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L},\mathcal{L}\backslash
b}_{M},$ (80)

where $\mathcal{L}\backslash b$ in the superscript of $\Psi$ indicates that
the level $b$ is blocked. If all levels are equal, so are the coefficients,
thus

$\psi_{\sigma}=C\frac{1}{\sqrt{\mathcal{L}}}\sum_{b=1}^{\mathcal{L}}\phi_{b\sigma}.$
(81)

Let us now fix the proportionality constant $C$:

$\langle\psi_{\sigma}|\psi_{\sigma}\rangle=\langle\Psi^{\mathcal{L}}_{M}|f_{\sigma}f^{\dagger}_{\sigma}|\Psi^{\mathcal{L}}_{M}\rangle=1-\langle\Psi^{\mathcal{L}}_{M}|f^{\dagger}_{\sigma}f_{\sigma}|\Psi^{\mathcal{L}}_{M}\rangle=C^{2}.$
(82)

We expand in the site basis:

$\begin{split}f^{\dagger}_{\sigma}f_{\sigma}&=\frac{1}{\mathcal{L}}\sum_{i,j}c^{\dagger}_{i,\sigma}c_{j,\sigma}\\\
&=\frac{1}{\mathcal{L}}\sum_{i}\hat{n}_{i,\sigma}+\frac{1}{\mathcal{L}}\sum_{i\neq
j}c^{\dagger}_{i,\sigma}c_{j,\sigma}\\\
&=\frac{\hat{N}_{\sigma}}{\mathcal{L}}+\frac{1}{\mathcal{L}}\sum_{i\neq
j}c^{\dagger}_{i,\sigma}c_{j,\sigma}.\end{split}$ (83)

The first term gives $M/\mathcal{L}$. The second gives zero, because when
applied to the ket it produces a linear combination of terms which all contain
singly occupied levels, and such a state is orthogonal to
$\Psi^{\mathcal{L}}_{M}$ in the bra. Finally

$C^{2}=1-\frac{M}{\mathcal{L}}=\frac{\mathcal{L}-M}{\mathcal{L}}.$ (84)

We conclude that

$\psi_{\sigma}=f^{\dagger}_{\sigma}\Psi^{\mathcal{L},\mathcal{L}}_{M}=\sqrt{\frac{\mathcal{L}-M}{\mathcal{L}}}\frac{1}{\sqrt{\mathcal{L}}}\sum_{b}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L},\mathcal{L}\backslash
b}_{M}.$ (85)

This expression is exact. For $M\to 0$, in the absence of any Cooper pairs,
the prefactor becomes equal to 1, since the electron can be added to any
level. For $M\to\mathcal{L}$, the prefactor tends to 0, since there are only a
few levels available to accommodate the additional electron. Finally, for half
filling the prefactor becomes $1/\sqrt{2}$, thus

$\psi_{\sigma}=f^{\dagger}_{\sigma}\Psi^{\mathcal{L},\mathcal{L}}_{\mathcal{L}/2}=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{\mathcal{L}}}\sum_{b}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L}\backslash
b}_{\mathcal{L}/2}.$ (86)

Note that this state $\psi_{\sigma}$ is not normalized to 1.

In the same way one derives

$\psi^{\prime}_{\sigma}=f_{\bar{\sigma}}\Psi^{\mathcal{L},\mathcal{L}}_{M}=(-1)^{\sigma}\sqrt{\frac{M}{\mathcal{L}}}\frac{1}{\sqrt{\mathcal{L}}}\sum_{b}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L},\mathcal{L}\backslash
b}_{M-1},$ (87)

where $(-1)^{\uparrow}=-1$ and $(-1)^{\downarrow}=1$. This expression is
exact. It simplifies at half filling to

$\psi^{\prime}_{\sigma}=(-1)^{\sigma}f_{\bar{\sigma}}\Psi^{\mathcal{L},\mathcal{L}}_{\mathcal{L}/2}=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{\mathcal{L}}}\sum_{b}c^{\dagger}_{b,\sigma}\Psi^{\mathcal{L},\mathcal{L}\backslash
b}_{\mathcal{L}/2-1},$ (88)

The sign factor is due to the order convention of spin-up and spin-down
operators in the hard-core boson operators $A_{i}$.

###  B.2 Two-quasiparticle state

The other state of interest is the spin-singlet state

$\begin{split}\psi_{2}&=f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\Psi^{\mathcal{L}}_{M}\\\
&=\frac{1}{\mathcal{L}}\sum_{i,j}c^{\dagger}_{i\uparrow}c^{\dagger}_{j\downarrow}\Psi^{\mathcal{L}}_{M}\\\
&=\frac{1}{\mathcal{L}}\sum_{i}A^{\dagger}_{i}\Psi^{\mathcal{L}}_{M}+\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}}_{M}\\\
&=\frac{1}{\sqrt{\mathcal{L}}}B^{\dagger}\Psi^{\mathcal{L}}_{M}+\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}}_{M}.\end{split}$
(89)

We have

$B^{\dagger}\Psi^{\mathcal{L}}_{M}=B^{\dagger}\mathcal{N}_{M}|M\rangle=\mathcal{N}_{M}|M+1\rangle=\frac{\mathcal{N}_{M}}{\mathcal{N}_{M+1}}\Psi^{\mathcal{L}}_{M+1}.$
(90)

The normalization factors $\mathcal{N}_{M}$ are evaluated in Appendix C. The
states with two blocked levels are spanned by the orthonormal set of
$\mathcal{L}^{2}-\mathcal{L}=\mathcal{L}(\mathcal{L}-1)$ states (for $b\neq
b^{\prime}$)

$\phi_{bb^{\prime}}=c^{\dagger}_{b,\uparrow}c^{\dagger}_{b^{\prime},\downarrow}\Psi^{\mathcal{L},\mathcal{L}\backslash
b,b^{\prime}}_{M}.$ (91)

We thus have

$\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}}_{M}=C\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\phi_{bb^{\prime}},$ (92)

with a proportionality constant $C$ that needs to be determined. For that
purpose, we compute the norm:

$C^{2}=\frac{1}{\mathcal{L}^{2}}\sum_{b\neq b^{\prime},c\neq
c^{\prime}}\langle\Psi^{\mathcal{L}}_{M}|c_{c^{\prime}\downarrow}c_{c\uparrow}c_{b\uparrow}^{\dagger}c_{b^{\prime}\downarrow}^{\dagger}|\Psi^{\mathcal{L}}_{M}\rangle.$
(93)

Only the terms with $c=b$ and $c^{\prime}=b^{\prime}$ contribute, otherwise
the states are orthogonal. Now

$c_{b^{\prime}\downarrow}c_{b\uparrow}c_{b\uparrow}^{\dagger}c_{b^{\prime}\downarrow}^{\dagger}=1-\hat{n}_{b\uparrow}-\hat{n}_{b^{\prime}\downarrow}+\hat{n}_{b\uparrow}\hat{n}_{b^{\prime}\downarrow}.$
(94)

We need to handle the restriction to $b\neq b^{\prime}$ carefully, noting that
$\sum_{b\neq b^{\prime}}1=\mathcal{L}(\mathcal{L}-1)$, $\sum_{b\neq
b^{\prime}}\hat{n}_{b\sigma}=N_{\sigma}(\mathcal{L}-1)$, and

$\sum_{b\neq
b^{\prime}}\hat{n}_{b\uparrow}\hat{n}_{b^{\prime}\downarrow}=\sum_{b,b^{\prime}}\hat{n}_{b\uparrow}\hat{n}_{b^{\prime}\downarrow}-\sum_{b}\hat{n}_{b\uparrow}\hat{n}_{b\downarrow}=\hat{N}_{\uparrow}\hat{N}_{\downarrow}-\sum_{b}\hat{P}_{2,b}.$
(95)

The operator $\hat{P}_{2}=\sum_{b}\hat{P}_{2,b}$ counts the number of pairs in
the system, therefore $\langle\Psi_{M}|\hat{P}_{2}|\Psi_{M}\rangle=M$.
Collecting all terms we then find

$C^{2}=\frac{1}{\mathcal{L}^{2}}\langle\Psi^{\mathcal{L}}_{M}|\left(\mathcal{L}(\mathcal{L}-1)-(\mathcal{L}-1)\hat{N}_{\uparrow}-(\mathcal{L}-1)\hat{N}_{\downarrow}+\hat{N}_{\uparrow}\hat{N}_{\downarrow}-\hat{P}_{2}\right)|\Psi^{\mathcal{L}}_{M}\rangle=\frac{(\mathcal{L}-M)(\mathcal{L}-M-1)}{\mathcal{L}^{2}}.$
(96)

For $M=0$, the factor $C^{2}=1-\frac{1}{\mathcal{L}}$, which goes to 1 in the
large-$\mathcal{L}$ limit. In the opposite limit of $M=\mathcal{L}$, it is
zero. Finally, for half filling, it tends to $C^{2}=1/2$ in the
large-$\mathcal{L}$ limit.

We write $\psi_{2}$ in the form

$\psi_{2}=f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\Psi^{\mathcal{L}}_{M}=\frac{\sqrt{(1+M)(\mathcal{L}-M)}}{\mathcal{L}}\Psi^{\mathcal{L}}_{M+1}+\frac{\sqrt{(\mathcal{L}-M)(\mathcal{L}-M-1)}}{\mathcal{L}}\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\phi_{bb^{\prime}}.$ (97)

This expression is exact. There is an overall $\sqrt{\mathcal{L}-M}$ factor.
The relative importance of the first and the second term scale as $\sqrt{M}$
and $\sqrt{\mathcal{L}-M}$, respectively, up to $1/\mathcal{L}$ corrections.

For $M=\mathcal{L}/2$, we finally find

$\psi_{2}=f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\Psi^{\mathcal{L}}_{M}=\frac{1}{2}\sqrt{\frac{2+\mathcal{L}}{\mathcal{L}}}\Psi^{\mathcal{L}}_{M+1}+\frac{1}{2}\sqrt{\frac{\mathcal{L}-2}{\mathcal{L}-1}}\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}c^{\dagger}_{b\uparrow}c^{\dagger}_{b^{\prime}\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M}.$ (98)

Note that for large $\mathcal{L}$ both square roots tend towards 1. The states
with two blocked levels are orthogonal to $\Psi^{\mathcal{L}}_{M+1}$ states.
Hence

$\langle\psi_{2}|\psi_{2}\rangle=\frac{(1+M)(\mathcal{L}-M)}{\mathcal{L}^{2}}+\frac{(\mathcal{L}-M)(\mathcal{L}-M-1)}{\mathcal{L}^{2}}=\frac{\mathcal{L}-M}{\mathcal{L}}.$
(99)

For $M=\mathcal{L}/2$, this gives

$\langle\psi_{2}|\psi_{2}\rangle=\frac{1}{2}.$ (100)

The state
$\psi_{2}=f^{\dagger}_{\uparrow}f^{\dagger}_{\downarrow}\Psi^{\mathcal{L}}_{M}$
is thus not normalized to 1.

In an analogous way one calculates

$\psi^{\prime}_{2}=f_{\downarrow}f_{\uparrow}\Psi^{\mathcal{L}}_{M}=\frac{\sqrt{M(1+\mathcal{L}-M)}}{\mathcal{L}}\Psi^{\mathcal{L}}_{M-1}-\frac{\sqrt{M(M-1)}}{\mathcal{L}}\frac{1}{\sqrt{\mathcal{L}(\mathcal{L}-1)}}\sum_{b\neq
b^{\prime}}\phi^{\prime}_{bb^{\prime}},$ (101)

with
$\phi^{\prime}_{bb^{\prime}}=c^{\dagger}_{b,\uparrow}c^{\dagger}_{b^{\prime},\downarrow}\Psi^{\mathcal{L}\backslash
b,b^{\prime}}_{M-2}$. This expression is exact. At half filling, it simplifies
to

$\psi^{\prime}_{2}=f_{\downarrow}f_{\uparrow}\Psi^{\mathcal{L}}_{\mathcal{L}/2}=\frac{1}{2}\sqrt{\frac{2+\mathcal{L}}{\mathcal{L}}}\Psi^{\mathcal{L}}_{\mathcal{L}/2-1}-\frac{1}{2}\sqrt{\frac{\mathcal{L}-2}{\mathcal{L}-1}}\frac{1}{\mathcal{L}}\sum_{b\neq
b^{\prime}}\phi^{\prime}_{bb^{\prime}},$ (102)

with the same prefactors as for $\psi_{2}$.

##  Appendix C Evaluation of normalization $\mathcal{N}_{M}$

We fix $\mathcal{N}_{M}$ from the condition

$\langle\Psi_{M}|\Psi_{M}\rangle=\mathcal{N}_{M}^{2}\langle M|M\rangle=1.$
(103)

We write

$d_{M}=\langle M|M\rangle=\langle 0|B^{M}(B^{\dagger})^{M}|0\rangle,$ (104)

so that $\mathcal{N}_{M}=(d_{M})^{-1/2}$. Clearly $d_{1}=1$ and
$\mathcal{N}_{1}=1$.

Recall Eq. (77), $B^{\dagger}B|M\rangle=c_{M}|M\rangle$. This can also be
expressed as $B|M\rangle=c_{M}|M-1\rangle$. Thus we find

$\begin{split}d_{M}&=\langle M|M\rangle\\\ &=\langle M-1|B|M\rangle\\\
&=c_{M}\langle M-1|M-1\rangle\\\ &=c_{M}d_{M-1}.\end{split}$ (105)

The solution of this recursion equation is

$d_{M}=(-\mathcal{U})^{-M}(1)_{M}(-\mathcal{U})_{M},$ (106)

where $(a)_{n}=\Gamma(a+n)/\Gamma(a)$ is the Pochhammer symbol. At half
filling $\mathcal{U}=\mathcal{L}$ and $M=\mathcal{L}/2$, we find

$d_{\mathcal{L}/2}=\frac{\mathcal{L}!}{\sqrt{\mathcal{L}^{\mathcal{L}}}}.$
(107)

Furthermore, for $\mathcal{U}=\mathcal{L}$ one has

$\frac{\mathcal{N}^{2}_{M}}{\mathcal{N}^{2}_{M+1}}=\frac{d_{M+1}}{d_{M}}=\frac{(1+M)(\mathcal{L}-M)}{\mathcal{L}}.$
(108)

##  Appendix D Reference calculation

The method presented in Sec. IV is exact, but for clarity we have presented in
that section only simplified results obtained in the $\mathcal{L}\to\infty$
limit at half-filling. Below are the exact expressions for
$H^{D}_{\mathrm{eff}}$ and $H^{S}_{\mathrm{eff}}$.

$H^{D}_{\mathrm{eff}}=\begin{pmatrix}E^{D}_{1}&0&v\sqrt{\frac{\mathcal{L}-M}{\mathcal{L}}}&-v\sqrt{\frac{M}{\mathcal{L}}}\\\
0&E^{D^{\prime}}_{1}&v\sqrt{\frac{M}{\mathcal{L}}}&v\sqrt{\frac{\mathcal{L}-M}{\mathcal{L}}}\\\
v\sqrt{\frac{\mathcal{L}-M}{\mathcal{L}}}&v\sqrt{\frac{M}{\mathcal{L}}}&E^{D}_{0}&0\\\
-v\sqrt{\frac{M}{\mathcal{L}}}&v\sqrt{\frac{\mathcal{L}-M}{\mathcal{L}}}&0&E^{D}_{2}\end{pmatrix}.$
(109)

$H^{S}_{\mathrm{eff}}=\begin{pmatrix}E^{S}_{1}&v\sqrt{\frac{2(M+1)}{\mathcal{L}}}&v\sqrt{\frac{2(\mathcal{L}-M-1)}{\mathcal{L}}}&v\sqrt{\frac{2(\mathcal{L}-M)}{\mathcal{L}}}&v\sqrt{\frac{2M}{\mathcal{L}}}\\\
v\sqrt{\frac{2(M+1)}{\mathcal{L}}}&E^{S}_{0}&0&0&0\\\
v\sqrt{\frac{2(\mathcal{L}-M-1)}{\mathcal{L}}}&0&{E^{S}_{0}}^{\prime}&0&0\\\
v\sqrt{\frac{2(\mathcal{L}-M)}{\mathcal{L}}}&0&0&E^{S}_{2}&0\\\
v\sqrt{\frac{2M}{\mathcal{L}}}&0&0&0&{E^{S}_{2}}^{\prime}\\\ \end{pmatrix}.$
(110)

We provide a Mathematica notebook with these definitions, as well as a
corresponding input file for the DMRG solver, as Supplemental material [52].

## References

  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Microscopic theory of superconductivity, Phys. Rev. 106, 162 (1957).
  * Micnas  _et al._ [1990] R. Micnas, J. Ranninger, and S. Robaszkiewicz, Superconductivity in narrow-band systems with local non-retarded interactions, Rev. Mod. Phys. 62, 113 (1990).
  * Tinkham [2004] M. Tinkham,  _Introduction to superconductivity_ , 2nd ed. (Dover Publications, 2004).
  * Meservey  _et al._ [1970] R. Meservey, P. M. Tedorow, and P. Fulde, Magnetic field splitting of the quasiparticle states in superconducting aluminum films, Phys. Rev. Lett. 25, 1270 (1970).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. T. chi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Nakagawa  _et al._ [2021] Y. Nakagawa, Y. Kasahara, T. Nomoto, R. Arita, T. Nojima, and Y. Iwasa, Gate-controlled BCS-BEC crossover in a two-dimensional superconductor, Science 372, 190 (2021).
  * Bohr  _et al._ [1958] A. Bohr, B. R. Mottelson, and D. Pines, Possible analogy between the excitation spectra of nuclei and those of the superconducting metallic state, Physical Review 110, 936 (1958).
  * Migdal [1959] A. Migdal, Superfluidity and the moments of inertia of nuclei, Nuclear Physics 13, 655 (1959).
  * Belyaev [1961] S. T. Belyaev, On the nature of the first excited states of even-even spherical nuclei, JETP Letters 12, 968 (1961).
  * Richardson and Sherman [1964] R. W. Richardson and N. Sherman, Exact eigenstates of the pairing-force Hamiltonian, Nuclear Physics 52, 221 (1964).
  * Zelevinsky and Volya [2003] V. Zelevinsky and A. Volya, Nuclear pairing: New perspectives, Physics of Atomic Nuclei 66, 1781 (2003).
  * Dean and Hjorth-Jensen [2003] D. J. Dean and M. Hjorth-Jensen, Pairing in nuclear systems: from neutron stars to finite nuclei, Reviews of Modern Physics 75, 607 (2003).
  * Nozières and Schmitt-Rink [1985] P. Nozières and S. Schmitt-Rink, Bose condensation in an attractive fermion gas: From weak to strong coupling superconductivity, Journal of Low Temperature Physics 59, 195 (1985).
  * Jochim  _et al._ [2003] S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. Riedl, C. Chin, J. H. Denschlag, and R. Grimm, Bose-Einstein condensation of molecules, Science 302, 2101 (2003).
  * Greiner  _et al._ [2003] M. Greiner, C. A. Regal, and D. S. Jin, Emergence of a molecular Bose–Einstein condensate from a fermi gas, Nature 426, 537 (2003).
  * Zwierlein  _et al._ [2003] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, S. Gupta, Z. Hadzibabic, and W. Ketterle, Observation of Bose-Einstein condensation of molecules, Physical Review Letters 91, 250401 (2003).
  * Bartenstein  _et al._ [2004] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H. Denschlag, and R. Grimm, Collective excitations of a degenerate gas at the BEC-BCS crossover, Physical Review Letters 92, 203201 (2004).
  * Ketterle and Zwierlein [2008] W. Ketterle and M. W. Zwierlein, Making, probing and understanding ultracold fermi gases, in  _Ultracold Fermi Gases, Proceedings of the International School of Physics ”Enrico Fermi”, Course CLXIV, Varenna, 20 - 30 June 2006_ , edited by M. Inguscio, W. Ketterle, and C. Salomon (IOS Press, Amsterdam, 2008).
  * Averin and Nazarov [1992] D. V. Averin and Y. V. Nazarov, Single-electron charging of a superconducting island, Physical Review Letters 69, 1993 (1992).
  * Tuominen  _et al._ [1992] M. T. Tuominen, J. M. Hergenrother, T. S. Tighe, and M. Tinkham, Experimental evidence for parity-based 2e periodicity in a superconducting single-electron tunneling transistor, Physical Review Letters 69, 1997 (1992).
  * Eiles  _et al._ [1993] T. M. Eiles, J. M. Martinis, and M. H. Devoret, Even-odd asymmetry of a superconductor revealed by the Coulomb blockade of Andreev reflection, Physical Review Letters 70, 1862 (1993).
  * von Delft  _et al._ [1996] J. von Delft, A. D. Zaikin, D. S. Golubev, and W. Tichy, Parity-affected superconductivity in ultrasmall metallic grains, Physical Review Letters 77, 3189 (1996).
  * von Delft and Ralph [2001] J. von Delft and D. Ralph, Spectroscopy of discrete energy levels in ultrasmall metallic grains, Physics Reports 345, 61 (2001).
  * Franceschi  _et al._ [2010] S. D. Franceschi, L. Kouwenhoven, C. Schönenberger, and W. Wernsdorfer, Hybrid superconductor-quantum dot devices, Nat. Nanotechnology 5, 703 (2010).
  * Krogstrup  _et al._ [2015] P. Krogstrup, N. L. B. Ziino, W. Chang, S. M. Albrecht, M. H. Madsen, E. Johnson, J. Nygård, C. M. Marcus, and T. S. Jespersen, Epitaxy of semiconductor–superconductor nanowires, Nature Materials 14, 400 (2015).
  * Randeria and Taylor [2014] M. Randeria and E. Taylor, Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the unitary Fermi gas, Annual Review of Condensed Matter Physics 5, 209 (2014).
  * Yu [1965] L. Yu, Bound state in superconductors with paramagnetic impurities, Acta Phys. Sin. 21, 75 (1965).
  * Shiba [1968] H. Shiba, Classical spins in superconductors, Prog. Theor. Phys. 40, 435 (1968).
  * Rusinov [1969] A. I. Rusinov, Superconductivity near a paramagnetic impurity, JETP Lett. 9, 85 (1969), zh. Eksp. Teor. Fiz. Pisma Red. 9, 146 (1968).
  * Shiba [1973] H. Shiba, A Hartree-Fock theory of transition-metal impurities in a superconductor, Prog. Theor. Phys. 50, 50 (1973).
  * Sakurai [1970] A. Sakurai, Comments on superconductors with magnetic impurities, Prog. Theor. Phys. 44, 1472 (1970).
  * Matsuura [1977] T. Matsuura, The effects of impurities on superconductors with Kondo effect, Prog. Theor. Phys. 57, 1823 (1977).
  * Yazdani  _et al._ [1997] 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).
  * Balatsky  _et al._ [2006] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Impurity-induced states in conventional and unconventional superconductors, Rev. Mod. Phys. 78, 373 (2006).
  * Alloul  _et al._ [2009] H. Alloul, J. Bobroff, M. Gabay, and P. J. Hirschfeld, Defects in correlated metals and superconductors, Rev. Mod. Phys. 81, 45 (2009).
  * Franke  _et al._ [2011] K. J. Franke, G. Schulze, and J. I. Pascual, Competition of superconductivity phenomena and Kondo screening at the nanoscale, Science 332, 940 (2011).
  * Lee  _et al._ [2017] E. J. H. Lee, X. Jiang, R. Žitko, R. Aguado, C. M. Lieber, and S. De Franceschi, Scaling of subgap excitations in a superconductor-semiconductor nanowire quantum dot, Phys. Rev. B 95, 180502(R) (2017).
  * Meden [2019] V. Meden, The Anderson-Josephson quantum dot – a theory perspective, J. Phys.: Condens. Matter 31, 163001 (2019).
  * Richardson [1963] R. W. Richardson, A restricted class of exact eigenstates of the pairing-force Hamiltonian, Physics Letters 3, 277 (1963).
  * Richardson [1966] R. W. Richardson, Numerical study of the 8-32-particle eigenstates of the pairing Hamiltonian, Physical Review 141, 949 (1966).
  * von Delft and Braun [1999] J. von Delft and F. Braun, Superconductivity in ultrasmall grains: Introduction to Richardson’s exact solution, cond-mat/9911058 (1999).
  * Dukelsky  _et al._ [2004] J. Dukelsky, S. Pittel, and G. Sierra, Colloquium: Exactly solvable Richardson-Gaudin models for many-body quantum systems, Reviews of Modern Physics 76, 643 (2004).
  * Glazman and Catelani [2021] L. Glazman and G. Catelani, Bogoliubov quasiparticles in superconducting qubits, SciPost Phys. Lect. Notes , 31 (2021).
  * Ortiz and Dukelsky [2005] G. Ortiz and J. Dukelsky, BCS-to-BEC crossover from the exact BCS solution, Physical Review A 72 (2005).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Crosswhite  _et al._ [2008] G. M. Crosswhite, A. C. Doherty, and G. Vidal, Applying matrix product operators to model systems with long-range interactions, Phys. Rev. B 78, 035116 (2008).
  * Pavešić  _et al._ [2021] L. Pavešić, D. Bauernfeind, and R. Žitko, Subgap states in superconducting islands, Phys. Rev. B 104, L241409 (2021).
  * Pavešić and Žitko [2022] L. Pavešić and R. Žitko, Qubit based on spin-singlet Yu-Shina-Rusinov states, Phys. Rev. B 105, 075129 (2022).
  * Saldaña  _et al._ [2021] J. C. E. Saldaña, A. Vekris, L. Pavešić, P. Krogstrup, R. Žitko, K. Grove-Rasmussen, and J. Nygård, Excitations in a superconducting Coulombic energy gap, arXiv: 2101.10794 (2021), submitted.
  * [52] See Supplemental Material at [url will be inserted by publisher] for Mathematica notebooks with proofs of the mathematical statements made in this work, DMRG solver input file for the case of flat-band superconductor, as well as the corresponding Mathematica notebook with an exact calculation of the ground state properties.
  * Žitko  _et al._ [2022] R. Žitko, L. Pavešić, and D. Bauernfeind, QD-SI DMRG solver (2022).
  * Matthias  _et al._ [1958] B. T. Matthias, H. Suhl, and E. Curenzwit, Spin exchange in superconductors, Phys. Rev. Lett. 1, 92 (1958).
  * Anderson [1961] P. W. Anderson, Localized magnetic states in metals, Phys. Rev. 124, 41 (1961).
  * Soda  _et al._ [1967] T. Soda, T. Matsuura, and Y. Nagaoka, s-d exchange interaction in a superconductor, Prog. Theor. Phys. 38, 551 (1967).
  * Satori  _et al._ [1992] K. Satori, H. Shiba, O. Sakai, and Y. Shimizu, Numerical renormalization group study of magnetic impurities in superconductors, J. Phys. Soc. Japan 61, 3239 (1992).
  * Sakai  _et al._ [1993] O. Sakai, Y. Shimizu, H. Shiba, and K. Satori, Numerical renormalization group study of magnetic impurities in superconductors. II. Dynamical excitations spectra and spatial variation of the order parameter, J. Phys. Soc. Japan 62, 3181 (1993).
  * Yoshioka and Ohashi [1998] T. Yoshioka and Y. Ohashi, Ground state properties and localized excited states around a magnetic impurity described by the anisotropic $s$-$d$ interaction in superconductivity, J. Phys. Soc. Japan 67, 1332 (1998).
  * Yoshioka and Ohashi [2000] T. Yoshioka and Y. Ohashi, Numerical renormalization group studies on single impurity Anderson model in superconductivity: a unified treatment of magnetic, nonmagnetic impurities, and resonance scattering, J. Phys. Soc. Japan 69, 1812 (2000).
  * Bauer  _et al._ [2007] J. Bauer, A. Oguri, and A. C. Hewson, Spectral properties of locally correlated electrons in a Bardeen-Cooper-Schrieffer superconductor, J. Phys.: Condens. Matter 19, 486211 (2007).
  * Pillet  _et al._ [2013] J.-D. Pillet, P. Joyez, R. Žitko, and M. F. Goffman, Tunneling spectroscopy of a single quantum dot coupled to a superconductor: From Kondo ridge to Andreev bound states, Phys. Rev. B 88, 045101 (2013).
  * Saldaña  _et al._ [2018] J. C. E. Saldaña, A. Vekris, G. Steffensen, R. Žitko, P. Krogstrup, J. Paaske, and K. G.-R. J. Nygård, Supercurrent in a double quantum dot, Phys. Rev. Lett. 121, 257701 (2018).
  * Kezilebieke  _et al._ [2019] S. Kezilebieke, R. Žitko, M. Dvorak, T. Ojanen, and P. Liljeroth, Observation of coexistence of Yu-Shiba-Rusinov states and spin-flip excitations, Nano Letters 19, 4614 (2019).
  * Rubio-Verdú  _et al._ [2021] C. Rubio-Verdú, J. Zaldívar, R. Žitko, and J. I. Pascual, Coupled Yu-Shiba-Rusinov states induced by a many-body molecular spin on a superconductor, Physical Review Letters 126, 017001 (2021).
  * Luitz and Assaad [2010] D. J. Luitz and F. F. Assaad, Weak-coupling continuous-time quantum Monte Carlo study of the single impurity and periodic Anderson models with s-wave superconducting baths, Physical Review B 81, 024509 (2010).
  * Žonda  _et al._ [2015] M. Žonda, V. Pokorný, V. Janiš, and T. Novotný, Perturbation theory of a superconducting 0-$\pi$ impurity quantum phase transition, Scientific Reports 5, 8821 (2015).
  * Pokorný and Žonda [2018] V. Pokorný and M. Žonda, Correlation effects in superconducting quantum dot systems, Physica B: Condensed Matter 536, 488 (2018).
  * Domański  _et al._ [2017] T. Domański, M. Žonda, V. Pokorný, G. Górski, V. Janiš, and T. Novotný, Josephson-phase-controlled interplay between correlation effects and electron pairing in a three-terminal nanostructure, Physical Review B 95, 045104 (2017).
  * Zalom  _et al._ [2021] P. Zalom, V. Pokorný, and T. Novotný, Spectral and transport properties of a half-filled anderson impurity coupled to phase-biased superconducting and metallic leads, Physical Review B 103, 10.1103/physrevb.103.035419 (2021).
  * Rombouts  _et al._ [2006] S. M. A. Rombouts, K. Van Houcke, and L. Pollet, Loop updates for quantum Monte Carlo simulations in the canonical ensemble, Physical Review Letters 96, 180603 (2006).
  * Yuzbashyan  _et al._ [2003] E. A. Yuzbashyan, A. A. Baytin, and B. L. Altshuler, Strong-coupling expansion for the pairing Hamiltonian for small superconducting metallic grains, Physical Review B 68, 214509 (2003).
  * Khodel’ and Shaginyan [1990] V. A. Khodel’ and V. R. Shaginyan, Superfluidity in system with fermion condensate, JETP Letters 51, 488 (1990).
  * Volovik [2018] G. E. Volovik, Graphite, graphene, and the flat band superconductivity, JETP Letters 107, 516 (2018).
  * Braun and von Delft [1999] F. Braun and J. von Delft, Superconductivity in ultrasmall metallic grains, Physical Review B 59, 9527 (1999).
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Relation between the Anderson and Kondo hamiltonians, Phys. Rev. 149, 491 (1966).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–Wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).
  * Hecht  _et al._ [2008] T. Hecht, A. Weichselbaum, J. von Delft, and R. Bulla, Numerical renormalization group calculation of near-gap peaks in spectral functions of the Anderson model with superconducting leads, J. Phys. Condens. Mat. 20, 275213 (2008).
  * Affleck  _et al._ [2000] I. Affleck, J.-S. Caux, and A. M. Zagoskin, Andreev scattering and josephson current in a one-dimensional electron liquid, Physical Review B 62, 1433 (2000).
  * Vecino  _et al._ [2003] E. Vecino, A. Martín-Rodero, and A. Yeyati, Josephson current through a correlated quantum level: Andreev states and pi junction behavior, Phys. Rev. B 68, 035105 (2003).
  * Bergeret  _et al._ [2007] F. Bergeret, A. Yeyati, and A. Martín-Rodero, Josephson effect through a quantum dot array, Phys. Rev. B 76, 174510 (2007).
  * Grove-Rasmussen  _et al._ [2018] K. Grove-Rasmussen, G. Steffensen, A. Jellinggaard, M. H. Madsen, R. Žitko, J. Paaske, and J. Nygård, Yu-Shiba-Rusinov screening of spins in double quantum dots, Nat. Commun. 9, 2376 (2018).
  * Yosida [1966] K. Yosida, Bound state due to the s-d exchange interaction, Phys. Rev. 147, 223 (1966).
  * Žitko and Fabrizio [2017] R. Žitko and M. Fabrizio, Non-Fermi-liquid behavior in quantum impurity models with superconducting channels, Physical Review B 95, 10.1103/physrevb.95.085121 (2017).
  * von Oppen and Franke [2021] F. von Oppen and K. J. Franke, Yu-Shiba-Rusinov states in real metals, Physical Review B 103, 10.1103/physrevb.103.205424 (2021).
  * Schmid  _et al._ [2022] H. Schmid, J. F. Steiner, K. J. Franke, and F. von Oppen, Quantum Yu-Shiba-Rusinov dimers, arXiv:2203.15011 (2022).