# Generalized transmon Hamiltonian for Andreev spin qubits

Luka Pavešić<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39, SI-1000
Ljubljana, Slovenia  Faculty of Mathematics and Physics, University of
Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia  Rok Žitko
<EMAIL_ADDRESS>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 solve the problem of an interacting quantum dot embedded in a Josephson
junction between two superconductors with finite charging energy described by
the transmon (Cooper pair box) Hamiltonian. The approach is based on the flat-
band approximation of the Richardson model, which reduces the Hilbert space to
the point where exact diagonalisation is possible while retaining all states
that are necessary to describe the low energy phenomena. The presented method
accounts for the physics of the quantum dot, the Josephson effect and the
Coulomb repulsion (charging energy) at the same level. In particular, it
captures the quantum fluctuations of the superconducting phase as well as the
coupling between the superconducting phase and the quantum dot (spin) degrees
of freedom. The method can be directly applied for modelling Andreev spin
qubits embedded in transmon circuits in all parameter regimes, for describing
time-dependent processes, and for the calculation of transition matrix
elements for microwave-driven transmon, spin-flip and mixed transitions that
involve coupling to charge or current degree of freedom.

##  I Introduction

Superconducting circuits are one of the leading platforms for the realization
of various quantum technological applications. Most implementations of
superconducting qubits are based on creating an anharmonic oscillator by
replacing the inductor in an LC-circuit with a Josephson junction (JJ), which
has non-linear inductance [1, 2, 3]. The realization where the charging energy
due to capacitance is small compared to the Josephson energy is called the
transmon qubit [4]. As all other superconducting devices, the transmons
utilize the macroscopic coherence of superconducting states to encode and
manipulate quantum information [5]. Their popularity is due to their
robustness with respect to charge fluctuations, which is one of the main
decoherence mechanisms in superconducting qubits.

In pursuit of further enhancing these devices, a novel approach has emerged,
combining the robust coherence of superconductors (SCs) with the
controllability of spin qubits built out of semiconducting quantum dots (QD).
The idea consists of embedding a QD into the JJ and storing quantum
information in the spin of the quasiparticle trapped in discrete subgap states
that emerge in the few-channel regime of the JJ. The architecture is called
the Andreev spin qubit (ASQ) [6, 7, 8, 9, 10]. The spin-orbit coupling (SOC)
permits manipulation of the spin degree of freedom using the supercurrent or
the electric field, as well as advanced readout based on circuit quantum
electrodynamics (cQED) techniques. If the ASQ is embedded in a transmon, such
setup can actually support two physical qubits, one defined in the QD spin and
the other in the standard computational subspace of transmon excitations [7].

Modelling of the transmon qubit typically relies on neglecting the
superconducting quasiparticles and only accounting for the dynamics of Cooper
pairs [4, 2]. The SC gap—the energy scale for creating quasiparticles—is
typically much larger than the energy of transmon excitations, so this is a
good approximation. However, the presence of an interacting QD in the JJ
induces Cooper pair-breaking processes, and an accurate description of the
physics requires solving the full electronic problem. This is particularly
important for modelling the ASQ, where a QD with a large charging energy is
favored as it ensures a ground state with a single localized spin trapped in
the QD.

Two non-interacting SC leads coupled to an interacting QD constitute a quantum
impurity problem, which is numerically solvable with standard impurity
solvers, such as the numerical renormalization group [11, 12]. The numerical
procedure hinges on the separation of high and low energy scales, a feature of
non-interacting leads. However, simultaneously including also the charging
energy of the SC islands—a critical feature of the transmon—makes the leads
interacting. As far as we are aware, a method for attacking such problems does
not exist.

In this work, we present a method for solving the microscopic QD-transmon
model. It is based on the flat-band approximation for the SC islands, which
are described by a single active orbital coupled to a condensate of Cooper
pairs. This approach exponentially reduces the Hilbert space to the point
where exact numerical diagonalisation is possible for a system of thousands of
electrons, while retaining the key subspace for capturing the low-energy
phenomena.

The paper is formatted as follows: In Section II we describe the basic
properties of the QD–transmon system we aim to reproduce. In Sec. III we
introduce the model and the reduced basis. Section IV contains benchmark
results that validate the method. We reproduce known behaviour and point out
parameter regimes where the standard transmon approximation is not adequate.
Section V presents three specific examples of possible applications of our
approach to problems that cannot be solved in any other way: 1) the case where
the QD and SC degrees of freedom are strongly coupled, 2) the case of time-
dependent perturbations, such as a microwave pulse, 3) the calculation of
transition matrix elements involving both qubit degrees of freedom.

##  II Basic properties of the quantum dot Josephson junctions

In this section we review the basic properties of quantum dot Josephson
junctions to set the stage for the discussion of the full model. The focus is
on the interplay of electronic processes that determine the nature of low-
energy states in the superconducting gap.

###  II.1 Cooper pair tunneling

In a conventional JJ without an embedded QD quasiparticles play no role at
temperatures much lower than the superconducting gap $\Delta$. An approximate
description in terms of Cooper pair hopping is adequate to capture the
Josephson effect. Microscopically, it is shown to arise form the coherent
transfer of electron pairs across the junction [13, 14, 15].

This can be expressed by writing the Hamiltonian in the charge basis
$|m\rangle$, with $m=0,1,-1,\dots$ the difference in the number of Cooper
pairs occupying the two SCs. We allow processes where a single Cooper pair
hops across the junction. These terms couple $|m\rangle$ to the neighbouring
$|m\pm 1\rangle$, and the Hamiltonian in this basis is tridiagonal [14]. As
the different $|m\rangle$ states are equivalent, this corresponds to a tight-
binding chain in $m$-space. It is diagonalised by Fourier transforming
$|m\rangle$ into states labelled by the dual quantity—this is the emergent
superconducting phase difference $\phi$. The eigenstates are Bloch waves,
superpositions of $|m\rangle$ states,
$|\phi\rangle\propto\sum_{m}e^{im\phi}|m\rangle$. Their dispersion is given by

$E(\phi)=-E_{J}\cos\phi,$ (1)

where the emergent energy scale, the Josephson energy $E_{J}$, is given by
twice the hopping matrix element between $|m\rangle$ and $|m\pm 1\rangle$,
i.e., the energy associated with transferring a single Cooper pair across the
junction.

###  II.2 Charging energy

In nanoscopic structures, the small size of the SCs that form the junction
implies large Coulomb repulsion between the electrons occupying the device
[16, 17]. This is encapsulated in the charging energy $E_{c}$, which also
includes other capacitance effects (such as the capacitive shunting of the JJ
in the transmon [4, 1]).

The problem is described by the transmon Hamiltonian (also known as the Cooper
pair box Hamiltonian), which combines pair hopping and charging terms. It was
first introduced to describe the Cooper pair box [18, 19] and transmon [4]
qubits. Today, it is ubiquitous in modelling superconducting circuits [13, 2,
3]. The Hamiltonian can be expressed in the mixed, phase, or charge basis:

$\begin{split}H_{T}&=4E_{c}\hat{m}^{2}-E_{J}\cos\hat{\phi}\\\
&=-4E_{c}\partial_{\phi}^{2}-E_{J}\cos\phi\\\
&=4E_{c}\sum_{m}m^{2}|m\rangle\langle
m|-\frac{E_{J}}{2}\sum_{m}|m\rangle\langle m+1|+\mathrm{H.c.}\end{split}$ (2)

The Hamiltonian is analytically solvable using Mathieu functions [4, 20]. This
description applies to transmons of different types, including tunables ones
(gatemons) [21, 22, 23].

The second line of Eq. (2) describes a particle in the $\phi$-space with an
effective mass $m\propto 1/E_{c}$, trapped in a potential
$V(\phi)=-E_{J}\cos\phi$ [24]. For $E_{c}\ll E_{J}$, corresponding to a very
massive particle, the ground state is localized at the bottom of the potential
well at $\phi=0$ (for $E_{J}>0$, i.e., in a 0-junction), while for $E_{c}\gg
E_{J}$, corresponding to a very low mass, the particle uniformly occupies the
entire $\phi$-range. It is thus localized in the dual $m$-space: the charge
difference across the junction is well defined and there are large
fluctuations of the superconducting phase. In general, the wavefunction has a
finite width in either space, controlled by the ratio $E_{c}/E_{J}$.

The wavefunctions of the excited states are similarly localized, but with an
increasing number of peaks and nodes. This is a general property of a particle
trapped in a potential well.

###  II.3 Embedded QD

The presence of a QD in the junction introduces pair-breaking processes and
lowers the energy of quasiparticle states by binding them to the QD spin [25,
26, 27, 28]. Therefore, the approximate treatment with a model formulated in
terms of Cooper pairs alone is no longer adequate and the full electron
dynamics has to be considered. In order to accurately incorporate the QD
physics into the transmon equation, the basis has to be extended with the QD
degrees of freedom and (this is the key point) one must allow for the presence
of quasiparticles in the SCs.

The QD can alter the sinusoidal Josephson potential in complex ways, and in
turn influence the properties of the transmon excitations. The most direct
example is the $\pi$-junction [29, 30, 31, 32, 33, 34]: if the QD is occupied
by a single electron, transferring a Cooper pair across the junction requires
permuting it over the QD electron. This process produces a fermionic minus
sign, which effectively flips the sinusoidal potential from $-\cos\phi$ to
$+\cos\phi$, so that its minimum is at $\phi=\pi$.

For $E_{c}\ll E_{J}$, a good approximate approach consists of dividing the
problem into two steps. First, one solves the impurity problem of a QD coupled
to SC leads at fixed $\phi$ values to obtain the effective Josephson
potentials $V(\phi)$, one for each eigensolution of the QD problem, and then
uses these potentials as an input for the transmon equation where it replaces
the $\cos\phi$ term [35, 36, 7, 37]. However, this adiabatic approximation
does not capture the possible dynamic coupling between the QD degrees of
freedom contained in $V(\phi)$ and the transmon excitations. In particular,
the approximation is expected to break down when two eigenstates approach each
other or even cross (similar to the breakdown of the Born-Oppenheimer
approximation due to vibronic coupling in molecules).

###  II.4 Reference Josephson junction

For various applications, it is important to have the ability to impose the
phase difference in the ground state [38]. The simplest way to achieve this is
to embed the QD JJ into a larger superconducting circuit, with a second
standard Josephson junction with larger Josephson energy connecting the two
SCs [39, 40, 4, 6, 7]. Here we refer to it as the reference JJ. The phase
difference is then controlled by piercing the resulting loop with a tunable
magnetic flux. We will neglect the loop inductance and the charging energy of
the reference junction, but they could be included in the model if so
required.

The system containing two Josephson junctions has an effective potential [41,
42]

$V(\phi)=V_{\mathrm{QD}}(\phi)-E_{J}^{\mathrm{ref}}\cos(\phi-\phi_{\mathrm{ext}}),$
(3)

with $E_{J}^{\mathrm{ref}}$ the Josephson energy of the reference JJ and
$\phi_{\mathrm{ext}}=2\pi\Phi/\Phi_{0}$ the enforced external phase due to the
magnetic flux $\Phi$, with $\Phi_{0}=h/(2e)$ the magnetic flux quantum. The
effective potential of the QD-junction, $V_{\mathrm{QD}}$, depends on the QD
parameters and is state-dependent. Enforcing the phase is possible if the
Josephson energy of the reference junction is much larger than that of the QD-
junction. Then, the minimum of $V(\phi)$, $\phi_{\mathrm{min}}$, is close to
$\phi_{\mathrm{ext}}$, and thus the ground state $\phi$ tends to approximately
follow $\phi_{\mathrm{ext}}$.

###  II.5 Spin-orbit coupling

Spin-orbit terms couple spin to the supercurrent and break spin degeneracy
even in the absence of the external magnetic field (except at $\phi=0$ and
$\phi=\pi$, where anti-unitary symmetry leads to Kramers degeneracy). In ASQs,
the matrix-valued potential energy for the doublet is given by

$V(\phi)=E_{0}\cos\phi-
E_{\mathrm{SO}}\bm{\sigma}\cdot\mathbf{n}\sin\phi+g\mu_{B}\frac{1}{2}\bm{\sigma}\cdot\mathbf{B},$
(4)

where $\bm{\sigma}$ is the spin operator, $\mathbf{n}$ is a unit vector along
the spin-polarization direction of SOC, $E_{\mathrm{SO}}$ and $E_{0}$ are
spin-dependent and spin-independent Cooper pair tunneling rates [43], $g$ is
the g-factor, $\mu_{B}$ the Bohr magneton, and $B$ the magnetic field. The
$g$-factor is itself $\phi$-dependent due to the impurity Knight shift [44].

##  III Model formulation

In the following we first present a model of a QD embedded between two
superconducting islands, as sketched in Fig. 1. We introduce the flat-band
approximation which is necessary to reduce the Hilbert space and present the
procedure to generate the reduced basis.

Figure 1:  Model sketch. The QD is a single energy level embedded between two
superconductors in the flat-band limit of the Richardson model. The QD is
coupled to active superconducting orbitals $f_{L}$ and $f_{R}$ via single-
electron hopping $v$. Far from the QD the two superconductors form another
auxiliary junction, modelled by the pair-hopping
$t_{p}e^{i\phi_{\mathrm{ext}}}$ terms and represented by the doublet black
arrow. This SQUID-like geometry enables control over the phase difference
$\phi$. The black dot represents a possible quasiparticle occupying the
$f$-orbital.

###  III.1 Model

Our approach is based on describing the superconducting contacts by the
Richardson model [45, 46, 47], a charge conserving pairing Hamiltonian. It is
equivalent to the BCS mean-field theory in the thermodynamic limit [48, 49]
and encapsulates the finite-size effects observed in nanoscopic metallic
grains [50]. Importantly, charge conservation enables a trivial implementation
of the SC charging energy terms and the modeling of coupled QDs [47].

The system consists of two SCs $(\beta=L,R)$, coupled through an embedded QD.
The full model Hamiltonian is the sum of

