# Exchange interaction between two quantum dots coupled through a
superconducting island

Ádám Bácsi<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39, Ljubljana
SI-1000, Slovenia  Department of Mathematics and Computational Sciences,
Széchenyi István University, 9026 Győr, Hungary  Luka Pavešić  Jožef Stefan
Institute, Jamova 39, Ljubljana SI-1000, Slovenia  Faculty of Mathematics and
Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
Rok Žitko  Jožef Stefan Institute, Jamova 39, Ljubljana SI-1000, Slovenia
Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19,
SI-1000 Ljubljana, Slovenia

###### Abstract

We present a theoretical study of a system consisting of a superconducting
island and two quantum dots, a possible platform for building qubits and
Cooper pair splitters, in the regime where each dot hosts a single electron
and, hence, carries a magnetic moment. We focus on the case where the dots are
coupled to overlapping superconductor states and we study whether the spins
are ferromagnetically or antiferromagnetically aligned. We show that if the
total number of particles is even, the spins align antiferromagnetically in
the flatband limit, i.e., when the bandwidth of the superconductor is
negligibly small, but ferromagnetically if the bandwidth is finite and above
some value. If the total number of particles is odd, the alignment is
ferromagnetic independently from the bandwidth. This implies that the results
of the flatband limit are applicable only within restricted parameter regime
for realistic superconducting qubit systems and that some care is required in
applying simplified models to devices such as Cooper pair splitters.

##  I Introduction

Superconducting devices are attracting strong interest due to their
applications in quantum computing [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11].
Building upon the ideas of semiconductor spin qubits [12, 13], the combination
of spin carrying quantum dots and superconducting islands is also a possible
qubit realization [14, 15, 16, 17, 18] which remains to be explored more in
depth.Such devices are particularly promising due to the presence of tunable
long-lived subgap states. They could find use in a wide range of applications,
such as Cooper-pair splitters [19], for linking spatially separated
superconducting qubits [20, 21, 22], or building Andreev molecules [23].

Quantum dots (QDs) can be designed to possess a single relevant electron state
well separated from other levels [24]. In this case, the QD may be modeled as
a single Anderson impurity with two characteristic energy scales: the energy
of the QD level, $\varepsilon$, and the on-site interaction energy, $U$. If
$U$ is large and $\varepsilon$ is deep below the Fermi level, the QD hosts a
single electron and carries a spin just like magnetic impurities.

Two magnetic impurities in the same host material can be coupled in different
ways. In a metallic host, they couple via the Ruderman-Kittel-Kasuya-Yosida
(RKKY) interaction which, depending on the inter-impurity distance, may result
in a ferromagnetic (FM) or antiferromagnetic (AFM) spin alignment. For very
short distances, when the distance is on the order of the lattice constant,
the atomistic details play a crucial role and the coupling is rather described
as superexchange with the coupling sign determined by the electron filling of
the intermediary atoms. Superexchange is also the essential mechanism
governing the interaction between non-neighboring QDs of a chain [25].

In the case of superconducting mediator, the RKKY interaction is shifted
toward AFM alignment [26] due to additional processes involving the subgap
states. The effect of superconductivity on superexchange, however, has not
been studied much so far.

In this paper, we study the system with two QDs coupled to a superconducting
island (SI) as sketched in Fig. 1. Our goal is to elucidate the sign and the
strength of the inter-dot exchange coupling in the case where the two dots are
coupled to the same superconductor state: this corresponds to the limit of
small separation between the dots.

Figure 1: Schematic of a superconducting island bridging two quantum dots.
Each dot carries a spin.

##  II Model of the superconductor

Conventional superconductors are successfully described by the BCS mean-field
theory [27]. An alternative approach to superconductivity consists of solving
the pairing model without performing the mean-field decoupling in the
canonical ensemble. The simplified version of the pairing Hamiltonian is also
known as the Richardson or the picket-fence model [28, 29, 30]. In the
original work of Richardson and Sherman [31], several energy levels
$\epsilon_{i}$ are taken into account which, in the context of solid state
physics, can be thought of as the states in the vicinity of the Fermi energy
that participate in the attractive pairing interaction. The pairing
interaction consists of quartic terms [see Eq. (1)] leading to the formation
of Cooper pairs.

The Richardson model is exactly solvable for an isolated superconductor (SC).
However, if the SC is coupled to QDs, exact analytical solutions are available
only in the flatband limit. In this limit, the energy levels are all
degenerate, $\epsilon_{i}\equiv\text{const}$. Despite of this simplification,
the flatband limit was proven to provide a proper description for systems of
one QD coupled to the SI [16]. In this paper, we extend these studies to the
case of two QDs coupled to the island.

Let us consider a system with $N$ fermionic single-particle levels and we
write the pairing Hamiltonian of the Richardson model as

$\displaystyle
H_{SC}=-\frac{g}{N}\sum_{i,j}^{N}c_{i\uparrow}^{+}c_{i\downarrow}^{+}c_{j\downarrow}c_{j\uparrow},$
(1)

where $c_{i\sigma}$ is the fermionic annihilation operator corresponding to
the energy level $i$ and spin $\sigma$. In contrast to the BCS Hamiltonian,
the number of fermions is conserved. For $N$ superconducting levels and $M$
Cooper pairs, the ground state is given by
$|\Psi_{M}^{N}\rangle=\sum_{conf}|\mathrm{conf}\rangle/\sqrt{\mathcal{N}_{M}^{N}}$
which is the completely symmetric combination of all possible configurations
of $M$ Cooper pairs occupying the $N$ levels. The number of configurations is
$\mathcal{N}_{M}^{N}=\binom{N}{M}$. The ground state (GS) energy is expressed
as[30, 16]

$\displaystyle E_{0}=-\frac{g}{N}M(N+1-M).$ (2)

A great advantage of the flatband limit is that the model is invariant under
any unitary transformations as long as the opposite spin levels are
transformed with complex conjugate matrices. To be more precise, for any
unitary matrix $U_{ij}$, the Hamiltonian is invariant under

$\begin{split}f_{j\uparrow}&=\sum_{i}U_{ji}c_{i\uparrow},\\\
f_{j\downarrow}&=\sum_{i}U_{ji}^{*}c_{i\downarrow}.\end{split}$ (3)

The physical significance of the complex conjugation is that it maintains the
time-reversal invariance of the Hamiltonian.

At this point, we have complete freedom how to choose the unitary
transformation. We will use this freedom to transform into a basis in which
the hybridization between the dots and the superconductor has a convenient
form. In general, the hybridization operators are defined as

$\displaystyle H_{hyb}=H_{hyb,L}+H_{hyb,R},$ (4)

where

$\displaystyle
H_{hyb,L}=v_{L}\sum_{i}\left(\gamma_{Li}d_{L\uparrow}^{+}c_{i\uparrow}+\gamma_{Li}^{*}d_{L\downarrow}^{+}c_{i\downarrow}+h.c.\right)$
(5)

describes the coupling between the left dot and the superconductor.
$d_{L\sigma}$ is the fermionic annihilation operator corresponding to the left
dot level. The strength of the coupling is determined by $v_{L}$ while the
coefficients $\gamma_{Li}$ describe the distribution of the couplings obeying
$\sum_{i}|\gamma_{Li}|^{2}=1$ (i.e., they define the wave function of the
superconductor state that hybridizes with the QD level). Note that for
different spin orientations, we use complex conjugate coupling constants as a
consequence of the time-reversal symmetry. The coupling between the right dot
and the superconductor is given by the same formula as (5), but with $v_{R}$
and $\gamma_{Ri}$. The coefficients are also normalized as
$\sum_{i}|\gamma_{Ri}|^{2}=1$.

The overlap between the distribution of $\gamma_{Li}$ and $\gamma_{Ri}$ is
given by [18]

$\alpha=\sum_{i}\gamma_{Li}^{*}\gamma_{Ri}$ (6)

whose absolute value ranges from 0 to 1. Within the framework of the
Richardson model, $\alpha$ mimics the distance between the dots. Namely,
$\alpha=0$ means that there is no overlap between the distributions and the
dots are far from each other. On the other hand, $|\alpha|=1$ implies that the
dots are coupled to the same superconducting state and are, hence, close to
each other. In the present paper, we focus on the latter limit, i.e., when the
distributions $\gamma_{Li}$ and $\gamma_{Ri}$ are basically the same apart
from an overall phase factor. This limit is valid as long as the Yu-Shiba-
Rusinov states of the dots have a large overlap. Recent experiments show that
the spatial extension of these states can reach even the range of 50 nm [32].
We also note that $\alpha$ fully characterizes the relation between the
coefficients $\gamma_{Li}$ and $\gamma_{Ri}$ in the flatband limit, while for
finite bandwidth their detailed values (energy dependence) also matter.

As shown in the Appendix A, there exists a unitary transformation such that
$f_{0\uparrow}=\sum_{i}\gamma_{L/Ri}c_{i\uparrow}$ is one of the new basis
states which forms an orthonormal system together with the other new basis
states. The state corresponding to $f_{0\sigma}$ will be referred to as the
distinguished state of the superconductor. The unitary transformation leaves
the pairing Hamiltonian invariant in form so that

