# Supplementary information for Non-topological zero bias peaks in full-shell
nanowires induced by flux tunable Andreev states

Marco Valentini1∗, Fernando Peñaranda2, Andrea Hofmann1, Matthias Brauns1,
Robert Hauschild1, Peter Krogstrup3, Pablo San-Jose2, Elsa Prada4, Ramón
Aguado2∗, Georgios Katsaros1<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
1Institute of Science and Technology Austria, Am Campus 1, 3400
Klosterneuburg, Austria.
2Instituto de Ciencia de Materiales de Madrid (ICMM), Consejo Superior de
Investigaciones Científicas (CSIC), Sor Juana Inés de la Cruz 3, 28049 Madrid,
Spain. Research Platform on Quantum Technologies (CSIC).
3Microsoft Quantum Materials Lab and Center for Quantum Devices, Niels Bohr
Institute, University of Copenhagen, Kanalvej 7, 2800 Kongens Lyngby, Denmark.
4Departamento de Física de la Materia Condensada, Condensed Matter Physics
Center (IFIMAC) and Instituto Nicolás Cabrera, Universidad Autónoma de Madrid,
E-28049 Madrid, Spain

##  I Model

The theoretical model used to describe the experimental results of this work
is presented in what follows. It consists of a single impurity Anderson model
(SIAM) coupled to a superconducting reservoir, the so-called superconducting
Anderson model, in which we take into account the Little-Parks (LP) modulation
of the gap with the applied flux, $\phi$ (see the details below). The impurity
is also coupled to a normal reservoir, hence forming a superconductor-quantum
dot-normal (SC-QD-N) junction, like in the experimental samples. Despite its
simplicity, this model correctly captures the interplay between two competing
mechanisms, namely: Coulomb blockade in the QD region and Andreev processes at
the superconducting interface. Each of them favors a different ground state
(GS) configuration: the former favors a doublet phase, $|D\rangle$, with spin
$\frac{1}{2}$, whereas the latter a singlet phase, $|S\rangle$, with spin 0.

The Hamiltonian of the system reads:
$H=H_{QD}+H_{N}+H_{S}+H_{T}^{N}+H_{T}^{S}$. Here, $H_{QD}$ refers to the
decoupled dot Hamiltonian given by

$H_{QD}=\sum_{\sigma}\epsilon_{\sigma}d^{\dagger}_{\sigma}d_{\sigma}+Un_{\uparrow}n_{\downarrow},$
(1)

where $d_{\sigma}$ ($d_{\sigma}^{\dagger}$) destroys (creates) an electron
with spin $\sigma=\\{\uparrow,\downarrow\\}$ and energy
$\epsilon_{\sigma}=\epsilon_{0}\mp V_{Z}$, where $V_{Z}=g\mu_{B}B/2$ refers to
the Zeeman energy owing to the external magnetic field $B$ (with $g$ and
$\mu_{B}$ being, respectively, the gyromagnetic factor and the Bohr’s
magneton). Electron-electron correlations in the dot occupations
$n_{\sigma}=d^{\dagger}_{\sigma}d_{\sigma}$ are included by means of the
charging energy $U$. The Hamiltonian of the normal reservoir is given by

$H_{N}=\sum_{k_{N}\sigma}\varepsilon_{k_{N}\sigma}c^{\dagger}_{k_{N}\sigma}c_{k_{N}\sigma},$
(2)

where $c^{\dagger}_{k_{N}\sigma}$ ($c_{k_{N}\sigma}$) creates (destroys) an
electron with momentum $k_{N}$, spin $\sigma$ and single-particle energy
$\varepsilon_{k_{N}\sigma}$ in the normal reservoir.

Since the superconducting lead (full-shell proximitized nanowire) to which the
dot is coupled has a doubly connected geometry and is longitudinally threaded
by a magnetic flux, $\phi=\pi R^{2}B$ (with $R$ being the nanowire radius),
the LP effect will manifest in the form of a modulation of the
Bardeen–Cooper–Schrieffer (BCS) gap, $\Delta(\phi)$ vaitieknas2019anomalous ;
schwiete_2010 ; shah_2007 . Such modulation may result in a full suppression
of superconductivity in finite flux regions for cylinders characterized by
small radius and thickness $d$, as compared to the zero temperature
superconducting coherence length $\xi$; namely $d,R\ll\xi$ schwiete_2010 ,
which is the relevant experimental situation. The resulting Hamiltonian is
given by

$H_{S}=\sum_{k_{S}\sigma}\varepsilon_{k_{S}\sigma}c^{\dagger}_{k_{S}\sigma}c_{k_{S}\sigma}+\sum_{k_{S}}\Delta(\phi)(c^{\dagger}_{k_{S}\uparrow}c^{\dagger}_{-k_{S}\downarrow}+h.c.),$
(3)

where $c^{\dagger}_{k_{S}\sigma}$ ($c_{k_{S}\sigma}$) denotes the creation
(destruction) operator in the superconducting lead with momentum $k_{S}$, spin
$\sigma$ and energy $\varepsilon_{k_{S}\sigma}$. All energies are measured
with respect to the chemical potential of the SC.

As was shown in schwiete_2010 , a complete analogy can be established between
pair-breaking mechanisms due to flux-induced circulating supercurrents in the
shell and that of paramagnetic impurities in a disordered SC. Thus, the time-
dependent Ginzburg-Landau formalism used to describe the latter can be applied
to the former, giving rise to the following expression for the critical
temperature $T_{C}(\phi)$ (see larkin  for a comprehensive derivation):

$\displaystyle
0=\ln\left(\frac{T_{C}}{T_{C}^{0}}\right)+\psi\left(\frac{1}{2}+\frac{\Lambda_{n}(\phi)}{2\pi
T_{C}({\phi})}\right)-\psi\left(\frac{1}{2}\right),$ (4)

where $T_{C}^{0}$ is the critical temperature at zero flux, $n$ is the fluxoid
winding number, $\psi$ is the digamma function, and $\Lambda_{n}$ is the pair-
breaking term corresponding to a hollow superconducting cylinder. This pair-
breaking term can be derived solving the Ginzburg-Landau equations in the
presence of impurities Groff:PR , and yields

$\Lambda_{n}(\phi)=\frac{T_{C}^{0}}{\pi}\frac{\xi^{2}}{R^{2}}\Big{[}4\left(n-\frac{\phi}{\phi_{0}}\right)^{2}+\frac{d^{2}}{R^{2}}\left(\frac{\phi^{2}}{\phi_{0}^{2}}+\left(\frac{1}{3}+\frac{d^{2}}{20R^{2}}\right)n^{2}\right)\Big{]},$
(5)

where $\phi_{0}$ denotes the superconducting flux quantum, $\phi_{0}=h/2e$.
Once $T_{C}(\phi)$ is known, in small devices ($d,R\ll\xi$) the field-
dependent critical current $I_{C}(\phi)$ can be expressed as
$I_{C}(\phi)\propto T_{C}(\phi)^{3/2}$, thanks to a Ginzburg-Landau
approximated relation in the presence of a magnetic field Bardeen1962 . In the
ballistic regime, $\Delta(\phi)$ is proportional to $I_{C}(\phi)$ and thus