$\begin{split}H_{\mathrm{QD}}&=\epsilon
n_{\mathrm{QD}}+Un_{\mathrm{QD}\uparrow}n_{\mathrm{QD}\downarrow}=\frac{U}{2}(n_{\mathrm{QD}}-\nu)^{2}+\mathrm{const.},\\\
H_{\mathrm{SC}}^{(\beta)}&=\sum_{i}\epsilon_{i}n_{\beta
i}+g\frac{1}{N}\sum_{ij}c^{\dagger}_{\beta i\uparrow}c^{\dagger}_{\beta
i\downarrow}c_{\beta j\downarrow}c_{\beta j\uparrow}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad+E_{c}^{(\beta)}\left(n_{\beta}-n_{0}^{(\beta)}\right)^{2},\\\
H_{\mathrm{hyb}}^{(\beta)}&=\frac{v_{\beta}}{\sqrt{N}}\sum_{i,\sigma}\left(d^{\dagger}_{\sigma}c_{\beta
i\sigma}+c_{\beta i\sigma}^{\dagger}d_{\sigma}\right),\\\
H_{\mathrm{ref}}&=t_{p}e^{i\phi_{\mathrm{ext}}}\frac{1}{N}\sum_{i}c_{Li\uparrow}^{\dagger}c_{Li\downarrow}^{\dagger}\sum_{j}c_{Rj\downarrow}c_{Rj\uparrow}+\mathrm{H.c.}\end{split}$
(5)

Here $n_{\mathrm{QD}\sigma}=d^{\dagger}_{\sigma}d_{\sigma}$ is the QD number
operator and
$n_{\mathrm{QD}}=n_{\mathrm{QD}\uparrow}+n_{\mathrm{QD}\downarrow}$.
$\epsilon$ is the QD level and $U$ the on-site Coulomb repulsion. In the
alternative formulation of $H_{\mathrm{QD}}$, $U/2$ plays the role of the QD
charging energy and $\nu=\frac{1}{2}-\frac{\epsilon}{U}$ is its optimal
occupation in the units of particle number (gate charge).

$\epsilon_{i}$ is the set of $N$ equally spaced energy levels that represent
the SC, with $c_{\beta i\sigma}$ the corresponding annihilation operators, and
$n_{\beta i}$ the number operators. $g$ is the superconducting pairing
strength [51]. $E_{c}^{(\beta)}$ is the SC charging energy with
$n_{\beta}=\sum_{i}n_{\beta i}$ the total SC charge operator and
$n_{0}^{(\beta)}$ the optimal occupation in units of particle number. The QD
and SCs are coupled by single particle hopping with strength $v_{\beta}$.

The $H_{\mathrm{ref}}$ term describes the reference JJ. Because it is far from
the QD, it can be safely assumed that it does not contain quasiparticles, and
can thus be described by the pair-hopping processes. $t_{p}$ is the pair
hopping strength and $\phi_{\mathrm{ext}}$ the phase difference externally
imposed by the magnetic flux that pierces the superconducting loop.

We consider symmetric parameters in both superconductors and thus drop the
subscript $\beta$ in the parameter symbols.

The Hamiltonian has a compact matrix product operator (MPO) representation and
can in principle be solved with the density matrix renormalization group
(DMRG) [52]. However, excitations in the regime we are interested in here,
with $E_{c}$ much smaller than the superconducting gap $\Delta$, present
certain challenges.

The reason is hidden in the nature of the lowest excitations of Eq. (5).
Consider the uncoupled limit with $v=0$, and $E_{c}\rightarrow 0$. The lowest
excitations are due to charge redistribution: if the ground state
configuration (electron occupations of constituent parts) is
$(n_{L},n_{\mathrm{QD}},n_{R})$, there exists an otherwise equivalent
configuration with a single shifted Cooper pair,
$(n_{L}-2,n_{\mathrm{QD}},n_{R}+2)$. Its excitation energy ($2d$ for moving a
Cooper pair to a higher single-particle level) is only due to the finite level
spacing, which is a consequence of the finite size of the system. In the
thermodynamic limit ($N\rightarrow\infty$, $d\rightarrow 0$) such states are
degenerate and at finite $v$ they form equal superpositions with a well-
defined phase difference $\phi$. In numerical calculations with a finite
system, however, the finite excitation energy precludes a clear emergence of
$\phi$. (See Sec. II D of Ref. 52 for a detailed description of the problems
and limitations of the full model.)

The effect of the level spacing $d$ on the Josephson effect in the Richardson
model was studied in the context of nanoscopic metallic grains in Ref. 14.
Interestingly, they found that as $d$ is increased from zero, the effective
$E_{J}$ at first decreases as finite $d$ acts like a charging energy term,
incurring a small energy penalty for each charge transferred. Additionally,
increasing $d$ increases the pair-hopping matrix elements. At values of $d$
comparable to the superconducting gap $\Delta$ the second effect dominates,
and $E_{J}$ actually increases beyond the BCS value found at $d\rightarrow 0$.

However, the goal of this work is to model superconductors in the regime where
the level spacing is negligible and the charging energy comes from other
effects, such as the capacitance of the constituent parts of the device.

###  III.2 Flat-band approximation

The desired behavior—the emergence of a well-defined $\phi$ from an equal
superposition of degenerate $|m\rangle$ states—can be recovered by
disregarding the kinetic energy of the SC levels by setting all $\epsilon_{i}$
to zero. This is the flat-band approximation. We have shown in Ref. 51 that
employing this approximation does not importantly change the low-energy
physics of a QD-SC system, so that the results remain qualitatively correct.
Simultaneously, the model simplifies to the point where certain limits are
analytically solvable. Here, this treatment is extended to the SC-QD-SC setup.
In the following we present simplified expressions valid in the limit of a
large number of SC levels for half-filled bands. 111For detailed derivations
and expressions at general filling see Ref. 51.

Making all SC levels nominally equivalent allows us to define a single active
orbital in each SC:

$f_{\beta\sigma}=\frac{1}{\sqrt{N}}\sum_{i}c_{\beta i\sigma}.$ (6)

This greatly simplifies the hybridisation term, which now only involves the
$d$ and $f$ orbitals:

$H_{\mathrm{hyb}}^{(\beta)}=v_{\beta}\sum_{\sigma}\left(d^{\dagger}_{\sigma}f_{\beta\sigma}+f^{\dagger}_{\beta\sigma}d_{\sigma}\right).$
(7)

Because the kinetic energy term is zero, the $H_{\mathrm{SC}}$ terms simplify
as well. After subtracting the condensation energy of $m_{\beta}$ pairs
($\frac{g}{2}m_{\beta}$), they read:

$\begin{split}H_{\mathrm{SC}}^{(\beta)}=&\frac{g}{2}\sum_{\sigma}f^{\dagger}_{\beta\sigma}f_{\beta\sigma}+\\\
&\quad
E_{c}^{(\beta)}\left[\sum_{\sigma}f^{\dagger}_{\beta\sigma}f_{\beta\sigma}+2m_{\beta}-n_{0}\right]^{2}.\end{split}$
(8)

The pairing term simply counts the number of quasiparticles occupying the
$f$-orbitals, each contributing $\frac{g}{2}$ to total energy. The SC gap
$\Delta$ is thus proportional to pairing strength, which is a well-known
feature of the flat-band systems [54, 55, 56]. Note that the approximation
retains the information about the number of Cooper pairs in each SC,
$m_{\beta}$. A state of the SC island in the flat-band approximation is thus
fully determined by the state of the $f$-orbital (0, $\uparrow$, $\downarrow$,
$2=\downarrow\uparrow$) and $m_{\beta}$, the number of Cooper pairs in the
condensate.

As all SC levels are equivalent, the condensate is spread equally across all
of them. The ground state of $H_{\mathrm{SC}}^{\beta}$ containing $m_{\beta}$
pairs is:

$|m_{\beta},0\rangle=\mathcal{N}\left(\sum_{i=1}^{N}c^{\dagger}_{\beta
i\downarrow}c^{\dagger}_{\beta i\uparrow}\right)^{m_{\beta}}|0\rangle,$ (9)

with normalization $\mathcal{N}$ [51] and electronic vacuum $|0\rangle$. The
second label in the ket of the left hand side denotes the state of the
$f$-orbital.

Excitations containing quasiparticles are

$|m_{\beta},\sigma\rangle=\frac{1}{\sqrt{N}}\sum_{b=1}^{N}c^{\dagger}_{\beta
b\sigma}\mathcal{N}\left(\sum_{i\neq b}^{N}c^{\dagger}_{\beta
i\downarrow}c^{\dagger}_{\beta i\uparrow}\right)^{m_{\beta}}|0\rangle.$ (10)

where the quasiparticle blocks the $b$-th level from participating in pairing
[50], while the rest contain the pair condensate. Similarly

$|m_{\beta},2\rangle=\frac{1}{N}\sum_{bb^{\prime}}^{N}c^{\dagger}_{\beta
b\downarrow}c^{\dagger}_{\beta b^{\prime}\uparrow}\mathcal{N}\left(\sum_{i\neq
b,b^{\prime}}^{N}c^{\dagger}_{\beta i\downarrow}c^{\dagger}_{\beta
i\uparrow}\right)^{m_{\beta}}|0\rangle,$ (11)

with $b$ and $b^{\prime}$ blocked. The superconducting nature of the
quasiparticles occupying the $f$-orbitals is contained in their algebra, given
by [51]:

$\begin{split}f^{\dagger}_{\beta\sigma}|m_{\beta},0\rangle&=\frac{1}{\sqrt{2}}|m_{\beta},\sigma\rangle,\\\
f^{\dagger}_{\beta\downarrow}f^{\dagger}_{\beta\uparrow}|m_{\beta},0\rangle&=\frac{1}{2}\left(|m_{\beta},2\rangle+|m_{\beta}+1,0\rangle\right).\end{split}$
(12)

Importantly, the second equation shows the possibility of recombination of two
quasiparticles in the orbital $f$ into a Cooper pair. This is precisely the
process involved in the Cooper pair transfers in a JJ.

###  III.3 Full active orbital set

The full basis is generated as the tensor product of eigenstates of each
subsystem for zero coupling. The basis states are determined by the three
active orbitals $d,f_{L},f_{R}$ and quantum numbers $m_{L}$ and $m_{R}$.
Additionally, after fixing the total charge in the system to $n$, the
$(m_{L},m_{R})$ pair can be replaced by the difference, $m=m_{L}-m_{R}$. We
thus denote the basis states by

$d\leavevmode\nobreak\ |m,f_{L},f_{R}\rangle$ (13)

with an operator
$d=\mathds{1}_{d},d^{\dagger}_{\uparrow},d^{\dagger}_{\downarrow},d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}$,
and labels $f_{L/R}=0,\uparrow,\downarrow,2$. These contain all states where
the quasiparticles are coupled to the QD [51]. This step beyond basis sets
that only involve Cooper pairs [57] is necessary in the presence of magnetic-
impurity-induced pair-breaking processes.

Furthermore, our model conserves the total spin $S$. We consider the singlet
($n=\mathrm{even}$, $S=0$) and the doublet ($n=\mathrm{odd}$, $S=1/2$,
$S_{z}=+1/2$) sectors. Because Cooper pairs are singlets, $S$ is solely
determined by the state of the three active orbitals. For each $m$, there are
14 possible states in each sector. (See App. A for their definitions.)

We generate the matrix representation of the Hamiltonian by computing symbolic
expressions for all matrix elements for general $[m_{L},m_{R}]$ using symbolic
algebra software [58] and then selecting the allowed $m$ configurations for a
given total charge $n$. The basis size grows linearly with $n$, which enables
efficient exact diagonalization for $n$ in the thousands. The computer code
for performing numerical calculations is available in a public repository
[59].

###  III.4 Fourier transform to the phase-space

The phase difference $\phi$ and the difference in Cooper pair number $m$ are
conjugate quantities [13], i.e., $[\phi,m]=i$, and thus related via Fourier
transform. When performed on the $m$-states, we obtain a basis in the
$\phi$-space:

$|\phi,f_{L},f_{R}\rangle=\frac{1}{\sqrt{m_{\mathrm{tot}}+1}}\sum_{m}e^{i\phi\left(m+m_{\mathrm{tot}}\right)/2}|m,f_{L},f_{R}\rangle,$
(14)

with $m_{\mathrm{tot}}$ the total number of Cooper pairs in the system. $m$
runs across all possible configurations,
$m=-m_{\mathrm{tot}},-m_{\mathrm{tot}}+2,\dots,m_{\mathrm{tot}}-2,m_{\mathrm{tot}}$;
there are $m_{\mathrm{tot}}+1$ such configurations. The multiplicative factor
of $1/2$ in the exponent is necessary because of the steps of 2. Furthermore,
we shift $m$ by $m_{\mathrm{tot}}$ for convenience. The phase $\phi$ is
defined in the interval $[0,2\pi)$:

$\phi=2\pi\frac{l}{m_{\mathrm{tot}}+1},$ (15)

with $l=0,\ldots,m_{\mathrm{tot}}$. In the thermodynamic limit, $N\to\infty$,
$\phi$ becomes a periodic continuous variable [60].

For $E_{c}=0$ and neglecting finite-size effects, $\phi$ is a good quantum
number. Therefore, in this case the Hamiltonian—a block tridiagonal matrix in
the $m$-basis—has a block-diagonal shape when expressed in the $\phi$-basis,
with the matrix elements inside the block describing the internal degrees of
freedom associated with the full active orbital set (unpaired particles), see
Fig. 2. These are $14\times 14$ matrices, with the phase variable $\phi$
appearing in the blue blocks of hopping matrix elements, similar to what
happens in the Bloch state basis in the tight-binding description of electrons
on a lattice. We remind the reader that the matrices are different in the
singlet and doublet sectors. By diagonalizing them we obtain the eigenstates
in the $\phi$-basis.

Figure 2:  Block diagonalizing the Hamiltonian at $E_{c}=0$ with the Fourier
transform. In the analogy with the tight-binding chain, the $t_{2}$ block
denotes coupling within each chain site, while $t_{1}$ is hopping between the
sites.

###  III.5 Charge vs. phase basis: Similarities and differences compared to
the zero-bandwidth BCS

The model resembles the zero-bandwidth BCS approximation (ZBW BCS), which is a
popular tool for qualitatively describing the physics of hybrid
superconducting QD devices [61, 31, 62, 63, 64, 65, 66]. Indeed, in the
simplest case of a QD coupled to a single SC, the two models are actually
mathematically equivalent [51]. However, we argue that the flat-band
Richardson formulation presented here is better suited for modelling extended
quantum devices with possibly complex geometries, and where charging energies
of the constituent parts play an important role.