$\displaystyle
H_{SC}=-\frac{g}{N}\sum_{ij}f_{i\uparrow}^{+}f_{i\downarrow}^{+}f_{j\downarrow}f_{j\uparrow}$
(7)

due to the flatband limit. In the presence of a non-zero bandwidth, the
superconductor levels would have non-uniform $\epsilon_{i}$ which, under the
unitary transformation, would generate inter-level hopping terms.

With the distinguished level $f_{0\sigma}$, let us define the many-body
superconductor states

$\displaystyle|\Psi_{0}\rangle=|0,\Psi_{M}^{N-1}\rangle\quad\mbox{and}\quad|\Psi_{2}\rangle=|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,$
(8)

where $\Psi_{0}$ ($\Psi_{2}$) describes the state where the $f_{0}$-level is
empty (occupied by a Cooper pair) and the remaining $N-1$ superconductor
levels are filled with $M$ (with $M-1$) Cooper pairs in a completely symmetric
fashion. The pairing Hamiltonian is represented on the basis of these states
as

$\displaystyle\mathbf{H}=E_{0}\mathbf{I}+\mathbf{H}_{SC},$
$\displaystyle\mathbf{H}_{SC}=\left[\begin{array}[]{cc}g\nu&\beta\\\
\beta&g(1-\nu)\end{array}\right],$ (11)

where

$\nu=M/N$ (12)

is the filling factor and

$\beta=-g\sqrt{\nu(1-\nu)}.$ (13)

The eigenvalues of the matrix are $0$ and $g$ corresponding to the GS and the
excited state (ES), respectively. This indicates that the gap of the model is
related to the pairing energy as $2\Delta=g$. Of course, there are many other
states with the excitation energy of $g$, but this is the only one mixing with
the GS in the presence of hybridization. The minimal model of the
superconductor, Eq. (11), is sufficient to study the coupling of two dots with
the overlap parameter $\alpha=1$ in the flatband limit. Note that if
$\alpha\neq 1$, there would exist two distinguished levels (see Appendix A)
and instead of Eq. (8) we would need to define a more complex set of many-body
basis states.

##  III Two quantum dots coupled to the superconductor

We study the effect of hybridization on the electronic structure in the case
of $\alpha=1$, when the two dots are coupled to the same superconductor level.
Furthermore, it is assumed that the system has parity symmetry, i.e.,
$v_{L}=v_{R}=v$, $U_{L}=U_{R}=U$ and
$\varepsilon_{L}=\varepsilon_{R}=\varepsilon$.

The total Hamiltonian of the system is expressed as

$\displaystyle H=H_{SC}+H_{QD}+H_{hyb},$ (14)

where $H_{SC}$ is given by Eq. (7),

$\displaystyle
H_{QD}=\varepsilon\sum_{\sigma}\left(\hat{n}_{L\sigma}+\hat{n}_{R\sigma}\right)+U\left(\hat{n}_{L\uparrow}\hat{n}_{L\downarrow}+\hat{n}_{R\uparrow}\hat{n}_{R\downarrow}\right),$
(15)

and the coupling between the dots and the superconductor is determined by

$\displaystyle
H_{hyb}=\sum_{\sigma}\left[v\left(d_{L\sigma}^{+}+d_{R\sigma}^{+}\right)f_{0\sigma}+h.c.\right],$
(16)

where $d_{L\sigma}$ ($d_{R\sigma}$) is the annihilation operator of the left
(right) quantum dot.

The total Hamiltonian conserves both the number of particles and the total
spin. In the followings, we determine the lowest-energy states in different
charge and spin sectors by using analytical and numerical methods.

In the analytical calculation, we apply the approximation that the on-site
interaction $U$ is infinitely large and, hence, none of the dots can be doubly
occupied. Furthermore, we take the thermodynamic limit of macroscopically
large number of states in the SI.

Numerical calculations are performed with Density Matrix Renormalization Group
(DMRG) technique [33, 34] in which $U$ has a finite value and the system is
finite, $N=80$. For this purpose, we extended the approach described in Refs.
35, 36 and implemented a matrix product operator expression of the Hamiltonian
for the QD-SI-QD system; see Appendix B for details.

We start the presentation of the analytical model by introducing four
different sectors which are not coupled to each other because of the
conservation of particle number and spin. For even number of particles, the
total spin is either $0$ (singlet) or $1$ (triplet) since the superconductor
is an overall singlet, but the electrons hosted by the dots can have arbitrary
spin orientations. For odd number of particles, the total spin of the system
is always $1/2$ (doublet), but based on the parity symmetry of states, a
symmetric and an anti-symmetric sector can be identified.

###  III.1 Analytical results for even number of particles

In the sector of even number of particles and $S=0$, the relevant states are
given by

$\begin{split}|\phi_{0}^{S}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|0,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{2}^{S}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{0X}^{S}\rangle&=\frac{1}{2}\Big{(}(|\uparrow,0\rangle+|0,\uparrow\rangle)\otimes|\downarrow,\Psi_{M}^{N-1}\rangle\\\
&\quad-(|0,\downarrow\rangle+|\downarrow,0\rangle)\otimes|\uparrow,\Psi_{M}^{N-1}\rangle\Big{)},\\\
|\phi_{0D}^{S}\rangle&=|0,0\rangle\otimes|\uparrow\downarrow,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{-2D}^{S}\rangle&=|0,0\rangle\otimes|0,\Psi_{M+1}^{N-1}\rangle.\end{split}$
(17)

The first part always encodes the state of the dots while the second part
after $\otimes$ describes the state of the superconductor. It is also worth to
note that all states in this sector are symmetric under the spatial reflection
$L\leftrightarrow R$. On this basis, the Hamiltonian is represented by the
block matrix

$\displaystyle\mathbf{H}_{S}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c|c}2\varepsilon+\mathbf{H}_{SC}&\begin{array}[]{c}\sqrt{2}v\\\
0\end{array}&\text{\large 0}\\\
\hline\cr\begin{array}[]{cc}\sqrt{2}v&0\end{array}&\varepsilon+g\nu&\begin{array}[]{cc}2v&0\end{array}\\\
\hline\cr\text{\large 0}&\begin{array}[]{c}2v\\\
0\end{array}&\mathbf{H}_{SC}+g(2\nu-1)\end{array}\right],$ (27)

where $\mathbf{H}_{SC}$ is the $2\times 2$ matrix describing the internal
structure of the superconductor as defined in Eq. (11). We note that for
finite system size, the matrix elements of the Hamiltonian would include
further terms proportional to $1/N$ which, however, are neglected in the
thermodynamic limit that we take here.

In the S=1 sector, we consider the $S_{z}=1$ states, other $S_{z}$ subspaces
being equivalent. The basis vectors are given as

$\begin{split}|\phi_{0}^{T}\rangle&=|\uparrow,\uparrow\rangle\otimes|0,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{2}^{T}\rangle&=|\uparrow,\uparrow\rangle\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{0X}^{T}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle-|0,\uparrow\rangle\right)\otimes|\uparrow,\Psi_{M}^{N-1}\rangle,\end{split}$
(28)

which are all parity anti-symmetric states. On this subspace, the Hamiltonian
is represented by

$\displaystyle\mathbf{H}_{T}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c}2\varepsilon+\mathbf{H}_{SC}&\begin{array}[]{c}\sqrt{2}v\\\
0\end{array}\\\
\hline\cr\begin{array}[]{cc}\sqrt{2}v&0\end{array}&\varepsilon+g\nu\end{array}\right]\,.$
(34)

To compare the two sectors with an even number of particles, we compute the
ground state of the Hamiltonian in each sector as a function of $\varepsilon$.
The results shown in Fig. 2 indicate that the $S=0$ sector always has lower
minimal energy. In fact, one can prove using perturbation theory that for
$v\ll|\varepsilon|$ the difference of minimal energies

$E_{T}-E_{S}=\frac{2v^{4}}{|\varepsilon|^{3}}$ (35)

is always positive. The scaling of $v^{4}$ implies that the energy difference
stems from a fourth order process which is characteristic to both
superexchange and to RKKY interaction. The fourth order processes involve
depopulation of the dots in the singlet sector through the states of
$\psi_{0D}^{S}$ and $\phi_{-2D}^{S}$. However, in the triplet sector, the dots
cannot be completely depopulated due to Pauli principle and because single
particles cannot move from the distinguished level to other SC levels in the
flatband limit.

Figure 2: Comparison of DMRG (denoted by symbols $\times$ and
$\blacktriangle$) and analytical results (solid and dashed lines). Different
colors correspond to different hybridization strengths, $v/g=0.2$, $v/g=0.5$
and $v/g=1$.

###  III.2 Analytical results for odd number of particles

If the total number of particles is odd, the system has an overall spin of
$S=1/2$ describing a doublet state. Here, we take $S_{z}=1/2$. This sector can
be split into a symmetric and an anti-symmetric subsector based on parity
symmetry.

In the symmetric subsector, the relevant states are given by

$\begin{split}|\phi_{1}^{Dsym}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|\uparrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{0}^{Dsym}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle+|0,\uparrow\rangle\right)\otimes|0,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{2}^{Dsym}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle+|0,\uparrow\rangle\right)\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{-1}^{Dsym}\rangle&=|0,0\rangle\otimes|\uparrow,\Psi_{M}^{N-1}\rangle,\end{split}$
(36)

and the Hamiltonian is represented by

