Evidence for spin-polarized bound states in
semiconductor-superconductor-ferromagnetic insulator islands

InAs wires with hexagonal cross section were grown to a length of $\sim 10~\mu$m and diameter of $\sim 120$ nm using molecular beam epitaxy Krogstrup et al. (2015). Partly-overlapping, two-faceted EuS (as grown thickness 8 nm) and Al (as grown thickness 6 nm) shells were grown in situ using electron beam evaporation, as shown in the main-text Fig. 1(a) Liu et al. (2020c, a). Devices were fabricated on a Si substrate with 200 nm SiOx capping. Coulomb islands were formed by wet-etching (MF-321 photoresist developer, 30 s, room temperature) and contacting the exposed InAs/EuS with Ti/Au (5/150 nm) after in situ Ar milling (15 W, 7 min). Devices were then coated with a HfO₂ (8 nm) dielectric layer, followed by the deposition of Ti/Au (5/150 nm) gates patterned using electron beam lithography. Additional details can be found in Ref. Vaitiekėnas et al. (2021).

The samples were cooled to base temperature of the dilution refrigerator at zero applied magnetic field. Following cooldown, an external magnetic field was applied $\mu_{0}H_{\parallel}=150$ mT along the wire axis then returned to zero. Unless otherwise noted, the measurements were carried out at zero applied field. Differential conductance, $G=$ d$I/$d$V$, was measured by sourcing voltage bias, $V$, through one of the outer leads, floating the opposite end of the wire, and draining the current, $I$, through the common lead between the two islands, as shown for the 400 nm island in the main-text Fig. 1(b) and for the 800 nm island in Fig. S2(a). The islands were tuned into the Coulomb blockade regime using negative voltages applied to the cutter gates, $V_{\rm C1}$ and $V_{\rm C2}$, and back gate, $V_{\rm BG}$ [labeled only in Fig. S2(a)]. Over a range much larger than the Coulomb peak spacing, $V_{\rm BG}$ also tuned the chemical potential of the island. The upper-gate voltage, $V_{\rm U}$, on the side coated with the Al shell, was predominantly used to tune the occupancy of the hybrid island, whereas the lower-gate voltage, $V_{\rm L}$, on the EuS shell side, was used to tune the charge carrier density in the semiconductor.

We consider a superconducting island hosting a subgap state that can be spin-slit (for example, due to exchange-coupling to the ferromagnetic insulator) by an energy $E_{\rm Z}$. For simplicity, we describe the continuum of states as a spin-degenerate quasiparticle state at energy $\Delta$, using the so-called zero bandwidth model Grove-Rasmussen et al. (2018). We take the charging energy of the island, $E_{\rm C}$, to be the largest and the coupling to the leads to be the smallest energy scales in the system, allowing us to treat the electron tunneling in a perturbative way.

The Hamiltonian of the system is given by

where the leads are described by

with the energy, $\xi_{\nu k\sigma}$, and the annihilation operator, $c_{\nu k\sigma}$, of an electron in lead $\nu\in\{L,R\}$ with momentum $k$ and spin $\sigma\in\{\uparrow,\downarrow\}$. We assume that each lead remains in internal equilibrium described by a Fermi-Dirac distribution with chemical potential $\mu_{\nu}$.

The superconducting island is described by

where $j\in\{0,\Delta\}$ is the state index, $\vartheta_{j\sigma}$ is the (spin-dependent) energy of subgap and continuum states, $\gamma_{j\sigma}$ is the Bogoliubov quasiparticle annihilation operator, and $E_{\rm el}$ is the electrostatic repulsion term that depends on the number of electrons in the island, $N=2N_{\rm C}+n_{\rm qp}$, where $N_{\rm C}$ is the number of Cooper pairs and $n_{\rm qp}$ is the total charge in the quasiparticle states. For the subgap state, we take $\vartheta_{0\downarrow}=\varepsilon$ and $\vartheta_{0\uparrow}=\varepsilon+E_{\rm Z}$. For simplicity, we take spin-degenerate continuum states at $\vartheta_{1\downarrow}=\vartheta_{1\uparrow}=\Delta$. Similar results are found for spin-split continuum states. The annihilation operator for an electron in the island, $d_{j\sigma}$, is related to $\gamma_{j\sigma}$ by