An important aspect is that the degrees of freedom associated with the Cooper
pair condensate do not appear explicitly in the ZBW BCS, while in the flat-
band Richardson model we keep track of them with the labels $m_{\beta}$. While
this is not important for a single SC, the condensate plays a key role when
two (or more) SCs are coupled; as discussed just above, the transfers of
Cooper pairs (changes of $m=m_{L}-m_{R}$ by $\pm 2$) generate the Josephson
effect and lead to the emergence of the phase difference $\phi$. This is a
statement equivalent to the notion that the phase of a single isolated SC is
arbitrary (gauge freedom), but the phase difference between two coupled SCs is
physically measurable and gauge invariant.

These details are not immediately obvious in the BCS model which is formulated
in the phase basis, with $\phi$ taking a fixed value due to the gauge symmetry
breaking. This becomes an issue, however, when considering non-zero charging
energy. In that case $\phi$ is not a good quantum number, and thus cannot be
represented by a constant value, but needs to be promoted to an operator
quantity in order to allow for the phase fluctuations in the model. The
problem is naturally resolved in the charge conserving formalism of the
Richardson model, where $\phi$ is an emergent global degree of freedom,
arising as a property of a superposition of many $|m\rangle$ states. This
means that the model is able to accurately capture its quantum fluctuations at
finite $E_{c}$.

Alternatively, it is possible to describe the problem in the $\phi$-basis,
while also retaining information about the pair condensate. Indeed, this is
achieved by the Fourier transform in Eq. (14). In that case the charging
energy terms appear as off-diagonal matrix elements which couple different
$|\phi\rangle$-states.

###  III.6 Magnetic field and spin-orbit coupling

The method can be easily extended to include additional effects. For example,
the magnetic field on the QD along the $z$ axis is included by substituting
$\epsilon\rightarrow\epsilon_{\sigma}=\epsilon+\sigma E_{Z}/2$, where $E_{Z}$
is the Zeeman splitting of the QD level, $E_{Z}=g\mu_{B}B$, and $\sigma=\pm
1$, depending on the electron spin. The $x$ and $y$ components of the magnetic
field appear as off-diagonal elements coupling the $S_{z}=1/2$ and
$S_{z}=-1/2$ basis states.

The spin-orbit coupling can be included in similar way as in Ref. 43, by
adding two terms: a spin-flip SC-QD hopping $v_{\uparrow\downarrow}$ and a SC-
SC single-electron hopping $t_{\mathrm{sc}}$ (see also the Supplementary
material of Ref. 43):

$\begin{split}\sum_{\sigma}\left(iv_{\uparrow\downarrow}d^{\dagger}_{\sigma}f_{L\bar{\sigma}}+iv_{\uparrow\downarrow}f^{\dagger}_{R\sigma}d_{\sigma}+\mathrm{H.c.}\right)\\\
+\sum_{\sigma}\left(t_{\mathrm{sc}}f_{L\sigma}^{\dagger}f_{R\sigma}+\mathrm{H.c.}\right).\end{split}$
(16)

Here $\bar{\sigma}$ indicates the reversed spin index. With this definition,
the spin-polarization direction of SOC is along the $x$-axis. The field along
the $x$ axis is said to be “parallel”, while fields along the $y$ and $z$ axes
are said to be “perpendicular”.

##  IV Benchmark results

In this section, we demonstrate that the presented model reproduces the
expected results in appropriate limits. We also point out the existence of
parameter regimes where the generalized formalism reveals less familiar
behavior.

###  IV.1 Role of quasiparticles

In our formalism, four electron hopping events are necessary for a Cooper pair
to pass through the QD, and thus we expect $E_{J}$, the half-width of the
cosine-like excitation band, to increase proportionally to $v^{4}$. However,
$\phi$-dependent processes of second order in $v$, when present, will dominate
the dynamics of the junction in the perturbative regime. Here we show that
this happens in the presence of a superconducting quasiparticle and identify
the parameter ranges where such contributions become large due to the presence
of the QD. In such cases the treatment with an effective pair-hopping
Hamiltonian is not adequate, and solving the QD-JJ Hamiltonian on the level of
single-electron processes is necessary.

####  IV.1.1 Subgap spectrum of the SC-QD-SC junction

Figure 3:  Spectra of the SC-QD-SC junction for $E_{c}=0$, $t_{p}=0$. (a, b)
$\phi$ dependence of the lowest singlet (red) and the lowest doublet (blue)
states for (a) the resonant case of $\epsilon=0$ and (b) for the generic case
of $\epsilon=0.1\Delta$. Here $U=0$. (c) $v$ dependence of the two lowest
singlets and the lowest doublet for $\phi=0$ (dashed) and $\phi=\pi$ (solid),
for $U=0$. (d) $U$ dependence of the spectra at $\phi=0$ close to the
particle-hole symmetric point, $\epsilon=-0.8U/2$. The energy zero coincides
with the energy of the lowest doublet. Here we use
$v/\Delta=0.2(1+\sqrt{U/\Delta})$, because the main goal is to study the
effect of the $U/\Delta$ ratio, while keeping the spin exchange interaction
$v^{2}/U$ constant in the large-$U$ regime.

We begin by reproducing the basic properties of a JJ with an embedded QD,
starting with $t_{p}=0$, i.e., without the reference junction, and $E_{c}=0$,
i.e., without any charging energy terms. In this simplest case the phase
difference $\phi$ is a conserved scalar quantity (good quantum number) that
labels the eigenstates of the system.

We start by showing the $\phi$ dependence of the two lowest singlets and the
lowest doublet for a non-interacting ($U=0$) QD in Fig. 3(a,b). The SC leads
induce superconducting pairing on the QD level, which results in a
$\phi$-dependent splitting of the states according to the occupation of the
QD. This is known as the proximity effect. In the doublet state the QD is
occupied by a single electron, while in the singlet states it contains a
Cooper pair or two Bogoliubov quasiparticles (broken Cooper pair),
respectively. In the electron occupation basis, these are represented as
orthogonal equal superpositions of $|0\rangle$ and $|2\rangle$ at the QD. (For
$\phi=0$, we have $|0\rangle+|2\rangle$ for a Cooper pair, and
$|0\rangle-|2\rangle$ for a broken pair.)

The resonant situation ($\epsilon=0$), Fig. 3(a), shows some particularities.
The singlets are split through a second order $v^{2}$ process where a virtual
quasiparticle hops from the QD to a SC orbital and back. This SC orbital is
equally distributed across the two SCs (assuming symmetric $v=v_{L}=v_{R}$),
with the phase between the contributions of the two leads corresponding to the
$\phi$-dependent symmetry of the system [67]. This results in a
$\phi$-dependent correction of energy $(2v^{2}/\Delta)\cos(\phi/2)$ [31]. The
resulting $4\pi$ periodicity of the energy dispersion in the singlet state and
the quadratic dependence on $v$ are the main signatures of resonant behavior.

In Fig. 3(b) we show the generic (non-resonant, particle-hole asymmetric) case
with non-zero $\epsilon$. This regime corresponds to a standard JJ with subgap
Andreev bound states (ABS) [68, 13]. In comparison to the resonant case shown
in Fig. 3(a), here the particle-hole asymmetry splits the singlets at
$\phi=\pi$. The dispersion is then $2\pi$-periodic. The ground state has a
larger contribution of the empty QD state $|0\rangle$ (for $\epsilon>0$),
while the excited singlet state has a mostly doubly occupied QD, state
$|2\rangle$. The previously discussed second order process is energetically
suppressed, for it requires the formation of Cooper pairs – equal
contributions of $|0\rangle$ and $|2\rangle$. The $\phi$-dependence of singlet
states is thus dominated by the next higher order ($\propto v^{4}\cos\phi$)
contribution, coming from coherent transport of two electrons—a Cooper
pair—across the junction. These are the processes that generate the Josephson
supercurrent, and we find the Josephson energy $E_{J}\propto v^{4}$.

Figure 3(c) shows the $v$-dependence of the ABS, illustrating the increase of
$E_{J}$ that is quantified by the spread between the dashed ($\phi=0$) and
solid ($\phi=\pi$) lines. Note the significant difference between the doublet
(blue) and singlet (red) states.

In Fig. 3(d) we show the evolution of the spectra with increasing $U$ while
staying close to the particle-hole symmetric point by setting
$\epsilon=-0.8U/2$, so that the average QD occupancy is approximately one. In
the lowest doublet state the QD contains a single spin already at $U=0$. This
state remains largely unperturbed for increasing $U$. In contrast, the singlet
states strongly depend on the interaction strength. The lowest singlet state
gradually transforms from the ABS at $U\rightarrow 0$ to a Yu-Shiba-Rusinov
(YSR) state for $U\gg\Delta$, where the electron in the singly-occupied QD
forms a bound state with a Bogoliubov quasiparticle in the superconductor
through the exchange interaction. The first excited singlet retains large
contributions of the $|0\rangle$ and $|2\rangle$ QD states and is thus quickly
pushed to higher energies with increasing charge repulsion $U$ [52].

####  IV.1.2 Effective Josephson energy

In analogy with a standard JJ, we define the effective Josephson energy as the
half-width of the energy band obtained by varying $\phi$,

$E_{J}^{\mathrm{eff}}=\frac{1}{2}\left(\max_{\phi}E(\phi)-\min_{\phi}E(\phi)\right),$
(17)

where $E(\phi)$ is the energy of the lowest state in the singlet or doublet
subspace (the result is state-dependent). The extreme values are taken at
$\phi=0$ and $\phi=\pi$ when $t_{p}=0$.

Fig. 4(a) shows a log-log plot of the $v$ dependence of $E_{J}^{\mathrm{eff}}$
at small $\epsilon=10^{-3}$, with $U=0$. For the doublet ground state (blue)
we find $E_{J}^{\mathrm{eff}}\propto v^{4}$ for all $v$, as the spin in the QD
is not coupled to the superconducting leads and the transfer of Cooper pairs
across the junction remains the leading contribution for $v<\Delta$. (The
presence of the QD spin does however change the prefactor of the Josephson
current, a phenomenon known as the $\pi$-junction.)

However, in the singlet ground state (red) we find two regimes, with
$E_{J}^{\mathrm{eff}}\propto v^{2}$ for large $v$ and
$E_{J}^{\mathrm{eff}}\propto v^{4}$ for small $v$. This is a consequence of
the competition between $\epsilon$ splitting the $|0\rangle$ and $|2\rangle$
QD states and the proximity effect hybridizing them into Cooper pairs. The
crossover occurs at $\epsilon=2v^{2}/\Delta$ (vertical black line). The right
hand side of the equality corresponds to the hybridisation coming from the
proximity effect. The regime where $\epsilon<\frac{2v^{2}}{\Delta}$
qualitatively corresponds to Fig. 3(a), while an example of the
$\epsilon>\frac{2v^{2}}{\Delta}$ case is shown in Fig. 3(b)

Fig. 4(b) shows the $\epsilon$ dependence of $E_{J}^{\mathrm{eff}}$. For the
singlet state (red) we see the evolution from the $\epsilon$-independent
$v^{2}$ regime to the $v^{4}$ dependence as discussed above. Interestingly, by
increasing $\epsilon$ beyond $\Delta$, we observe $\propto v^{2}$ behavior in
the doublet state (blue) as well. This is a regime where it becomes favorable
for the unpaired electron from the QD to enter the SC as a quasiparticle.
Similarly to the second order process that causes the proximity effect in the
singlet sector, here the quasiparticle occupies an orbital symmetrically
distributed across the two SCs, and the second order $\phi$-dependent
correction comes from processes where the quasiparticle traverses the junction
(or conversely, from the $\phi$-dependent shape of the orbital).

We have thus established a generic property of the system that can only be
captured by accounting for single electron hopping processes: the presence of
a SC quasiparticle always leads to a dominant second-order contribution to
$E_{J}^{\mathrm{eff}}$. This is further corroborated in Fig. 4(c,d). In the
singlet ground state, Fig. 4(c), a quasiparticle is induced by increasing $U$
beyond $\Delta$, transforming the nature of the ground state from ABS to YSR.
Correspondingly, $E_{J}^{\mathrm{eff}}(v)$ changes from $v^{4}$ for $U<\Delta$
to $v^{2}$ for $U>\Delta$. Similarly, a quasiparticle is favored in the
doublet ground state if $\epsilon>\Delta$ (for $\epsilon>0$), or
$\epsilon>\Delta+2\epsilon+U$ ($\epsilon<0$). Fig. 4(d) shows the negative
$\epsilon$ case, where the crossover from $v^{4}$ to $v^{2}$ behavior occurs
when $|\epsilon|>\Delta+U$.

Figure 4:  Effective Josephson energy $E_{J}^{\mathrm{eff}}$, defined as the
difference in the energy of the lowest state at $\phi=0$ and $\phi=\pi$. (a)
and (b) $U=0$ results. Full lines: singlet (red) and doublet (blue); black
dashed lines: expected perturbative corrections. (c) $E_{J}^{\mathrm{eff}}$
for the singlet ground state for different $U$ with $\epsilon=-0.8(U/2)$,
sightly shifted away from the particle-hole symmetric point. (d)
$E_{J}^{\mathrm{eff}}$ for the doublet ground state with increasingly negative
$\epsilon$, at $U/\Delta=0.2$. $-\epsilon=\Delta+U$ is the point where a
quasiparticle becomes present in the superconductor.

###  IV.2 Reference Josephson junction: $t_{p}$ dependence

Next, we investigate the effect of the reference Josephson junction in the
SQUID geometry as a means of enforcing the phase bias
$\phi=\phi_{\mathrm{ext}}$ across the QD JJ using external magnetic flux [42].
The Josephson energy of the reference junction is given by
$E_{J}^{\mathrm{ref}}=2t_{p}$. The effective Josephson energy of the isolated
QD junction is denoted $E_{J}^{\mathrm{eff}}(0)$ and corresponds to
$E_{J}^{\mathrm{eff}}$, defined in Eq. (17), at $t_{p}=0$. The total potential
is the sum of the QD-induced potential and the
$-E_{J}^{\mathrm{ref}}\cos(\phi-\phi_{\mathrm{ext}})$ term, as in Eq. (3). We
note that this setup corresponds to the ASQ embedded in a transmon that was
experimentally explored in Refs. 36, 43, 7.