$\displaystyle\mathbf{H}^{sym}_{D}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c|c}2\varepsilon+g(1-\nu)&\begin{array}[]{cc}0&-v\end{array}&0\\\
\hline\cr\begin{array}[]{c}0\\\
-v\end{array}&\varepsilon+\mathbf{H}_{SC}&\begin{array}[]{c}\sqrt{2}v\\\
0\end{array}\\\ \hline\cr
0&\begin{array}[]{cc}\sqrt{2}v&0\end{array}&g\nu\end{array}\right]\,.$ (46)

In the anti-symmetric subsector, the relevant states are

$\begin{split}|\phi_{1}^{Dasym}\rangle&=\frac{1}{\sqrt{6}}\left(2|\uparrow,\uparrow\rangle\otimes|\downarrow,\Psi_{M-1}^{N-1}\rangle\right.-\\\
&-\left.\left(|\uparrow,\downarrow\rangle+|\downarrow,\uparrow\rangle\right)\otimes|\uparrow,\Psi_{M-1}^{N-1}\rangle\right),\\\
|\phi_{0}^{Dasym}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle-|0,\uparrow\rangle\right)\otimes|0,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{2}^{Dasym}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle-|0,\uparrow\rangle\right)\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\end{split}$
(47)

and the Hamiltonian is given by

$\displaystyle\mathbf{H}^{asym}_{D}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c}2\varepsilon+g(1-\nu)&\begin{array}[]{cc}0&\sqrt{3}v\end{array}\\\
\hline\cr\begin{array}[]{c}0\\\
\sqrt{3}v\end{array}&\varepsilon+\mathbf{H}_{SC}\end{array}\right]\,.$ (53)

The minimal energies are numerically obtained and shown in Fig. 2. The results
show that the anti-symmetric sector has lower minimal energies which is
further confirmed by the perturbative expansion for $|\varepsilon|\gg v$
leading to

$E_{D}^{sym}-E_{D}^{asym}=\frac{2v^{2}}{|\varepsilon|}.$ (54)

Note that the anti-symmetric sector, particularly the state
$|\phi_{1}\rangle$, features FM alignment of the spins sitting on the dots
even though the total system is a spin-1/2 state. The FM coupling between the
spins differs essentially from what has been found in the case of even number
of particles, where the spins show AFM alignment in the singlet sector.

###  III.3 Discussion

The difference in behavior for odd and even electron number can be understood
by recalling the properties of superexchange [37, 38]. If the superexchange
between two magnetic moments is realized through a single orbital level, it
can be shown that the nature of coupling depends on the occupancy of the
mediating orbital, see Refs. [37, 39] and the Appendix C. If it is singly
occupied, the magnetic moments are preferentially FM aligned, while for empty
or doubly occupied mid-state, the coupling is AFM. For an odd total number of
particles, the energetically most favourable configuration is that the extra
particle occupies the distinguished level of the superconductor ensuring
single occupancy in the orbital that bridges both quantum dot levels.

For even number of particles, however, all electrons of the superconductor
tend to form Cooper pairs and, hence, the distinguished level is more likely
to be empty or doubly occupied. We note that the probability of having a
single particle on the distinguished level is zero only in the strict flatband
limit. For finite bandwidth, single particle hopping takes place among the SC
levels which indicates that the probability of single occupancy is not zero.
This effect is further analyzed in Sec. IV.

###  III.4 Comparison with DMRG results

We compare the analytical results with the results of a DMRG calculation
carried out with $N=80$ superconducting levels and $n=82$, 83 particles. The
numerical simulation has been performed with a uniform distribution of
couplings, $\gamma_{Li}=\gamma_{Ri}=1/\sqrt{N}$ which corresponds to
$\alpha=1$ as explained after Eq. 5. The on-site interaction has been set to
$U/g=40$.

The comparison in Fig. 2 confirms that the analytical approach produces
correct results as regards the spin configuration. Namely, numerical
simulations also indicate AFM coupling for even and FM coupling for odd number
of particles. However, some quantitative differences are observed. These are
the consequence of finite $U$ and of the finite system size in DMRG. The dot
energy level has been swept between $-20g$ and $-5g$ to cover the regime where
the probability of single occupancy of the QDs is high. We do not study
$\varepsilon<-20g$ because charge fluctuations through doubly occupied states,
which are not included in our analytical model, would play a more important
role. The fingerprints of these charge fluctuations are captured by the DMRG
close to $\varepsilon=-U/2=-20g$ but these have no effect on the spin
configuration.

##  IV Effects of finite bandwidth

The flatband limit has been shown to represent an adequate model for the
situation of a single QD coupled to the SI [16]. There, the effects of finite
bandwidth lead to some quantitative modifications but no change in the general
structure has been observed.

We now investigate whether the nature of spin couplings of the double dot
system uncovered in the flatband limit is preserved for finite bandwidth.
Since the Richardson model with finite bandwidth and attached QDs cannot be
handled analytically, we mostly rely on DMRG calculations (below we will also
introduce a minimal analytical model which qualitatively describes the main
new feature). For simplicity, we will consider the case of a uniform
distribution of energy levels $\epsilon_{i}$ in the SI spanning the interval
$[-W:W]$. The numerical simulations are performed with several hybridization
strengths, $v/g=0.2$, $v/g=0.5$ and $v/g=1$. The QD parameters are $U/g=40$
and $\varepsilon/g=-10$.

We find that for an odd number of particles the finite bandwidth does not
change the spin configuration in the ground state. In the studied regime, the
spins are always FM aligned.

However, for the even number sector, we do find a difference: whereas the
flatband limit exhibits AFM alignment, we observe FM alignment at large enough
values of the bandwidth. Fig. 3 shows that above some limit value of the
bandwidth, which is in the range of $g$, the triplet sector has lower energy.

Figure 3: Effects of finite bandwidth for even number of particles at
different values of hybridization strength, DMRG results. The energy
difference $E_{T}-E_{S}$ (solid lines, left vertical axis) changes sign with
increasing bandwidth indicating that the ground state becomes a triplet state.
This feature is further verified by calculating the spin-spin correlation
between the dots (dashed lines, right vertical axis).

We now present an analytical model which recovers this feature. The finite
bandwidth is the consequence of the kinetic energy of the particles in the
superconductor. In the original basis determined by the operators $c_{i}$, the
kinetic energy is represented by a diagonal matrix with the entries of
$\epsilon_{i}$. After applying the unitary transformation $c_{i}\rightarrow
f_{j}$ transforming the hybridization into a more convenient form, the kinetic
energy term becomes non-diagonal. In principle all-to-all couplings could
appear but we only consider those in which the distinguished state $f_{0}$
participates. Therefore, the kinetic energy term is written as

$\displaystyle H_{kin}=\frac{t}{\sqrt{N}}\sum_{i\neq
0,\sigma}\left(f_{i\sigma}^{+}f_{0\sigma}+h.c.\right)\,.$ (55)

This hopping process breaks Cooper pairs inside the superconductor. Therefore,
the basis set of $|\Psi_{0}\rangle$ and $|\Psi_{2}\rangle$ as defined in Eq.
(8) must be complemented with the Cooper-pair-broken state

$\displaystyle|\Psi_{1}\rangle=\frac{1}{\sqrt{2}}\Big{(}|\uparrow,\Psi_{M-1,+\downarrow}^{N-1}\rangle-|\downarrow,\Psi_{M-1,+\uparrow}^{N-1}\rangle\Big{)},$
(56)

where

$\displaystyle|\uparrow,\Psi_{M-1,+\downarrow}^{N-1}\rangle=\frac{1}{\sqrt{N-1}}\sum_{i=1}^{N-1}f_{0\uparrow}^{+}f_{i\downarrow}^{+}|0,\Psi_{M-1}^{N-1}\rangle$
(57)

describes a state with a single particle occupying the distinguished site and
its partner residing deeper inside the superconductor. Note that this state
characterizes one broken Cooper pair. Due to the coupling to the quantum dots,
states with additional broken Cooper pairs could also occur. Since breaking
further pairs always costs an energy proportional to $g$, we neglect states
with more than one broken pair.

On the subspace spanned by $|\Psi_{0}\rangle$, $|\Psi_{1}\rangle$ and
$|\Psi_{2}\rangle$, the Hamiltonian of the superconductor in isolation is now
represented by

$\displaystyle\mathbf{H}_{SC+kin}=\left[\begin{array}[]{ccc}g\nu&\sqrt{2\nu}\,t&\beta\\\
\sqrt{2\nu}\,t&g&\sqrt{2(1-\nu)}\,t\\\
\beta&\sqrt{2(1-\nu)}\,t&g(1-\nu)\end{array}\right],$ (61)