where $u_{j}$ and $v_{j}$ are the Bogoliubov coefficients, $\rho_{\uparrow/\downarrow}=\pm 1$ depending on the spin, $\phi$ is the superconducting phase operator such that $e^{i\phi}$ creates a Cooper pair on the island, and $\bar{\sigma}$ is the spin opposite to $\sigma$. The electrostatic repulsion term is given by

where $n_{\rm G}$ is the dimensionless gate-induced charge offset.

The tunneling Hamiltonian is described by

where $t$ denotes the tunneling amplitudes between the island and leads.

The transport properties through the superconducting island are calculated using the T-matrix formalism Bruus and Flensberg (2004). The transition rate between initial, $\left|i\right\rangle$, and final, $\left|f\right\rangle$, states can be computed by

where $W_{if}$ weights the rate through thermal distributions of initial and final states at energies $E_{i}$ and $E_{f}$, respectively, $\delta$ is the Dirac delta function and

which can be truncated at the desired order. Here the term linear in $H_{\rm T}$ describes the sequential tunneling; the higher order terms describe the cotunneling contributions, which become progressively more important when the tunneling coupling between the island and the leads increases.

The quantum state of the island is described by $\left|a\right\rangle$=$\left|N,N_{\rm C},\textbf{n}\right\rangle$, where n is a vector representing the occupation of the subgap and continuum excited states. The time derivative of the occupation probability of a given state can be written as

which describes the stationary probability of occupation, $P^{\rm{stat}}_{b}$, by imposing $\dot{P}^{\rm{stat}}_{b}=0$ and $\sum_{b}P^{\rm{stat}}_{b}=1$. The current through the device can be computed using the stationary distribution of probabilities and the rates. The tunneling current is given by

where $s=+1$ for rightwards sequential tunneling and $-1$ for the opposite direction, $\delta n(L\to R)$ is the net charge transferred from left to right in a cotunneling process, and the sum runs over all the possible rates connecting the island state $\left|b\right\rangle$ with any state $\left|a\right\rangle$.

The sequential rates are given by

where $n_{\rm F}$ is the Fermi-Dirac distribution function and $\mu_{\nu}$ is the chemical potential of the lead $\nu$. Here, we have used the wideband approximation, where the tunneling rates are energy independent and $\Gamma_{\nu}=2\pi\rho_{F}|t_{\nu}|^{2}$, with the lead density of states at the Fermi level $\rho_{F}$. Note that for these rates n and $\textbf{n}^{\prime}$ differ by the occupation of one state.

We consider the cotunneling rates transferring an electron from one lead to the other, which are the dominant ones in the limit $E_{\rm C}\gg\max\left(\Delta,\,v\equiv\mu_{\rm L}-\mu_{\rm R}\right)$, where $v$ is the source-drain voltage bias used in calculations. These processes, which conserve the charge on the island, can be expressed by

where $H_{\rm T}(\nu)$ describes the tunneling between $\nu$ lead and the island, $\bar{\nu}$ denotes a lead opposite to $\nu$, and $m_{1,2}$ are virtual intermediate states. To derive this expression we have imposed energy conservation, which leads to a function dependent on the energy of the tunneling electron from/to one of the leads, $\omega_{1}$. The cotunneling rate can be written as

This expression for the cotunneling rates is divergent due to the appearing sequential tunneling. To avoid the divergent behaviour, we regularize the divergences as explained in Ref. Koller et al. (2010). The resulting integral can be formally solved analytically Koch et al. (2004), which leads to a complicated expression involving special functions. In the limit of $T/E_{\rm C}\gtrsim 10^{-3}$, it turns out to be more computationally efficient to simplify it by expanding $n_{\rm F}$ into a sum of complex Matsubara-Ozaki frequencies Ozaki (2007)

where $\beta_{\alpha}$ and $r_{\alpha}$ are the approximated Matsubara frequencies and residues, respectively. Finally, Eq. (13) can be evaluated using the residue theorem yielding

This sum can be truncated at $\alpha\approx 100$ Matsubara frequencies for the parameters used in the calculations.