We set $U=3\Delta$ and remain close to the particle-hole symmetric point, a
parameter regime relevant for the ASQ, where the ground state is a doublet
with a singly occupied QD, $E_{J}^{\mathrm{eff}}\propto v^{4}$, and $V(\phi)$
has a minimum at $\phi_{\mathrm{min}}=\pi$. Figs. 5(a-c) show the spectra with
increasing $t_{p}$ for $\phi_{\mathrm{ext}}=0.4\pi$. The phase difference in
the ground state, $\phi_{\mathrm{min}}$, depends on the competition between
the reference and the QD junction, which we quantify by the ratio

$r=\frac{E_{J}^{\mathrm{ref}}}{E_{J}^{\mathrm{eff}}(0)}.$ (18)

Fig. 5(d) shows the $r$ dependence of $\phi_{\mathrm{min}}$ as $t_{p}$ is
increased, and the systems transitions from a QD-induced $\pi$-junction to the
one with enforced $\phi_{\mathrm{min}}=\phi_{\mathrm{ext}}$. There is weak
$\phi_{\mathrm{ext}}$ dependence of $\phi_{\mathrm{min}}$ at finite $r$: the
phase difference is not equally strongly enforced for all
$\phi_{\mathrm{ext}}$.

This is explicitly shown in Fig. 5(e), where we plot the location of the
potential minimum $\phi_{\mathrm{min}}$ with changing $\phi_{\mathrm{ext}}$
for various values of $r$. For experimentally relevant $r\sim 5-10$ we find
relatively small deviation from the optimal
$\phi_{\mathrm{ext}}=\phi_{\mathrm{min}}$ situation (black dashed line).

The depth of the effective total potential $V(\phi)$, defined in Eq. (3),
depends on $\phi_{\mathrm{ext}}$. This is shown in Fig. 5(f), where we plot
the $r$ dependence of the effective Josephson energy $E_{J}^{\mathrm{eff}}$
from Eq. (17). $V(\phi)$ is a sum of two sinusoidal terms [see Eq. (3)], one
coming from the QD with a state-dependent phase factor ($0$ for singlet, $\pi$
for doublet) and the other from the auxiliary junction with the phase
$\phi_{\mathrm{ext}}$. When the two phases are aligned at
$\phi_{\mathrm{ext}}=\pi$, the amplitude of $V(\phi)$ is the sum of the two
amplitudes (upper black dashed line). On the other hand, when the two
junctions are out of phase at $\phi_{\mathrm{ext}}=0$ their contributions
subtract (lower black dashed line). At $r=1$ this leads to a completely flat
$V(\phi)$ (zero $E_{J}^{\mathrm{eff}}$) despite the Josephson energies of the
two junctions possibly being very large. Any finite value of $E_{c}$ close to
this point would cause large fluctuations of $\phi$.

Figure 5:  Effect of the tunneling amplitude of Cooper pairs across the
reference junction, $t_{p}$. (a-c) $\phi$-dependence of the spectra for three
increasing values of $t_{p}$, at fixed $\phi_{\mathrm{ext}}=0.4\pi$ (indicated
with vertical lines). The black dashed curves show the Josephson potential of
the auxiliary Josephson junction, $-2t_{p}\cos(\phi-\phi_{\mathrm{ext}})$. (d)
Location of the minimum of the ground state for different
$\phi_{\mathrm{ext}}$ as a function of the ratio of the Josephson energy of
the auxiliary and the QD-junction
$r=E_{J}^{\mathrm{ref}}/E_{J}^{\mathrm{eff}}(0)$. Dashed black lines indicate
the corresponding $\phi_{\mathrm{ext}}$. (e) Position of $\phi_{\mathrm{min}}$
vs. $\phi_{\mathrm{ext}}$ for different values of $r$. (f) The black dashed
line corresponds to $\phi_{\mathrm{min}}=\phi_{\mathrm{ext}}$.
$E_{J}^{\mathrm{eff}}$ vs. $r$, with $\phi_{\mathrm{ext}}$ ranging from $0$
(orange) to $\pi$ (blue) in steps of $0.2\pi$. The top black dashed line
corresponds to the sum of contributions from both junctions,
$E_{J}^{\mathrm{eff}}(0)+E_{J}^{\mathrm{ref}}$, while the bottom one is the
absolute difference
$\left|E_{J}^{\mathrm{eff}}(0)-E_{J}^{\mathrm{ref}}\right|$. Other parameters:
$U/\Delta=3$, $\epsilon=-0.9\frac{U}{2}$, $v/\Delta=0.4$, $E_{c}=0$.

###  IV.3 QD-transmon: charging energy effects

After demonstrating that the model reproduces a number of standard properties
of a QD JJ, we turn our attention to the case of finite charging energy.
$E_{c}$ induces mixing between the states with a different $\phi$, and it is
no longer possible to find eigenstates in each $\phi$ subspace separately. We
thus use the full basis as defined in Sec. III.3 and App. A.

The effect of $E_{c}$ on the spectrum is shown in Fig. 6(a). With increasing
$E_{c}$ the eigenstates lose their well-defined phase and gradually transform
into states with a well-defined $m$. This is reflected in the disassociation
of the band of width $2E_{J}^{\mathrm{eff}}$ into discrete states with
excitation energy growing as $4E_{c}m^{2}$. In agreement with the effective
transmon model (Eq. 2), the crossover occurs over a protracted range of
$E_{c}$ values, centered at $E_{c}/E_{J}^{\mathrm{eff}}\sim 10^{-3}$ .

The nature of the eigenstates can be gauged by tracing over the active orbital
degrees of freedom and plotting the amplitudes in the $m$-basis ($\alpha_{m}$)
and in the Fourier transformed $\phi$-basis ($\alpha_{\phi}$), see Fig. 6(c).
In the $E_{c}\rightarrow 0$ limit the eigenvectors tend to a $\delta$ peak in
the $\phi$-space, corresponding to a wide distribution in the dual $m$-space.
The opposite is found in the $E_{c}\to\infty$ regime. In the crossover regime
the distribution of $\alpha$ exhibits substantial width in both basis spaces.
This indicates that neither $\phi$ nor $m$ are well defined, and there is no
privileged basis for the description of eigenstates.

In Fig. 6(b), we show $|\alpha_{\phi}|^{2}$ for the doublet states in the
transmon regime, exhibiting an increasing number of peaks and nodes in higher-
lying excitations. These curves can be interpreted as effective wavefunctions
of the transmon degree of freedom, exhibiting the expected characteristics of
a particle trapped in a potential well.

Figure 6:  Effect of the charge repulsion $E_{c}$. (a) Evolution of the
excitation spectrum with increasing $E_{c}$ in the doublet (blue) and the
singlet (red) sector. (b) Absolute values of the amplitudes in the
$\phi$-basis for the ground state and first two excitations in the doublet
sector in the intermediate regime, $E_{c}/(E_{J}^{\mathrm{eff}}/2)=10^{-3}$.
(c) Absolute values of amplitudes in the $\phi$-basis ($\alpha_{\phi}$) and
$m$-basis ($\alpha_{m}$) for the doublet ground state for different values of
$E_{c}$. (d) Variance of $m$, $\mu=\langle m^{2}\rangle-\langle m\rangle^{2}$
in the doublet ground state vs. $E_{c}/E_{J}^{\mathrm{eff}}$ at
$\phi_{\mathrm{ext}}=\pi$. The black vertical line corresponds to
$8E_{c}=E_{J}^{\mathrm{eff}}$. $E_{J}^{\mathrm{eff}}(0)/\Delta=0.027$,
$r=1.1$. Other parameters are $U/\Delta=3$, $\epsilon=-0.8\frac{U}{2}$,
$v/\Delta=0.5$, $N=601$ and $n_{0}=300$. In (a-c) $t_{p}/\Delta=0.1$,
$\phi_{\mathrm{ext}}=\pi$.

The quantum fluctuations of $m$ reflect the competition between the charging
energy $E_{c}$, favoring states with a well defined $m$, and the effective
potential $V(\phi)$, favoring a state with a well defined $\phi$. In Fig. 6(d)
we plot the variance of $m$ in the doublet ground state,

$\mu=\langle m^{2}\rangle-\langle m\rangle^{2},$ (19)

for different values of $U$. $\mu$ decreases with increasing $E_{c}$, as a
certain charge distribution becomes favored and the probabilities of other $m$
values are reduced. Plotting in log-log scales uncovers two regimes. At small
$E_{c}$ we find roughly $\mu\propto 1/\sqrt{E_{c}}$, independent of the value
of $U$. In this regime the charging energy provides a large effective mass
$m^{*}$ in the harmonic-oscillator picture, $m^{*}\propto 1/E_{c}$, which then
leads directly to $\mu\propto m^{*}\omega\propto\sqrt{m^{*}}\propto
1/\sqrt{E_{c}}$.

The inflection point is at $8E_{c}=E_{J}^{\mathrm{eff}}$ (dashed vertical
line). Here $8E_{c}$ corresponds to the charging energy penalty of the excited
state with $m=\pm 2$, obtained by transferring one Cooper pair across the
junction; the sum of the $E_{c}n^{2}$ terms in each SC gives
$E_{c}\left(2^{2}+(-2)^{2}\right)=8E_{c}$. Notably, at this point $\mu\sim 1$,
indicating a close to equal superposition of the $|m=0\rangle$ and $|m=\pm
2\rangle$ states.

At larger $E_{c}$, only the $|m=0\rangle$ state remains populated.
Furthermore, the electrostatic effects in the QD start to play a role, as seen
by the dependence of $\mu$ on $U$. For large $U$, the fluctuations of charge
on the QD are strongly prohibited, and thus $\mu$ decreases much faster with
increasing $E_{c}$ than for smaller values of $U$ [47].

###  IV.4 Spin-orbit splitting of the doublet states

Figure 7:  Spin-orbit splitting of the doublet states at zero magnetic field.
We plot (a) the eigenenergies $E$ of the states $|0,\sigma\rangle$ and
$|1,\sigma\rangle$, where $\sigma$ is the helicity index, and (b) the
corresponding transition frequencies $\omega$ between these states. The line
styles and colors for transitions resemble those from Figs. 4(a,b) in Ref. 7.
Blue and purple: pure transmon transitions. Yellow and orange: pure spin-flip
transitions. Dashed and dotted: mixed transitions involving both spin and
transmon degrees of freedom. Parameters: $U/\Delta=3$,
$\epsilon=-\frac{U}{2}$, $v/\Delta=0.5$, $v_{\uparrow\downarrow}/\Delta=0.2$,
$v_{\mathrm{sc}}/\Delta=0.2$, $t_{p}/\Delta=0.1$, $E_{c}/\Delta=0.02$.

If the QD is in the doublet state, the transmon excitations (labelled by an
integer $i=0,1,2,\ldots$) split into $|i,\uparrow\rangle$ and
$|i,\downarrow\rangle$ states due to the Zeeman splitting of the QD level in
magnetic field, or due to SOC-induced spin splitting for general values of
$\phi_{\mathrm{ext}}$. Figure 7 shows the low-energy doublet states
$|i,\sigma\rangle$ for the two lowest transmon levels $i=0,1$ in the presence
of the SOC at zero magnetic field, as well as all six transitions between
these four states. This figure can be compared with Fig. 4b in Ref. 7, showing
the two-tone spectroscopy of the joint two-qubit (transmon and ASQ) system. In
addition to the spin-conserving transmon transitions (blue and purple), there
are two “mixed” transition lines involving both spin and transmon degrees of
freedom in the presence of coherent coupling between them [7] (dashed and
dotted). The pure spin-flip transitions (orange and yellow) reflect the
different spin-orbit splitting, which is here found to be much stronger for
the $i=0$ transmon ground state (yellow) than the $i=1$ transmon excited state
(orange). The general parameter dependence of the amplitude of SOC splitting
in the different $i$ levels of the transmon ASQs is an interesting open
question for future work.

##  V Applications

Having established that the presented model reproduces all key features of the
QD physics, as well as phenomena arising from the interplay of the Josephson
effect and charging, we now show that the method is able to address more
challenging problems such as those of strong-coupling between the QD and
transmon degrees of freedom, time-dependent problems, as well as the
calculation of matrix elements that quantify the possibility of driving
various transitions.

###  V.1 Strong coupling of transmon and spin degrees of freedom

In the absence of SOC $S_{z}$ is a conserved quantum number, so with
increasing magnetic field the $|i=0,\uparrow\rangle$ and
$|i=1,\downarrow\rangle$ states cross at a certain field strength. In the
presence of the SOC, however, the states instead mix which leads to an avoided
crossing. This effect can be utilized to implement coupling between two
qubits, one encoded in the first two transmon states and the other in the spin
of the QD electron, see Ref. 7 for an experimental realization.

Here we present a minimal example of such strong mixing effects. Fig. 8(a)
shows the splitting of the doublet states with increasing Zeeman splitting on
the QD, $E_{Z}$. Dashed lines correspond to the case with no SOC, while solid
lines exhibit the avoided crossing for finite SOC. Our method enables direct
access to the wavefunction, which in general is not a product state of the two
subsystems, i.e., the two degrees of freedom become entangled. In Fig. 8(b) we
plot $|\alpha_{\phi}|^{2}$ for the two states for a range of Zeeman energies
$E_{Z}$ through the avoided crossing, gauging the mixing in the transmon
subspace. The $|0,\uparrow\rangle$ state evolves from a single-peak shape,
characteristic of the transmon ground state, to two peaks at the crossing
point $E_{Z}/\Delta=0.25$. The single peak is reinstated as the energy
difference between the states increases at large $E_{Z}$. We find that this
behavior is ubiquitous and largely independent of the values of the
parameters.

Such avoided crossing has been recently observed experimentally in an ASQ-
transmon device, see Fig. 4(c,d) in Ref. 7. We return to this phenomenon in
Sec. V.3.

Figure 8: SOC-induced level repulsion. (a) Spectra of the doublet manifold
with increasing Zeeman splitting of the QD level with (dashed) no spin-orbit
coupling and with (solid) $t_{sc}=0.2\Delta$,
$v_{\uparrow\downarrow}=0.2\Delta$. In both cases, the ground state at
$E_{Z}=0$ is set to zero energy. (b) Amplitudes in the $\phi$-space of the
$|0,\uparrow\rangle$ and $|1,\downarrow\rangle$ eigenstates through the
avoided crossing. Other parameters: $U/\Delta=3$, $\epsilon=-\frac{U}{2}$,
$v/\Delta=0.5$, $t_{p}/\Delta=0.1$, $\phi_{\mathrm{ext}}=\pi$, producing
$E_{J}^{\mathrm{eff}}\approx 0.22\Delta$. $E_{c}=0.02\Delta\approx
0.1E_{J}^{\mathrm{eff}}$. $N=601$ and $n_{0}=300$.