$\frac{\Delta(\phi)}{\Delta}\approx\left(\frac{T_{C}(\phi)}{T_{C}^{0}}\right)^{3/2},$
(6)

where $\Delta$ is the gap at $\phi=0$. Departures from a power-law dependence
of $\Delta(\phi)$ on $T_{C}(\phi)$ in the presence of pair-breaking mechanisms
are predicted to happen as we go beyond the Ginzburg-Landau/BCS description
abrikosov1961 ; gennes66 ; Skalski64 . However, Eq. (6) suffices to accurately
capture the LP modulation found in the experiment.

The remaining two terms in the Hamiltonian:

$H_{T}^{S}=\sum_{k_{S}\sigma}\left(V_{k_{S}}c^{\dagger}_{k_{S}\sigma}d_{\sigma}+h.c.\right),$
(7)

and

$H_{T}^{N}=\sum_{k_{N}\sigma}\left(V_{k_{N}}c^{\dagger}_{k_{N}\sigma}d_{\sigma}+h.c.\right),$
(8)

stand for the couplings between the dot and the SC lead and the normal
reservoir, respectively. These two tunneling processes define the two
following tunneling rates:
$\Gamma_{N,S}=\pi\sum_{k_{N,S}}|V_{k_{N,S}}|^{2}\delta(\omega-\epsilon_{k_{N,S}})=\pi|V_{{N,S}}|^{2}\rho_{{N,S}}$,
with $\rho_{{N,S}}$ being the density of states of the normal/superconducting
lead evaluated at the Fermi energy.

In order to compute spectral observables, the relevant quantity that has to be
determined is the equilibrium retarded Green’s function in the dot, namely
$\hat{\mathcal{G}}_{\sigma}(\omega)$. In a Nambu $2\times 2$ basis spanned by
$\Psi_{\sigma}=\\{d_{\sigma},d^{\dagger}_{-\sigma}\\}^{T}$, the Dyson equation
reads

$\hat{\mathcal{G}}_{\sigma}^{-1}(\omega,\phi)=\hat{\mathcal{G}}_{0\sigma}^{-1}(\omega,\phi)-\hat{\Sigma}_{\sigma}(\omega,\phi).$
(9)

Here, $\hat{\mathcal{G}}_{0\sigma}(\omega)$ refers to the decoupled and non-
interacting dot retarded Green’s function and
$\hat{\Sigma}_{\sigma}=\hat{\Sigma}_{\sigma}^{QD}+\hat{\Sigma}_{\sigma}^{N}+\hat{\Sigma}_{\sigma}^{S}$
is the self-energy that takes into account the coupling to the leads,
$\hat{\Sigma}_{\sigma}^{N,S}$, and the Coulomb interaction in the dot,
$\hat{\Sigma}_{\sigma}^{QD}$. While a complete description of the problem
needs to fully take into account electron-electron correlations, the Coulomb
blockade regime is well captured by considering a diagrammatic expansion of
the dot self-energy up to the lowest order in $U$:

$\hat{\Sigma}^{QD,HF}_{\sigma}=U\begin{pmatrix}\langle
n_{-\sigma}\rangle&\langle d_{\uparrow}d_{\downarrow}\rangle\\\ \langle
d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}\rangle&-\langle
n_{\sigma}\rangle\end{pmatrix}.$ (10)

This mean-field approximation, which neglects high-order correlations, is
known as the Hartree-Fock-Bogoliubov (HF) approximation PhysRevB.68.035105 ;
Rodero_2012 ; Lee2014 . This level of approximation is enough to achieve an
excellent agreement with the experimental results (away from the Kondo
regime). Within this mean field picture, the retarded dot Green’s function
reads

$\hat{\mathcal{G}}_{\sigma}^{r,HF}(\omega,\phi)=\frac{1}{D(\omega,\phi)}\begin{pmatrix}\omega+\epsilon_{-\sigma}+i\Gamma_{N}+\Gamma_{S}\frac{\omega}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}+U\langle
n_{\sigma}\rangle&\Gamma_{S}\frac{\Delta(\phi)}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}+U\langle
d_{\uparrow}d_{\downarrow}\rangle\\\
\Gamma_{S}\frac{\Delta(\phi)}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}+U\langle
d^{\dagger}_{\downarrow}d^{\dagger}_{\uparrow}\rangle&\omega-\epsilon_{\sigma}+i\Gamma_{N}+\Gamma_{S}\frac{\omega}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}-U\langle
n_{-\sigma}\rangle\end{pmatrix},$ (11)

where $D(\omega,\phi)$ is given by the determinant of
$[\hat{\mathcal{G}}_{\sigma}^{r,HF}]^{-1}(\omega,\phi)$. The Andreev bound
states (ABSs) are given by the zeroes of the determinant and can be written in
a compact form as:

$\displaystyle\left[\omega+i\Gamma_{N}+\Gamma_{S}\frac{\omega}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}\pm\tilde{V}_{Z}\right]^{2}=\tilde{\epsilon}^{2}+\left|\frac{\Gamma_{S}\Delta(\phi)}{\sqrt{\Delta(\phi)^{2}-\omega^{2}}}+\tilde{\Delta}\right|^{2},$
(12)

with $\tilde{V}_{Z}=V_{Z}+\frac{U}{2}(\langle n_{\uparrow}\rangle-\langle
n_{\downarrow}\rangle)$, $\tilde{\epsilon}=\epsilon_{0}+\frac{U}{2}(\langle
n_{\uparrow}\rangle+\langle n_{\downarrow}\rangle)$ and
$\tilde{\Delta}=U\langle d_{\uparrow}d_{\downarrow}\rangle$. Physically,
$\tilde{V}_{Z}$ and $\tilde{\epsilon}$ can be physically interpreted as an
effective exchange field and a level shift due to correlations. When
$\Delta\rightarrow 0$ they define two Coulomb resonances at
$\tilde{\epsilon}\pm\tilde{V}_{Z}$ and broadened by $\Gamma_{N}+\Gamma_{S}$.
When $\tilde{\epsilon}=0$, these two resonances are located at
$\pm(V_{Z}+U/2)$. In the presence of superconducting correlations, these
resonances become ABSs whose location and physical character depend on the
ratios of different parameters [$\Delta(\phi)/U$, $\Gamma_{S}/U$, etc.] that
give rise to different physical regimes, as we explain below.