where $\beta=-g\sqrt{\nu(1-\nu)}$ as it was defined in Eq. (13). It is
instructive to study the probability of single occupancy of the distinguished
level as a function of $t$. If $t=0$, the eigenvalues are $0$, $g$ and $g$;
the lowest energy of 0 corresponds to a mixture of $|\Psi_{0}\rangle$ and
$|\Psi_{2}\rangle$ indicating that the probability of single occupancy is
zero. If the hopping parameter has a small but non-zero value, the
$|\Psi_{1}\rangle$ state is admixed to the ground state, indicating that the
ground state contains broken Cooper pairs. Second order perturbation theory
indicates that the ground state energy is shifted to $-2t^{2}/g$ in the
presence of a small $t$. In the opposite limit, $t\rightarrow\infty$, the
eigenvalues are $-\sqrt{2}t$, $0$ and $\sqrt{2}t$. The ground state is
$\big{(}|\Psi_{0}\rangle-\sqrt{2}|\Psi_{1}\rangle+|\Psi_{2}\rangle\big{)}/2$.
In this state, the probability of having a single particle on the
distinguished level is $1/2$. Between the two limiting cases, the probability
increases smoothly from $0$ to $1/2$.

It should be emphasized that in reality singly occupied levels do not appear
in the GS of superconductors in isolation. Blocked levels in the Richardson
model are only present in the excited states [28] or due to the Cooper-pair
breaking action of magnetic impurities. The results are thus physically
meaningful only when the QD is coupled to the superconductor and the local
magnetic moment breaks pairs even in the GS.

Taking also the dots into account, one has to add further states with one
broken Cooper pair to all sectors investigated in the flatband limit. The
calculation of the Hamiltonian on these extended Hilbert spaces is tedious but
straightforward and can be found in the Appendix D.

The minimal energies of the obtained Hamiltonians are computed numerically and
the difference between the minimal energies are shown in Fig. 4. The results
show that for small values of $t$, the ground state forms a singlet, while by
increasing the strength of the kinetic energy, the ground state becomes a
triplet. This is further confirmed by calculating the spin-spin correlation
between the two quantum dots which suddenly switches from approximately -3/4
to close to 1/4 as the levels of singlet and triplet states cross each other.

Figure 4: Comparison of minimal energies in the singlet ($E_{S}$) and the
triplet ($E_{T}$) sector at different values of hybridization strength,
analytical results. Solid lines illustrate the energy difference while dashed
lines show the spin-spin correlation between the dots in the ground state. By
exceeding some limit value of the bandwidth, the ground state changes from
singlet to triplet. The transition value is around $0.3g$ and depends only
slightly on the hybridization strength.

The minimal model for even number of particles implies that if the hopping
processes between superconducting energy levels are large enough, i.e., the
kinetic energy term is significant, then it is energetically favourable to
break a Cooper pair inside the superconductor and allow the dot spins to form
a triplet through the mediation of one member of the broken pair. This is also
the ground state for a normal-state bath in the double quantum dot impurity
problem [40, 41, 42, 43, 44, 45, 46].

##  V Charging energy

In the case of a small SI with a large charging energy $E_{C}$, we add
$E_{C}(\hat{N}_{SC}-n_{0})^{2}$ to the Hamiltonian, where $\hat{N}_{SC}$ is
the SI occupancy operator and $n_{0}$ is the target number of particles
controlled in experimental setups by tuning the gate voltage applied on the
SI.

We study by using DMRG whether the spin configuration of the dots changes when
tuning $n_{0}$. The parameters were chosen to fit typical experimental setups,
$U=6\Delta$, $\varepsilon=-3\Delta$, $E_{C}=2\Delta$. We use $g=0.8W$, which
corresponds to $\Delta=0.16W$. For the hybridization, $v=0.25W$ is set which
is equivalent to $\Gamma/U=0.102$ with $\Gamma=\pi v^{2}/(2W)$ the
hybridization strength. The number of sites is $N=80$ and the the number of
total particles are $n=82$, 83 and 84. Fig. 5 shows the $n_{0}$-dependence of
the energy of both the ground and the first excited state.For odd number of
particles, the solid line indicates the parity anti-symmetric sector with FM
alignment while the excited state is the symmetric, AFM aligned state. It is
readily seen that the charging energy only shifts these energy states.

Figure 5: Energy of the ground and first excited state as a function of
$n_{0}$ for states with a different total number of particles $n$. For $n=83$,
the ground state is the parity-antisymmetric state with FM alignment, and the
excited state is symmetric and AFM aligned. For even number of particles,
solid lines represent the triplet configuration with FM alignment, dashed
lines correspond to singlet states.

For even number of particles, however, one can observe level-crossing. The
ground state turns from an FM aligned triplet state to an AFM singlet when
$n_{0}$ differs from $n$ the actual number of particles. The kink in the solid
lines occurs where the energy of the subgap state reaches the gap edge.

The results indicate that the charging energy can also have an effect on the
spin configuration. The detailed behavior depends both on the ratios of
various energy scales $E_{C}$, $U$, $g$, $W$, and also on the gate voltage
through $n_{0}$. Exploring these dependencies is beyond the scope of the
present paper.

##  VI Conclusion

In this paper, the system consisting of a superconducting island and two
quantum dots is modelled within the framework of the charge-conserving
Richardson model. We have focused on the regime when both dots are singly
occupied with a spin-carrying particle and studied how the two spins are
effectively coupled to each other. For odd number of particles, the spins
align ferromagnetically. We find this result both in the flatband limit and in
the full model with finite bandwidth. For even number of particles, the spins
align antiferromagnetically in the flatband limit and ferromagnetically in the
wideband limit, with a transition between the two cases at some intermediate
ratio of bandwidth to pairing strength. This significant change in the nature
of the ground state implies that the flatband limit results are not fully
generic and cannot be extended to any finite bandwidth case, not even at
qualitative level. This is in contrast with what has been found for a single
QD case [16].

On the technical side, we introduced an extension of the DMRG method for
solving the general problem of two QDs coupled through a SI for arbitrary
parameters. We furthermore presented a minimal analytical model that provides
qualitatively adequate description of the spin configurations even for a
finite bandwidth.

Finally, we note that our results are also relevant for applications aiming to
split Cooper pairs by coupling a superconductor to two quantum dots. We have
shown that the singlet configuration (AF alignment) of the QD spins is the
ground state only in a limited part of the parameter space and can be realized
only when the pairing interaction is strong.

###### Acknowledgements.

We acknowledge the support of the Slovenian Research Agency (ARRS) under
P1-0416 and J1-3008.

##  Appendix A Löwdin orthogonalization of superconductor levels

We present a unitary transformation of the SC levels leading to a more
convenient form of the hybridization. In general, the hybridizations are given
by

$\displaystyle
H_{hyb,L}=v_{L}\sum_{i}\left(\gamma_{Li}d_{L\uparrow}^{+}c_{i\uparrow}+\gamma_{Li}^{*}d_{L\downarrow}^{+}c_{i\downarrow}+h.c.\right),$
$\displaystyle
H_{hyb,R}=v_{R}\sum_{i}\left(\gamma_{Ri}d_{R\uparrow}^{+}c_{i\uparrow}+\gamma_{Ri}^{*}d_{R\downarrow}^{+}c_{i\downarrow}+h.c.\right)\,.$
(62)

The operators $\sum_{i}\gamma_{Li}c_{i\uparrow}$ and
$\sum_{i}\gamma_{Ri}c_{i\uparrow}$ describe one (two) distinguished state(s)
if the overlap parameter fulfills $|\alpha|=1$ ($|\alpha|<1$). Here,
$\alpha=\sum_{i}\gamma_{Li}^{*}\gamma_{Ri}$.

First, we study the case of $|\alpha|<1$ with two distinguished states which
are not necessarily orthogonal to each other. We choose the unitary
transformation such that two of the new basis vectors lie in the subspace of
the two distinguished states, $\mathcal{D}$, while others form an orthonormal
system which is also orthogonal to $\mathcal{D}$. Inside $\mathcal{D}$, we
apply Löwdin orthogonalization [47, 48] to obtain a symmetric form of the
hybridization.

Let us introduce the operators

$\displaystyle f_{L\uparrow}^{\prime}=\sum_{i}\gamma_{Li}c_{i\uparrow},$
$\displaystyle f_{R\uparrow}^{\prime}=\sum_{i}\gamma_{Ri}c_{i\uparrow}$ (63)

which are not orthogonal for $\alpha\neq 0$. The Löwdin orthogonalization is
performed as

$\displaystyle
f_{L\uparrow}=rf_{L\uparrow}^{\prime}-p^{*}f_{R\uparrow}^{\prime}$
$\displaystyle f_{R\uparrow}=-pf_{L\uparrow}^{\prime}+rf_{R\uparrow}^{\prime}$
(64)

with

$\displaystyle
r=\frac{\sqrt{1+|\alpha|}+\sqrt{1-|\alpha|}}{2\sqrt{1-|\alpha|^{2}}},$
$\displaystyle
p=e^{i\delta}\frac{\sqrt{1+|\alpha|}-\sqrt{1-|\alpha|}}{2\sqrt{1-|\alpha|^{2}}}\,.$
(65)

The phase factor occurs if the overlap parameter is a complex number,
$\alpha=|\alpha|e^{i\delta}$. After the transformation, the hybridizations can
be rewritten as

$\displaystyle
H_{hyb,L}=v_{L}\sqrt{1-|\alpha|^{2}}\left(rd_{L\uparrow}^{+}f_{L\uparrow}+p^{*}d_{L\uparrow}^{+}f_{R\uparrow}+rd_{L\downarrow}^{+}f_{L\downarrow}+pd_{L\downarrow}^{+}f_{R\downarrow}+h.c.\right),$
$\displaystyle
H_{hyb,R}=v_{R}\sqrt{1-|\alpha|^{2}}\left(pd_{R\uparrow}^{+}f_{L\uparrow}+rd_{R\uparrow}^{+}f_{R\uparrow}+p^{*}d_{R\downarrow}^{+}f_{L\downarrow}+rd_{R\downarrow}^{+}f_{R\downarrow}+h.c.\right)\,.$
(66)