###  V.2 Exact time evolution

Because of the exponential reduction of the Hilbert space achieved with the
flat-band approximation, by obtaining the full eigendecomposition of the
Hamiltonian via exact diagonalisation it is straightforward to compute the
time evolution and calculate various non-equilibrium properties by numerically
integrating the time-dependent Schrödinger equation. Such calculations are out
of reach of other methods.

A minimal example of Rabi oscillations is shown in Fig. 9 [39, 40]. We turn on
a small Zeeman splitting $E_{Z}/\Delta=0.1$ to split the two doublet states,
and apply a Gaussian spin-flip pulse (for simplicity we assume that we are at
resonance, that we work in the interaction picture, and that we have dropped
the counter-rotating terms in the rotating-wave approximation [40]):

$H^{\prime}(t)=Ae^{-(t-t_{0})^{2}/2\sigma^{2}}\hat{S}_{x}.$ (20)

The time evolution of the expectation value of $S_{z}$ at the QD is shown in
Fig. 9(a). For the chosen values of $A$ and $\sigma$ the pulse causes a single
oscillation of the spin, which then settles at $S_{z}\sim-1/2$. Fig. 9(b)
shows the final value of $\langle S_{z}\rangle$ long after the pulse with
varying $A$ and $\sigma$. Typical Rabi fringes are observed.

An obvious utility of such calculations is in predicting optimal pulses for
controlling the qubits [69]. This could be done by interfacing our code with a
quantum optimal control library where the required input is the system
Hamiltonian, eg. Ref. 70. For example, one could predict and optimize two-
qubit gates in the ASQ-transmon devices [7].

Figure 9:  Coherent evolution in response to spin-rotating pulses. (a) Time
evolution of the expectation value of $S_{z}$ of the QD spin when applying a
Gaussian $\hat{S}_{x}$ pulse (shown on top). (b) Expectation value of $S_{z}$
at long times, $t\gg\sigma$, for a pulse of height $A$ and width $\sigma$.
Other parameters: $U=3\Delta$, $\epsilon=-\frac{U}{2}$, $v/\Delta=0.5$,
$E_{Z}/\Delta=0.1$, $t_{p}=0$, $E_{c}=0$, no SOC.

###  V.3 Transition matrix elements

Figure 10:  Matrix elements for transitions: overview of the
$\phi_{\mathrm{ext}}$ dependence. a) Charge operator, $|\langle
i|\hat{n}_{\mathrm{QD}}|j\rangle|^{2}$, b) current operator, $|\langle
i|J|j\rangle|^{2}$, with $J=\partial\hat{H}/\partial\phi_{\mathrm{ext}}$. c)
dipole operator, $|\langle
i|\hat{n}_{\mathrm{L}}-\hat{n}_{\mathrm{R}}|j\rangle|^{2}$. Vertical positions
indicate the transition frequency. Symbol areas are proportional to the
squared absolute values of the matrix elements; the number indicated in the
legend corresponds to the maximum value attained (the two charge operators are
dimensionless, while the “current” operator is given in the energy units of
$\Delta$). The numbers are not directly comparable, because the overall
transition probability also depends on different prefactors. Parameters:
$U/\Delta=3$, $\epsilon=-\frac{U}{2}$, $v/\Delta=0.5$, $t_{sc}=0.2\Delta$,
$v_{\uparrow\downarrow}=0.2\Delta$, $t_{p}/\Delta=0.1$,
$\phi_{\mathrm{ext}}=\pi$, producing $E_{J}^{\mathrm{eff}}\approx 0.22\Delta$.
$E_{c}=0.02\Delta\approx 0.1E_{J}^{\mathrm{eff}}$ (as in Fig. 8), for a
perpendicular magnetic field of $E_{Z}/\Delta=0.02$. $N=101$ and $n_{0}=50$.
Figure 11:  Magnetic field dependence: parallel field $E_{x}$. First column:
$\phi_{\mathrm{ext}}=0$, second column: $\phi_{\mathrm{ext}}=\pi/4$, third
column: $\phi_{\mathrm{ext}}=\pi/2$, last column: $\phi_{\mathrm{ext}}=\pi$.
a,e,i,m) Eigenenergies. b,f,j,n) Transition frequencies. c,g,k,o) Matrix
elements for the charge operator. d,h,l,p) Matrix elements for the current
operator. Color scheme for transitions as in Fig. 7. Parameters: $U/\Delta=3$,
$\epsilon=-\frac{U}{2}$, $v/\Delta=0.5$, $t_{sc}=0.2\Delta$,
$v_{\uparrow\downarrow}=0.2\Delta$, $t_{p}/\Delta=0.1$, $E_{c}=0.02\Delta$.
$N=101$ and $n_{0}=50$.  Figure 12:  Magnetic field dependence: perpendicular
field $E_{z}$ (indicated in units of $\Delta$). For easier comparison of
matrix elements for different types of transitions on the same vertical scale,
we multiply all matrix elements for pure spin-flip transition matrix elements
(yellow, orange) by 100, and all matrix elements for mixed transitions
(dashed, dotted) by 10; the matrix elements for pure transmon transitions
(blue, purple) are not scaled. At the upper range of magnetic field we start
to observe anomalies due to the anticrossing between the
$|1,\downarrow\rangle$ and $|2,\uparrow\rangle$ states that occurs at
$E_{z}/\Delta\approx 0.095$. Parameters as in Fig. 11.  Figure 13:  Asymptotic
low-field behavior of current matrix elements for spin-flip transitions at
$\phi_{\mathrm{ext}}\approx 0$. Parameters as in Fig. 12.

The Hamiltonian includes an effective low-energy description of the QD and its
coupling to the rest of the circuit, at the level of single-orbital Anderson
impurity model with extensions for the SOC. While approximate, this
description is nonetheless microscopically realistic, it contains all terms
which are believed to be relevant (in the renormalisation group sense, as well
as considering the symmetry constraints), and it is most likely adequate to
address many questions about the interplay between the different degrees of
freedom. In particular, one can investigate various coupling mechanisms that
permit the qubit control. In the experiments of Ref. [7], the transitions were
driven by applying microwave pulses to the QD gate electrode, which has a
dominant capacitive coupling to the QD orbital. All types of transitions shown
in Fig. 7 were observed, both pure transmon, pure spin-flip, as well as mixed
transitions. Due to the SOC, the electric field modulates not only the orbital
wavefunction in real space but indirectly also the spin degrees of freedom, a
phenomenon known as the electric dipole spin resonance (EDSR) [71, 72, 73, 74,
75, 7]. To shed light on the physical mechanisms of state transitions we now
investigate the off-diagonal matrix elements of some key operators, in
particular the charge operator, the electric dipole $\sum
q_{i}\mathbf{r}_{i}$, and the current operator. For simplicity, we use the
following operating definitions: the charge operator is taken to be the QD
occupancy operator $n_{\mathrm{QD}}$, the dipole is defined as $n_{L}-n_{R}$
(which does not involve the QD directly, but contains information about the
charge asymmetry in the device), and the current is defined as $J=\partial
H/\partial\phi_{\mathrm{ext}}$. The current matrix element is calculated using
the off-diagonal form of the Hellmann-Feynman theorem:

$J_{ij}=\langle i|\frac{\partial
H}{\partial\phi_{\mathrm{ext}}}|j\rangle=(E_{i}-E_{j})\left(\frac{\partial\langle
i|}{\partial\phi_{\mathrm{ext}}}\right)|j\rangle,$ (21)

where $E_{i}$, $E_{j}$ are the eigenenergies of the states $|i\rangle$ and
$|j\rangle$. To calculate the derivative we use finite differences, taking
$\Delta\phi_{\mathrm{ext}}=0.01\pi$.

We first consider the $\phi$-dependence of these operators for an
experimentally realistic parameter set, Fig. 10, at finite magnetic field
applied in the direction that is perpendicular to the SOC spin-polarization
direction (“perpendicular field”). For all three operators considered, we
generally observe that the matrix elements are the largest for the transmon
transitions, the weakest for the pure spin-flip transitions, and intermediate
for the mixed transitions. The operators considered couple to the spin degree
of freedom only indirectly, via the SOC, which naturally explains this
hierarchy. Importantly, we also observe different flux dependencies, depending
on the operator considered and the transition type. The matrix elements for
the charge operator, Fig. 10a), have a nearly sinusoidal dependence for the
pure transmon transition; they reach a maximum value at
$\phi_{\mathrm{ext}}\approx\pm\pi/2$, and they go through a zero at
$\phi_{\mathrm{ext}}=0$ and $\phi_{\mathrm{ext}}=\pm\pi$. The pure spin-flip
transitions are anharmonic, with a maximum close to
$\phi_{\mathrm{ext}}=\pm\pi/4$, and go through a zero at
$\phi_{\mathrm{ext}}=0$ and $\phi_{\mathrm{ext}}=\pm\pi$. The mixed
transitions have a maximum at $\phi_{\mathrm{ext}}\approx\pm\pi/3$ and remain
finite at both $\phi_{\mathrm{ext}}=0$ and $\phi_{\mathrm{ext}}=\pm\pi$. The
matrix elements for the current operator, Fig. 10b), have a weak
$\phi_{\mathrm{ext}}$ dependence for the transmon transitions. Spin-flip and
mixed transitions are strongly $\phi_{\mathrm{ext}}$-dependent and they are
out-of-phase (not easily seen in Fig. 10 due to the chosen scale): the pure
spin-flip matrix elements are maximal at $\phi_{\mathrm{ext}}=0$ and go
through zero for $\phi_{\mathrm{ext}}=\pi/2$; the mixed transition matrix
element are, conversely, large at $\phi_{\mathrm{ext}}=\pi/2$ and go through
zero at $\phi_{\mathrm{ext}}=0$. The matrix elements for the dipole moment,
Fig. 10c), have a similar $\phi$ dependence as those for the current, although
they are not exactly proportional. Because of this similarity, in the
following we only consider the charge and current operators, assuming that
these two operators represent the two important classes of the
$\phi_{\mathrm{ext}}$ dependence. (We have also considered the quadrupole
operator with an operational definition $n_{\mathrm{L}}+n_{\mathrm{R}}$; in
our model it brings no new information because it is trivially related to the
charge operator $n_{\mathrm{QD}}$.)

The experiments suggest that driving becomes easier with the application of
the magnetic field [7]. A magnetic field that is parallel to the SOC spin-
polarization direction (“parallel field”) has little effect, see Fig. 11. The
most striking feature is the fact that only transmon transitions show finite
matrix elements in this case, spin-flip and mixed transitions are altogether
absent. From this we immediately conclude that within our model a
perpendicular component of magnetic field is necessary to enable transitions
involving the spin degree of freedom if the coupling is only via charge or
current.

We next investigate the case of a perpendicular magnetic field, Fig. 12. The
pure transmon transition matrix elements (blue, violet) are not affected
strongly by the field, yet they do show some non-trivial field dependence. The
pure spin-flip and mixed transitions have even more complex field variation
which furthermore depends on the value of $\phi_{\mathrm{ext}}$. The pure
spin-flip matrix elements (orange, yellow) remain zero for the charge operator
at $\phi_{\mathrm{ext}}=0$, and for the current operator at
$\phi_{\mathrm{ext}}=\pi/2$; for general $\phi_{\mathrm{ext}}$, they increase
with the field, although there is a sizable dependence on the transmon level.
Mixed transitions (dashes, dotted) are found to have non-monotonic field
dependence in most cases, with a maximum value attained for some intermediate
field strength. We note that the cases of $\phi_{\mathrm{ext}}\approx 0$ and
$\phi_{\mathrm{ext}}\approx\pi$ appear to be anomalous, and the $E_{z}\to 0$
dependence of current matrix elements is misleading: at very low field values
the matrix elements for transitions involving spin-flips actually drop to zero
rather the saturate, see Fig. 13. Nevertheless, even a small perpendicular
magnetic field is apparently sufficient to enable the spin-flip transitions at
these $\phi_{\mathrm{ext}}$ settings.

It turns out that the role of perpendicular field can be explained via a
simplified zero-bandwidth-approximation (ZBA) calculation, see App. B: the
low-field behavior is indeed found to be linear in $E_{z}$. Such calculations
also reproduce some trends observed in the $\phi$-dependence (see also App.
C).

For completeness, we now study the variation of eigenenergies and matrix
elements upon variation of other key model parameters. We consider the case of
small perpendicular field $E_{z}=0.02\Delta$ at $\phi_{\mathrm{ext}}=\pi/2$.
The variation of the QD level $\epsilon$, see Fig. 14, has a weak effect on
the energies, but we observe a large effect on the charge matrix elements,
which is monotonic for transmon and mixed transitions, while for pure spin-
flip transitions it changes sign not far from the particle-hole symmetric
point at $\epsilon=-U/2=-1.5\Delta$. This behavior is also predicted by the
ZBA calculation in App. B. Furthermore, we note that the relative magnitude of
pure spin-flip and mixed transition matrix elements is a strong function of
$\epsilon$. For the current operators, the dependence on $\epsilon$ is very
weak for the transmon and mixed transitions; for the pure spin-flip transition
there is again a sign change, but these matrix elements are anyhow quite low.
Rather similar results are obtained upon variation of the electron-electron
repulsion $U$, see Fig. 15.

Figure 14:  Quantum dot level $\epsilon$ dependence. Parameters: $U/\Delta=3$,
$v/\Delta=0.5$, $t_{sc}=0.2\Delta$, $v_{\uparrow\downarrow}=0.2\Delta$,
$t_{p}/\Delta=0.1$, $E_{c}=0.02\Delta$, $\phi_{\mathrm{ext}}=\pi/2$,
$E_{z}/\Delta=0.02$. $N=101$ and $n_{0}=50$.  Figure 15:  Quantum dot
electron-electron repulsion $U$ dependence. Parameters:
$\epsilon/\Delta=-1.5$, $v/\Delta=0.5$, $t_{sc}=0.2\Delta$,
$v_{\uparrow\downarrow}=0.2\Delta$, $t_{p}/\Delta=0.1$, $E_{c}/\Delta=0.02$,
$\phi_{\mathrm{ext}}=\pi/2$, $E_{z}/\Delta=0.02$.