Table S1:  Parameters used in the numerical simulations of the main text. $\Delta$ values were determined from experimental evidence. $R$ and $d$ have been extracted by using them as free parameters in Eqs. (5) and (6) and finding the best agreement with the experimental LP data. Going from $\phi$ to B we assumed a circular effective section of radius $R+10\rm{nm}$ which captures the finite penetration length of B in the SC shell and is required to achieve experimental agreement. | $U[\rm{meV}]$ | $\epsilon_{0}/U$ | $g$ | $\Delta[\rm{meV}]$ | $\Gamma_{S}/U$ | $\Gamma_{S}/\Delta$ | $\Gamma_{N}/\Delta$ |  $R$ [nm] |  $d$ [nm]
---|---|---|---|---|---|---|---|---|---
Fig. 4(E1) | 2.5 | -0.97 | 13 | 0.18 | 0.036 | 0.5 | 0.07 | 65 | 25
Fig. 4(E2) | 2.5 | -0.95 | 13 | 0.18 | 0.036 | 0.5 | 0.07 | 65 | 25
Fig. 4(E3) | 2.5 | -0.93 | 13 | 0.18 | 0.036 | 0.5 | 0.07 | 65 | 25
Fig. 5(G-H) | 1.5 | -0.93 | 13 | 0.14 | 0.075 | 0.8 | 0.1 | 64 | 23
Fig. S5(D) | 2.5 | -0.5 | 13 | 0.18 | 0.18 | 2.5 | 0.3 | 65 | 25
Fig. S9 | 1.5 | -0.91 | 14.3 | 0.18 | 0.12 | 1 | 0.16 | 50 | 35
Fig. S10(B) | 1.5 | -0.93 | 17 | 0.13 | 0.062 | 0.72 | 0.25 | 64 | 23

The average values concerning the dot level occupations in Eqs. (11) and (12)
are determined self-consistently in equilibrium. Once
$\hat{\mathcal{G}}_{\sigma}^{r,HF}(\omega)$ is known, we can easily calculate
the density of states (DOS)
$\rho(\omega,\phi)=-\frac{1}{\pi}{\rm{Im\,Tr}}[\mathcal{G}_{\sigma}^{r,HF}(\omega,\phi)]$.
At small bias voltage $V$ and temperature $T$, the differential conductance
$dI/dV$ in the tunneling limit provides an approximate measure of the QD DOS
at energy $eV$. The sets of parameters for each simulation in the main text
are collected in Table S2.

While the full solution of the zeroes in Eq. (12) is needed to find the ABSs
at arbitrary energy, one may find some useful analytic boundaries in different
parameter regimes, as we discuss now.

###  I.1 Large $U$ approximations

The experimental regime discussed in the main text corresponds to small
$\Delta(\phi)/U$ (large-$U$ limit). In this case, singlets are Yu-Shiba-
Rusinov (YSR)-like superpositions between the singly occupied state in the QD
and quasiparticles below the gap in the SC. The ABS position in this limit can
be much lower than the original gap and, assuming that $\tilde{U}$ is the
largest scale in the problem, namely that the system is in a doublet GS, we
can find the poles of Eq. (12) as ($\Gamma_{N}=0$):

$\displaystyle\omega\approx\pm\Delta(\phi)\frac{\tilde{V}^{2}_{Z}-\tilde{\Gamma}^{2}_{S}-\tilde{\epsilon}^{2}}{\sqrt{[\tilde{\Gamma}^{2}_{S}+\tilde{\epsilon}^{2}-\tilde{V}^{2}_{Z}]^{2}+4(\tilde{\Gamma}_{S}-\tilde{\Delta})^{2}\tilde{V}^{2}_{Z}}},$
(13)

which is a generalization of the classical expression for YSR states in a SC
Yu:APS65 ; Shiba:PTP68 ; Rusinov:SPJ69 ; RevModPhys.78.373  (see also Ref.
PhysRevB.92.235422  for a different derivation using an effective Kondo
model). Note that the poles cross zero energy when PhysRevB.91.045441

$\displaystyle\tilde{V}^{2}_{Z}=\tilde{\Gamma}^{2}_{S}+\tilde{\epsilon}^{2}$
(14)