The phase of $\alpha$ can be scaled out with the gauge transformation

$\displaystyle d_{L\uparrow}=\tilde{d}_{L\uparrow}e^{-i\frac{\delta}{2}}\qquad
d_{R\uparrow}=\tilde{d}_{R\uparrow}e^{i\frac{\delta}{2}}$ $\displaystyle
f_{L\uparrow}=\tilde{f}_{L\uparrow}e^{-i\frac{\delta}{2}}\qquad
f_{R\uparrow}=\tilde{f}_{R\uparrow}e^{i\frac{\delta}{2}}$ $\displaystyle
d_{L\downarrow}=\tilde{d}_{L\downarrow}e^{i\frac{\delta}{2}}\qquad
d_{R\downarrow}=\tilde{d}_{R\downarrow}e^{-i\frac{\delta}{2}}$ $\displaystyle
f_{L\downarrow}=\tilde{f}_{L\downarrow}e^{i\frac{\delta}{2}}\qquad
f_{R\downarrow}=\tilde{f}_{R\downarrow}e^{-i\frac{\delta}{2}}\,.$ (67)

Henceforth, we will drop the notation of $\sim$ and consider that both
$\alpha$ and $p$ are equal to their absolute value. By introducing the angle
$\varphi$ as $\sin(2\varphi)=\alpha$, we obtain

$\displaystyle
H_{hyb,L}=v_{L}\sum_{\sigma}\left(\cos(\varphi)d_{L\sigma}^{+}f_{L\sigma}+\sin(\varphi)d_{L\sigma}^{+}f_{R\sigma}+h.c.\right),$
$\displaystyle
H_{hyb,R}=v_{R}\sum_{\sigma}\left(\sin(\varphi)d_{R\sigma}^{+}f_{L\sigma}+\cos(\varphi)d_{R\sigma}^{+}f_{R\sigma}+h.c.\right)\,.$
(68)

This hybridization is formally the same as that considered in Ref. [18].

The great advantage of this Hamiltonian is that in the special case of
$U_{L}=U_{R}$, $\varepsilon_{L}=\varepsilon_{R}$ and $v_{L}=v_{R}$, it has a
parity symmetry. The symmetry ensures that the system decouples to a symmetric
and an anti-symmetric sectors.

In the case of $\alpha=1$, there is no need for Löwdin orthogonalization. The
unitary transformation must be chosen such that one of the new basis vectors
is described by $f_{0\sigma}=\sum_{i}\gamma_{L/Ri}c_{i\sigma}$, while the
other new basis vectors must only fulfill orthonormality.

##  Appendix B Matrix product operator representation of the Hamiltonian

The matrix product representation of the Hamiltonian is given as
$H=\prod_{j=0}^{N+1}W_{j}$. The first element, $W_{0}$, is a vector and
defines the site with the left quantum dot:

$W_{0}=\setcounter{MaxMatrixCols}{13}\begin{pmatrix}I&[\epsilon\hat{n}_{L}+U\hat{n}_{L\uparrow}\hat{n}_{L\downarrow}]&0&0&0&0&d_{L\uparrow}&d_{L\downarrow}&d^{\dagger}_{L\uparrow}&d^{\dagger}_{L\downarrow}&0&0&0\end{pmatrix}.$
(69)

$d_{L\sigma}$ are the left quantum dot operators,
$\hat{n}_{L\sigma}=d^{\dagger}_{L\sigma}d_{L\sigma}$, and
$\hat{n}_{L}=\hat{n}_{L\uparrow}+\hat{n}_{L\downarrow}$. The superconducting
sites $j=1,\ldots,N$ are represented by matrices $W_{j}$:

$W_{j}=\setcounter{MaxMatrixCols}{13}\begin{pmatrix}1&X_{j}&-v_{R}\gamma_{Rj}c_{j\uparrow}F_{j}&-v_{R}\gamma_{Rj}^{*}c_{j\downarrow}F_{j}&v_{R}\gamma_{Rj}^{*}c^{\dagger}_{j\uparrow}F_{j}&v_{R}\gamma_{Rj}c^{\dagger}_{j\downarrow}F_{j}&0&0&0&0&gy_{j}c_{j\downarrow}c_{j\uparrow}&gy_{j}c^{\dagger}_{j\uparrow}c^{\dagger}_{j\downarrow}&2E_{c}n_{j}\\\
0&I&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&F_{j}&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&F_{j}&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&F_{j}&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&F_{j}&0&0&0&0&0&0&0\\\
0&v_{L}\gamma_{Lj}^{*}c^{\dagger}_{j\uparrow}&0&0&0&0&F_{j}&0&0&0&0&0&0\\\
0&v_{L}\gamma_{Lj}c^{\dagger}_{j\downarrow}&0&0&0&0&0&F_{j}&0&0&0&0&0\\\
0&v_{L}\gamma_{Lj}c_{j\uparrow}&0&0&0&0&0&0&F_{j}&0&0&0&0\\\
0&v_{L}\gamma_{Lj}^{*}c_{j\downarrow}&0&0&0&0&0&0&0&F_{j}&0&0&0\\\
0&y_{j}c^{\dagger}_{j\uparrow}c^{\dagger}_{j\downarrow}&0&0&0&0&0&0&0&0&1&0&0\\\
0&y_{j}c_{j\downarrow}c_{j\uparrow}&0&0&0&0&0&0&0&0&0&1&0\\\
0&n_{j}&0&0&0&0&0&0&0&0&0&0&1\end{pmatrix},$ (70)

where the on-site potential term is

$X_{j}=[\epsilon_{j}+E_{c}(1-2n_{0})]\hat{n}_{j}+(y_{j}^{2}g+2E_{c})\hat{n}_{j\uparrow}\hat{n}_{j\downarrow}.$
(71)

$F_{j}=(-1)^{\hat{n}_{j}}$ is the local fermionic parity operator. It takes
care of the fermionic anticommutation rules by giving a factor of $-1$ if the
$j$-th level is occupied by a single particle. Finally, the right quantum dot
is represented by a vector $W_{N+1}$:

$W_{N+1}=\setcounter{MaxMatrixCols}{13}\begin{pmatrix}[\epsilon_{R}+Un_{R\uparrow}n_{R\downarrow})n_{R}]&1&d^{\dagger}_{R\uparrow}&d^{\dagger}_{R\downarrow}&d_{R\uparrow}&d_{R\downarrow}&0&0&0&0&0&0&0\end{pmatrix}.$
(72)

##  Appendix C Superexchange between two dots

In this section, we consider a single level mediating exchange interaction
between two QDs. Each dot hosts a single particle and we assume symmetric
configuration of parameters. The interaction then depends on the occupancy of
the intermediate level.

The Hamiltonian is given as

$\begin{split}H=&\varepsilon\sum_{\sigma}\left(d^{+}_{L\sigma}d_{L\sigma}+d^{+}_{R\sigma}d_{R\sigma}\right)+U\sum_{I=L,R}n_{I\uparrow}n_{I\downarrow}\\\
&+\varepsilon_{f}\sum_{\sigma}f^{+}_{\sigma}f_{\sigma}+v\sum_{\sigma}\left(d^{+}_{L\sigma}f_{\sigma}+d^{+}_{R\sigma}f_{\sigma}+h.c.\right),\end{split}$
(73)

where the on-site interaction $U$ may be finite. Note that in the pairing
model of the main text, there is no on-site potential term; that situation
corresponds to $\epsilon_{f}=0$. In what follows, we systematically study all
cases depending on the occupancy of the intermediary level and compare the AF
and FM configurations.

Zero occupancy:

For FM aligned case, a representative state is given as
$|\uparrow,\uparrow\rangle\otimes|0\rangle$, where the first part
$|\uparrow,\uparrow\rangle$ describes the state of the two dots, while the
second part is the state of the intermediate level. The basis vectors in this
spin sector are

$\displaystyle|\uparrow,\uparrow\rangle\otimes|0\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle-|0,\uparrow\rangle\right)\otimes|\uparrow\rangle\,.$
(74)

The Hamiltonian is represented by

$\displaystyle\mathbf{H}_{T_{0}}=\left[\begin{array}[]{cc}2\varepsilon&\sqrt{2}v\\\
\sqrt{2}v&\varepsilon+\varepsilon_{f}\end{array}\right],$ (77)

where the index $T_{0}$ indicates that the system is an overall triplet state.

For AF aligned case, the basis vectors are

$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|0\rangle,$
$\displaystyle\frac{|\uparrow,0\rangle+|0,\uparrow\rangle}{2}\otimes|\downarrow\rangle-\frac{|\downarrow,0\rangle+|0,\downarrow\rangle}{2}\otimes|\uparrow\rangle,$
$\displaystyle|0,0\rangle\otimes|\uparrow\downarrow\rangle,$ (78)
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow,0\rangle+|0,\uparrow\downarrow\rangle\right)\otimes|0\rangle,$
(79)