The charging energy $E_{c}$, see Fig. 16, strongly affects the transmon
transition frequency, as expected. More interestingly, increasing $E_{c}$
suppresses the pure spin-flip frequency, i.e., it renormalizes the SOC
splitting, the effect being particularly strong for the excited $i=1$ level.
For both charge and current matrix elements, the strongest $E_{c}$ dependence
is observed in the mixed transitions, while pure transmon and pure spin-flip
transitions (esp. for $i=0$) appear to be less sensitive to the value of
$E_{c}$. We note that the proposed model is the minimal Hamiltonian that
allows to study how the QD properties are renormalized by the phase
fluctuations in the superconducting circuit. One could reexamine, for example,
the impurity Knight shift explored in Ref. 44 and determine how a finite
$E_{c}$ affects the renormalization of the Zeeman splitting by the coupling to
the superconducting leads.

Figure 16:  Charging energy $E_{c}$ dependence. Parameters: $U/\Delta=3$,
$\epsilon=-U/2$, $v/\Delta=0.5$, $t_{sc}=0.2\Delta$,
$v_{\uparrow\downarrow}=0.2\Delta$, $t_{p}/\Delta=0.1$,
$\phi_{\mathrm{ext}}=\pi/2$, $E_{z}/\Delta=0.02$. $N=101$ and $n_{0}=50$.

The strength of SOC is proportional to the spin-flip hybridisation
$v_{\uparrow\downarrow}$ (in fact, it is proportional to the product of
$vv_{\uparrow\downarrow}t_{\mathrm{sc}}$, see Ref. [43]). Indeed, the SOC
level splittings are directly reflected in the approximately linear trends of
transition frequencies, see Fig. 17. As a general trend (with some exceptions
for weak SOC), the pure spin-flip matrix elements tend to increase with
$v_{\uparrow\downarrow}$, while the mixed transition tend to decrease; this is
the case for both charge and current operators. The transitions involving the
transmon degree of freedom show a sign change of matrix elements for the
charge operator at some large value of SOC (this is also reflected in mixed
transitions). As for the spin-flip transitions, we observe a sign change for
the current operator. This is, in fact, predicted by the ZBA calculation, see
App. B, which predict a proportionality of transition matrix elements to
$v^{2}-v_{\uparrow\downarrow}^{2}$. Clearly, this is renormalized by the
coupling to transmon degrees of freedom, so that the cancellation point
becomes transmon level dependent.

Figure 17:  SOC (spin-flip hybridization $v_{\uparrow\downarrow}$) dependence.
Parameters: $U/\Delta=3$, $\epsilon=-U/2$, $v/\Delta=0.5$, $t_{sc}=0.2\Delta$,
$t_{p}/\Delta=0.1$, $E_{c}/\Delta=0.02$, $\phi_{\mathrm{ext}}=\pi/2$,
$E_{z}/\Delta=0.02$. $N=101$ and $n_{0}=50$.

We now return to the case of strong ASQ-transmon coupling, already introduced
in Sec. V.1. We thus extend the plots similar to that in Fig. 12 to a wider
range of magnetic fields to fully reveal the variation of the transition
matrix elements in the vicinity of anticrossings, of which there are three in
the range considered: between $|1,\downarrow\rangle$ and $|2,\uparrow\rangle$
at $E_{z}/\Delta\approx 0.095$, between $|0,\downarrow\rangle$ and
$|1,\uparrow\rangle$ at $E_{z}/\Delta\approx 0.25$ (the case discussed in
relation to Fig. 18 and the one experimentally explored in Fig. 3c,d) in Ref.
7), and between $|1,\downarrow\rangle$ and $|1,\uparrow\rangle$ at
$E_{z}/\Delta\approx 0.34$. It should be noted that we here use labelling that
adiabatically continues the zero-field labels, although it is clear, for
example, that the third anticrossing occurs between the states that actually
correspond to the $i=2$ and $i=0$ transmon levels. Let us focus on the
anticrossing at $E_{z}/\Delta\approx 0.25$. The anticrossing between
$|0,\downarrow\rangle$ and $|1,\uparrow\rangle$ (panel a) is indirectly
revealed through the behavior of transition frequencies as an “anticrossing”
between a pure transmon transition and a mixed (spin-flipping transmon)
transition (yellow and purple lines involving the $|0,\uparrow\rangle$ ground
state in panel b, but also orange and blue lines involving the
$|1,\downarrow\rangle$ excited state). Across the anticrossing range the
matrix elements are changing between two asymptotic value, which is a
signature of the changing nature of the states. Importantly, at the point of
closest approach all matrix elements are sizable, even if the asymptotic value
on one side of the transition is very small.

Figure 18:  Perpendicular magnetic field dependence. Parameters as in Fig. 12,
except $\phi_{\mathrm{ext}}=(3/4)\pi$. The values of matrix elements are not
rescaled in any way.

The examples in this subsection clearly illustrate that the proposed method
can provide guidance in choosing system parameters that lead to the optimal
operating point, both in terms of level energies as well as for the transition
matrix elements. We have observed that the variation of the matrix elements is
in general not trivial, it can easily span multiple decades in amplitude, it
is often non-monotonic, and there are sign changes. While some general trends
for pure spin-flip and pure transmon transition can be extracted for
simplified ZBA calculation that do not involve the dynamical superconductor
(transmon) degrees of freedom (App. B and C), the full model is required to
capture all the details and to account for the mixed transitions.

##  VI Conclusion

We have investigated a system of a QD embedded in a junction between two
superconductors from the viewpoint of the charge-conserving Richardson model
of superconductivity. We argue that this approach enables clearer
understanding of various properties of the system, especially those related to
the emergence of the superconducting phase difference and its dynamics. Most
importantly, it is able to access parameters regimes that BCS-type
approximations cannot, most notably the case of finite charging energy in the
superconducting islands. The model treats on the same level all important
physical processes—the impurity physics of the interacting QD, the charging
energy and pairing terms in the superconductors, and the Josephson effect.
Furthermore, by including spin-orbit coupling this becomes a model for Andreev
spin qubits in transmon devices.

The central step of our solution is the flat-band approximation, which
consists of neglecting the kinetic energy in the superconductors. This leads
to an exponential reduction of the Hilbert space while retaining all states
through which the QD is coupled to the SCs and are important for the low-
energy physics of the system.

The obtained Hamiltonian is similar, in spirit, to the zero-bandwidth BCS
approximation for the superconducting Anderson impurity model, however it
retains information about the distribution of the Cooper pairs between the two
superconductors, as in the conventional transmon Hamiltonian, thereby
combining aspects of both approaches. It is formulated in the conjugate space
of well-defined charge, which enables natural implementation of the charging
energy terms in the SCs and allows us to capture the quantum fluctuations of
the phase difference.

We have demonstrated that the model reproduces all features expected of a
quantum dot Josephson junction in the limit of zero charging energy, and we
have singled out the parameter regimes where standard approximations fail,
pointing out the importance of quasiparticles in the superconductor induced by
the pair-breaking effects of the “magnetic” quantum dot. Next, we investigated
the effect of the charging energy and its interplay with the Josephson effect.
We observed and quantified the evolution between the limits of well-defined
phase difference and well-defined charge difference on the level of eigenstate
properties and quantum fluctuations of phase. Finally, we showcased the
applications of the method with three experimentally relevant problems: 1) the
case where the transmon and the QD degrees of freedom mix, 2) the calculation
of QD spin flipping induced by a pulse, 3) the investigation of transition
matrix elements.

The model can be applied to other states as well. While in this work we have
focused on the doublet subspace, as appropriate for Andreev spin qubits, it is
equally possible to study the singlet subspace which is relevant for Andreev
level qubits based on the superpositions of empty and doubly occupied quantum
dot orbital.

Finally, noting that the Hilbert space is only moderately large, it seems
quite feasible to include relaxation and dephasing effects in the model
through appropriate Lindblad operators and to study decoherence in the
associated open quantum system.

###### Acknowledgements.

We thank Marta Pita Vidal and Arno Bargerbos for suggestions concerning the
transition matrix elements and other applications of the presented method. We
acknowledge the support of the Slovenian Research and Innovation Agency (ARIS)
under P1-0416 and J1-3008.

##  Appendix A Basis states and Hamiltonian generation

Here we list the basis states that span the singlet ($S=0$, $S_{z}=0$) and the
doublet ($S=1/2$, $S_{z}=+1/2$) subspaces, and describe how the Hamiltonian
matrix is generated.

The code, including the matrices, is available on GitHub [59].

###  A.1 Singlet states

The singlet manifold has dimension 14, consisting of one state with no
unpaired particles:

  * •

$|m_{L},0,m_{R},0\rangle$,

six states with two unpaired particles:

  * •

$|m_{L},0,m_{R},2\rangle$,

  * •

$|m_{L},2,m_{R},0\rangle$,

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},0,m_{R},0\rangle$,

  * •

$\left(d^{\dagger}_{\uparrow}|m_{L},\downarrow,m_{R},0\rangle-d^{\dagger}_{\downarrow}|m_{L},\uparrow,m_{R},0\rangle\right)/\sqrt{2}$,

  * •

$\left(d^{\dagger}_{\uparrow}|m_{L},0,m_{R},\downarrow\rangle-d^{\dagger}_{\downarrow}|m_{L},0,m_{R},\uparrow\rangle\right)/\sqrt{2}$,

  * •

$\left(|m_{L},\uparrow,m_{R},\downarrow\rangle-|m_{L},\downarrow,m_{R},\uparrow\rangle\right)/\sqrt{2}$,

six states with four unpaired particles:

  * •

$|m_{L},2,m_{R},2\rangle$,

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},2,m_{R},0\rangle$,

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},0,m_{R},2\rangle$,

  * •

$\left(d^{\dagger}_{\uparrow}|m_{L},\downarrow,m_{R},2\rangle-d^{\dagger}_{\downarrow}|m_{L},\uparrow,m_{R},2\rangle\right)/\sqrt{2}$,

  * •

$\left(d^{\dagger}_{\uparrow}|m_{L},2,m_{R},\downarrow\rangle-d^{\dagger}_{\downarrow}|m_{L},2,m_{R},\uparrow\rangle\right)/\sqrt{2}$,

  * •

$\left(d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},\uparrow,m_{R},\downarrow\rangle-d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},\downarrow,m_{R},\uparrow\rangle\right)/\sqrt{2}$,

and one state with six unpaired particles:

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},2,m_{R},2\rangle$.

###  A.2 Doublet states

The doublet manifold also has dimension 14. It contains three states with one
unpaired particle:

  * •

$d^{\dagger}_{\uparrow}|m_{L},0,m_{R},0\rangle$,

  * •

$|m_{L},\uparrow,m_{R},0\rangle$,

  * •

$|m_{L},0,m_{R},\uparrow\rangle$,

eight states with three unpaired particles:

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},\uparrow,m_{R},0\rangle$,

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},0,m_{R},\uparrow\rangle$,

  * •

$d^{\dagger}_{\uparrow}|m_{L},2,m_{R},0\rangle$,

  * •

$d^{\dagger}_{\uparrow}|m_{L},0,m_{R},2\rangle$,

  * •

$|m_{L},\uparrow,m_{R},2\rangle$,

  * •

$|m_{L},2,m_{R},\uparrow\rangle$,

  * •

$\left(d^{\dagger}_{\downarrow}|m_{L},\uparrow,m_{R},\uparrow\rangle-d^{\dagger}_{\uparrow}|m_{L},\downarrow,m_{R},\uparrow\rangle\right)/\sqrt{2}$,

  * •

$\Bigl{(}d^{\dagger}_{\downarrow}|m_{L},\uparrow,m_{R},\uparrow\rangle+d^{\dagger}_{\uparrow}|m_{L},\downarrow,m_{R},\uparrow\rangle\\\
\hskip
56.9055pt-2d^{\dagger}_{\uparrow}|m_{L},\uparrow,m_{R},\downarrow\rangle\Bigr{)}/\sqrt{6}$,

and three states with five unpaired particles:

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},2,m_{R},\uparrow\rangle$,

  * •

$d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}|m_{L},\uparrow,m_{R},2\rangle$,

  * •

$d^{\dagger}_{\uparrow}|m_{L},2,m_{R},2\rangle$.

These are the states with $S_{z}=+1/2$. In cases of broken time-reversal
symmetry the $S_{z}=\pm 1/2$ degeneracy is broken and the basis should be
extended with the corresponding $S_{z}=-1/2$ states.

###  A.3 Hamiltonian generation

The diagonal part of the Hamiltonian is

$\begin{split}H_{\mathrm{diag}}=&\epsilon
n_{\mathrm{QD}}+Un_{\mathrm{QD}\uparrow}n_{\mathrm{QD}\downarrow}+\frac{g}{2}\sum_{\beta\sigma}f^{\dagger}_{\beta\sigma}f_{\beta\sigma}\\\
&+\sum_{\beta}E_{c}^{\beta}\left(\sum_{\sigma}f^{\dagger}_{\beta\sigma}f_{\beta\sigma}+2m_{\beta}-n_{0}\right)^{2}\\\
\end{split}$ (22)

and can be read out directly for each basis state.

The off-diagonal part is generated by first calculating analytical expressions
for all matrix elements

$\langle
d^{\prime},m^{\prime}_{L},f^{\prime}_{L},m^{\prime}_{R},f^{\prime}_{R}|H_{\mathrm{hyb}}|d,m_{L},f_{L},m_{R},f_{R}\rangle.$
(23)

Only states where a single electron is transferred from the QD to a SC are
coupled, so the matrix is sparse, with elements of the type

$\begin{split}&\langle
m^{\prime}_{L},\uparrow,m^{\prime}_{R},0|H_{\mathrm{hyb}}d^{\dagger}_{\uparrow}|m_{L},0,m_{R},0\rangle\\\
&\hskip
28.45274pt=v_{L}\delta_{m_{L},m^{\prime}_{L}+1}\delta_{m_{R},m^{\prime}_{R}}.\end{split}$
(24)

The matrices of $H_{\mathrm{hyb}}$ in the singlet (Eq. (25)) and doublet (Eq.
(26)) are given below. To shorten the notation, we denote
$\delta^{\beta}_{x}=\delta_{m_{\beta},m^{\prime}_{\beta+x}}$.