We note that this criterion for parity crossings of YSR states in a QD is
essentially the same as the criterion for topological superconductivity and
Majorana zero modes in a proximitized nanowire (by just substituting
$\tilde{\Gamma}_{S}\rightarrow\Delta^{NW}$,
$\tilde{\epsilon}\rightarrow\mu^{NW}$ and $\tilde{V}_{Z}\rightarrow
V_{Z}^{NW}$ lutchyn_majorana_2010 ; oreg_helical_2010 . This is yet another
example that emphasises the difficulty of making a clear distinction between
parity crossings in QDs and Majorana zero modes in nanowires. The connection
between Eq. (13) and the standard expression for YSR states Yu:APS65 ;
Shiba:PTP68 ; Rusinov:SPJ69 ; RevModPhys.78.373  is evident in the limit
$\tilde{\epsilon}\rightarrow 0$, $\tilde{\Delta}\rightarrow 0$, where we
recover

$\omega=\pm\Delta(\phi)\frac{1-\left(\frac{\Gamma_{S}}{\tilde{V}_{Z}}\right)^{2}}{1+\left(\frac{\Gamma_{S}}{\tilde{V}_{Z}}\right)^{2}}=\pm\Delta(\phi)\frac{1-J^{2}}{1+J^{2}}.$
(15)

Using $\langle n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle=1$ in this
limit, we can recover $J=\pi jS\rho_{S}/2$ with $S=1/2$ and
$j=8|V_{S}|^{2}/(U+2V_{Z})$.

Typically, $U$ can be experimentally extracted from the size of Coulomb
blockade diamonds and we can use Eq. (15) to find the value of $\Gamma_{S}$ in
the doublet regime by just substituting the ABS energy position at $B=0$,
namely:

$\Gamma_{S}=\frac{U}{2}\sqrt{\frac{(1-\omega_{+}/\Delta)}{(1+\omega_{+}/\Delta)}}.$
(16)

Using this procedure, we find an excellent agreement between the analytical
YSR expression of Eq. (15) and the experiments in the full magnetic field
range [see Fig. 4(d) of the main text]. Further improvement in the agreement
can be found by performing full HF numerics directly from Eq. (12) [typically
using slightly smaller $\Gamma_{S}$ values as the ones directly extracted from
Eq. (16)]. Note that the way the magnetic field enters in Eqs. (13) and (15)
(through $\Delta(\phi)$ and $V_{Z}$) strongly deviates from a standard linear
Zeeman dispersion, which explains why the ABSs in the experiments show almost
no dispersion against magnetic field.

These YSR singlets eventually become Kondo singlets in the destructive LP
regions when $\Delta(\phi)\rightarrow 0$. One can estimate the flux at which
Kondo physics dominates over superconductivity by using the expression for the
$V_{Z}=0$ Kondo temperature at the symmetry point $T_{K}\approx
0.3\sqrt{\Gamma_{S}U}e^{-\frac{\pi U}{8\Gamma_{S}}}\approx
0.2U\sqrt{J}e^{-\frac{\pi}{4J}}$. For example, using the parameters of Fig.
4(g) of the main text, we can estimate that a Kondo GS develops when
$\Delta(\phi)\lesssim T_{K}\approx 0.044\;\rm{meV}\approx 0.24\Delta$, which
first happens at $\phi\approx 0.35\phi_{0}$, in very good agreement with the
experimental observation of a Kondo effect in the first destructive LP region
at $B\approx 50$ mT. Nevertheless, we emphasise that a full calculation of the
flux at which the Kondo effect develops constitutes a rather involved problem,
well beyond the scope of this theoretical analysis, and should be the subject
of future studies.

When approaching a singlet-doublet crossing near a charge-degeneracy point as
a function of $\epsilon_{0}$, the self-consistent HF method is known to
exhibit small discontinuities between non-magnetic and magnetic solutions
Rodero_2012 ; Lee2014 . In these cases [Fig. 5(h) of the main text], and since
these discontinuities can result in unphysical energy shifts of the ABSs near
the crossing, we find the mean-field values in Eq. (12) that provide the best
agreement with the experiment.

###  I.2 Large $\Delta$ approximations

For completeness, we also discuss the opposite large gap limit, where the
coupling to the SC $\Gamma_{S}$ mainly induces local superconducting
correlations in the QD. This can be easily seen by taking the so-called
$\Delta\rightarrow\infty$ atomic limit PhysRevB.91.045441 ; Bauer_2007 ;
PhysRevB.79.224521 , where superconducting singlets are just bonding and
antibonding superpositions of empty and doubly-occupied states in the QD,
namely $|S_{\pm}\rangle=|0\rangle\pm|2\rangle$. In this atomic limit, the
singlet-doublet boundary (parity crossing) can be easily obtained by taking
the $\omega\rightarrow 0$ limit in Eq. (12) ($V_{Z}=0$ and $\Gamma_{N}=0$):

$\left(\frac{U}{2}\right)^{2}=\Gamma_{S}^{2}+\left(\epsilon_{0}+\frac{U}{2}\right)^{2}.$
(17)

Expanding the determinant further to order $\omega^{2}$ allows to capture
finite $\Delta$ corrections to the above atomic limit boundary. This is called
the generalized atomic limit (GAL) GAL1 ; GAL2 , whose singlet-doublet
boundary now reads

$\frac{U^{2}}{4(1+\frac{\Gamma_{S}}{\Delta})^{2}}=\Gamma_{S}^{2}+\left(\epsilon_{0}+\frac{U}{2}\right)^{2}.$
(18)

It can be used to demonstrate that HF provides a very good approximate
solution of the problem, particularly near the symmetry point
$\epsilon_{0}=-U/2$, when neglecting the above-gap continuum contribution to
the anomalous propagator. Further improvement can be reached by performing
finite-$\Delta$ self-consistent perturbation theory around the atomic limit
both for $V_{Z}=0$ PhysRevB.79.224521  and $V_{Z}\neq 0$ PhysRevB.94.085151 ;
PhysRevResearch.2.012065 . The resulting boundaries are, however, no longer
analytic.

##  II Phase diagrams in different coupling regimes

Figure S1:  Gap modulation of the phase diagram through the Little-Parks
effect. (a) LP modulation of the superconducting gap in the destructive regime
vs normalized flux $\phi/\phi_{0}$ . (b-c) Full phase diagram of $\langle
n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle$ vs $\Gamma_{S}/U$ and
$\epsilon_{0}/U$ for $\phi=0$ in the small and large $\Delta/U$ regimes,
respectively. Only the LP effect is included here so we take $g=0$. Hartree-
Fock (HF) solution: Inside the dark blue dome the superconductor-quantum dot
(SC-QD) system is in a doublet ground state (GS), whereas it is in a singlet
GS in the white region. Grey, red and cyan dashed lines mark the analytic
boundaries in the $\Delta/U\rightarrow\infty$ (atomic), generalized atomic
(GAL) and $\Delta/U\rightarrow 0$ (Anderson) approximations, respectively.
(d-e) $\langle n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle$ vs
$\phi/\phi_{0}$ for fixed $\epsilon_{0}=-U/2$ and increasing values of
$\Gamma_{S}/U$ from 0.06 (black) to 0.36 (d) and 0.45 (e) (light grey). The
orange lines mark the destructive LP regime. Note that for the larger
couplings, the singlet phase occurs before the LP destructive region.
Similarly, the system is back to a doublet well outside the LP destructive
region. (f-g) $\langle n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle$ vs
$\Gamma_{S}/U$ for fixed $\epsilon_{0}=-U/2$ and increasing fluxes marked with
colours in (a) for the two regimes of (b) and (c), respectively, also in the
absence of Zeeman effect. The LP effect results in an oscillatory modulation
of the dome. Parameters: $g=0$, $\Gamma_{N}\approx 0$, $\Delta=0.2$ meV,
$R=50$ nm, $d=35$ nm, and $U=0.5$ meV ($U=0.05$ meV) in (b), (d) and (f) [in
(c), (e) and (g)].  Figure S2:  Gap plus Zeeman modulation of the phase
diagram owing to external flux.  (a,b) Full phase diagram of $\langle
n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle$ vs $\Gamma_{S}/U$ and
$\epsilon_{0}/U$ at $\phi/\phi_{0}=0.9$ ($\phi/\phi_{0}=0.3$) in the small
(large) $\Delta/U$ regime. Hartree-Fock (HF) solution: Inside the dark blue
dome the SC-QD system is in a doublet ground state (GS), whereas it is in a
singlet GS in the white region. Grey, red and cyan dashed lines mark the
analytic boundaries in the $\Delta/U\rightarrow\infty$ (atomic), generalized
atomic (GAL) and $\Delta/U\rightarrow 0$ (Anderson) approximations,
respectively. (c,d) $\langle n_{\uparrow}\rangle-\langle
n_{\downarrow}\rangle$ vs $\Gamma_{S}/U$ for fixed $\epsilon_{0}=-U/2$ in the
small and large $\Delta/U$ regimes, respectively, and increasing fluxes marked
with colors in (a). The dashed cyan line in (c) corresponds to the $B=0$,
$\Delta/U\rightarrow 0$ limit, whereas the dashed grey line in (d) corresponds
to the $B=0$, $\Delta/U\rightarrow\infty$ limit. In this figure, both the
Zeeman effect, through a non-zero $g=11$, and the LP effect are included, the
rest of parameters are chosen as those of Fig. S1. The doublet dome exhibits a
non trivial evolution with the applied flux due to two different mechanisms:
the Zeeman effect that expands the doublet dome in both $x$ and $y$ directions
and the LP effect which leads to an oscillatory modulation.

Whether the SC-QD GS is a singlet or a doublet is governed by the magnetic
field evolution of the subgap states, with a change in GS parity whenever they
cross zero energy. This in turn results from the competition of two
mechanisms: LP effect that modulates the gap and the Zeeman effect that shifts
the dot level energy. The first effect is illustrated in Fig. S1(a), where we
show the oscillatory modulation of the superconducting gap vs the applied
normalized flux $\phi/\phi_{0}$. The gap exhibits the so-called destructive LP
effect, whereby the gap decreases in the zeroth lobe from its maximum value at
$\phi=0$, disappears in the first destructive region, reemerges in the fist
lobe, and so on. In Figs. S1(b,c) we analyze the role of LP modulation of the
gap on the phase diagram of the problem, disregarding for the moment the
Zeeman effect (i.e., taking $g=0$). We plot the difference of spin-polarised
occupations $\langle n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle$ in the
self-consistent HF approximation versus $\Gamma_{S}/U$ and $\epsilon_{0}/U$.
Both panels show the phase diagram for the largest gap of the problem,
$\Delta$ (i.e., the gap at $\phi=0$), but for different charging $U$,
resulting in small, Fig. S1(b), and large, Fig. S1(c), $\Delta/U$ ratios.
Already at this level, it it important to stress how the different ratios
influence the singlet-doublet boundary. The gap evolution across LP
oscillations (without Zeeman) modifies the singlet-doublet boundary most
strongly at the symmetry point $\epsilon_{0}=-U/2$. This LP modulation of the
boundary is clearly seen in Figs. S1(d,e), where we show how the spin
polarisation varies with flux at the symmetry point $\epsilon_{0}=-U/2$ for
increasing $\Gamma_{S}$ (from top to bottom curves). Both panels show that
doublet-singlet-doublet transitions can be induced by LP modulation. Note that
while our HF does not provide a quantitative result at doublet phases (the
spin should be fully screened in this phase, $\langle
n_{\uparrow}\rangle-\langle n_{\downarrow}\rangle\rightarrow 0$, owing to
Kondo correlations), the overall trend is perfectly captured. These non-
monotonic doublet-singlet-doublet transitions are, in turn, responsible for
the phenomenology discussed in Fig. 4 of the main text. In this regime, the
measured excitations in the experiments clearly show a YSR singlet-Kondo (co-
tunneling) excitation-YSR singlet transition mediated by flux. Note that for
the largest coupling $\Gamma_{S}$ shown in the plot, the boundaries occur well
before/after the destructive LP region (delimited by orange lines). In this
latter case, the singlet phase extends well beyond the destructive LP region.
This could result in a Kondo effect and a singlet-doublet parity crossing in
flux regions well inside the first lobe where Majoranas are also predicted.
Apart from a non-monotonic dependence of the boundary on $\phi$, owing to the
non-monotonic $\Delta(\phi)$, we clearly observe how the two analytic
boundaries of the problem, the normal Anderson model boundary at
$\Gamma_{S}/U=1/\pi$ in the $\Delta\rightarrow 0$ (cyan dashed line) and the
$\Delta/U\rightarrow\infty$ boundary at $\Gamma_{S}/U=0.5$ (grey dashed line),
are approached. The first limit, in particular, is approached in the
destructive LP phases where the superconducting gap fully closes. As we
mentioned above, a Kondo effect (not captured here) should develop when
$T_{K}\gtrsim\Delta(\phi)$ in these Anderson regions doi:10.1143/JPSJ.69.1812
; Bauer_2007 ; PhysRevB.95.180502 .

Having discussed the non-trivial role that LP modulation has on the singlet-
doublet boundary, we now take into account also the Zeeman effect. The full
phenomenology showing the interplay of LP and Zeeman effects induced by $\phi$
is shown in Fig. S2. The upper plots show the full phase diagram for small
[Fig. S2(a)] and large [Fig. S2(b)] $\Delta/U$ ratios, respectively. The role
of Zeeman is now evident with growing doublet regions as $\phi$ increases. As
expected, for the same $g$-factor (which we take as $g=11$ here) the degree of
spin-polarisation of the GS also depends on whether $U$ is much larger [Fig.
S2(a)] or comparable [Fig. S2(b)] to the typical Zeeman energy scales. The
latter limit, in particular, results in a greatly enhanced doublet dome even
for moderate fluxes. This can be easily understood already from the
$\Delta\rightarrow\infty$ atomic limit which now reads

$\left(V_{Z}+\frac{U}{2}\right)^{2}=\Gamma_{S}^{2}+\left(\epsilon_{0}+\frac{U}{2}\right)^{2}.$
(19)

The lower panels, Figs. S2(c,d), show typical cuts for increasing flux at the
symmetry point $\epsilon_{0}=-U/2$. Important for our discussion is that the
singlet-doublet boundary for $\phi=0$ in the small $\Delta/U$ regime, Fig.
S2(c), moves to considerably larger $\Gamma_{S}/U$ ratios (see e.g. the
$\phi=0.9\phi_{0}$ boundary) due to the reentrant LP gap combined with the
increased doublet phase owing to Zeeman. This overall increase of the doublet
phase near $\phi\approx\phi_{0}$ greatly favours singlet-doublet parity
crossings in the first LP lobe as we discuss now.

Figure S3:  Suppressed Zeeman splitting in the zeroth lobe and ZBAs induced by
singlet-doublet parity crossings. (a,b) Poles of Eq. (11) shown in red
superimposed to the numerical DOS as a function of the bias voltage $V$ and
$\epsilon_{0}/U$ for two different normalized fluxes $\phi/\phi_{0}$ within
the zeroth LP lobe, namely 0 and 0.2, respectively. (c) Poles and numerical
DOS versus $V$ and $\phi/\phi_{0}$ at $\epsilon_{0}=-U/2$. $U=1$ meV,
$\Gamma_{S}=0.43$ meV, and $\Gamma_{N}=0.2$ meV. The rest of parameters are
the same as in Fig. 4(H1) (see Table S1). Within the zeroth lobe, the LP
modulation of the gap together with a small Zeeman energy (even for a sizeable
$g=13$) lead to a Zeeman splitting that cannot be experimentally resolved, as
shown in (b). Using the above parameters $V_{Z}/\Gamma_{N}=0.45$.

##  III Parity crossings versus flux $\phi$

As it is clear from the previous discussion, the fate of the GS character with
external flux depends on the initial $\phi=0$ position of the system
parameters in the $\epsilon_{0}/U$-$\Gamma_{S}/U$ plane. If these parameters
are far away from the singlet-doublet boundary, the system will remain in the
same GS for moderate fluxes. This is the case of Figs. 5(a),(b) and (d) of the
main text where the system remains in a singlet GS for fluxes deep into the
first lobe. Since the GS is a singlet, a further effect, apart from LP
modulation, is the Zeeman splitting of the subgap excitations. However, such
Zeeman splitting is not clearly seen in the experiment across the zeroth lobe.
Figure S3 shows theoretical calculations in this regime. In agreement with the
experiment, the Zeeman splitting in the zeroth lobe cannot be resolved. This
can be explained by the combined effect of the LP gap closing (whimch strongly
modifies the linear Zeeman dispersion at low $B$ fields), the small Zeeman
energy (even for a sizeable $g=13$) and the tunneling broadening of the lines
(compare poles with the full DOS). In Fig. S3(b) $V_{Z}/\Gamma_{N}=0.45$.

When the initial point is relatively close to the dome boundary, one or more
singlet-doublet transitions will probably occur at moderate fluxes. In
particular, Fig. 6(i) in the main text discusses how we approach the singlet-
doublet transition that we present in Figs. 6(a,c,e) (the corresponding point
in parameter space for said figure is marked by a yellow circle). The system
starts in a singlet GS for $\phi=0$. As the applied flux increases, the dome
approaches the yellow point which crosses the singlet-doublet line at
$\phi\approx\phi_{0}$, remaining in the doublet phase afterwards. We note
that, if we had a system with an initial singlet state slightly different from
the one of the yellow point, but same $g$-factor, it may cross the dome at
very different values of $\phi$. For example, experimental points in this
phase diagram with the same $\Gamma_{S}$ but slightly different gate voltage
outside the region $\epsilon_{0}/U\in[-1,0]$ (red circle), or for the same
gate voltage ($\epsilon_{0}$) but larger coupling to the SC $\Gamma_{S}$
(light blue circle), are still in the singlet phase for the largest $\phi$
shown here. This exemplifies how QD physics can lead to $\phi$-induced parity
crossings, and concomitant zero bias anomalies in transport, at very different
points in the parameter space of this complex interplay of different physical
phenomena.

The simulation for the specific parity crossing of Fig. 6(a) is shown in Figs.
6(g,h), where we plot the poles of Eq. (11), shown in red, superimposed to the
numerical DOS as a function of bias voltage $V$ and magnetic flux.

##  IV Additional experimental data

###  IV.1 Materials and methods

The hexagonal InAs NWs, with a diameter of about
$110\,$\mathrm{n}\mathrm{m}$$, were grown via the VLS technique and the
$30\,$\mathrm{n}\mathrm{m}$$ Al shell was epitaxially grown in situ, giving a
highly transparent semiconductor-SC interface chang_hard_2015 ; krogstrup_2015
. The NWs were deposited on a heavily doped silicon substrate covered with
$285\,$\mathrm{n}\mathrm{m}$$ of silicon oxide. The Al shell is completely
removed from one side of the NW via transene etching and the bare InAs is
employed as tunnel barrier. The Al shell, $L$, is longer than
$1.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ for all the measured devices.
Hence, effects of overlapping potential MZMs could be neglected as $L$ is much
larger than the coherence length sarma2012splitting ;
albrecht_exponential_2016 . The gates and contacts consist of normal-metal
($5\,$\mathrm{n}\mathrm{m}$$/$180\,$\mathrm{n}\mathrm{m}$$ Ti/Au bilayer). The
differential conductance, $dI/dV$, was measured via a lock-in technique in a
dilution refrigerator with a base temperature of $20\,$\mathrm{m}\mathrm{K}$$.

Figure S4:  Devices with a short junction. In the right-upper corner of each plot, a letter indicates the device. (A), (C), (E), (G), (I) , (K), (M) and (O) $dI/dV$ as a function of $V$ and $V_{\textrm{bg}}$ for devices with a short junction in absence of magnetic fields. No ABSs are present for these devices because of the short junction. More concretely, device F has $X\approx 95\,$\mathrm{n}\mathrm{m}$$, device G has $X\approx 60\,$\mathrm{n}\mathrm{m}$$, device H has $X\approx 85\,$\mathrm{n}\mathrm{m}$$, device I has $X\approx 70\,$\mathrm{n}\mathrm{m}$$, device J has $X\approx 90\,$\mathrm{n}\mathrm{m}$$, device K has $X\approx 60\,$\mathrm{n}\mathrm{m}$$, device L has $X\approx 80\,$\mathrm{n}\mathrm{m}$$ and device M has $X\approx 40\,$\mathrm{n}\mathrm{m}$$. (B), (D), (F), (H), (J), (L), (N) and (P) $dI/dV$ as a function of $V$ and $B$ at a fixed value of $V_{\textrm{bg}}$. The magnetic field scans have been taken at $V_{\textrm{bg}}=-12\,$\mathrm{V}$$ for device f, at $V_{\textrm{bg}}=-32\,$\mathrm{V}$$ for device G, at $V_{\textrm{bg}}=-12\,$\mathrm{V}$$ for device H, at $V_{\textrm{bg}}=-7.5\,$\mathrm{V}$$ for device I, at $V_{\textrm{bg}}=-19\,$\mathrm{V}$$ for device J, at $V_{\textrm{bg}}=-15.5\,$\mathrm{V}$$ for device K, at $V_{\textrm{bg}}=-5\,$\mathrm{V}$$ for device L and at $V_{\textrm{bg}}=-30\,$\mathrm{V}$$ for device M. The different shapes and positions of the lobes are due to the slightly different NW and Al shell dimensions. However, neither ZBPs nor other subgap features are observed. Table S2: Dimensions of short junction devices. Dimension extracted from the lobes using an hexagonal shape vaitieknas2019anomalous . In particular, $R$ is the distance between the center of the NW and the semiconductor facet, as indicated in Fig. 1. Name of the device |  $R$ [nm] |  $d$ [nm] |  $\xi$ [nm]
---|---|---|---
A | 61 | 23 | 190
F | 59 | 32 | 190
G | 71 | 20 | 155
H | 62 | 26 | 159
I | 61 | 26 | 174
J | 60 | 27 | 180
K | 62 | 24 | 160
L | 61 | 28 | 169
M | 63 | 24 | 163

###  IV.2 Kondo effect and g-factor

As mentioned, YSR singlets are the superconducting counterparts of Kondo
singlets and a full screening of the quantum impurity spin is expected when
the exchange between the QD and the SC is large enough, such that
$k_{B}T_{K}\approx\Delta$. While previous works have studied this doublet-YSR
singlet quantum phase transition (QPT) by a controlled tuning of $\Gamma_{S}$
Lee:PRB17, we here can reach a Kondo regime owing to the destructive LP
regions where the gap is suppressed. From Eq. 1 in the main text, we can
extract an effective exchange coupling $J=2\Gamma_{S}/U$ (at zero magnetic
field), which enters the expression of the Kondo temperature
$k_{B}T_{K}\approx 0.3\sqrt{U\Gamma_{S}}e^{-\frac{\pi U}{8\Gamma_{S}}}\approx
0.2U\sqrt{J}e^{-\frac{\pi}{4J}}$ at the symmetry point [the dashed line in
Fig. 3A]. Using the experimental values for $U$ and $\Gamma_{S}$, we get
$k_{B}T_{K}\approx 0.044\;\rm{meV}\approx 0.24\Delta(\phi=0)$. Owing to the LP
reduction of the gap, a full Kondo effect is expected to develop when
$k_{B}T_{K}\approx\Delta(\phi)$ doi:10.1143/JPSJ.69.1812 , which first happens
at $\phi\approx 0.35\phi_{0}$. Indeed, we find a symmetric pair of excitations
that appear in the first destructive LP regime, starting at $\approx 50$mT
which corresponds to $0.33\phi_{0}$ [cyan trace in Fig. S5C], that we
interpret as cotunneling processes through a spin-split Kondo resonance,
supporting the previous picture. At larger fluxes, the gap reopens, and the
system goes back to a doublet GS with YSR singlet subgap excitations. This
remarkable modulation from YSR singlets to Kondo and back is only possible due
to the destructive LP effect. This YSR-Kondo QPT mediated by the LP effect
occurs for typical QD parameters around the same flux regions where Majoranas
are predicted to occur in full-shell nanowires Vaitiekenaseaav3392 ;
fern2019evenodd , which adds further complexity to the problem. This Kondo
interpretation is supported by the decreasing height of the conductance peaks
for increasing temperatures [Fig. S5E]. As expected, the hard superconducting
gap precludes the formation of a Kondo resonance, which is absent in the
zeroth and first lobe. On the contrary, the soft gap in the second lobe allows
the formation of an above-gap split Kondo resonance, which, in this regime,
coexists with the YSR states Lee:PRL12; ramon2015  [see line traces in Fig.
S5C]. The progression of the Kondo peaks in magnetic field allows for an
estimation of the NW’s $g$-factor. The extracted values are $g=6.4\pm 0.8$ in
the first destructive region, $g=9.2\pm 0.3$ in the second destructive one and
$g=10.5\pm 0.1$ from the end of the second lobe to
$B=900\,$\mathrm{m}\mathrm{T}$$ [Fig. S6E]. The data indicate a non linear
behaviour of $g$ with the magnetic field. This apparent non-linearity relates
to the LP effect as the flux which can penetrate the core of the NW depends on
the magnetic field value. Nonlinearities of the effective $g$-factor owing to
Kondo interactions in a superconducting system are also expected ramon2015 .

[figure]capposition=bottom

Figure S5:  Kondo effect in the Little-Parks regime with a doublet GS from
device B. (A) Same plot as Fig. 3D but without the analytical solution
superimposed. (B) and (C) Line traces extracted from (A) for different values
of $B$. The line traces in (B) show that ABSs in the doublet GS, indicated
with 1, increase their energy with magnetic field. In the absence of magnetic
field, $\zeta\approx\pm 63\,$\mathrm{\SIUnitSymbolMicro}\mathrm{V}$$ (deep
blue trace), while at $B=110\,$\mathrm{m}\mathrm{T}$$ (green trace)
$\zeta\approx\pm 80\,$\mathrm{\SIUnitSymbolMicro}\mathrm{V}$$, even though the
superconducting gap is smaller at finite $B$. (C) Line traces showing the
presence of a split Kondo peak (K) in the destructive regimes, namely at
$B=60\,$\mathrm{m}\mathrm{T}$$ (cyan trace), at
$B=170\,$\mathrm{m}\mathrm{T}$$ (purple trace) and at
$B=280\,$\mathrm{m}\mathrm{T}$$ (red trace). Furthermore, the Kondo peaks
coexist with superconductivity and ABSs in the 2L (orange trace taken at
$B=230\,$\mathrm{m}\mathrm{T}$$). (D) Numerical simulation of the $dI/dV$
versus bias voltage $V$ and magnetic flux $\phi/\phi_{0}$ in a QD-S (modelled
as a superconducting Anderson model in the Hartree-Fock approximation), with
the SC in the destructive LP regime and the QD in a doublet GS (parameters can
be found in Table S1 of the SI). (E) $dI/dV$ as a function of $V$ at
$B=60\,$\mathrm{m}\mathrm{T}$$ for different temperatures. The split Kondo
peaks broaden and decrease their conductance forming a ZBP at
$T=240\,$\mathrm{m}\mathrm{K}$$, while they completely disappear for
$T=2817\,$\mathrm{m}\mathrm{K}$$. (F) Current, $I$, as a function of $V$ and
$V_{\textrm{bg}}$ for an out-of-plane magnetic field of
$50\,$\mathrm{m}\mathrm{T}$$ that destroys superconductivity revealing Coulomb
diamonds. The charging energy is $\approx
2.5\,$\mathrm{m}\mathrm{e}\mathrm{V}$$.

[figure]capposition=beside,capbesideposition=center,right

Figure S6:  Kondo data for higher magnetic fields, from device B. (A) $dI/dV$
as a function of $V$ and $B$ for $V_{\textrm{bg}}=-2.62\,$\mathrm{V}$$, hence
in the doublet GS. (B) Line traces extracted from (A) for
$B=550\,$\mathrm{m}\mathrm{T}$$ (blue trace) and for
$B=800\,$\mathrm{m}\mathrm{T}$$ (green trace).

[figure]capposition=bottom

Figure S7:  ZBP in the second lobe due to ABSs in device B. (A) Same plot as
Fig. 4D but withcoloured markers. (B) ABSs splitting in the 1L at positive
bias versus $B$, (upper panel) and effective g-factor $g_{\rm{eff}}$ extracted
from their splitting (lower panel). (C) Line traces extracted from (A) for the
0L, 1L and 2L (from left to right) illustrating how the ABSs merge at zero
bias at the end of the 1L. The lower trace of the 0L is taken at $B=0$ and the
upper one is taken at $B=42\,$\mathrm{m}\mathrm{T}$$. They are shifted by
$0.3{e^{2}}/{h}$ with respect to each other. For the 1L, the line traces are
shifted by 0.12 $e^{2}/{h}$ and they range from $B=82\,$\mathrm{m}\mathrm{T}$$
(lower one) to $B=152\,$\mathrm{m}\mathrm{T}$$ (upper one). Finally, for the
2L, the line traces are shifted by $0.04{e^{2}}/{h}$ and range from
$B=200\,$\mathrm{m}\mathrm{T}$$ (lower one) to $B=256\,$\mathrm{m}\mathrm{T}$$
(upper one).

[figure]capposition=bottom

Figure S8:  ZBP in the first lobe due to ABSs in device B. Line traces
extracted in the first lobe from Fig.4E showing $dI/dV$ as a function of $V$
from $B=74\,$\mathrm{m}\mathrm{T}$$ (lower one) to
$B=158\,$\mathrm{m}\mathrm{T}$$ (upper one); they are shifted by
$0.01\frac{e^{2}}{h}$ with respect to each other.

[figure]capposition=beside,capbesideposition=center,right

Figure S9:  Theory plot for device C. Numerical simulation of the $dI/dV$
versus $V$ and $\phi/\phi_{0}$ for Fig. 5A in a QD-S (modelled as a
superconducting Anderson model in the Hartree-Fock approximation), with the SC
in the destructive LP regime and the QD in a doublet GS (parameters can be
found in Table S1 of the SI).

[figure]capposition=bottom

Figure S10:  Dependence of ZBP on the backgate voltage and additional data for
device D. (A) $dI/dV$ as a function of $V$ and $V_{\textrm{bg}}$ at in absence
of magnetic field. Multiple ABSs are clearly visible. (B) Numerical simulation
of the $dI/dV$ versus $V$ and $\phi/\phi_{0}$ for Fig.6c in a QD-S (modelled
as a superconducting Anderson model in the Hartree-Fock approximation), with
the SC in the destructive LP regime and the QD in a doublet GS (parameters can
be found in Table S1 of the SI). (C) $dI/dV$ as a function of $V$ and
$V_{\textrm{bg}}$ at $B=108\,$\mathrm{m}\mathrm{T}$$, i.e., in the first lobe.
The presence of a ZBP depends on $V_{\textrm{bg}}$, proving that it has no
topological origin. (D) Line traces extracted from (C); the blue curve
corresponds to $V_{\textrm{bg}}=-1.25\,$\mathrm{V}$$ while the green one to
$V_{\textrm{bg}}=-1.34\,$\mathrm{V}$$.

[figure]capposition=beside,capbesideposition=center,right

Figure S11:  Zoom of Fig. 5C. $dI/dV$ as a function of $B$ and $V$, same as
Fig. 5C but with a finer resolution.

[figure]capposition=bottom

## References

  * (1) Vaitiekėnas, S., Krogstrup, P. & Marcus, C. M.  Anomalous metallic phase in tunable destructive superconductors.  _Phys. Rev. B_ 101, 060507 (2020).
  * (2) Schwiete, G. & Oreg, Y.  Fluctuation persistent current in small superconducting rings.  _Phys. Rev. B_ 82, 214514 (2010).
  * (3) Shah, N. & Lopatin, A.  Microscopic analysis of the superconducting quantum critical point: Finite-temperature crossovers in transport near a pair-breaking quantum phase transition.  _Phys. Rev. B_ 76, 094511 (2007).
  * (4) Larkin, A. & Varlamov, A.  _Theory of Fluctuations in Superconductors_ (Oxford University Press, 2005).
  * (5) Groff, R. P. & Parks, R. D.  Fluxoid quantization and field-induced depairing in a hollow superconducting microcylinder.  _Phys. Rev._ 176, 567–580 (1968).
  * (6) Bardeen, J.  Critical fields and currents in superconductors.  _Rev. Mod. Phys._ 34, 667–681 (1962).
  * (7) Abrikosov, A.  On the problem of the knight shift in superconductors.  _Soviet Physics JETP_ 12, 337–339 (1961).
  * (8) de Gennes, P. G.  _Superconductivity in Metals and Alloys_ (W. A. Benjamin, Inc., New York, 1966).
  * (9) Skalski, S., Betbeder-Matibet, O. & Weiss, P. R.  Properties of superconducting alloys containing paramagnetic impurities.  _Phys. Rev._ 136, A1500–A1518 (1964).
  * (10) Vecino, E., Martín-Rodero, A. & Yeyati, A. L.  Josephson current through a correlated quantum level: Andreev states and $\pi$ junction behavior.  _Phys. Rev. B_ 68, 035105 (2003).
  * (11) Martín-Rodero, A. & Yeyati, A. L.  The andreev states of a superconducting quantum dot: mean field versus exact numerical results.  _Journal of Physics: Condensed Matter_ 24, 385303 (2012).
  * (12) Lee, E. J. H. _et al._ Spin-resolved andreev levels and parity crossings in hybrid superconductor-semiconductor nanostructures.  _Nat. Nanotechnol._ 9, 79–84 (2014).
  * (13) Yu, L.  Bound state in superconductors with paramagnetic impurities.  _Acta. Phys. Sin._ 21, 75–9 (1965).
  * (14) Shiba, H.  Classical spins in superconductors.  _Prog. Theor. Phys._ 40, 435–451 (1968).
  * (15) Rusinov, A.  Superconductivity near a paramagnetic impurity.  _Sov. Phys. JETP_ 9, 85 (1969).
  * (16) Balatsky, A. V., Vekhter, I. & Zhu, J.-X.  Impurity-induced states in conventional and unconventional superconductors.  _Rev. Mod. Phys._ 78, 373–433 (2006).
  * (17) Kiršanskas, G., Goldstein, M., Flensberg, K., Glazman, L. I. & Paaske, J.  Yu-shiba-rusinov states in phase-biased superconductor–quantum dot–superconductor junctions.  _Phys. Rev. B_ 92, 235422 (2015).
  * (18) Žitko, R., Lim, J. S., López, R. & Aguado, R.  Shiba states and zero-bias anomalies in the hybrid normal-superconductor anderson model.  _Phys. Rev. B_ 91, 045441 (2015).
  * (19) Lutchyn, R. M., Sau, J. D. & Das Sarma, S.  Majorana Fermions and a Topological Phase Transition in Semiconductor-Superconductor Heterostructures.  _Phys. Rev. Lett._ 105, 077001 (2010).
  * (20) Oreg, Y., Refael, G. & von Oppen, F.  Helical Liquids and Majorana Bound States in Quantum Wires.  _Phys. Rev. Lett._ 105, 177002 (2010).
  * (21) Bauer, J., Oguri, A. & Hewson, A. C.  Spectral properties of locally correlated electrons in a bardeen–cooper–schrieffer superconductor.  _Journal of Physics: Condensed Matter_ 19, 486211 (2007).
  * (22) Meng, T., Florens, S. & Simon, P.  Self-consistent description of andreev bound states in josephson quantum dot devices.  _Phys. Rev. B_ 79, 224521 (2009).
  * (23) Žonda, M., Pokorný, V., Janiš, V. & Novotný, T.  Perturbation theory of a superconducting 0-$\pi$ impurity quantum phase transition.  _Scientific Reports_ 5, 8821 (2015).
  * (24) Žonda, M., Pokorný, V., Janiš, V. & Novotný, T.  Perturbation theory for an anderson quantum dot asymmetrically attached to two superconducting leads.  _Phys. Rev. B_ 93, 024523 (2016).
  * (25) Wentzell, N., Florens, S., Meng, T., Meden, V. & Andergassen, S.  Magnetoelectric spectroscopy of andreev bound states in josephson quantum dots.  _Phys. Rev. B_ 94, 085151 (2016).
  * (26) García Corral, A. _et al._ Magnetic-field-induced transition in a quantum dot coupled to a superconductor.  _Phys. Rev. Research_ 2, 012065 (2020).
  * (27) Yoshioka, T. & Ohashi, Y.  Numerical renormalization group studies on single impurity anderson model in superconductivity: A unified treatment of magnetic, nonmagnetic impurities, and resonance scattering.  _J. Phys. Soc. Jpn_ 69, 1812–1823 (2000).
  * (28) Lee, E. J. H. _et al._ Scaling of subgap excitations in a superconductor-semiconductor nanowire quantum dot.  _Phys. Rev. B_ 95, 180502 (2017).
  * (29) Chang, W. _et al._ Hard gap in epitaxial semiconductor–superconductor nanowires.  _Nat. Nanotechnol._ 10, 232–236 (2015).
  * (30) Krogstrup, P. _et al._ Epitaxy of semiconductor–superconductor nanowires.  _Nat. Mater._ 14, 400–406 (2015).
  * (31) Sarma, S. D., Sau, J. D. & Stanescu, T. D.  Splitting of the zero-bias conductance peak as smoking gun evidence for the existence of the majorana mode in a superconductor-semiconductor nanowire.  _Phys. Rev. B_ 86, 220506 (2012).
  * (32) Albrecht, S. M. _et al._ Exponential protection of zero modes in Majorana islands.  _Nature_ 531, 206–209 (2016).
  * (33) Vaitiekėnas, S. _et al._ Flux-induced topological superconductivity in full-shell nanowires.  _Science_ 367, eaav3392 (2020).
  * (34) Peñaranda, F., Aguado, R., San-Jose, P. & Prada, E.  Even-odd effect and majorana states in full-shell nanowires.  _Phys. Rev. Research_ 2, 023171 (2020).
  * (35) Žitko, R., Lim, J. S., López, R. & Aguado, R.  Shiba states and zero-bias anomalies in the hybrid normal-superconductor anderson model.  _Phys. Rev. B_ 91, 045441 (2015).