which describe overall singlet states. On this subspace, the Hamiltonian is
represented by

$\displaystyle\mathbf{H}_{S_{0}}=\left[\begin{array}[]{cccc}2\varepsilon&\sqrt{2}v&0&0\\\
\sqrt{2}v&\varepsilon+\varepsilon_{f}&2v&\sqrt{2}v\\\
0&2v&2\varepsilon_{f}&0\\\ 0&\sqrt{2}v&0&2\varepsilon+U\end{array}\right].$
(84)

By calculating the minimal energy for both $T_{0}$ and $S_{0}$, we find that
the singlet sector has lower minimal energy, see the blue lines in Fig. 6.

Single occupancy:

In this case, the $f$ level hosts one particle with spin-down. The total spin
is 1/2 and, hence, the sector describes doublet states but we distinguish
between parity symmetric and anti-symmetric subsectors exhibiting AF and FM
alignment on the dots, respectively. In the parity-anti-symmetric sector, the
basis vectors are given by

$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,0\rangle-|0,\uparrow\rangle\right)\otimes|\uparrow\downarrow\rangle,$
$\displaystyle\frac{1}{\sqrt{6}}\left(2|\uparrow,\uparrow\rangle\otimes|\downarrow\rangle-\left(|\uparrow,\downarrow\rangle+|\downarrow,\uparrow\rangle\right)\otimes|\uparrow\rangle\right),$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,\uparrow\downarrow\rangle-|\uparrow\downarrow,\uparrow\rangle\right)\otimes|0\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|0,\uparrow\downarrow\rangle-|\uparrow\downarrow,0\rangle\right)\otimes|\uparrow\rangle,$
(85)

and the Hamiltonian is represented by

$\displaystyle\mathbf{H}_{D,asym}=\left[\begin{array}[]{cccc}\varepsilon+2\varepsilon_{f}&\sqrt{3}v&0&v\\\
\sqrt{3}v&2\varepsilon+\varepsilon_{f}&\sqrt{3}v&0\\\
0&\sqrt{3}v&3\varepsilon+U&v\\\
v&0&v&2\varepsilon+U+\varepsilon_{f}\end{array}\right]\,.$ (90)

In the parity symmetric sector, the basis vectors are defined by

$\displaystyle\frac{1}{\sqrt{2}}\left(|\downarrow,0\rangle+|0,\downarrow\rangle\right)\otimes|\uparrow\downarrow\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|\downarrow\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow,\downarrow\rangle+|\downarrow,\uparrow\downarrow\rangle\right)\otimes|0\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow,0\rangle+|0,\uparrow\downarrow\rangle\right)\otimes|\downarrow\rangle,$
(91)

and the Hamiltonian reads

$\displaystyle\mathbf{H}_{D,sym}=\left[\begin{array}[]{cccc}\varepsilon+2\varepsilon_{f}&-v&0&-v\\\
-v&2\varepsilon+\varepsilon_{f}&-v&0\\\ 0&-v&3\varepsilon+U&v\\\
-v&0&v&2\varepsilon+U+\varepsilon_{f}\end{array}\right]\,.$ (96)

By calculating the minimal energy for both $D_{sym}$ and $D_{asym}$, we find
that the sector $D_{asym}$ featuring ferromagnetic alignment has lower minimal
energy, see the red lines in Fig. 6.

We note that the quadruplet states (e.g. $S_{z}=3/2$) are always higher in
energy.

Double occupancy:

This sector is analogous to the case of zero occupancy. The $f$ level now
hosts two particles. In the triplet sector, the basis states are

$\displaystyle|\uparrow,\uparrow\rangle\otimes|\uparrow\downarrow\rangle,$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow,\uparrow\rangle-|\uparrow,\uparrow\downarrow\rangle\right)\otimes|\uparrow\rangle,$
(97)

with the Hamiltonian

$\displaystyle\mathbf{H}_{T_{\uparrow\downarrow}}=\left[\begin{array}[]{cc}2\varepsilon+2\varepsilon_{f}&\sqrt{2}v\\\
\sqrt{2}v&3\varepsilon+U+\varepsilon_{f}\end{array}\right]\,.$ (100)

In the singlet sector, we have the following states.

$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|\uparrow\downarrow\rangle,$
$\displaystyle\frac{1}{2}\left((|\uparrow\downarrow,\downarrow\rangle+|\downarrow,\uparrow\downarrow\rangle)\otimes|\uparrow\rangle-(|\uparrow\downarrow,\uparrow\rangle+|\uparrow,\uparrow\downarrow\rangle)\otimes|\downarrow\rangle\right),$
$\displaystyle\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow,0\rangle+|0,\uparrow\downarrow\rangle\right)\otimes|\uparrow\downarrow\rangle,$
$\displaystyle|\uparrow\downarrow,\uparrow\downarrow\rangle\otimes|0\rangle.$
(101)

The Hamiltonian is represented by the matrix

$\displaystyle\mathbf{H}_{S_{\uparrow\downarrow}}=\left[\begin{array}[]{cccc}2\varepsilon+2\varepsilon_{f}&\sqrt{2}v&0&0\\\
\sqrt{2}v&3\varepsilon+U+\varepsilon_{f}&-\sqrt{2}v&-2v\\\
0&-\sqrt{2}v&2\varepsilon+U+2\varepsilon_{f}&0\\\
0&-2v&0&4\varepsilon+2U\end{array}\right]\,.$ (106)

As in the zero occupancy sector, the singlet has lower minimal energy, see the
green lines in Fig. 6.

Figure 6: Comparison of minimal energies of anti-ferromagnetically (solid
line) and ferromagnetically (dashed line) aligned configurations in case of
different occupancies of the intermediate level. For the plot,
$\varepsilon_{f}=2$, $U=10$ and $\varepsilon=-5$ were chosen in arbitrary
units.

So far, we have seen that for zero/single/double occupancy of the
intermediating level, the sector $S_{0}$/$D_{asym}$/$S_{\uparrow\downarrow}$
has the lowest minimal energy. The phase diagram in Fig. 7 shows the sector
with the lowest overall energy as a function of $v$ and $\varepsilon_{f}$.

Figure 7: Phase diagram of the ground state at $\varepsilon=-5$ and $U=10$.

##  Appendix D Minimal analytical model of kinetic energy to capture

The kinetic energy can be taken into account with the hopping process between
the distinguished level and all other levels. The corresponding Hamiltonian is
defined as

$\displaystyle H_{kin}=\frac{t}{\sqrt{N}}\sum_{j\neq
0,\sigma}\left(f_{0\sigma}^{+}f_{j\sigma}+h.c.\right),$ (107)

which is added to the Hamiltonian in the flat-band limit as given in Eq. (7).
We note that, in general, one could consider hopping terms between all
remaining levels but that would break analytical solvability. We shall see
that this minimal model is sufficient to capture the singlet-triplet
transition as a function of bandwidth.

In this Appendix, we study only the case of an even number of particles
because this is the only sector where the DMRG shows level-crossing.

First, we consider the states in the singlet sector. The set of basis states
of Eq. (17) must be extended as follows:

$\begin{split}|\phi_{0}^{S}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|0,\Psi^{N-1}_{M}\rangle,\\\
|\phi_{1}^{S}\rangle&=\frac{1}{2}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes\left(|\uparrow,\Psi^{N-1}_{M-1+\downarrow}\rangle-|\downarrow,\Psi^{N-1}_{M-1+\uparrow}\rangle\right),\\\
|\phi_{2}^{S}\rangle&=\frac{1}{\sqrt{2}}\left(|\uparrow,\downarrow\rangle-|\downarrow,\uparrow\rangle\right)\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{-1X}^{S}\rangle&=\frac{1}{2}\Big{(}\left(|\uparrow,0\rangle+|0,\uparrow\rangle\right)\otimes|0,\Psi^{N-1}_{M+\downarrow}\rangle-\left(|\downarrow,0\rangle+|0,\downarrow\rangle\right)\otimes|0,\Psi^{N-1}_{M+\uparrow}\rangle\Big{)},\\\
|\phi_{0X}^{S}\rangle&=\frac{1}{2}\Big{(}(|\uparrow,0\rangle+|0,\uparrow\rangle)\otimes|\downarrow,\Psi_{M}^{N-1}\rangle-(|0,\downarrow\rangle+|\downarrow,0\rangle)\otimes|\uparrow,\Psi_{M}^{N-1}\rangle\Big{)},\\\
|\phi_{1X}^{S}\rangle&=\frac{1}{2}\Big{(}\left(|\uparrow,0\rangle+|0,\uparrow\rangle\right)\otimes|\uparrow\downarrow,\Psi^{N-1}_{M-1+\downarrow}\rangle-\left(|\downarrow,0\rangle+|0,\downarrow\rangle\right)\otimes|\uparrow\downarrow,\Psi^{N-1}_{M-1+\uparrow}\rangle\Big{)},\\\
|\phi_{-2D}^{S}\rangle&=|0,0\rangle\otimes|0,\Psi_{M+1}^{N-1}\rangle,\\\
|\phi_{-1D}^{S}\rangle&=\frac{1}{\sqrt{2}}|0,0\rangle\otimes\left(|\uparrow,\Psi_{M+\downarrow}^{N-1}\rangle-|\downarrow,\Psi_{M+\uparrow}^{N-1}\rangle\right),\\\
|\phi_{0D}^{S}\rangle&=|0,0\rangle\otimes|\uparrow\downarrow,\Psi_{M}^{N-1}\rangle.\end{split}$
(108)