[0 -vLδ0Lδ0R2-vRδ0Lδ0R2vLδ1Lδ0R2vRδ0Lδ1R20 0 0 0 0 0 0 0 0 -vLδ0Lδ0R20 0 0 0
vLδ1Lδ0R20 0 0 \- vRδ0Lδ1R2\- 3vRδ0Lδ1R20 0 0 -vRδ0Lδ0R20 0 0 0 0 0 0
vRδ0Lδ1R2vLδ1Lδ0R0 0 0 0 vLδ-1Lδ0R20 0 0 0 vLδ0Lδ0R20 0 0 \- vRδ0Lδ0R2\-
3vRδ0Lδ0R20 0 0 vRδ0Lδ-1R20 0 0 0 0 0 0 vRδ0Lδ0R2vLδ0Lδ0R0 0 0 0 0 vLδ-1Lδ0R20
vLδ0Lδ0R20 0 -vRδ0Lδ0R20 0 0 0 vRδ0Lδ1R20 0 0 0 0 0 0 -vRδ0Lδ0R20 0 0
vLδ0Lδ0R0 0 0 vRδ0Lδ1R20 0 0 0 0 0 0 0 -vLδ0Lδ0R2\- vRδ0Lδ0R2\- 3vRδ0Lδ0R20 0
vLδ1Lδ0R20 0 vRδ0Lδ-1R20 vRδ0Lδ0R20 0 -vLδ0Lδ0R20 0 0 0 vLδ1Lδ0R20 \-
vRδ0Lδ-1R2vLδ-1Lδ0R\- vRδ0Lδ0R2vLδ0Lδ0R0 vLδ0Lδ0R\- vRδ0Lδ0R20 0 0
-vLδ1Lδ0RvRδ0Lδ1R20 0 \- 3vRδ0Lδ-1R20 \- 3vRδ0Lδ0R20 0 0 \- 3vRδ0Lδ0R20 0 0 0
3vRδ0Lδ1R20 0 0 0 0 vRδ0Lδ-1R20 0 0 -vLδ-1Lδ0R0 0 0 vRδ0Lδ0R20 0 0 0 0 0 0 0
vLδ-1Lδ0R212vRδ0Lδ-1R123vRδ0Lδ-1R0 0 vLδ0Lδ0R20 0 0 0 0 0
vRδ0Lδ-1R2vLδ-1Lδ0R20 0 0 vRδ0Lδ0R2vLδ0Lδ0R20 ]

(25)

[0 0 0 0 -vLδ1Lδ0R-vRδ0Lδ1R0 0 0 0 0 0 0 0 0 0 0 0 0 -vRδ0Lδ0R0 0 0 0
vLδ1Lδ0R0 0 0 0 0 0 0 -vLδ0Lδ0R0 0 0 0 0 0 vRδ0Lδ1R0 0 0 0 0 0
-vLδ0Lδ0R-vRδ0Lδ0R0 0 0 0 0 0 0 0 -vLδ-1Lδ0R0 -vLδ0Lδ0R-vLδ0Lδ0R0 0 vRδ0Lδ0R20
0 vLδ1Lδ0R0 0 vRδ0Lδ1R20 -vRδ0Lδ-1R-vRδ0Lδ0R0 -vRδ0Lδ0R0 0 vLδ0Lδ0R20
vRδ0Lδ1R0 0 0 vLδ1Lδ0R20 0 0 0 0 vRδ0Lδ0R2vLδ0Lδ0R20 0 0 0 vRδ0Lδ1R2vLδ1Lδ0R20
0 0 0 0 0 0 0 0 0 0 0 vLδ0Lδ0RvRδ0Lδ0R0 0 0 0 0 0 0 vRδ0Lδ-1R0 0 0 0 vLδ0Lδ0R0
0 0 0 0 0 0 vLδ-1Lδ0R0 0 0 0 0 0 vRδ0Lδ0R0 0 0 vLδ-1Lδ0R0 0 0 0
vRδ0Lδ-1R2vLδ0Lδ0RvLδ0Lδ0R0 0 0 -vRδ0Lδ0R2-vLδ1Lδ0R0 0 vRδ0Lδ-1R0 0 0
vLδ-1Lδ0R2vRδ0Lδ0R0 vRδ0Lδ0R0 0 -vLδ0Lδ0R2-vRδ0Lδ1R0 0 0 0
vRδ0Lδ-1R2vLδ-1Lδ0R20 0 0 0 -vRδ0Lδ0R2-vLδ0Lδ0R20 0 0 0 0 0 0 0 0 0 0 0
-vLδ-1Lδ0R-vRδ0Lδ-1R0 0 ]

(26)

##  Appendix B Zero-bandwidth-approximation calculation of spin-flip
transition matrix elements

The transition matrix elements for pure spin-flip transitions can be obtained
using a simplified zero-bandwidth approximation (ZBA) for the superconductors,
neglecting all dynamical effects (order-parameter fluctuations). In this
appendix we present the results of such a calculation and comment on the
differences that originate from coupling to the transmon degrees of freedom in
the full problem.

The Hamiltonian is
$H=H_{\mathrm{SC}}^{(L)}+H_{\mathrm{SC}}^{(R)}+H_{\mathrm{QD}}+H_{\mathrm{hop}}+H_{\mathrm{LR}}$,
with

$\begin{split}H_{\mathrm{SC}}^{(\beta)}&=\sum_{k\sigma}\epsilon_{k}c^{\dagger}_{\beta
k\sigma}c_{\beta k\sigma}+\Delta\sum_{k}e^{i\phi_{\beta}}c^{\dagger}_{\beta
k\downarrow}c^{\dagger}_{\beta k\uparrow}+\text{H.c.},\\\
H_{\mathrm{QD}}&=\epsilon
n_{\mathrm{QD}}+Un_{\mathrm{QD}\uparrow}n_{\mathrm{QD}\downarrow}+E_{z}/2(n_{\mathrm{QD}\uparrow}-n_{\mathrm{QD}\downarrow}),\\\
H_{\mathrm{hop}}&=\frac{v}{\sqrt{N}}\sum_{\beta=L,R}\sum_{k\sigma}d^{\dagger}_{\sigma}c_{\beta
k\sigma}+\mathrm{H.c.}\\\
&+\frac{iv_{\uparrow\downarrow}}{\sqrt{N}}\sum_{k\sigma}\left(d^{\dagger}_{\sigma}c_{L,k\bar{\sigma}}+c^{\dagger}_{R,k\bar{\sigma}}d_{\sigma}+\mathrm{H.c.}\right),\\\
H_{\mathrm{LR}}&=\frac{t_{\mathrm{sc}}}{N}\sum_{k,k^{\prime},\sigma}c^{\dagger}_{L,k\sigma}c_{R,k^{\prime}\sigma}+\text{H.c.}\end{split}$

We set $\phi_{L}=-\phi/2$ and $\phi_{R}=\phi/2$. We apply the gauge
transformation $c_{\beta k\sigma}\to e^{i\phi_{\beta}/2}c_{\beta k\sigma}$
[33, 32] to remove the phase from the pairing terms and transfer it to the
hybridisation part. We introduce the combinations of states
$f_{\beta\sigma}=(1/\sqrt{N})\sum_{k}c_{\beta k\sigma}$, and drop all other
modes from consideration. This leads to

$\begin{split}H_{\mathrm{SC}}^{(\beta)}&=\Delta
f^{\dagger}_{\beta\downarrow}f^{\dagger}_{\beta\uparrow}+\text{H.c.},\\\
H_{\mathrm{hop}}&=\sum_{\sigma}\Bigl{[}\sum_{\beta=L,R}ve^{i\phi_{\beta}}d^{\dagger}_{\sigma}f_{\beta\sigma}+\text{H.c.}\\\
&+iv_{\uparrow\downarrow}\left(e^{i\phi_{L}}d^{\dagger}_{\sigma}f_{L\bar{\sigma}}+e^{-i\phi_{R}}f^{\dagger}_{R\sigma}d_{\bar{\sigma}}+\text{H.c.}\right)\Bigr{]},\\\
H_{\mathrm{LR}}&=t_{\mathrm{sc}}\sum_{\sigma}e^{-i\phi_{L}+i\phi_{R}}f^{\dagger}_{L\sigma}f_{R\sigma}+\text{H.c.}\end{split}$

Finally, we perform the Bogoliubov transformation on both superconductors:

$\begin{split}f_{\beta,\uparrow}&=(b_{\beta,\uparrow}+b^{\dagger}_{\beta,\downarrow})/\sqrt{2},\\\
f^{\dagger}_{\beta,\downarrow}&=(b^{\dagger}_{\beta,\downarrow}-b_{\beta,\uparrow})/\sqrt{2},\end{split}$
(27)

and

$\begin{split}f_{\beta,\downarrow}&=(b_{\beta,\downarrow}-b^{\dagger}_{\beta,\uparrow})/\sqrt{2},\\\
f^{\dagger}_{\beta,\uparrow}&=(b^{\dagger}_{\beta,\uparrow}+b_{\beta,\downarrow})/\sqrt{2}.\end{split}$
(28)

In this basis the Hamiltonian is diagonal dominant: the only out-of-diagonal
matrix elements are couplings ($v$, $v_{\uparrow\downarrow}$,
$t_{\mathrm{sc}}$).

For a decoupled system ($v=v_{\uparrow\downarrow}=t_{\mathrm{sc}}=0$), the
low-energy subsector is spanned by the states
$|\sigma\rangle=d^{\dagger}_{\sigma}|\mathrm{BCS}\rangle$, where
$|\mathrm{BCS}\rangle$ is the BCS ground states of both superconductors.
Excited states are generated by creating quasiparticles in superconductors or
adding/removing an electron on the QD site. There are in total 30 states in
the $S_{z}=\pm 1/2$ sector of the Hilbert space, two in the low-energy
subsector and 28 in the high-energy subsector. To study the effect of admixing
high-energy states into the ground state doublet, we make use of the recursive
Schrieffer-Wolff transformation (RSWT) [76]. The results in the lowest order
are sufficiently compact to reproduce them here.

The transition matrix element for the charge operator is

$\langle\sigma|n|\bar{\sigma}\rangle=8E_{z}(v^{2}+v_{\mathrm{\uparrow\downarrow}}^{2})\left(\frac{1}{(U+2\delta+2\Delta)^{3}}-\frac{1}{(U-2\delta+2\Delta)^{3}}\right)\approx-\frac{96E_{z}\delta(v^{2}+v_{\mathrm{\uparrow\downarrow}}^{2})}{(U+2\Delta)^{4}}.$
(29)

The transition matrix element for the current operator
$J=\mathrm{d}H/\mathrm{d}\phi$ is

$\langle\sigma|J|\bar{\sigma}\rangle=2E_{z}\delta
t_{\mathrm{sc}}(v^{2}-v_{\uparrow\downarrow}^{2})\frac{U^{3}-4U\delta^{2}+18U^{2}\Delta+8\delta^{2}\Delta+60U\Delta^{2}+56\Delta^{3}}{\Delta(U-2\delta+2\Delta)^{3}(U+2\delta+2\Delta)^{3}}\sin\phi\approx\frac{2E_{z}\delta
t_{\mathrm{sc}}(v^{2}-v_{\uparrow\downarrow}^{2})(U+14\Delta)}{\Delta(U+2\Delta)^{4}}\sin\phi.$
(30)

(We note the difference in definition of current operator: here it is defined
as the derivative with respect to the phase difference, for the full model it
is defined as the derivative with respect to the external magnetic flux. At
qualitative level these definitions should be equivalent, but there could be
subtle differences in detailed behavior.) As in the full model, a
perpendicular field is needed to enable the spin-flip transitions. Both matrix
elements go through zero at the particle-hole symmetric point $\delta=0$; in
the full model, the zeros are shifted to some transmon-level-dependent non-
zero value of $\delta$, see Figs. 14 and 15. The current matrix element also
goes to zero for a given ratio of spin-flip and non-spin-flip tunneling rates,
in the ZBA limit for $v_{\uparrow\downarrow}/v=1$, and in the full model for
some renormalized values, see Fig. 17. Such cancellation does not occur for
the charge matrix element, due to a different relative sign for the $v^{2}$
and $v_{\uparrow\downarrow}^{2}$ contributions. At this order of the
approximation, the charge matrix elements do not depend on $\phi$, while the
$\phi$-dependence of current matrix elements is sinusoidal. Neither of these
results agrees with the behavior in the full model, which shows anharmonic
$\phi$-dependence for $n$ and a phase shift of $\pi/2$ ($\cos\phi$ term) for
$J$. This appears to be a limitation of the lowest order RSWT, not of the
simplified ZBA description of the problem. At higher orders (or by performing
an exact diagonalisation of the ZBA Hamiltonian) we find that both matrix
elements are $\phi$-dependent (sinusoidal) with the same phase offsets as in
the full problem (charge matrix element maximal at $\pi/2,3\pi/2$ and passing
through zero at $0$ and $\pi$, current matrix elements maximal at $0$ and
$\pi$ and passing through zero at $\pi/2,3\pi/2$). Furthermore, we find that
zeros are shifted slightly away from the particle-hole symmetric point
$\delta=0$, which is expected given that the ZBA Hamiltonian in the presence
of SOC has no particular symmetry at $\delta=0$.

##  Appendix C Transmon transition matrix elements

The model proposed in the previous section can also provide some clues about
the $\phi$-dependence of the matrix elements for pure transmon transitions.
For example, it is clear that the QD occupancy (charge) affects the value of
the Josephson energy $E_{J}$, hence a modulation of dot filling can produce
transitions between the different transmon levels. The coupling between the
charge and SC phase degrees of freedom is quantified by the $\phi$-dependence
of the diagonal matrix elements (expectation values) of the charge operator
for the QD eigenstates. Observing that the $\phi$-dependence obtained in low-
order RSWT is not reliable, we investigated this problem numerically using
exact diagonalisation. We observe that the QD charge has a contribution to the
expectation value that is proportional to $\cos\phi$, while the Josephson
current has the expected contribution proportional to $\sin\phi$, see Fig. 19.
Based on this, one would expect that the transmon transitions would be
strongest for $\phi$ of strongest variation: close to $\pi/2$ for the charge
operator and close to $0,\pi$ for the current operator. The first prediction
is confirmed in the full model that has a sinusoidal $\phi_{\mathrm{ext}}$
dependence, see Fig. 10a), but not the second: Fig. 10b) shows that the matrix
element for $J$ actually remains fairly constant for all
$\phi_{\mathrm{ext}}$.