Calculated Coulomb spectra for a spin-degenerate ($E_{\rm Z}=0$) and weakly spin-split ($E_{\rm Z}<\Delta$) subgap states are shown in Figs. S3(a) and S3(g). The colored ticks mark the onset of inelastic cotunneling events that excite the system into a higher-energy state, resulting in a step-like increase in conductance. The transitions in the odd and even valleys are represented separately in panels (b) and (c) for the spin-degenerate, and panels (h) and (i) for the weakly spin-split cases. Some examples of the cotunneling transport mechanisms for the two cases are illustrated in panels (d)–(f) and (j)–(l), respectively. In general, there are four cotunneling lines in both spin-degenerate and weakly spin-split cases (except for the specific instances where $\varepsilon$ is fine tuned such that two transitions are degenerate in energy). In contrast, our experimental data show only two cotunneling steps (see, for example, Fig. 3 in the main text). A qualitative line-cut comparison between experimental data and spin-polarized as well as spin-degenerate cases is shown for the 400 nm island on wire 1 in Fig. S4. Aggregated line-cuts from several consecutive even and odd Coulomb valleys measured for the 400 nm island on wire 2 display two cotunneling steps, independent of lower-gate voltage and subgap-state energy, see Fig. S5.

The measured cotunneling features decrease in energy without splitting as an external magnetic field is applied, see Fig. S6(b). For comparison, we calculate conductance for a spin-polarized and two spin-degenerate cases—with and without magnetic hysteresis. As an input, we use approximated expressions for the field dependence of the superconducting gap and the subgap states. We describe the gap closing as $\Delta(H_{\parallel})=\Delta_{0}\left[1-\left(H_{\parallel}/H_{\rm C}\right)^{2}\right]$, with critical field $H_{\rm C}=60$ mT and the zero-field gap $\Delta_{0}=60~\mu$eV, see the black curves in Fig. S6(c), (e), and (g). For simplicity we assume that the subgap state depends linearly on the field, see the red and blue curves in Fig. S6(c), (e), and (g). To account for the asymmetry in the measured data, we include magnetic hysteresis of 10 mT (consistent with previous experiments Vaitiekėnas et al. (2021)) in the spin-polarized and one of the spin-degenerate cases. The magnetic field in the spin-degenerate case leads to the appearance of additional cotunneling features that are not observed experimentally, see Fig. S6(f) and (h). We therefore find that the measurements agree best with the spin-polarized spectrum shown in Fig. S6(d).

We note that the number of cotunneling steps increases in similar hybrid island devices without EuS shell (see the replotted data from Ref. Higginbotham et al. (2015) in Fig. S7) consistent with Fig. S3(a).

These observations together suggest that the investigated hybrid islands are in the strongly spin-split limit ($E_{Z}\gtrsim\Delta$), discussed in the main-text Fig. 2.

The relative height of the cotunneling steps in the even and odd Coulomb valleys depends on the relative phase, $\varphi$, between the Bogoliubov components of the subgap state, $u_{0}\equiv|u_{0}|\exp(i\varphi)$ and $v_{0}\equiv|v_{0}|$, as illustrated in Fig. S9 for $E_{\rm Z}\gg\Delta$ and $\varepsilon=0$. The effect can be explained by the interference between different cotunneling mechanisms that involve the same initial and final states. However, $\varphi$ cannot be determined unambiguously in our setup, as the strong spin splitting and gate-dependent $\varepsilon$ makes the global ground state unknown. This remains an open problem, relevant for distinguishing trivial and topological states.

In the experiment, the magnitude and, in some cases, the sign of $\varepsilon$ can be tuned using electrostatic gate electrodes (Fig. 4 in the main text). In our model, we account for this behavior by changing $\vartheta_{0\sigma}$ in Eq. (3). A change in the sign of $\varepsilon$ is equivalent to exchanging the Bogoliubov components of the corresponding subgap state. In Fig. S10 we show the calculated conductance in the strongly spin-split case for $\varepsilon$ greater than, equal to, and less than zero. For $\varepsilon=0$, the Coulomb blockade is 1$e$-periodic—the onsets of the lowest cotunneling steps in both odd and even valleys align at $v=\Delta$. For $\varepsilon$ away from zero, the spectrum is even-odd periodic with the inelastic onset at $v=\Delta-\varepsilon$ in the odd and $v=\Delta+\varepsilon$ in the even valleys. The sign of $\varepsilon$ determines the relative size of odd and even Coulomb diamonds.