On this space, the Hamiltonian is represented by the block matrix

$\displaystyle\mathbf{H}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c|c}2\varepsilon+\mathbf{H}_{SC+kin}&\begin{array}[]{ccc}0&\sqrt{2}v&0\\\
0&0&-v\\\ 0&0&0\end{array}&0\\\ \hline\cr\begin{array}[]{ccc}0&0&0\\\
\sqrt{2}v&0&0\\\
0&-v&0\end{array}&\varepsilon+\tilde{\mathbf{H}}_{SC+kin}&\begin{array}[]{ccc}0&\sqrt{2}v&0\\\
0&0&2v\\\ 0&0&0\end{array}\\\ \hline\cr 0&\begin{array}[]{ccc}0&0&0\\\
\sqrt{2}v&0&0\\\
0&2v&0\end{array}&g(2\nu-1)+\mathbf{H}_{SC+kin}\end{array}\right],$ (124)

where $\mathbf{H}_{SC+kin}$ is defined in Eq. (61) and

$\displaystyle\tilde{\mathbf{H}}_{SC+kin}=\left[\begin{array}[]{ccc}2g\nu&t\sqrt{1-\nu}&\beta\\\
t\sqrt{1-\nu}&g\nu&-t\sqrt{\nu}\\\ \beta&-t\sqrt{\nu}&g\end{array}\right]\,.$
(128)

We note that in principle the operator $H_{kin}$ leads out of this space and
we should also take into account states in which two Cooper pairs are broken.
However, as we mentioned before, these are very high in energy, therefore we
neglect them.

In the triplet sector, the basis states of Eq. (28) must be extended as
follows.

$\begin{split}|\phi_{0}^{T}\rangle&=|\uparrow,\uparrow\rangle\otimes|0,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{1}^{T}\rangle&=\frac{1}{\sqrt{2}}|\uparrow,\uparrow\rangle\otimes\Big{(}|\uparrow,\Psi_{M-1+\downarrow}^{N-1}\rangle-|\downarrow,\Psi_{M-1+\uparrow}^{N-1}\Big{)},\\\
|\phi_{2}^{T}\rangle&=|\uparrow,\uparrow\rangle\otimes|\uparrow\downarrow,\Psi_{M-1}^{N-1}\rangle,\\\
|\phi_{1}^{\prime
T}\rangle&=\frac{1}{2}\left[|\uparrow,\uparrow\rangle\otimes\Big{(}|\uparrow,\Psi_{M-1+\downarrow}^{N-1}\rangle+|\downarrow,\Psi_{M-1+\uparrow}^{N-1}\Big{)}-\Big{(}|\uparrow,\downarrow\rangle+|\downarrow,\uparrow\rangle\Big{)}\otimes|\uparrow,\Psi_{M-1+\uparrow}^{N-1}\rangle\right],\\\
|\phi_{-1X}^{T}\rangle&=\frac{1}{\sqrt{2}}\Big{(}|\uparrow,0\rangle-|0,\uparrow\rangle\Big{)}\otimes|0,\Psi_{M+\uparrow}^{N-1}\rangle,\\\
|\phi_{0X}^{T}\rangle&=\frac{1}{\sqrt{2}}\Big{(}|\uparrow,0\rangle-|0,\uparrow\rangle\Big{)}\otimes|\uparrow,\Psi_{M}^{N-1}\rangle,\\\
|\phi_{1X}^{T}\rangle&=\frac{1}{\sqrt{2}}\Big{(}|\uparrow,0\rangle-|0,\uparrow\rangle\Big{)}\otimes|\uparrow\downarrow,\Psi_{M-1+\uparrow}^{N-1}\rangle.\end{split}$
(129)

We note that these states all belong to the total $S_{z}=1$. The sectors
corresponding to $S_{z}=0$ and $S_{z}=-1$ would lead to the same Hamiltonian
matrix. On this space, the Hamiltonian is represented by

$\displaystyle\mathbf{H}=E_{0}\mathbf{I}+\left[\begin{array}[]{c|c|c}2\varepsilon+\mathbf{H}_{SC}&\begin{array}[]{c}0\\\
0\\\ 0\end{array}&\begin{array}[]{ccc}0&\sqrt{2}v&0\\\ 0&0&-v\\\
0&0&0\end{array}\\\
\hline\cr\begin{array}[]{ccc}0&0&0\end{array}&2\varepsilon+g&\begin{array}[]{ccc}0&0&\sqrt{2}v\end{array}\\\
\hline\cr\begin{array}[]{ccc}0&0&0\\\ \sqrt{2}v&0&0\\\
0&-v&0\end{array}&\begin{array}[]{c}0\\\ 0\\\
\sqrt{2}v\end{array}&\varepsilon+\tilde{\mathbf{H}}_{SC}\end{array}\right].$
(147)

To compare the singlet and triplet sectors, we obtain the lowest minimal
energies within the perturbation theory. For $t,v\ll g\ll\varepsilon$, the
energy difference is found to be

$E_{T}-E_{S}=\frac{2v^{4}}{|\varepsilon|^{3}}\left(1-\frac{3|\varepsilon|t^{2}}{2g^{3}}\right),$
(148)

indicating that the energy difference decreases with increasing $t$. The
decreasing trend is induced by the sixth order process which admixes the state
$|\phi_{1}^{\prime T}\rangle$, which has a relatively low energy of
$2\varepsilon+g$, to the ground state. The sixth order process involves
breaking a Cooper pair and a superexchange procedure between the states
$|\phi_{1}^{\prime T}\rangle$ and $|\phi_{1}^{T}\rangle$ through one member of
the broken pair. This process only exists in the triplet sector since there is
no singlet counterpart to the $\phi_{1}^{\prime T}$ state. As $t$ increases,
the sixth order process competes with the fourth order processes favouring AF
alignment and will eventually overtake it at
$t=t_{0}=\sqrt{2g^{3}/3|\epsilon|}$.

Finally, we provide the perturbation theory results for both $E_{S}$ and
$E_{T}$. We list only the most important sixth order terms, which involve
states with both dots occupied since these are the most relevant for large
$|\varepsilon|$. For small $v$ and $t$,

$\begin{split}E_{S}&=E_{0}+2\varepsilon+\frac{v^{2}}{\varepsilon-\frac{g}{2}}-\frac{2t^{2}}{g}+\frac{2t^{2}v^{2}}{g\left(\varepsilon-\frac{g}{2}\right)^{2}}-\frac{v^{4}}{g\left(\varepsilon-\frac{g}{2}\right)^{2}}+\frac{4t^{4}}{g^{3}}+\frac{v^{4}}{\varepsilon\left(\varepsilon-\frac{g}{2}\right)^{2}}+E_{S}^{(6)}\\\
E_{T}&=E_{0}+2\varepsilon+\frac{v^{2}}{\varepsilon-\frac{g}{2}}-\frac{2t^{2}}{g}+\frac{2t^{2}v^{2}}{g\left(\varepsilon-\frac{g}{2}\right)^{2}}-\frac{v^{4}}{g\left(\varepsilon-\frac{g}{2}\right)^{2}}+\frac{4t^{4}}{g^{3}}-\frac{v^{4}}{\left(\varepsilon-\frac{g}{2}\right)^{3}}+E_{T}^{(6)}\\\
E_{S}^{(6)}&=\frac{t^{2}v^{4}}{g^{2}(2\varepsilon-g)}\left(\frac{1}{\left(\varepsilon-\frac{g}{2}\right)^{2}}+\frac{1}{\left(\varepsilon-\frac{3g}{2}\right)^{2}}-\frac{1}{\left(\varepsilon-\frac{g}{2}\right)\left(\varepsilon-\frac{3g}{2}\right)}\right)\\\
E_{T}^{(6)}&=\frac{t^{2}v^{4}}{g^{2}(-g)}\left(\frac{1}{\left(\varepsilon-\frac{g}{2}\right)^{2}}+\frac{1}{\left(\varepsilon-\frac{3g}{2}\right)^{2}}+\frac{1}{\left(\varepsilon-\frac{g}{2}\right)\left(\varepsilon-\frac{3g}{2}\right)}\right)\end{split}$
(149)

and in the limit of $g\ll|\varepsilon|$, we obtain

$\begin{split}E_{S}&=E_{0}+2\varepsilon+\frac{v^{2}}{\varepsilon}-\frac{2t^{2}}{g}+\frac{2t^{2}v^{2}}{g^{2}\varepsilon}-\frac{v^{4}}{g\varepsilon^{2}}+\frac{4t^{4}}{g^{3}}+\frac{v^{4}}{\varepsilon^{3}}+\frac{t^{2}v^{4}}{2g^{2}\varepsilon^{3}}\\\
E_{T}&=E_{0}+2\varepsilon+\frac{v^{2}}{\varepsilon}-\frac{2t^{2}}{g}+\frac{2t^{2}v^{2}}{g^{2}\varepsilon}-\frac{v^{4}}{g\varepsilon^{2}}+\frac{4t^{4}}{g^{3}}-\frac{v^{4}}{\varepsilon^{3}}-\frac{3t^{2}v^{4}}{g^{3}\varepsilon^{2}}\end{split}$
(150)