Figure 19:  Ground state expectation value of charge and Josephson current
operators in the zero-bandwidth-approximation model of the JJ. Parameters are:
$U/\Delta=3$, $\epsilon/\Delta=-1$, $v/\Delta=0.2$,
$v_{\uparrow\downarrow}/\Delta=0.1$, $t_{\mathrm{sc}}/\Delta=0.1$,
$E_{z}/\Delta=0.1$.

## References

  * 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).
  * Devoret  _et al._ [2004] M. H. Devoret, A. Wallraff, and J. M. Martinis, Superconducting qubits: A short review (2004).
  * Devoret [1997] M. Devoret,  _Quantum Fluctuations: Volume 63_ , edited by S. Reynaud, E. Giacobino, and F. David, Les Houches (North-Holland, Oxford, England, 1997).
  * 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, 10.1103/physreva.76.042319 (2007).
  * Wendin [2017] G. Wendin, Reports on Progress in Physics 80, 106001 (2017).
  * Hays  _et al._ [2021] M. Hays, V. Fatemi, D. Bouman, J. Cerrillo, S. Diamond, K. Serniak, T. Connolly, P. Krogstrup, J. Nygård, A. L. Yeyati, A. Geresdi, and M. H. Devoret, Science 373, 430 (2021).
  * Pita-Vidal  _et al._ [2023a] M. Pita-Vidal, A. Bargerbos, R. Žitko, L. J. Splitthoff, L. Grünhaupt, J. J. Wesdorp, Y. Liu, L. P. Kouwenhoven, R. Aguado, B. van Heck, A. Kou, and C. K. Andersen, Nature Physics 10.1038/s41567-023-02071-x (2023a).
  * Pita-Vidal  _et al._ [2023b] M. Pita-Vidal, J. J. Wesdorp, L. J. Splitthoff, A. Bargerbos, Y. Liu, L. P. Kouwenhoven, and C. K. Andersen, Strong tunable coupling between two distant superconducting spin qubits (2023b), arXiv:2307.15654 [quant-ph] .
  * Chtchelkatchev and Nazarov [2003] N. M. Chtchelkatchev and Y. V. Nazarov, Phys. Rev. Lett. 90, 226806 (2003).
  * Zazunov  _et al._ [2003] A. Zazunov, V. S. Shumeiko, E. N. Bratus’, J. Lantz, and G. Wendin, Physical Review Letters 90, 10.1103/physrevlett.90.087003 (2003).
  * Wilson [1975] K. G. Wilson, Reviews of Modern Physics 47, 773 (1975).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Reviews of Modern Physics 80, 395 (2008).
  * Tinkham [2004] M. Tinkham,  _Introduction to superconductivity_ , 2nd ed. (Dover Publications, 2004).
  * Gobert  _et al._ [2004] D. Gobert, U. Schollwöck, and J. von Delft, The European Physical Journal B - Condensed Matter and Complex Systems 38, 501 (2004).
  * Glazman and Catelani [2021] L. Glazman and G. Catelani, SciPost Physics Lecture Notes 10.21468/scipostphyslectnotes.31 (2021).
  * Tuominen  _et al._ [1992] M. T. Tuominen, J. M. Hergenrother, T. S. Tighe, and M. Tinkham, Physical Review Letters 69, 1997 (1992).
  * Lafarge  _et al._ [1993] P. Lafarge, P. Joyez, D. Esteve, C. Urbina, and M. H. Devoret, Physical Review Letters 70, 994 (1993).
  * Bouchiat  _et al._ [1998] V. Bouchiat, D. Vion, P. Joyez, D. Esteve, and M. H. Devoret, Physica Scripta T76, 165 (1998).
  * Nakamura  _et al._ [1999] Y. Nakamura, Y. A. Pashkin, and J. S. Tsai, Nature 398, 786 (1999).
  * Cottet [2002] A. Cottet,  _Implementation of a quantum bit in a superconducting circuit_ , Ph.D. thesis, Ecole normale superieure de Paris (2002).
  * Larsen  _et al._ [2015] T. W. Larsen, K. D. Petersson, F. Kuemmeth, T. S. Jespersen, P. Krogstrup, J. Nygård, and C. M. Marcus, Physical Review Letters 115 (2015).
  * de Lange  _et al._ [2015] G. de Lange, B. van Heck, A. Bruno, D. J. van Woerkom, A. Geresdi, S. R. Plissard, E. P. A. M. Bakkers, A. R. Akhmerov, and L. DiCarlo, Physical Review Letters 115 (2015).
  * Zheng  _et al._ [2023] H. Zheng, L. Y. Cheung, N. Sangwan, A. Kononov, R. Haller, J. Ridderbos, C. Ciaccia, J. H. Ungerer, A. Li, E. P. A. M. Bakkers, A. Baumgartner, and C. Schönenberger, Coherent control of a few-channel hole type gatemon qubit (2023), arXiv:2312.06411 [cond-mat.mes-hall] .
  * Clarke  _et al._ [1988] J. Clarke, A. N. Cleland, M. H. Devoret, D. Esteve, and J. M. Martinis, Science 239, 992–997 (1988).
  * Yu [1965] L. Yu, Acta Phys. Sin. 21, 75 (1965).
  * Shiba and Soda [1969] H. Shiba and T. Soda, Progress of Theoretical Physics 41, 25 (1969).
  * Rusinov [1969] A. I. Rusinov, JETP Lett. 9, 85 (1969), zh. Eksp. Teor. Fiz. Pisma Red. 9, 146 (1968).
  * Hewson [1993] A. C. Hewson,  _The Kondo Problem to Heavy Fermions_ (Cambridge University Press, 1993).
  * Rozhkov and Arovas [1999] A. V. Rozhkov and D. P. Arovas, Physical Review Letters 82, 2788 (1999).
  * Clerk and Ambegaokar [2000] A. A. Clerk and V. Ambegaokar, Physical Review B 61, 9109 (2000).
  * Vecino  _et al._ [2003] E. Vecino, A. Martín-Rodero, and A. L. Yeyati, Physical Review B 68, 10.1103/physrevb.68.035105 (2003).
  * Choi  _et al._ [2004] M.-S. Choi, M. Lee, K. Kang, and W. Belzig, Physical Review B 70, 10.1103/physrevb.70.020502 (2004).
  * Oguri  _et al._ [2004] A. Oguri, Y. Tanaka, and A. C. Hewson, Journal of the Physical Society of Japan 73, 2494 (2004).
  * van Dam  _et al._ [2006] J. A. van Dam, Y. V. Nazarov, E. P. A. M. Bakkers, S. D. Franceschi, and L. P. Kouwenhoven, Nature 442, 667 (2006).
  * Žitko [2022] R. Žitko, Josephson potentials for single impurity Anderson impurity in a junction between two superconductors (2022).
  * Bargerbos  _et al._ [2022] A. Bargerbos, M. Pita-Vidal, R. Žitko, J. Ávila, L. J. Splitthoff, L. Grünhaupt, J. J. Wesdorp, C. K. Andersen, Y. Liu, L. P. Kouwenhoven, R. Aguado, A. Kou, and B. van Heck, PRX Quantum 3, 10.1103/prxquantum.3.030311 (2022).
  * Sahu  _et al._ [2023] M. R. Sahu, F. J. Matute-Cañadas, M. Benito, P. Krogstrup, J. Nygård, M. F. Goffman, C. Urbina, A. L. Yeyati, and H. Pothier, Ground state phase diagram and ”parity flipping” microwave transitions in a gate-tunable josephson junction (2023), arXiv:2312.12914 [cond-mat.mes-hall] .
  * Bulaevskii  _et al._ [1978] L. Bulaevskii, V. Kuzii, and A. Sobyanin, Solid State Communications 25, 1053 (1978).
  * Vion  _et al._ [2002] D. Vion, A. Aassime, A. Cottet, P. Joyez, H. Pothier, C. Urbina, D. Esteve, and M. H. Devoret, Physica Scripta T102, 162 (2002).
  * Vion  _et al._ [2003] D. Vion, A. Aassime, A. Cottet, P. Joyez, H. Pothier, C. Urbina, D. Esteve, and M. Devoret, Fortschritte der Physik 51, 462 (2003).
  * Della Rocca  _et al._ [2007] M. L. Della Rocca, M. Chauvin, B. Huard, H. Pothier, D. Esteve, and C. Urbina, Physical Review Letters 99, 127005 (2007).
  * Ginzburg  _et al._ [2018] L. V. Ginzburg, I. E. Batov, V. V. Bol’ginov, S. V. Egorov, V. I. Chichkov, A. E. Shchegolev, N. V. Klenov, I. I. Soloviev, S. V. Bakurskiy, and M. Y. Kupriyanov, JETP Letters 107, 48–54 (2018).
  * Bargerbos  _et al._ [2023] A. Bargerbos, M. Pita-Vidal, R. Žitko, L. J. Splitthoff, L. Grünhaupt, J. J. Wesdorp, Y. Liu, L. P. Kouwenhoven, R. Aguado, C. K. Andersen, A. Kou, and B. van Heck, Phys. Rev. Lett. 131, 097001 (2023).
  * Pavešić  _et al._ [2023] L. Pavešić, M. Pita-Vidal, A. Bargerbos, and R. Žitko, SciPost Physics 15, 2542 (2023).
  * Richardson and Sherman [1964] R. W. Richardson and N. Sherman, Nuclear Physics 52, 221 (1964).
  * Dukelsky and Sierra [2000] J. Dukelsky and G. Sierra, Physical Review B 61, 12302 (2000).
  * Pavešić  _et al._ [2021] L. Pavešić, D. Bauernfeind, and R. Žitko, Physical Review B 104, 10.1103/physrevb.104.l241409 (2021).
  * Román  _et al._ [2002] J. Román, G. Sierra, and J. Dukelsky, Nuclear Physics B 634, 483–510 (2002).
  * Dukelsky  _et al._ [2004] J. Dukelsky, S. Pittel, and G. Sierra, Reviews of Modern Physics 76, 643 (2004).
  * von Delft and Ralph [2001] J. von Delft and D. Ralph, Physics Reports 345, 61 (2001).
  * Žitko and Pavešić [2022] R. Žitko and L. Pavešić, Physical Review B 106, 10.1103/physrevb.106.024513 (2022).
  * Pavešić and Žitko [2022] L. Pavešić and R. Žitko, Physical Review B 105, 10.1103/physrevb.105.075129 (2022).
  * Note [1] For detailed derivations and expressions at general filling see Ref. flatband.
  * Belyaev [1961] S. T. Belyaev, JETP Letters 12, 968 (1961).
  * Khodel’ and Shaginyan [1990] V. A. Khodel’ and V. R. Shaginyan, JETP Letters 51, 488 (1990).
  * Volovik [2018] G. E. Volovik, JETP Letters 107, 516 (2018).
  * Matute-Cañadas  _et al._ [2023] F. J. Matute-Cañadas, L. Tosi, and A. L. Yeyati, Quantum circuits with multiterminal josephson-andreev junctions (2023), arXiv:2312.17305 [cond-mat.mes-hall] .
  * Wolfram Research, Inc. [2023] Wolfram Research, Inc., Mathematica, Version 13.3 (2023), Champaign, IL, 2023.
  * Pavešič [2023] L. Pavešič, Flat band two channel, https://github.com/PavesicL/flat_band_two_channel (2023).
  * Devoret [2021] M. H. Devoret, Journal of Superconductivity and Novel Magnetism 34, 1633–1642 (2021).
  * Affleck  _et al._ [2000] I. Affleck, J.-S. Caux, and A. M. Zagoskin, Physical Review B 62, 1433 (2000).
  * Bergeret  _et al._ [2007] F. S. Bergeret, A. L. Yeyati, and A. Martín-Rodero, Physical Review B 76, 10.1103/physrevb.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, Nature Communications 9, 10.1038/s41467-018-04683-x (2018).
  * Žonda  _et al._ [2022] M. Žonda, P. Zalom, T. Novotný, G. Loukeris, J. Bätge, and V. Pokorný, Generalized atomic limit of a double quantum dot coupled to superconducting leads (2022).
  * Hermansen  _et al._ [2022] C. Hermansen, A. L. Yeyati, and J. Paaske, Physical Review B 105, 10.1103/physrevb.105.054503 (2022).
  * Schmid  _et al._ [2022] H. Schmid, J. F. Steiner, K. J. Franke, and F. von Oppen, Physical Review B 105, 10.1103/physrevb.105.235406 (2022).
  * Pavešič  _et al._ [2023] L. Pavešič, R. Aguado, and R. Žitko, Strong-coupling theory of quantum dot josephson junctions: role of the residual quasiparticle (2023).
  * Andreev [1964] A. F. Andreev, Zh. Eksperim. i Teor. Fiz. 46,  (1964).
  * Koch  _et al._ [2022] C. P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, EPJ Quantum Technology 9, 10.1140/epjqt/s40507-022-00138-x (2022).
  * Rossignolo  _et al._ [2023] M. Rossignolo, T. Reisser, A. Marshall, P. Rembold, A. Pagano, P. J. Vetter, R. S. Said, M. M. Müller, F. Motzoi, T. Calarco, F. Jelezko, and S. Montangero, Computer Physics Communications 291, 108782 (2023).
  * Golovach  _et al._ [2006] V. N. Golovach, M. Borhani, and D. Loss, Physical Review B 74, https://doi.org/10.1103/PhysRevB.74.165319 (2006).
  * Nowack  _et al._ [2007] K. C. Nowack, F. H. L. Koppens, Y. V. Nazarov, and L. M. K. Vandersypen, Science 318, 1430 (2007).
  * Pioro-Ladrière  _et al._ [2008] M. Pioro-Ladrière, T. Obata, Y. Tokura, Y.-S. Shin, T. Kubo, K. Yoshida, T. Taniyama, and S. Tarucha, Nature Physics 4, 776 (2008).
  * Nadj-Perge  _et al._ [2010] S. Nadj-Perge, S. M. Frolov, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Nature 468, 1084 (2010).
  * van den Berg  _et al._ [2013] J. W. G. van den Berg, S. Nadj-Perge, V. S. Pribiag, S. R. Plissard, E. P. A. M. Bakkers, S. M. Frolov, and L. P. Kouwenhoven, Physical Review Letters 110, https://doi.org/10.1103/PhysRevLett.110.066806 (2013).
  * Li  _et al._ [2022] B. Li, T. Calarco, and F. Motzoi, PRX Quantum 3, 10.1103/prxquantum.3.030313 (2022).