which is reformulated in Eq. (148) in the limiting case of
$g\ll|\varepsilon|$.

## References

  * Yu et al. [2002] Y. Yu, S. Han, X. Chu, S.-I. Chu, and Z. Wang, Science 296, 889 (2002), eprint https://www.science.org/doi/pdf/10.1126/science.1069452, URL https://www.science.org/doi/abs/10.1126/science.1069452.
  * Koch et al. [2007] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Physical Review A 76 (2007), URL https://doi.org/10.1103/physreva.76.042319.
  * Devoret and Schoelkopf [2013] M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013), URL https://doi.org/10.1126/science.1231930.
  * Wendin [2017] G. Wendin, Reports on Progress in Physics 80, 106001 (2017), URL https://doi.org/10.1088/1361-6633/aa7e1a.
  * Arute et al. [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, et al., Nature 574, 505 (2019), URL https://doi.org/10.1038/s41586-019-1666-5.
  * Kjaergaard et al. [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Annual Review of Condensed Matter Physics 11, 369 (2020), URL https://doi.org/10.1146/annurev-conmatphys-031119-050605.
  * Clerk et al. [2020] A. A. Clerk, K. W. Lehnert, P. Bertet, J. R. Petta, and Y. Nakamura, Nature Physics 16, 257 (2020), URL https://doi.org/10.1038/s41567-020-0797-9.
  * Aguado [2020] R. Aguado, Applied Physics Letters 117, 240501 (2020).
  * Blais et al. [2021] A. Blais, A. L. Grimsmo, S. Girvin, and A. Wallraff, Reviews of Modern Physics 93 (2021), URL https://doi.org/10.1103/revmodphys.93.025005.
  * Siddiqi [2021] I. Siddiqi, Nature Reviews Materials 6, 875 (2021), URL https://doi.org/10.1038/s41578-021-00370-4.
  * Hendrickx et al. [2021] N. W. Hendrickx, W. I. L. Lawrie, M. Russ, F. van Riggelen, S. L. de Snoo, R. N. Schouten, A. Sammak, G. Scappucci, and M. Veldhorst, Nature 591, 580 (2021), URL https://doi.org/10.1038/s41586-021-03332-6.
  * Taylor et al. [2005] J. M. Taylor, H.-A. Engel, W. Dür, A. Yacoby, C. M. Marcus, P. Zoller, and M. D. Lukin, Nature Physics 1, 177 (2005), URL https://doi.org/10.1038/nphys174.
  * Burkard et al. [2021] G. Burkard, T. D. Ladd, J. M. Nichol, A. Pan, and J. R. Petta, _Semiconductor spin qubits_ (2021), eprint 2112.08863.
  * Karrasch et al. [2011] C. Karrasch, S. Andergassen, and V. Meden, Physical Review B 84, 134512 (2011), URL https://doi.org/10.1103/physrevb.84.134512.
  * Padurariu and Nazarov [2010] C. Padurariu and Y. V. Nazarov, Physical Review B 81, 144519 (2010).
  * Žitko and Pavešić [2022] R. Žitko and L. Pavešić, Phys. Rev. B 106, 024513 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.106.024513.
  * Pavešić and Žitko [2022] L. Pavešić and R. Žitko, Phys. Rev. B 105, 075129 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.075129.
  * Malinowski [2023] F. K. Malinowski, _Two anderson impurities coupled through a superconducting island: charge stability diagrams and double impurity qubit_ (2023), eprint 2303.14410.
  * Kürtössy et al. [2021] O. Kürtössy, Z. Scherübl, G. Fülöp, I. E. Lukács, T. Kanne, J. Nygård, P. Makk, and S. Csonka, Nano Letters 21, 7929 (2021), URL https://doi.org/10.1021/acs.nanolett.1c01956.
  * Petersson et al. [2012] K. D. Petersson, L. W. McFaul, M. D. Schroer, M. Jung, J. M. Taylor, A. A. Houck, and J. R. Petta, Nature 490, 380 (2012), URL https://doi.org/10.1038/nature11559.
  * Mi et al. [2017] X. Mi, J. V. Cady, D. M. Zajac, P. W. Deelman, and J. R. Petta, Science 355, 156 (2017), URL https://doi.org/10.1126/science.aal2469.
  * Borjans et al. [2019] F. Borjans, X. G. Croot, X. Mi, M. J. Gullans, and J. R. Petta, Nature 577, 195 (2019), URL https://doi.org/10.1038/s41586-019-1867-y.
  * Pillet et al. [2019] J.-D. Pillet, V. Benzoni, J. Griesmar, J.-L. Smirr, and Ç. O. Girit, Nano Letters 19, 7138 (2019), URL https://doi.org/10.1021/acs.nanolett.9b02686.
  * Zwanenburg et al. [2013] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Reviews of Modern Physics 85, 961 (2013), URL https://doi.org/10.1103/revmodphys.85.961.
  * Kandel et al. [2021] Y. P. Kandel, H. Qiao, and J. M. Nichol, Applied Physics Letters 119, 030501 (2021), URL https://doi.org/10.1063/5.0055908.
  * Yao et al. [2014] N. Y. Yao, L. I. Glazman, E. A. Demler, M. D. Lukin, and J. D. Sau, Phys. Rev. Lett. 113, 087202 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.087202.
  * Bardeen et al. [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 106, 162 (1957), URL https://link.aps.org/doi/10.1103/PhysRev.106.162.
  * von Delft and Ralph [2001] J. von Delft and D. Ralph, Physics Reports 345, 61 (2001), URL https://doi.org/10.1016/s0370-1573(00)00099-5.
  * Dean and Hjorth-Jensen [2003] D. J. Dean and M. Hjorth-Jensen, Reviews of Modern Physics 75, 607 (2003), URL https://doi.org/10.1103/revmodphys.75.607.
  * Tempere and Huybrechts [2019] J. Tempere and D. Huybrechts, The European Physical Journal B 92, 129 (2019), ISSN 1434-6036, URL https://doi.org/10.1140/epjb/e2019-90564-5.
  * Richardson and Sherman [1964] R. Richardson and N. Sherman, Nuclear Physics 52, 221 (1964), ISSN 0029-5582, URL https://www.sciencedirect.com/science/article/pii/002955826490687X.
  * Scherübl et al. [2020] Z. Scherübl, G. Fülöp, C. P. Moca, J. Gramich, A. Baumgartner, P. Makk, T. Elalaily, C. Schönenberger, J. Nygård, G. Zaránd, et al., Nature Communications 11 (2020), URL https://doi.org/10.1038/s41467-020-15322-9.
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2011] U. Schollwöck, Annals of Physics 326, 96 (2011), ISSN 0003-4916, january 2011 Special Issue, URL http://www.sciencedirect.com/science/article/pii/S0003491610001752.
  * Pavešć et al. [2021] L. Pavešć, D. Bauernfeind, and R. Žitko, Phys. Rev. B 104, L241409 (2021).
  * Saldaña et al. [2022] J. C. E. Saldaña, A. Vekris, L. Pavešić, P. Krogstrup, R. Žitko, K. Grove-Rasmussen, and J. Nygård, Nat. Commun. 13, 2243 (2022).
  * Anderson [1950] P. W. Anderson, Phys. Rev. 79, 350 (1950), URL https://link.aps.org/doi/10.1103/PhysRev.79.350.
  * Goodenough [1955] J. B. Goodenough, Phys. Rev. 100, 564 (1955), URL https://link.aps.org/doi/10.1103/PhysRev.100.564.
  * Anderson [1959] P. W. Anderson, Phys. Rev. 115, 2 (1959), URL https://link.aps.org/doi/10.1103/PhysRev.115.2.
  * Vojta et al. [2002] M. Vojta, R. Bulla, and W. Hofstetter, Phys. Rev. B 65, 140405(R) (2002).
  * Hofstetter and Schoeller [2002] W. Hofstetter and H. Schoeller, Phys. Rev. Lett. 88, 016803 (2002).
  * Žitko and Bonča [2006] R. Žitko and J. Bonča, Phys. Rev. B 74, 045312 (2006).
  * Pustilnik and Borda [2006] M. Pustilnik and L. Borda, Phys. Rev. B 73, 201301(R) (2006).
  * Žitko and Bonča [2007] R. Žitko and J. Bonča, Phys. Rev. B 76, 241305(R) (2007).
  * Logan et al. [2009] D. E. Logan, C. J. Wright, and M. R. Galpin, Phys. Rev. B 80, 125117 (2009).
  * Žitko et al. [2012] R. Žitko, J. Mravlje, and K. Haule, Phys. Rev. Lett. 108, 066602 (2012).
  * Löwdin [1950] P. Löwdin, The Journal of Chemical Physics 18, 365 (1950), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/18/3/365/7387148/365_1_online.pdf, URL https://doi.org/10.1063/1.1747632.
  * Piela [2014] L. Piela, in _Ideas of Quantum Chemistry (Second Edition)_ , edited by L. Piela (Elsevier, Oxford, 2014), pp. e99–e103, second edition ed., ISBN 978-0-444-59436-5, URL https://www.sciencedirect.com/science/article/pii/B9780444594365000301.