# Interference and parity blockade in transport through a Majorana box

Maximilian Nitsch  Division of Solid State Physics and NanoLund, Lund
University, S-22100 Lund, Sweden  Rubén Seoane Souto  Division of Solid State
Physics and NanoLund, Lund University, S-22100 Lund, Sweden  Center for
Quantum Devices, Niels Bohr Institute, University of Copenhagen, DK-2100
Copenhagen, Denmark  Martin Leijnse  Division of Solid State Physics and
NanoLund, Lund University, S-22100 Lund, Sweden  Center for Quantum Devices,
Niels Bohr Institute, University of Copenhagen, DK-2100 Copenhagen, Denmark

###### Abstract

A Majorana box – two topological superconducting nanowires coupled via a
trivial superconductor – is a building block in devices aiming to demonstrate
nonabelian physics, as well as for topological quantum computer architectures.
We theoretically investigate charge transport through a Majorana box and show
that current can be blocked when two Majoranas couple to the same lead, fixing
their parity. In direct analogy to Pauli spin blockade in spin qubits, this
parity blockade can be used for fast and high-fidelity qubit initialization
and readout, as well as for current-based measurements of decoherence times.
Furthermore, we demonstrate that transport can distinguish between a clean
Majorana box and a disordered box with additional unwanted Majorana or Andreev
bound states.

_Introduction._ Topological $p$-wave superconductors host Majorana bound
states (MBSs) [1, 2, 3, 4, 5, 6] at edges and defects, which have nonlocal and
nonabelian properties. Semiconductor nanowires are one of the most promising
systems for creating and detecting MBSs, where a combination of spin-orbit
coupling, proximity-induced superconductivity, and external magnetic field can
lead to $p$-wave superconductivity [7, 8]. By now, many experiments have
observed zero-bias conductance peaks, consistent with MBSs at the nanowire
ends (see Refs. [9, 10, 11, 12, 13, 14] for a few examples, similar results
have been obtained also in other MBS platforms). However, nontopological
states provide an alternative explanation for most of the experimental
observations [15, 16, 17, 18, 19, 20, 21].

A measurement of the nonabelian properties of MBSs is still missing, but would
provide definite evidence of a topological superconducting phase, constituting
at the same time a first step towards topological quantum computing. One
promising path towards a demonstration of nonabelian physics uses repeated
measurements of MBS pairs to perform topologically protected qubit operations
[22], with a possibility to move towards a scalable quantum computer platform
[23, 24, 25]. A simple building block for these technologies is the Majorana
box qubit [26], where a qubit is encoded in four MBSs with overall parity
fixed by a large charging energy. Qubit readout can be done by charge sensing
of a quantum dot coupled to two MBSs [27, 28, 29, 30], or by measuring the
interference of cotunneling currents when the box is connected to external
leads [31]. Furthermore, coupling the Majorana box to four leads enables
measurements of the topological Kondo effect [32, 33, 34], and networks of
coupled Majorana boxes exhibit additional interesting transport physics [35,
36, 37, 38, 39].

Figure 1: (a) Sketch of Majorana box, where two topological superconducting
wires (blue), connected by a nontopological superconductor (dark blue), host
four end MBSs (red crosses, operators $\gamma_{rm}$). Disorder might lead to
additional unwanted MBSs (orange crosses). The box is tunnel coupled
(amplitudes $t_{rm}$) to two normal leads ($L$ and $R$), subject to a voltage
bias $V_{b}$. Magnetic fluxes $\Phi_{L},\,\Phi_{R}$ are threaded through the
loops associated with leads $L,\,R$ and cause relative phase differences
$\phi_{L},\,\phi_{R}$ between $t_{Lu}$ and $t_{Ld}$, $t_{Ru}$ and $t_{Rd}$. A
gate voltage $V_{g}$ controls the equilibrium number of electrons on the
Majorana box. (b) Current $I$ through the Majorana box as a function of
$V_{b}$ and $V_{g}$ with the remaining parameters specified in the text. (c)
$I$ at $V_{g}=0$ V, $V_{b}=10\,T$ (red cross in (a)) as a function of
$\phi_{L}$ at $\phi_{R}=0$.

In this work, we develop and employ a quantum master equation approach to
investigate charge transport through a Majorana box where the source and drain
contacts couple to two MBSs each, see Fig. 1(a). We show that the same
mechanism that allows quantum-dot-based parity readout [27, 28, 30] induces a
parity blockade in our transport setup, where the current is quenched and the
qubit is stuck in a well-defined state. This is in close analogy to the Pauli
spin blockade in double quantum dot spin qubits [40, 41]. Just like the Pauli
spin blockade, parity blockade can simplify various important qubit
experiments. Fast and high-fidelity qubit initialization can be achieved by
driving a current through the Majorana box which quickly gets stuck in the
blocked state. The same principle can be used for readout, by applying a bias
voltage such that an electron tunnels if the system is not in the blocking
state. Single-shot readout can then be accomplished by charge detection on the
box. Alternatively, measuring the current resulting from repeated operations
provides an averaged readout. We solve the quantum master equation
analytically for the clean box with four end MBSs and numerically for a
disordered Majorana box with additional unwanted MBSs or topologically trivial
Andreev bound states (ABSs). We show that the qubit coherence time can be read
off from the remnant steady-state current in the blocking regime. This
measurement requires neither fast manipulation, nor fast readout, only a DC
transport measurement. Finally, we explain how to distinguish the clean
Majorana box from the disordered system with additional MBSs or ABSs inside
the box.

_Model and transport theory._ We consider the Majorana box transport setup
sketched in Fig. 1(a). Two topological superconducting nanowires are connected
by a conventional (nontopological) superconductor and are tunnel coupled to
electrically biased normal source and drain contacts. The Hamiltonian is
$H=H_{MB}+H_{res}+H_{T}$. The Majorana box is described by ($e=\hbar=k_{B}=1$)

$\displaystyle
H_{MB}=\sum_{m=u,d}\frac{i}{2}\varepsilon_{m}\gamma_{Lm}\gamma_{Rm}+E_{C}(N-n_{g})^{2}+H_{MB}^{\mathrm{dis}},$
(1)

where $\gamma_{rm}$ are MBS operators and $\varepsilon_{m}$ is the overlap
between MBSs in the same wire (our results remain qualitatively the same in
the presence of additional overlaps between MBSs in different wires), $E_{C}$
is the charging energy, $N$ counts the number of electrons (including Cooper
pairs) on the Majorana box, and $n_{g}$ is the background charge controlled by
the gate voltage $V_{g}$, $n_{g}=\alpha_{g}V_{g}$ with gate lever arm
$\alpha_{g}$. $H_{MB}^{\mathrm{dis}}$ describes a number of additional
unwanted MBSs [orange crosses in Fig. 1(a)] induced by disorder, which may
overlap with each other and with the edge MBSs $\gamma_{rm}$, see specific
examples below. In this low-energy Hamiltonian we neglect the quasiparticle
states above the superconducting gap.

The lead Hamiltonian is $H_{res}=\sum_{r}H_{r}$, with
$H_{r}=\sum_{k}\xi_{rk}c_{rk}^{\dagger}c_{rk}$, where the $c_{rk}^{\dagger}$
create spinless electrons in lead $r=L,\,R$ with energies $\xi_{rk}$. We
assume the leads to remain in thermal equilibrium at temperature $T$ and
chemical potential $\mu_{L,R}=\pm V_{b}/2$. The tunneling between leads and
Majorana box is described by

$\displaystyle
H_{T}=\sum_{rmk}\gamma_{rm}\left(t_{rm}c_{rk}-t_{rm}^{*}c_{rk}^{\dagger}\right)+H_{T}^{\mathrm{dis}},$
(2)

with tunnel amplitudes $t_{rm}$ which we take to be energy independent
(wideband limit). We include magnetic fluxes $\Phi_{L},\,\Phi_{R}$ threaded
through the loops formed by leads $L,\,R$ and the end MBSs [Fig. 1(a)] by
adding a phase $\phi_{r}=2\pi\Phi_{r}/\Phi_{0}$ to the upper tunnel amplitude
of the left and right leads, $t_{ru}=|t_{ru}|e^{i\phi_{r}}$,
$t_{rd}=|t_{rd}|$, where $\Phi_{0}$ is the flux quantum. The amplitude for a
tunneling-induced transition between two many-body eigenstates $a$ and $b$ of
the Majorana box is related to the tunnel matrix element
$T_{r}^{ab}=\sum_{m=u,d}\langle a|\gamma_{rm}|b\rangle$. The typical time-
scale of electron tunneling is then given by the tunnel rates
$\Gamma_{r}^{ab}=2\pi\nu_{r}|T_{r}^{ab}|^{2}$, where we take the density of
states $\nu_{r}$ to be energy-independent within the bandwidth chosen as
$D=100T$. Unless stated otherwise, we will throughout the paper consider all
tunnel amplitudes and densities of states to be equal, $t_{rm}=t$ and
$\nu_{r}=\nu$, and define $\Gamma=2\pi\nu|t|^{2}$. $H_{T}^{\mathrm{dis}}$
describes the tunnel coupling of the disorder-induced MBSs to the leads.

Let us comment on two model assumptions which will be important for the
results and which place some constraints on an experimental realization.
First, the way $\gamma_{ru}$ and $\gamma_{rd}$ couple to the same lead channel
in Eq. (2) is only strictly correct for an effectively 1D lead, but is a good
approximation whenever tunneling from $\gamma_{ru}$ and $\gamma_{rd}$ occur
into points of the lead separated by less than the Fermi wavelength. Second,
considering spinless lead electrons is valid either when the magnetic field
needed to induce the topological superconducting phase has fully spin-
polarized the lead electrons around the Fermi level, or when the spin
directions associated with allowed tunneling into $\gamma_{ru}$ and
$\gamma_{rd}$ are aligned [42] (which is the case for two identical wires).

_Quantum master equations._ We focus on the regime of weak tunneling,
$\Gamma\ll T$, but strong electron-electron interaction $E_{C}$. Then it is
appropriate to use a quantum master equation for the reduced density matrix
$\rho$ of the Majorana box:

$\partial_{t}\rho=-i[H_{MB},\,\rho]+W\rho.$ (3)

The quantum master equation consists of a unitary time evolution determined by
$H_{MB}$ and a dissipative part introduced by the attached leads. We
diagonalize the Majorana box Hamiltonian in Eq. (1) to obtain the many-body
eigenstates $|a_{i}\rangle$ and solve the master equation for the stationary
state reduced density matrix $\rho_{a_{1}a_{2}}$ and current $I$, where
tunneling is treated in leading-order perturbation theory. We emphasize that,
because of the near-degenerate ground state, it is important to solve for the
full nondiagonal density matrix. All results presented below are obtained
within a 1st order von Neumann quantum master equation [43] (equivalent to
real-time diagrammatics in 1st order [44, 45, 46]). We have cross-checked that
other approximations [47, 48] provide similar results, details are given in
the supplementary information (SI) [49].

Because of the large $E_{C}$ we consider only two charge states, arbitrarily
denoted by $N=0$ and $N=1$, corresponding to the total parity of the MBSs
being even or odd. We tune $V_{g}$ such that the two parity sectors are almost
degenerate. The density matrix is diagonal in total parity, and a term in the
master equation describing an electron tunneling onto or out of the Majorana
box connects the two parity sectors.

_Parity blockade._ We first consider the clean box with
$H^{\mathrm{dis}}_{MB}=0$ and $\varepsilon_{m}=0$. The current $I$ as a
function of $V_{b}$ and $V_{g}$ for $\phi_{L}=\phi_{R}=0$ [Fig. 1(b)] shows
the Coulomb blockade pattern characteristic of transport through quantum dots
[50]. The current is finite for $V_{b},V_{g}$ such that there are available
electrons in one contact that can tunnel into the Majorana box ($N\rightarrow
N+1$) and available empty states in the other contact that can accept
electrons tunneling out of the box ($N+1\rightarrow N$). Otherwise current is
suppressed by charging effects (Coulomb blockade).

For the remainder of the paper, we fix the voltages within the conducting
regime [at the point marked by the red cross in Fig. 1(b). We now vary
$\phi_{L}$, see Fig. 1(c), and find $I(\phi_{L})\propto\cos^{2}\phi_{L}$,
meaning that the current is blocked for
$\phi_{L}=(2n+1)\frac{\pi}{2},\,n\in\mathrm{Z}$. To understand the blockade,
we construct fermion operators using the two left and the two right MBSs:
$f_{L}=(\gamma_{Lu}+i\gamma_{Ld})/2,\,f_{R}=(\gamma_{Rd}+i\,\gamma_{Ru})/2$.
The eigenstates of the number operators $\hat{n}_{r}=f_{r}^{\dagger}f_{r}$,
$r=L,R$, $|n_{L}n_{R}\rangle$, are also eigenstates of $H_{MB}$ when
$\varepsilon_{m}=0$ and $H_{MB}^{\mathrm{dis}}=0$. Note that the even
eigenstates, $|0_{L}0_{R}\rangle$ and $|1_{L}1_{R}\rangle$, are degenerate,
and so are the odd eigenstates, $|0_{L}1_{R}\rangle$ and $|1_{L}0_{R}\rangle$.
With this choice of basis and in this simple limit, the density matrix is
diagonal. At the chosen voltages, electrons tunnel into the Majorana box from
lead $L$ and out to lead $R$. For a current to flow, the state of the Majorana
box must change according to
$|0_{L}0_{R}\rangle\rightarrow|1_{L}0_{R}\rangle\rightarrow|1_{L}1_{R}\rangle\rightarrow|0_{L}1_{R}\rangle\rightarrow|0_{L}0_{R}\rangle\rightarrow\ldots$
(electron tunnels in from the left, out to the right, in from the left, out to
the right, $\ldots$). Note that because the number states $n_{L},R$ are not
charge eigenstates, it is possible to, for example, switch from $n_{L}=1$ to
$n_{L}=0$ by an electron _entering_ the box from contact $L$. Taking the
tunneling term that adds an electron from the left lead in Eq. (2), and
writing it in terms of the left/right fermion operator, we obtain:

$\displaystyle H_{T,L}\rightarrow
t\sum_{k}[c_{k}(e^{i\phi_{L}}+i)f_{L}^{\dagger}+c_{k}(e^{i\phi_{L}}-i)f_{L}].$
(4)

For $\phi_{L}=\pi/2$ the second term in Eq. (4) vanishes, which results in the
transition $|1_{L}1_{R}\rangle\rightarrow|0_{L}1_{R}\rangle$ being suppressed.
Therefore, the system becomes trapped in the blocking state
$|1_{L}1_{R}\rangle$ and no current can flow. For $\phi_{L}=3\pi/2$ the
blocking state is instead $|0_{L}0_{R}\rangle$. Reversing $V_{b}$ or changing
$\phi_{R}$ at the right lead causes blocking instead in an odd state
($|0_{L}1_{R}\rangle$ or $|1_{L}0_{R}\rangle$). We note that, in direct
analogy with the Pauli spin blockade [40, 51], this parity blockade can be
used for fast and high-fidelity initialization of a Majorana box qubit in any
of the blocking states, as well as for readout in the corresponding basis.

We now move on to investigate how the blockade is lifted and how to read off
qubit lifetimes from the remnant current
$I_{rem}=\min_{\phi_{L}}[I(\phi_{L})]$ in the blocked regime. First, we note
that the blockade is lifted for asymmetric tunnel couplings to the upper/lower
MBSs. We will quantify this more explicitly below and for now assume
$t_{Lu}=t_{Ld}$. For now we keep the assumption $H_{MB}^{\mathrm{dis}}=0$ but
take $\varepsilon_{m}\neq 0$. Then the eigenstates are $|n_{u}n_{d}\rangle$
rather than $\ket{n_{L}n_{R}}$, associated with the up/down fermions with
operators $f_{m}=(\gamma_{Lm}+i\gamma_{Rm})/2$ for $m=u,d$. The eigenenergies
$E_{n_{u}n_{d}}$ within each parity sector are split by the MBS overlap,
$2\,\Delta_{e}=E_{11}-E_{00}=\varepsilon_{u}+\varepsilon_{d}$ and
$2\,\Delta_{o}=E_{01}-E_{10}=\varepsilon_{u}-\varepsilon_{d}$. Moreover, the
coupling to the leads introduces a Lamb-shift given by

$H_{LS}=\Gamma
I_{P}\begin{pmatrix}\sigma_{x}(\sin\phi_{L}+\sin\phi_{R})&0_{2}\\\
0_{2}&-\sigma_{x}(\sin\phi_{L}-\sin\phi_{R})\end{pmatrix}$ (5)

proportional to the principle value integrals $I_{P}$, see [52, 53, 49].

We can write the master equation in terms of the probability $p_{e/o}$ to be
in the even/odd sector, and a pseudospin $\vec{s}_{e/o}$ that describes the
density matrix within each sector, where we choose the $z$-axis to be along
$|n_{L}n_{R}\rangle$. In the SI [49], we derive Bloch-like equations for the
pseudospin and show that the current is given by

$I=2e\,\Gamma(p_{e}+\sin\phi_{L}\,s_{e}^{z}).$ (6)

Without MBS overlaps, $s_{e,o}^{z}$ are decoupled from $s_{e,o}^{x,y}$. Finite
overlaps correspond to a magnetic field of strength $\Delta_{e,o}$ along the
$x$-direction. In leading order perturbation theory, $\Delta_{e,o}$ induces an
additional loss term of magnitude $\Delta^{2}/\Gamma^{2}(1+I_{P}^{2})$ in the
master equation for $\partial_{t}s_{e}^{z}$ at $\phi_{L}=\pi/2$. The blocking
state corresponds to $p_{e}=1-\frac{1}{2}\Delta^{2}/\Gamma^{2}(1+I_{P}^{2})$,
$s_{e}^{z}=-1+\Delta^{2}/\Gamma^{2}(1+I_{P}^{2})$, resulting in a current
$I_{\mathrm{rem}}=e\Delta_{e}^{2}/\Gamma(1+I_{P}^{2})$. This result can be
generalized to any mechanism that allows parity to escape from the left
Majorana pair ($n_{L}=0\rightarrow n_{L}=1$) without changing the total charge
on the Majorana box. If the parity escape rate is $\tilde{\Delta}/\hbar$, the
resulting remnant current is
$I_{\mathrm{rem}}=e\tilde{\Delta}^{2}/\Gamma(1+I_{P}^{2})$ in the blocking
regime. Interestingly, a larger tunnel coupling to the leads enhances the
lifetime of the blocking state and suppresses current. Thus, even though a
measurement of the remnant current directly gives the inverse lifetime of the
blocking state $1/\tau=I_{rem}/e$, this is not the same as that for the
isolated Majorana box. The Lamb-shift is not experimentally accessible and
experiments can only extract the decay rate of the coupled Majorana box qubit
$\tilde{\Delta}_{\mathrm{coup}}\equiv\tilde{\Delta}/(1+I_{P}^{2})^{1/2}$ which
is always smaller but of the same order of magnitude as the decay rate
$\tilde{\Delta}$ of the isolated Majorana box. To measure
$\tilde{\Delta}_{\mathrm{coup}}$ one should first extract $\Gamma$ from the
current in the non-blocked regime, see Fig. 1(c), and then measure
$\tilde{\Delta}_{\mathrm{coup}}^{2}/\Gamma$ from the current in the blocked
regime.

Figure 2: Inset: Lifting of the blockade at $\phi_{L}=\frac{\pi}{2}$ for
$\varepsilon_{Lu}=2\cdot 10^{-4}\,T$ and $\delta_{t}=10^{-2}$. The main plot
shows $I_{rem}$ as a function of $\Gamma$ for
$\delta_{t}=\delta_{\phi}=10^{-2}$. After a sharp increase linear in $\Gamma$
until $\Gamma\approx\tilde{\Delta}_{\mathrm{coup}}$, $I_{rem}$ decreases as
$e\,\tilde{\Delta}_{\mathrm{coup}}^{2}/\Gamma$ (black dotted line) before it
increases again as $2e\,\delta^{2}\Gamma$ (black dashed line).

We illustrate this with a specific model for a disordered device containing
unwanted MBSs. We assume that these MBSs are uncoupled to the leads,
$H_{T}^{\mathrm{dis}}=0$, and

$H_{MB}^{\mathrm{dis}}=\frac{i}{2}\sum_{m=u,d}\sum_{r=L,R}\tilde{\varepsilon}_{rm}\gamma_{rm}\tilde{\gamma}_{rm}+\Omega_{m}\tilde{\gamma}_{Lm}\tilde{\gamma}_{Rm},$
(7)

where the $\tilde{\gamma}$s are four additional disorder-induced MBSs with
couplings $\Omega_{m}$ between each other and couplings
$\tilde{\varepsilon}_{rm}$ to the end MBSs. For
$\tilde{\varepsilon}_{Lu}\gg\tilde{\varepsilon}_{Ld},\,\Omega_{m}$, the
relevant parity escape rate is
$\tilde{\Delta}\approx(\tilde{\varepsilon}_{Lu}\pm\tilde{\varepsilon}_{Ld})/2\approx\tilde{\varepsilon}_{Lu}/2$.
Figure 2 shows the current as a function of $\Gamma$. For
$\Gamma\ll\tilde{\Delta}$ the current is proportional to $\Gamma$, it peaks at
$\Gamma\approx\tilde{\Delta}$ and then decays with larger $\Gamma$ as
$\tilde{\Delta}_{\mathrm{coup}}^{2}/\Gamma$ (dotted black line in Fig. 2). We
also introduce deviations from the ideal blocked situation
$\absolutevalue{t_{Ld}/t_{Lu}}=1-\delta_{t},\,\phi_{L}=\frac{\pi}{2}+\delta_{\phi},\,\delta=\sqrt{\delta_{t}^{2}+\delta_{\phi}^{2}}$.
These deviations lead to a contribution to the current that is linear in
$\Gamma$ (dashed black line) which dominates $I_{rem}$ for
$\Gamma>\tilde{\Delta}_{\mathrm{coup}}/2\delta$ [49]. In particular, the
different scaling in $\Gamma$ makes it possible to distinguish experimentally
between a remnant current caused by an escape rate ($I\propto 1/\Gamma$)
compared to one due to finite $\delta$ ($I\propto\Gamma$).

_Distinguishing clean from disordered box._ Now we move on to showing that the
phase dependence of the Lamb-shift offers the possibility to distinguish
between the clean box with only four MBSs in total, and the disordered box
with additional MBSs or ABSs, Fig. 3(a). We model each ABS as two closely
spaced MBSs which both couple to the leads by

$H_{T}^{\mathrm{dis}}=\sum_{rmk}\tilde{\gamma}_{rm}\left(t_{rm}c_{rk}-t_{rm}^{*}c_{rk}^{\dagger}\right),$
(8)

but with no overlaps with the MBSs on the other side of the box,
$\Omega_{m},\epsilon_{U}=0$. For each of these three cases, we block the
current from the left lead with $\phi_{L}=\frac{\pi}{2}$ and investigate the
dependence of $I_{rem}$ on $\phi_{R}$. Figure 3(b) shows that the result is
qualitatively different for the clean box (blue lines) compared with the
disordered box (red lines) and ABS box (green lines), and this difference is
robust to various parameter choices (different line styles).

Figure 3: (a) Escape rate mechanisms leading to $I_{rem}$ for the clean box
(blue), the disordered box (red) and the ABS box (green). (b) $I_{rem}$ as a
function of $\phi_{R}$ at $\phi_{L}=\frac{\pi}{2}$ and $\Gamma=10^{-2}\,T$
renormalized by $I_{rem}(\phi_{R}=0)$ for all three box models in line colours
that match colours in (a). The overlaps are chosen as
$\varepsilon_{U},\varepsilon_{Lu}=10^{-5}\,T$ For each model we plot $I_{rem}$
for the perfectly fine-tuned setting (solid lines), deviation in the fine-
tuned tunneling $\delta_{t}=10^{-4}$ (dotted lines), mismatch in tunneling to
left/right lead $\Gamma_{R}=\Gamma_{L}/2$ (dash-dotted lines) and a very large
overlap $\varepsilon_{U},\varepsilon_{Lu}=10^{-3}\,T$ (dashed lines).

In the case of the clean box the Lamb-shifts introduced by both leads either
add up ($\phi_{R}=\frac{\pi}{2}$) or subtract ($\phi_{R}=\frac{3\pi}{2}$), see
Eq. (5). This leads to a decrease respectively increase in $I_{rem}$. This
qualitative dependence is stable under all investigated parameter settings.
For a large overlap, there appears a peak at $\phi_{R}=\frac{\pi}{2}$, which
corresponds to additional blockade at the right lead interfering with the
blockade at the left lead.

For the disordered box with additional MBSs or ABSs, the line-shape is
qualitatively different, with only very narrow dips (width
$\propto\tilde{\Delta}_{\mathrm{coup}}/\Gamma$) around
$\phi_{R}=\frac{\pi}{2},\frac{3\pi}{2}$. They mark a transition from a
blockade at the left lead to a blockade at the right lead with smaller
overlap/escape rate.

The qualitative difference of $I_{rem}(\phi_{R})$ is due to the different
escape mechanism from the parity blocked state to the next pair of MBSs, see
Fig. 3(a). In the clean box this nearest pair of MBSs is connected to the
right lead. Accordingly, the adjustment of the Lamb-shift via $\phi_{R}$
affects the remnant current. But for the disordered and ABS boxes the nearest
pair is located in the inner part/at the left end of the setup. There is no
connection to the right lead and accordingly no dependence on the phase
$\phi_{R}$.

_Conclusions._ We have used a quantum master equation approach to investigate
transport through a Majorana box coupled to normal leads. There is a blocking
regime, where the Majorana box becomes trapped in a well-defined state and the
current is suppressed. In analogy with the Pauli spin blockade, this parity
blockade can be used for qubit initialization and readout, as well as for
measuring qubit coherence times from DC transport. We believe that this can
become a key enabling technique for a first generation of MBS qubit
experiments, where single shot readout might be challenging due to limited
control or short qubit coherence times. Furthermore, the proposed setup makes
it possible to experimentally distinguish between a clean Majorana box and a
box with additional disorder-induced MBSs or ABSs. In our model, the parity
lifetime is limited by MBS overlaps, but we expect that if quasiparticle
poisoning is the dominant relaxation mechanism, the proposed measurement of
the remnant current will instead reveal the poisoning time.

_Acknowledgments._ We acknowledge stimulating discussions with Karsten
Flensberg, Michele Burello and Jens Schulenborg and funding from Nanolund, the
Swedish Research Council (VR) and the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme under
Grant Agreement No. 856526.

## References

  * Kitaev [2001] A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Alicea [2012] J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Leijnse and Flensberg [2012] M. Leijnse and K. Flensberg,  27, 124003 (2012).
  * Aguado [2017] R. Aguado, La Rivista del Nuovo Cimento 40, 523 (2017).
  * Beenakker [2020] C. W. J. Beenakker, SciPost Phys. Lect. Notes , 15 (2020).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
  * Deng  _et al._ [2012] M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Nano Letters 12, 6414 (2012), pMID: 23181691, https://doi.org/10.1021/nl303758w .
  * Finck  _et al._ [2013] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Phys. Rev. Lett. 110, 126406 (2013).
  * Deng  _et al._ [2016] M. T. Deng, S. Vaitiekėnas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus, Science 354, 1557 (2016), https://www.science.org/doi/pdf/10.1126/science.aaf3961 .
  * Nichele  _et al._ [2017] F. Nichele, A. C. C. Drachmann, A. M. Whiticar, E. C. T. O’Farrell, H. J. Suominen, A. Fornieri, T. Wang, G. C. Gardner, C. Thomas, A. T. Hatke, P. Krogstrup, M. J. Manfra, K. Flensberg, and C. M. Marcus, Phys. Rev. Lett. 119, 136803 (2017).
  * Lutchyn  _et al._ [2018] R. M. Lutchyn, E. P. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nature Reviews Materials 3, 52 (2018).
  * Prada  _et al._ [2012] E. Prada, P. San-Jose, and R. Aguado, Phys. Rev. B 86, 180503 (2012).
  * Kells  _et al._ [2012] G. Kells, D. Meidan, and P. W. Brouwer, Phys. Rev. B 86, 100503 (2012).
  * Moore  _et al._ [2018] C. Moore, C. Zeng, T. D. Stanescu, and S. Tewari, Phys. Rev. B 98, 155314 (2018).
  * Awoga  _et al._ [2019] O. A. Awoga, J. Cayao, and A. M. Black-Schaffer, Phys. Rev. Lett. 123, 117001 (2019).
  * Vuik  _et al._ [2019] A. Vuik, B. Nijholt, A. R. Akhmerov, and M. Wimmer, SciPost Phys. 7, 61 (2019).
  * Pan and Das Sarma [2020] H. Pan and S. Das Sarma, Phys. Rev. Research 2, 013377 (2020).
  * Hess  _et al._ [2021] R. Hess, H. F. Legg, D. Loss, and J. Klinovaja, Phys. Rev. B 104, 075405 (2021).
  * Bonderson  _et al._ [2008] P. Bonderson, M. Freedman, and C. Nayak, Phys. Rev. Lett. 101, 010501 (2008).
  * Vijay  _et al._ [2015] S. Vijay, T. H. Hsieh, and L. Fu, Phys. Rev. X 5, 041038 (2015).
  * Plugge  _et al._ [2016] S. Plugge, L. A. Landau, E. Sela, A. Altland, K. Flensberg, and R. Egger, Phys. Rev. B 94, 174514 (2016).
  * Karzig  _et al._ [2017] T. Karzig, C. Knapp, R. M. Lutchyn, P. Bonderson, M. B. Hastings, C. Nayak, J. Alicea, K. Flensberg, S. Plugge, Y. Oreg, C. M. Marcus, and M. H. Freedman, Phys. Rev. B 95, 235305 (2017).
  * Plugge  _et al._ [2017] S. Plugge, A. Rasmussen, R. Egger, and K. Flensberg, New Journal of Physics 19, 012001 (2017).
  * Munk  _et al._ [2020] M. I. K. Munk, J. Schulenborg, R. Egger, and K. Flensberg, Phys. Rev. Research 2, 033254 (2020).
  * Steiner and von Oppen [2020] J. F. Steiner and F. von Oppen, Phys. Rev. Research 2, 033255 (2020).
  * Smith  _et al._ [2020] T. B. Smith, M. C. Cassidy, D. J. Reilly, S. D. Bartlett, and A. L. Grimsmo, PRX Quantum 1, 020313 (2020).
  * Schulenborg  _et al._ [2021] J. Schulenborg, M. Burrello, M. Leijnse, and K. Flensberg, Phys. Rev. B 103, 245407 (2021).
  * Fu [2010] L. Fu, Phys. Rev. Lett. 104, 056402 (2010).
  * Béri and Cooper [2012] B. Béri and N. R. Cooper, Phys. Rev. Lett. 109, 156803 (2012).
  * Galpin  _et al._ [2014] M. R. Galpin, A. K. Mitchell, J. Temaismithi, D. E. Logan, B. Béri, and N. R. Cooper, Phys. Rev. B 89, 045143 (2014).
  * Buccheri and Egger [2020] F. Buccheri and R. Egger, in  _Strongly Coupled Field Theories for Condensed Matter and Quantum Information Theory_ , edited by A. Ferraz, K. S. Gupta, G. W. Semenoff, and P. Sodano (Springer International Publishing, Cham, 2020) pp. 131–153.
  * Béri [2013] B. Béri, Phys. Rev. Lett. 110, 216803 (2013).
  * Altland and Egger [2013] A. Altland and R. Egger, Phys. Rev. Lett. 110, 196401 (2013).
  * Herviou  _et al._ [2016] L. Herviou, K. Le Hur, and C. Mora, Phys. Rev. B 94, 235102 (2016).
  * Michaeli  _et al._ [2017] K. Michaeli, L. A. Landau, E. Sela, and L. Fu, Phys. Rev. B 96, 205403 (2017).
  * Väyrynen  _et al._ [2020] J. I. Väyrynen, A. E. Feiguin, and R. M. Lutchyn, Phys. Rev. Research 2, 043228 (2020).
  * Ono  _et al._ [2002] K. Ono, D. G. Austing, Y. Tokura, and S. Tarucha, Science 297, 1313 (2002), https://www.science.org/doi/pdf/10.1126/science.1070958 .
  * Hanson  _et al._ [2007] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Kjaergaard  _et al._ [2012] M. Kjaergaard, K. Wölms, and K. Flensberg, Phys. Rev. B 85, 020503 (2012).
  * Kiršanskas  _et al._ [2017] G. Kiršanskas, J. N. Pedersen, O. Karlström, M. Leijnse, and A. Wacker, Computer Physics Communications 221, 317 (2017).
  * König  _et al._ [1997] J. König, H. Schoeller, and G. Schön, Phys. Rev. Lett. 78, 4482 (1997).
  * Leijnse and Wegewijs [2008] M. Leijnse and M. R. Wegewijs, Phys. Rev. B 78, 235424 (2008).
  * Schoeller [2009] H. Schoeller, The European Physical Journal Special Topics 168, 179 (2009).
  * Nathan and Rudner [2020] F. Nathan and M. S. Rudner, Phys. Rev. B 102, 115109 (2020).
  * Kiršanskas  _et al._ [2018] G. Kiršanskas, M. Franckié, and A. Wacker, Phys. Rev. B 97, 035432 (2018).
  * [49] See Supplementary Information in (URL included by the Editor).
  * Kouwenhoven  _et al._ [2001] L. P. Kouwenhoven, D. G. Austing, and S. Tarucha, Reports on Progress in Physics 64, 701 (2001).
  * Danon and Nazarov [2009] J. Danon and Y. V. Nazarov, Phys. Rev. B 80, 041301 (2009).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Ptaszyński and Esposito [2019] K. Ptaszyński and M. Esposito, Phys. Rev. Lett. 122, 150603 (2019).
  * Evans [1977] D. E. Evans, Communications in Mathematical Physics 54, 293 (1977).
  * Evans and Hanche-Olsen [1979] D. E. Evans and H. Hanche-Olsen, Journal of Functional Analysis 32, 207 (1979).
  * Manzano [2020] D. Manzano, AIP Advances 10, 025106 (2020).

Supplementary information for interference and parity blockade in transport
through a Majorana box

In this Supplementary information (SI) we analytically derive the conditions
for parity blockade for the clean box in a first order quantum master equation
approach. We start by rewriting the quantum master equation into Lindblad
form. Afterwards, we explain how the parity blockade is established without
overlaps between the MBSs, $\varepsilon_{u}=\varepsilon_{d}=0,$ and with
perfectly fine tuned tunnel couplings. Finally, we discuss how finite overlaps
and deviation of the tunnel amplitudes lift this blockade.

##  .1 Transforming the first order von Neumann approach into Lindblad form

The starting point of the derivation is the formulation of electron transport
through a Majorana box in the 1st order von Neumann approach (1vN), see [43].
It describes the time evolution of the quantum dot density matrix $\rho$

$\displaystyle\partial_{t}\rho_{bb^{\prime}}=$
$\displaystyle-i(E_{b}-E_{b^{\prime}})\rho_{bb^{\prime}}$ (9)
$\displaystyle-i\sum_{b^{\prime\prime},r}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}\Gamma_{b^{\prime\prime}a,ab^{\prime}}^{r}I_{ba}^{r-}-\sum_{c}\Gamma_{b^{\prime\prime}c,cb^{\prime}}^{r}I_{cb}^{r+*}\Big{]}$
$\displaystyle-i\sum_{b^{\prime\prime},r}\rho_{b^{\prime\prime}b^{\prime}}\Big{[}\sum_{c}\Gamma_{bc,cb^{\prime\prime}}^{r}I_{cb^{\prime}}^{r+}-\sum_{a}\Gamma_{ba,ab^{\prime\prime}}^{r}I_{b^{\prime}a}^{r-*}\Big{]}$
$\displaystyle-i\sum_{aa^{\prime},r}\rho_{aa^{\prime}}\Gamma_{ba,a^{\prime}b^{\prime}}^{r}\Big{[}I_{b^{\prime}a}^{r+*}-I_{ba^{\prime}}^{r+}\Big{]}$
$\displaystyle-i\sum_{cc^{\prime},r}\rho_{cc^{\prime}}\Gamma_{bc,c^{\prime}b^{\prime}}^{r}\Big{[}I_{c^{\prime}b}^{r-*}-I_{cb^{\prime}}^{r-}\Big{]},$

where the first term describes the unitary time evolution of the quantum dot
and the remaining ones describe the dissipative dynamics induced by the leads
labeled by $r=L,R$. Depending on $N_{b}$ the indices $a,c$ sum over states
with lower or higher total charge:

$\displaystyle N_{a}=N_{b}-1$ $\displaystyle N_{c}=N_{b}+1$ (10)

The tunneling between leads and the Majorana box is described by the tunneling
rate matrix $\Gamma$

$\Gamma_{ba,a^{\prime}b^{\prime}}^{r}=2\pi\nu_{F}\,T_{ba}^{r}\,T_{a^{\prime}b^{\prime}}^{r}.$
(11)

$I^{r\pm}$ contains the lead contribution

$I_{cb}^{r\pm}\equiv\frac{1}{2\pi}\mathcal{P}\int_{-D_{r}}^{D_{r}}\frac{f\left(\pm\frac{E-\mu_{r}}{T_{r}}\right)}{E-E_{cb}}dE-\frac{i}{2}f(\pm
x_{cb}^{r})\theta(D_{r}-\absolutevalue{E_{cb}}),$ (12)

with the lead potential $\mu_{r}$, temperature $T_{r}$, bandwidth $D_{r}$, and
the eigenenergies of the system $E_{b}$. We also define

$\displaystyle E_{cb}\equiv E_{c}-E_{b},$ $\displaystyle
x_{cb}^{r}\equiv\frac{E_{cb}-\mu_{r}}{T_{r}}.$ (13)

Furthermore, we introduce the Fermi function $f(x)$ and the Heaviside step
function $\theta(x)$

$\displaystyle f(x)=[\exp(x)+1]^{-1},$
$\displaystyle\theta(x)=\begin{cases}0,\,x<0\\\ 1,\,x\geq 0\end{cases},$ (14)

as well as the principle value integral $\mathcal{P}\int_{-D_{r}}^{D_{r}}dE$.
Throughout the SI, we focus on the red cross in Fig. 1(b), were gate and bias
voltage are chosen such that current can flu in principle. Furthermore, we
assume the system and the leads to be in the limit

$E_{cb}\ll T_{r}\ll\absolutevalue{\mu_{r}}\ll D_{r},$ (15)

with the symmetries

$\displaystyle T_{L}=T_{R}=T,$ $\displaystyle-\mu_{R}=\mu_{L}=\mu>0,$
$\displaystyle D_{L}=D_{R}=D.$ (16)

Due to the chosen limit and symmetries we neglect the eigenenergies compared
to temperature, chemical potential and bandwidth by approximating

$\frac{+E_{cb}-\mu}{T}\approx-\frac{\mu}{T}$ (17)

within the principle value integrals, i.e. we drop the indices

$\displaystyle I^{r\pm}\approx$
$\displaystyle\frac{1}{2\pi}\mathcal{P}\int_{-D}^{D}\frac{f\left(\pm\frac{E-\mu_{r}}{T_{r}}\right)}{E}dE-\frac{i}{2}f(\pm
x^{r}),$ (18) $\displaystyle=$
$\displaystyle\frac{1}{2}\,I_{P}-\frac{i}{2}f(\pm x^{r}),$ $\displaystyle
I_{P}\equiv$
$\displaystyle\frac{1}{\pi}\mathcal{P}\int_{-D}^{D}\frac{f\left(\frac{-E+\mu}{T}\right)}{E}dE,\hskip
28.45274ptx^{r}\equiv-\frac{\mu_{r}}{T_{r}},$

and split $I^{r\pm}$ into the principle value integral $I_{P}$ and the Fermi
function. Note that the lead symmetries, Eq. (16), allow us to drop the index
$r$ on the principle value integral. The approximation enables us to
reformulate Eq. (9) as

$\displaystyle\partial_{t}\rho_{bb^{\prime}}=$
$\displaystyle-i(E_{b}-E_{b^{\prime}})\rho_{bb^{\prime}}$ (19)
$\displaystyle-i\sum_{b^{\prime\prime},r=L,R}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}\Gamma_{b^{\prime\prime}a,ab^{\prime}}^{r}I^{r-}-\sum_{c}\Gamma_{b^{\prime\prime}c,cb^{\prime}}^{r}I^{r+*}\Big{]}$
$\displaystyle-i\sum_{b^{\prime\prime},r=L,R}\rho_{b^{\prime\prime}b^{\prime}}\Big{[}\sum_{c}\Gamma_{bc,cb^{\prime\prime}}^{r}I^{r+}-\sum_{a}\Gamma_{ba,ab^{\prime\prime}}^{r}I^{r-*}\Big{]}$
$\displaystyle-i\sum_{aa^{\prime},r=L,R}\rho_{aa^{\prime}}\Gamma_{ba,a^{\prime}b^{\prime}}^{r}\Big{[}I^{r+*}-I^{r+}\Big{]}$
$\displaystyle-i\sum_{cc^{\prime},r=L,R}\rho_{cc^{\prime}}\Gamma_{bc,c^{\prime}b^{\prime}}^{r}\Big{[}I^{r-*}-I^{r-}\Big{]}$

These approximations allow us to reformulate the 1vN quantum master equation
into a Lindblad form in the next section. We will start with defining matrices
that describe the jump operators of the Lindblad master equation and the Lamb-
shift.

###  .1.1 Definition of jump operators and Lamb-shift

To obtain a Linblad form corresponding to Refs. [43, 47]

$\partial_{t}\rho=-i[H,\rho]+L\rho
L^{\dagger}-\frac{1}{2}\\{L^{\dagger}L,\rho\\}$ (20)

we introduce the jump operators

$L_{kl}^{r}=T_{kl}^{r}\frac{1}{t_{r}}\begin{cases}\sqrt{f(-x^{r})}\text{ if
}N_{k}<N_{l}\\\ \sqrt{f(+x^{r})}\text{ if }N_{k}>N_{l}.\end{cases}$ (21)

The operators absorb the Fermi function into the tunnel matrix. Furthermore,
we use this definition to identify the tunneling amplitude $t_{r}$, defining
the tunneling rate

$\Gamma_{r}\equiv 2\pi\nu_{F}\,\absolutevalue{t_{r}}^{2},$ (22)

which is assumed to be symmetric $\Gamma_{R}=\Gamma_{L}=\Gamma$ and smaller
than the temperature but larger than the energies $E_{ab}\ll\Gamma_{r}\ll T$.
We introduce the Lamb-shift as

$\displaystyle H_{LS,kl}=$
$\displaystyle\pi\nu_{F}I_{P}\sum_{r,h}T_{k,h}^{r}T_{h,l}^{r}\cdot\begin{cases}+1\text{
if }N_{k}<N_{h}\\\ -1\text{ if }N_{k}>N_{h}\\\ \end{cases},$ (23)

to absorb the unitary evolution introduced by the coupled leads. From the
hermiticity of the tunnel matrices $T^{r}$, it is straightforward to show that
$H$ is hermitian

$H_{LS,kl}^{*}=H_{LS,lk}.$ (24)

In the following section we use these definitions to rewrite the dissipation
terms of Eq. (19).

###  .1.2 Redefinition of the dissipation terms of the quantum master equation

The neglected indices allow us to simplify

$I^{r\pm*}-I^{r\pm}=-2i\imaginary(I^{r\pm})=if(\pm x^{r})$ (25)

which lifts the principle value integrals within the last two summands in Eq.
(19)

$\displaystyle-i$
$\displaystyle\sum_{aa^{\prime},r}\rho_{aa^{\prime}}\Gamma_{ba,a^{\prime}b^{\prime}}^{r}\Big{[}I^{r+*}-I^{r+}\Big{]}+$
(26) $\displaystyle-i$
$\displaystyle\sum_{cc^{\prime},r}\rho_{cc^{\prime}}\Gamma_{bc,c^{\prime}b^{\prime}}^{r}\Big{[}I^{r-*}-I^{r-}\Big{]}$
$\displaystyle\overset{\ref{eq:noPrinciple}}{=}\sum_{aa^{\prime},r}\rho_{aa^{\prime}}\Gamma_{ba,a^{\prime}b^{\prime}}^{r}f(+x^{r})+\sum_{cc^{\prime},r}\rho_{cc^{\prime}}\Gamma_{bc,c^{\prime}b^{\prime}}^{r}f(-x^{r})$
$\displaystyle\overset{\ref{eq:gammaMatrix}}{=}2\pi\nu_{F}\bigg{[}\sum_{aa^{\prime},r}\rho_{aa^{\prime}}T_{ba}^{r}\,T_{a^{\prime}b^{\prime}}^{r}f(+x^{r})+$
$\displaystyle\hskip
42.67912pt+\sum_{cc^{\prime},r}\rho_{cc^{\prime}}T_{bc}^{r}\,T_{c^{\prime}b^{\prime}}^{r}f(-x^{r})\bigg{]}$
$\displaystyle\overset{\ref{eq:jumpOp}}{=}2\pi\nu_{F}\bigg{[}\sum_{aa^{\prime},r}\absolutevalue{t_{r}}^{2}L_{ba}^{r}\,\rho_{aa^{\prime}}\,(L^{r\dagger})_{a^{\prime}b^{\prime}}+$
$\displaystyle\hskip
42.67912pt+\sum_{cc^{\prime},r}L_{bc}^{r}\,\rho_{cc^{\prime}}\,(L^{r\dagger})_{c^{\prime}b^{\prime}}\bigg{]}$
$\displaystyle=2\pi\nu_{F}\cdot\sum_{r}\absolutevalue{t_{r}}^{2}(L^{r}\,\rho\,L^{r\dagger})_{bb^{\prime}}$
$\displaystyle\overset{\ref{eq:gammaScalar}}{=}\sum_{r}\Gamma_{r}\cdot(L^{r}\,\rho\,L^{r\dagger})_{bb^{\prime}}.$

Next is the first dissipation term for which the Fermi function and the
principle value integral in $I^{r\pm}$ both yield a non-zero contribution. We
find

$\displaystyle-i\sum_{b^{\prime\prime},r}$
$\displaystyle\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}\Gamma_{b^{\prime\prime}a,ab^{\prime}}^{r}I^{r-}-\sum_{c}\Gamma_{b^{\prime\prime}c,cb^{\prime}}^{r}I^{r+*}\Big{]}$
(27) $\displaystyle\overset{\ref{eq:gammaMatrix}}{=}$
$\displaystyle-i\,2\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}T_{b^{\prime\prime}a}^{r}\,T_{ab^{\prime}}^{r}I^{r-}-\sum_{c}T_{b^{\prime\prime}c}^{r}\,T_{cb^{\prime}}^{r}I^{r+*}\Big{]}$
$\displaystyle\overset{\ref{eq:simplified}}{=}$
$\displaystyle-i\,\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}T_{b^{\prime\prime}a}^{r}\,T_{ab^{\prime}}^{r}I_{P}-\sum_{c}T_{b^{\prime\prime}c}^{r}\,T_{cb^{\prime}}^{r}I_{P}\Big{]}-\frac{1}{2}\,2\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}T_{b^{\prime\prime}a}^{r}\,T_{ab^{\prime}}^{r}f(-x^{r})+\sum_{c}T_{b^{\prime\prime}c}^{r}\,T_{cb^{\prime}}^{r}f(x^{r})\Big{]}$
$\displaystyle\overset{\ref{eq:jumpOp}}{=}$
$\displaystyle-i\,\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}T_{b^{\prime\prime}a}^{r}\,T_{ab^{\prime}}^{r}I_{P}-\sum_{c}T_{b^{\prime\prime}c}^{r}\,T_{cb^{\prime}}^{r}I_{P}\Big{]}-\frac{1}{2}\,2\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\absolutevalue{t_{r}}^{2}\,\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}(L^{r\dagger})_{b^{\prime\prime}a}\,L_{ab^{\prime}}^{r}+\sum_{c}(L^{r\dagger})_{b^{\prime\prime}c}\,L_{cb^{\prime}}^{r}\Big{]}$
$\displaystyle\overset{\ref{eq:lambShiftDef}}{=}$
$\displaystyle+i\sum_{b^{\prime\prime}}\rho_{bb^{\prime\prime}}H_{LS,b^{\prime\prime}b^{\prime}}-\frac{1}{2}\,2\pi\nu_{F}\,\sum_{b^{\prime\prime},r}\absolutevalue{t_{r}}^{2}\,\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}(L^{r\dagger})_{b^{\prime\prime}a}\,L_{ab^{\prime}}^{r}+\sum_{c}(L^{r\dagger})_{b^{\prime\prime}c}\,L_{cb^{\prime}}^{r}\Big{]}$
$\displaystyle\overset{\ref{eq:gammaScalar}}{=}$
$\displaystyle+i\,\sum_{b^{\prime\prime}}\rho_{bb^{\prime\prime}}H_{LS,b^{\prime\prime}b^{\prime}}-\frac{1}{2}\,\sum_{b^{\prime\prime},r}\Gamma_{r}\,\rho_{bb^{\prime\prime}}\Big{[}\sum_{a}(L^{r\dagger})_{b^{\prime\prime}a}\,L_{ab^{\prime}}^{r}+\sum_{c}(L^{r\dagger})_{b^{\prime\prime}c}\,L_{cb^{\prime}}^{r}\Big{]}$
$\displaystyle=$ $\displaystyle+i\,(\rho
H)_{b,b^{\prime}}-\frac{1}{2}\,\sum_{r}\Gamma_{r}\,(\rho\,L^{r\dagger}\,L^{r})_{b,b^{\prime}}$

With the same steps as in Eq. (27), we find for the last remaining summand of
Eq. (19)

$\displaystyle-i$
$\displaystyle\sum_{b^{\prime\prime},r}\rho_{b^{\prime\prime}b^{\prime}}\Big{[}\sum_{c}\Gamma_{bc,cb^{\prime\prime}}^{r}I^{r+}-\sum_{a}\Gamma_{ba,ab^{\prime\prime}}^{r}I^{r-*}\Big{]}$
(28) $\displaystyle=$
$\displaystyle-i\,(H\,\rho)_{b,b^{\prime}}-\frac{1}{2}\,\sum_{r}\Gamma_{r}\,(L^{r\dagger}\,L^{r}\,\rho)_{b,b^{\prime}}$

Finally, we substitute Eqs. (26), (27) and (28) into Eq. (19) and obtain

$\displaystyle\partial_{t}\rho_{bb^{\prime}}=$
$\displaystyle-i(E_{b}-E_{b^{\prime}})\rho_{bb^{\prime}}+i(H\,\rho)_{bb^{\prime}}-i\,(\rho\,H)_{b^{\prime\prime},b^{\prime}}$
(29)
$\displaystyle+\sum_{r}\Gamma_{r}\cdot(L^{r}\,\rho\,L^{r\dagger})_{bb^{\prime}}$
$\displaystyle-\frac{1}{2}\Gamma_{r}\,(\rho\,L^{r\dagger}\,L^{r})_{b,b^{\prime}}-\frac{1}{2}\,\sum_{r}\Gamma_{r}\,(L^{r\dagger}\,L^{r}\,\rho)_{b,b^{\prime}}$
$\displaystyle=$
$\displaystyle-i\left[H_{MB}+H_{LS},\,\rho\right]_{bb^{\prime}}$
$\displaystyle+\Gamma\sum_{r}\left(L^{r}\,\rho\,L^{r\dagger}-\frac{1}{2}\\{\rho,\,L^{r\dagger}\,L^{r}\\}\right)_{bb^{\prime}}.$

In the next section follows an explicit representation of the jump operators
and the Lamb-shift for the clean Majorana box.

##  .2 Lindblad master equation of the clean box

Without any overlaps, we are free to choose the basis in which we combine the
MBSs without obtaining a non-diagonal Hamiltonian for the Majorana box. For
simplicity, we combine the left and the right MBSs

$\displaystyle f_{L}=\frac{1}{2}(\gamma_{Lu}+i\,\gamma_{Ld}),$ $\displaystyle
f_{R}=\frac{1}{2}(\gamma_{Rd}+i\,\gamma_{Ru}).$ (30)

This defines the Fock states of the system as

$\displaystyle\ket{n_{L},n_{R}},$
$\displaystyle\mathcal{F}=\\{\ket{00},\,\ket{11},\,\ket{10},\,\ket{01}\\},$
(31)

where the choice $\ket{11}=f_{L}^{\dagger}f_{R}^{\dagger}\ket{00}$ fixes the
phase. Next, we use the fermionic operators $f_{L},\,f_{R}$ to rewrite the
tunnel Hamiltonian, Eq. (2), into

$\displaystyle H_{T}=t\sum_{k}\bigg{(}$ $\displaystyle
c_{L,k}\,[(e^{i\phi_{L}}+i)f_{L}^{\dagger}+(e^{i\phi_{L}}-i)f_{L}]$ (32)
$\displaystyle
c_{R,k}\,[(e^{i\phi_{R}}+i)f_{R}^{\dagger}+(e^{i\phi_{R}}-i)f_{R}]\bigg{)}+h.c.$

From this we obtain the tunnel matrices defined in Eq. (33).

###  .2.1 Evaluation of jump operators

Using the tunnel matrices, we evaluate the jump operators and the Lamb-shift.
Because of our assumption that the chemical potential is far bigger than the
temperature of the leads, $\mu\gg k_{B}T$, we approximate the Fermi functions
of Eq. (21) to 0/1 depending on the sign of $x^{r}$. The physical
interpretation for this is that the electrons are able to tunnel from the left
lead into the box and from the box into the right lead but not the other way
around. We obtain

$\displaystyle T_{L}=t\begin{pmatrix}\mbox{\Large
0}_{2}&&\begin{matrix}e^{-i\phi_{L}}-i&0\\\ 0&e^{-i\phi_{L}}+i\end{matrix}\\\
\begin{matrix}e^{i\phi_{L}}+i&0\\\ 0&e^{i\phi_{L}}-i\end{matrix}&&\mbox{\Large
0}_{2}\end{pmatrix},$ $\displaystyle T_{R}=t\begin{pmatrix}\mbox{\Large
0}_{2}&&\begin{matrix}0&-ie^{-i\phi_{R}}+1\\\
ie^{-i\phi_{R}}+1&0\end{matrix}\\\ \begin{matrix}0&-ie^{i\phi_{R}}+1\\\
ie^{i\phi_{R}}+1&0\end{matrix}&&\mbox{\Large 0}_{2}\end{pmatrix}.$ (33)
$\displaystyle L_{L}$ $\displaystyle=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\
e^{+i\phi_{L}}+i&0&0&0\\\ 0&e^{+i\phi_{L}}-i&0&0\\\ \end{pmatrix},$ (34)
$\displaystyle L_{R}$
$\displaystyle=\begin{pmatrix}0&0&0&-ie^{-i\phi_{R}}+1\\\
0&0&ie^{-i\phi_{R}}+1&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ \end{pmatrix}.$

Apart from the dynamics, we also need the jump operators to evaluate the
current

$I=\mathrm{Tr}\big{(}L_{L}\,\rho\,L_{L}^{\dagger}\big{)},$ (35)

see [43].

###  .2.2 Evaluation of the Lamb-shift

Inserting Eq. (33) into Eq. (23) yields

$\displaystyle H_{LS}=$ $\displaystyle\Gamma
I_{P}\bigg{[}2\begin{pmatrix}-I_{2}&0_{2}\\\ 0_{2}&I_{2}\end{pmatrix}$ (36)
$\displaystyle+\begin{pmatrix}\sigma_{z}(\sin\phi_{L}+\sin\phi_{R})&0_{2}\\\
0_{2}&-\sigma_{z}(\sin\phi_{L}-\sin\phi_{R})\end{pmatrix}\bigg{]}$

The first term of Eq. (36) only introduces a constant energy splitting between
the even and odd parities. It is of the order $\Gamma\ll\mu$, so we can safely
neglect its effect on the dynamics. The second term of the Lamb-shift, Eq.
(36), is more significant. It introduces a $\phi_{L},\phi_{R}$ dependent
rotation around the $z$-axis within each parity sector. Interestingly, this
second term takes the form of additional overlaps between the left (right)
pair of MBSs mediated via the leads. We can write it as

$H_{LS}=\frac{i}{2}\varepsilon_{L}\,\gamma_{Lu}\gamma_{Ld}+\frac{i}{2}\varepsilon_{R}\,\gamma_{Ru}\gamma_{Rd}$
(37)

for the effective overlaps

$\displaystyle\varepsilon_{L}\equiv-2\Gamma I_{P}\sin\phi_{L},$
$\displaystyle\varepsilon_{R}\equiv 2\Gamma I_{P}\sin\phi_{R}.$ (38)

In the next section, we insert these expressions into the Lindblad form and
rewrite the dynamics into Bloch equations.

##  .3 Representation in Bloch equations

Up to now, it was important to include both phases to understand how the
contributions from the left and right leads to the Lamb-shift add up or
subtract. In order to understand the blockade it is sufficient to consider one
phase, so from now on we choose $\phi_{R}=0$ such that $\varepsilon_{R}=0$. We
insert Eq. (34) and Eq. (36) into the Lindblad form in Eq. (29) and, after
some matrix multiplication, we obtain a differential equation of the form

$\displaystyle\partial_{t}\rho=-i[H_{LS},\rho]+2\,\Gamma(-\rho+D),$ (39)

where the matrices $H_{LS},D$ are defind in Eq. (44). Due to the total parity
being a good quantum number the density matrix is block-diagonal

$\displaystyle\rho=\begin{pmatrix}\rho_{e}&0_{2}\\\
0_{2}&\rho_{o}\end{pmatrix}.$ (40)

We separate it into sectors of even ($\rho_{e}$) and odd ($\rho_{o}$) total
parity. The last step of the derivation is to define the Bloch vectors for the
even/odd parity sector $s_{e}/s_{o}$. They are chosen such that they fulfill

$\displaystyle\rho_{e}=\frac{p_{e}\,I_{2}+\vec{s}_{e}\,\vec{\sigma}}{2},$
$\displaystyle\rho_{o}=\frac{p_{o}\,I_{2}+\vec{s}_{o}\,\vec{\sigma}}{2}.$ (41)

In this notation $I_{2}$ is the identity and $\vec{\sigma}$ a vector of Pauli
matrices. $p_{e}/p_{o}$ are the probabilities to measure the system in
even/odd parity. Probability normalization gives

$1=\mathrm{Tr}(\rho)=p_{e}+p_{o}.$ (42)

In total we have seven independent variables. The differential equation for
the even probability and the $z$ component of the Bloch vectors in both
sectors read

$\displaystyle\partial_{t}p_{e}=2\,\Gamma(-p_{e}+p_{o}-\sin(\phi_{L})\,s_{e}^{z}),$
(43)
$\displaystyle\partial_{t}s_{e}^{z}=2\,\Gamma(-s_{e}^{z}-s_{o}^{z}-\sin(\phi_{L})\,p_{e}),$
$\displaystyle\partial_{t}s_{o}^{z}=2\,\Gamma(-s_{o}^{z}+s_{e}^{z}+\sin(\phi_{L})\,p_{e}).$
$\displaystyle H_{LS}=\frac{1}{2}\begin{pmatrix}-\varepsilon_{L}&0&0&0\\\
0&\varepsilon_{L}&0&0\\\ 0&0&\varepsilon_{L}&0\\\
0&0&0&-\varepsilon_{L}\end{pmatrix},$ $\displaystyle
D=\begin{pmatrix}\rho_{33}-\sin(\phi_{L})\rho_{00}&-i\rho_{32}&0&0\\\
i\rho_{23}&\rho_{22}+\sin(\phi_{L})\rho_{11}&0&0\\\
0&0&(1+\sin(\phi_{L}))\rho_{00}&i\cos(\phi_{L})\rho_{01}\\\
0&0&-i\cos(\phi_{L})\rho_{10}&(1-\sin(\phi_{L}))\rho_{11}\\\ \end{pmatrix}$
(44)

The equations are decoupled from the ones for the $x$ and $y$ components. Note
that the $\sigma_{z}$ rotation of the Lamb-shift describes a unitary evolution
within the $x-y$-plane. Accordingly, there is no contribution from it to Eq.
(43). Before we go on to the solution, we also need to express the current,
specified in Eq. (35), in the Bloch representation as

$I=2\,\Gamma(p_{e}+\sin\phi_{L}\,s_{e}^{z}),$ (45)

which is also the version used in the main paper. Motivated by Evans theorem,
which guarantees us always at least one zero eigenvalue for a Lindblad
equation [54, 55, 56], we will now determine the stationary solution
($\partial_{t}\rho=0$) for large times.

Let’s first consider the case $\phi_{L}=0$. It is easy to show that the only
possibility for a stationary solution is given if both $z$ components are zero
and the probability to find the system with even or odd parity are equally
high

$\displaystyle\partial_{t}\rho=0$ $\displaystyle\Leftrightarrow$
$\displaystyle s_{e}^{z}=s_{o}^{z}=0,\,p_{e}=p_{o}=\frac{1}{2},$ (46)

leading to a current

$I=\Gamma.$ (47)

This is the unblocked case of electrons tunneling with rate $\Gamma$ through
the device. We now turn to the more interesting blocking case and choose
$\phi_{L}=\frac{\pi}{2}$ for which we find that

$\displaystyle\partial_{t}\rho=0$ $\displaystyle\Leftrightarrow$
$\displaystyle p_{e}=1,\,s_{e}^{z}=-1,\,p_{o}=s_{o}^{z}=0.$ (48)

In this setting the only possibility for a stationary solution is that the
even parity Bloch vector is aligned anti-parallel to the $z$-axis and the
system is projected on the even parity sector. This results in $I=0$. Indeed,
a blockade is established. In the next section we will solve the system for
finite overlaps between the MBSs and quantify the lifting of the blockade.

##  .4 Lifting of the parity blockade due to finite overlap

We introduce small but finite overlaps between the MBSs

$\displaystyle 0<\varepsilon_{u},\,\varepsilon_{d}\ll\Gamma.$ (49)

From this we obtain the Hamiltonian of the system

$\displaystyle H_{MB}=\begin{pmatrix}-\Delta_{e}\,\sigma_{x}&0_{2}\\\
0_{2}&\Delta_{o}\,\sigma_{x}\end{pmatrix},$ (50)

it introduces a rotation around the $x$-axis with the frequency

$\displaystyle\Delta_{e}\equiv\frac{\varepsilon_{u}+\varepsilon_{d}}{2},$
$\displaystyle\Delta_{o}\equiv\frac{\varepsilon_{u}-\varepsilon_{d}}{2},$ (51)

for the even and odd parity sector respectively. These rotations couple the
$z$-axis to the (previously uncoupled) $x-y$-plane. The corresponding unitary
evolution reads

$\displaystyle-i[H_{MB},\rho]=$
$\displaystyle-i\begin{pmatrix}-\Delta_{e}\,[\sigma_{x},\rho_{e}]&0_{2}\\\
0_{2}&\Delta_{o}\,[\sigma_{x},\rho_{o}]\end{pmatrix}$ (52) $\displaystyle=$
$\displaystyle\begin{pmatrix}\Delta_{e}\,(s_{e}^{z}\sigma_{y}-s_{e}^{y}\sigma_{z})&0_{2}\\\
0_{2}&-\Delta_{o}\,(s_{o}^{z}\sigma_{y}-s_{o}^{y}\sigma_{z})\end{pmatrix}.$

Because of the overlaps we need to include all seven degrees of freedom to
solve for the dynamics

$\displaystyle\partial_{t}p_{e}$
$\displaystyle=2\Gamma(1-2p_{e}-\sin\phi_{L}\,s_{e}^{z})$ (53)
$\displaystyle\partial_{t}s_{e}^{z}$
$\displaystyle=2\Gamma(-s_{e}^{z}-s_{o}^{z}-\sin\phi_{L}\,p_{e}-\xi\,s_{e}^{y})$
$\displaystyle\partial_{t}s_{o}^{z}$
$\displaystyle=2\Gamma(-s_{o}^{z}+s_{e}^{z}+\sin\phi_{L}\,p_{e}+\eta\,s_{o}^{y})$
$\displaystyle\partial_{t}s_{e}^{x}$
$\displaystyle=2\Gamma(-s_{e}^{x}+s_{o}^{y}+\varepsilon_{L}/2\Gamma\,s_{e}^{y})$
$\displaystyle\partial_{t}s_{e}^{y}$
$\displaystyle=2\Gamma(-s_{e}^{y}+s_{o}^{x}-\varepsilon_{L}/2\Gamma\,s_{e}^{x}+\xi\,s_{e}^{z})$
$\displaystyle\partial_{t}s_{o}^{x}$
$\displaystyle=2\Gamma(-s_{o}^{x}+\cos\phi_{L}\,s_{e}^{y}-\varepsilon_{L}/2\Gamma\,s_{o}^{y})$
$\displaystyle\partial_{t}s_{o}^{y}$
$\displaystyle=2\Gamma(-s_{o}^{y}-\cos\phi_{L}\,s_{e}^{x}+\varepsilon_{L}/2\Gamma\,s_{o}^{x}-\eta\,s_{o}^{z}),$

with the small parameters

$\displaystyle\xi\equiv\frac{\Delta_{e}}{\Gamma},$
$\displaystyle\eta\equiv\frac{\Delta_{o}}{\Gamma}.$ (54)

As we are interested in the blockade lifting, a general solution to this
equation is not needed. Instead, we solve this by expanding the solution for
$\phi_{L}=\frac{\pi}{2}$ to leading order in $\eta,\,\xi$

$\displaystyle s_{o}^{y},s_{o}^{x},s_{e}^{y},s_{e}^{x}=0+O(\xi,\eta),$
$\displaystyle p_{e}=+1+O(\xi,\eta),$ (55) $\displaystyle
s_{o}^{z}=0+O(\xi,\eta),$ $\displaystyle s_{e}^{z}=-1+O(\xi,\eta).$ (56)

Setting $\phi_{L}=\frac{\pi}{2}$ also sets the Lamb-shift to
$\varepsilon_{L}=-2\Gamma I_{P}$ and we obtain

$\displaystyle\partial_{t}p_{e}$ $\displaystyle=2\Gamma(1-2p_{e}-s_{e}^{z})$
(57) $\displaystyle\partial_{t}s_{e}^{z}$
$\displaystyle=2\Gamma(-s_{e}^{z}-s_{o}^{z}-p_{e}-\xi\,s_{e}^{y})$
$\displaystyle\partial_{t}s_{o}^{z}$
$\displaystyle=2\Gamma(-s_{o}^{z}+s_{e}^{z}+p_{e}+\eta\,s_{o}^{y})$
$\displaystyle\partial_{t}s_{e}^{x}$
$\displaystyle=2\Gamma(-s_{e}^{x}+s_{o}^{y}-I_{P}s_{e}^{y})$
$\displaystyle\partial_{t}s_{e}^{y}$
$\displaystyle=2\Gamma(-s_{e}^{y}+s_{o}^{x}+I_{P}s_{e}^{x}+\xi\,s_{e}^{z})$
$\displaystyle\partial_{t}s_{o}^{x}$
$\displaystyle=2\Gamma(-s_{o}^{x}+I_{P}s_{o}^{y})$
$\displaystyle\partial_{t}s_{o}^{y}$
$\displaystyle=2\Gamma(-s_{o}^{y}-I_{P}s_{o}^{x}-\eta\,s_{o}^{z}).$

Inserting $s_{o}^{z}=0+O(\xi,\eta)$ into the last two lines of Eq. (57) we
find that

$\displaystyle s_{o}^{x},s_{o}^{y}=O(\eta^{2},\xi\,\eta).$ (58)

Therefore, we can neglect their contribution on the remaining set of
equations. We find for the 5th and 6th line of Eq. (57)

$\displaystyle\partial_{t}s_{e}^{x}$
$\displaystyle=2\Gamma(-s_{e}^{x}-I_{P}s_{e}^{y})$ (59)
$\displaystyle\partial_{t}s_{e}^{y}$
$\displaystyle=2\Gamma(-s_{e}^{y}+I_{P}s_{e}^{x}+\xi\,s_{e}^{z}),$ (60)

which yields by inserting $s_{e}^{z}=-1+O(\xi,\eta)$

$\displaystyle s_{e}^{y}=\frac{-\xi}{1+I_{P}^{2}},$ $\displaystyle
s_{e}^{x}=\frac{-2I_{P}\xi}{1+I_{P}^{2}}.$ (61)

In turn we insert this into the 1st, 2nd and 3rd line of Eq. (57) and get

$\displaystyle\partial_{t}p_{e}$ $\displaystyle=2\Gamma(1-2p_{e}-s_{e}^{z}),$
(62) $\displaystyle\partial_{t}s_{e}^{z}$
$\displaystyle=2\Gamma(-s_{e}^{z}-s_{o}^{z}-p_{e}+\frac{\xi^{2}}{1+I_{P}^{2}}),$
$\displaystyle\partial_{t}s_{o}^{z}$
$\displaystyle=2\Gamma(-s_{o}^{z}+s_{e}^{z}+p_{e}).$

These equations are the same ones written in Eq. (43) for the fully blocked
system, with the effect of the overlaps included as an additional loss term
$\frac{\xi^{2}}{1+I_{P}^{2}}$. We finally obtain the stationary state solution

$\displaystyle p_{e}=1-\frac{1}{2}\frac{\xi^{2}}{1+I_{P}^{2}},$ $\displaystyle
s_{e}^{z}=-1+\frac{\xi^{2}}{1+I_{P}^{2}},$ $\displaystyle
s_{o}^{z}=\frac{1}{2}\frac{\xi^{2}}{1+I_{P}^{2}},$ (63)

which resembles the previous solution with a small misalignment
$\propto\xi^{2}$ from the south pole of the Bloch sphere. Finally, we insert
this into Eq. (45) and find for the remnant current

$\displaystyle
I_{rem}=\Gamma\,\frac{\xi^{2}}{1+I_{P}^{2}}\overset{\xi=\Delta_{e}/\Gamma}{=}\frac{\Delta_{e}^{2}}{1+I_{P}^{2}}\,\frac{1}{\Gamma},$
(64)

which is counter-intuitively proportional to $1/\Gamma$.

##  .5 Lifting of the parity blockade due to deviations in the tunnel coupling

In this last section, we provide a short argument about the effect of a small
deviation from the perfect blockade on the stationary current. For this we
introduce the small deviation parameters $\delta_{t},\delta_{\phi}\ll 1$ as

$\displaystyle t_{Ld}=(1-\delta_{t})\,t,$
$\displaystyle\phi_{L}=\frac{\pi}{2}+\delta_{\phi},$
$\displaystyle\delta=\sqrt{\delta_{t}^{2}+\delta_{\phi}^{2}}.$ (65)

In a fine-tuned setting the blockade occurs because one of the tunneling
matrix elements evaluates to zero

$\displaystyle
t_{Lu}+t_{Ld}=t\,(e^{i\phi_{L}}+i)\overset{\delta_{\phi}=\delta_{t}=0}{=}0.$
(66)

If we include finite deviations we obtain

$\displaystyle t_{Lu}+t_{Ld}$
$\displaystyle=t\,e^{i\frac{\pi}{2}+\delta_{\phi}}+it\,(1-\delta_{t})$ (67)
$\displaystyle\approx t(-\delta_{\phi}+i\delta_{t}),$

in leading order. We view this finite matrix element as a new tunnel coupling,
which is added to the blocked dynamics. It defines the deviation tunneling
rate as

$\displaystyle\Gamma_{dev}=\absolutevalue{t(-\delta_{\phi}+i\delta_{t})}^{2}=t^{2}\delta^{2}.$
(68)

As we are only interested in the leading order, we can neglect the effect of
this additional tunnel coupling on the stationary state. We insert
$\Gamma_{dev}$ into Eq. (45) for the stationary state from the perfect
blockade Eq. (48) to obtain

$\displaystyle I_{dev}=2\,\Gamma_{dev}.$ (69)

This leading order contribution also holds if we include finite overlaps
between the MBSs.


# Extreme value theory and the St. Petersburg paradox in the failure
statistics of wires

Alessandro Taloni1, Stefano Zapperi2,3 1 CNR-ISC, Institute for Complex
Systems, Via dei Taurini 19, 00185 Roma, Italy  2 Center for Complexity and
Biosystems and Department of Physics, University of Milano, Via Celoria 16,
20133 Milano, Italy  3 CNR-ICMATE, Institute of Condensed Matter Chemistry and
Technologies for Energies, Via Roberto Cozzi 53, 20125 Milano, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The fracture stress of materials typically depends on the sample size and is
traditionally explained in terms of extreme value statistics. A recent work
reported results on the carrying capacity of long polyamide and polyester
wires and interpret the results in terms of a probabilistic argument known as
the St. Petersburg paradox. Here, we show that the same results can be better
explained in terms of extreme value statistics. We also discuss the relevance
of rate dependent effects.

††: J. Phys. G: Nucl. Part. Phys.

## Introduction

The failure of macroscopic objects is determined both by intrinsic properties,
such as its crystalline structure, and by extrinsic factors as the grain
boundaries, vacancies, dislocations or microcracks. In particular, the
strength of brittle materials is largely influenced by the presence of flaws,
as it was first demonstrated by Griffith in his seminal paper 100 years ago
[Griffith(1921)]. The Griffith’s criterion is the foundation stone of fracture
mechanics, showing, by simple linear elasticity arguments, that a material
fails at a lower nominal applied stress when a crack is present.

The Griffith theory is also relevant to describe fracture size effects or the
observation that larger samples are generally weaker. In fact, the larger the
sample, the higher is the probability to find microcraks that might become
unstable when a threshold stress is reached. According to the Griffith’s
theory this threshold is inversely proportional to the square root of defect
size. Thus the strength of a material corresponds to the stress needed to
trigger the first crack instability event, namely the smallest threshold
corresponding to the largest flaw. In this sequence of arguments lies the
connection between fracture mechanics , materials failure and extreme value
statistics (for a review see [Taloni  _et al._(2018)Taloni, Vodret,
Costantini, and Zapperi]).

In a recent letter [Ball(2020), Fontana and Palffy-Muhoray(2020)], the authors
proposed a connection between materials failure statistics and the St.
Petersburg paradox by linking the average load carrying capacity of a wire to
its length. The result, however, was derived assuming that ”the force required
to fracture the fiber is a linear function of the defect size” [Fontana and
Palffy-Muhoray(2020)], which is in glaring contrast with fracture mechanics.
Here we address the problem combining extreme value theory (EVT) [Taloni  _et
al._(2018)Taloni, Vodret, Costantini, and Zapperi] with the Griffith’s
stability crack criterion [Griffith(1921)]. We also show that in the
asymptotic limit, the wire strength follows the Gumbel’s distribution, in full
agreement with the data reported in [Fontana and Palffy-Muhoray(2020)], as we
demonstrate using the maximum likelihood method. We thus conclude that the
load carrying capacity of the wires studied in [Fontana and Palffy-
Muhoray(2020)] follows EVT, in agreement with previous observations for
different materials [Taloni  _et al._(2018)Taloni, Vodret, Costantini, and
Zapperi].

##  1 Theory

Let us consider a wire of length $L$ which we divide in $N=L/L_{0}$
independent elements of size $L_{0}$. Following Ref. [Fontana and Palffy-
Muhoray(2020)], we want to relate the statistics of the micro-cracks present
in the wire with its failure strength.

Defining $P(n)$ as the probability density function (pdf) of micro-cracks of
length $w\equiv nL_{0}$, we introduce

$F(z)=\int_{0}^{z}dnP(n)$ (1)

as the probability that no micro-crack larger than z$L_{0}$ will be found in
the wire. Next, we look for the statistics of the largest micro-crack in the
wire $w_{nmax}=n_{max}L_{0}$ which satisfies the only constraint
$w_{max}=n_{max}L_{0}\ll L$ [Duxbury  _et al._(1987)Duxbury, Leath, and
Beale]. If $\rho_{N}(n_{max})$ is the pdf for the largest micro-crack, then

$F_{N}(z)=\int_{0}^{z}dn_{max}\,\rho_{N}(n_{max})=\left[F(z)\right]^{N}.$ (2)

The Fisher-Gnedenko-Tippet theorem ensures that $F_{N}(z)\to G(z)$ for large
$N$, where $G$ belongs to one of three families only: Weibull, Fréchet or
Gumbel [Fisher and Tippett(1928), Gnedenko(1943)]. The convergence to either
one of these universal distributions depends on the asymptotic properties of
$P(n)$ [Leadbetter  _et al._(2012)Leadbetter, Lindgren, and Rootzén, Manzato
_et al._(2012)Manzato, Shekhawat, Nukala, Alava, Sethna, and Zapperi]. If the
distribution of micro-cracks has an exponential tail [Kunz and
Souillard(1978), Manzato  _et al._(2012)Manzato, Shekhawat, Nukala, Alava,
Sethna, and Zapperi, Duxbury  _et al._(1987)Duxbury, Leath, and Beale], $G(z)$
converges asymptotically to the Gumbel distribution [Gumbel(2004)]:

$F_{N}(z)\to G(z)=e^{-Ne^{-\frac{z-\mu}{\beta}}}$ (3)

To derive the fracture strength, the authors of Ref. [Fontana and Palffy-
Muhoray(2020)] assume that it is linearly dependent on the size of the largest
defect ($n_{max}$), obtained through the analogy with the St. Petersburg
model. However, this assumption is not justified by fracture mechanics. A
relation between crack length and fracture strength in an elastic medium is
provided by the Griffith’s stability criterion. It states that a crack of
length $w$ subject to a normal stress $\sigma$ is stable as long as
$\sigma<K_{1C}/Y\,w^{-1/2}$[Griffith(1921)], where $K_{1C}$ is the critical
stress intensity factor and $Y$ is a geometric factor. In our context, the
wire should break when the largest micro-crack becomes unstable, i.e. when it
is reached the stress

$\sigma_{break}=K_{1C}/Y\,w_{max}^{-1/2}.$ (4)

This relation corresponds to the weakest link hypothesis in elastic mechanics
[Taloni  _et al._(2018)Taloni, Vodret, Costantini, and Zapperi]. Thanks to 4,
we can formally write down the probability that a wire of length $L$ brakes at
stress $\sigma_{break}$ as

$\rho_{L}(\sigma_{break})=-\frac{2K_{1C}^{2}}{Y^{2}L_{0}\sigma_{break}^{3}}\rho_{N}\left(n_{max}=\frac{K_{1C}^{2}}{Y^{2}L_{0}\sigma_{break}^{2}}\right),$
(5)

and, by means of 2, introduce the probability that the wire does not fail up
to a stress $\sigma$

$\Sigma_{L}(\sigma)=F_{N}(z)$ (6)

with $z=\frac{K_{1C}^{2}}{Y^{2}L_{0}\sigma^{2}}$. For large samples, the
convergence theorem 3 yields

$\Sigma_{L}(\sigma)\sim
e^{-\frac{L}{L_{0}}\,e^{-\left(\frac{\sigma_{0}}{\sigma}\right)^{2}}},$ (7)

with $\sigma_{0}=\frac{K_{1C}}{Y\sqrt{L_{0}\beta}}$. This is the Duxbury-
Leath-Beale distribution [Duxbury  _et al._(1987)Duxbury, Leath, and Beale],
which was shown to converge to the Gumbel distribution as $L\gg L_{0}$ , i.e.

$\Sigma_{L}(\sigma)\to e^{-\frac{L}{L_{0}}\,e^{\frac{\sigma-\mu}{\beta}}}$ (8)

[Manzato  _et al._(2012)Manzato, Shekhawat, Nukala, Alava, Sethna, and
Zapperi]. Hence, in analogy with the expression 2, we can define the survival
probability for a single crack as

$\Sigma_{0}(\sigma)=e^{-e^{\frac{\sigma-\mu}{\beta}}}.$ (9)

Finally, the average breaking stress of a wire of length $L$ is given by

$\langle\sigma\rangle=\int_{0}^{\infty}d\sigma\Sigma_{L}(\sigma)=\beta\left[\gamma-\ln\left(L/L_{0}\right)\right]+\mu,$
(10)

which recovers Eq.(1) of [Fontana and Palffy-Muhoray(2020)] and $\gamma$ is
the Euler-Mascheroni constant.

##  2 Comparison with the experiments

In Ref.[Fontana and Palffy-Muhoray(2020)] tensile experiments were performed
on polyester and polyamide wires. In a first series of experiments the
fracture strengths of samples ranging from 1 mm to 1 km were reported. A
second sequence of experiments showed a small strain rate dependence of the
fracture forces for wires of fixed lengths, in both materials.

To capture and reinterpret the experimental results within the framework of
the classical fracture mechanics, we have used the Gumbel form of the
generalized EVT, which is purposely designed to account for strain and thermal
effects [Taloni  _et al._(2018)Taloni, Vodret, Costantini, and Zapperi,
Sellerio  _et al._(2015)Sellerio, Taloni, and Zapperi]. In this context we can
introduce the survival distribution function as the probability that a wire of
volume $V_{i}$ is still intact up to a stress $\sigma_{i}$, if the tensile
experiment is conducted at a strain rate $\dot{\varepsilon_{i}}$ and
temperature $T$,

$\Sigma_{V_{i}}(\sigma_{i};\dot{\varepsilon_{i}},T)=\left[\int_{\sigma_{i}/E}^{\infty}d\varepsilon_{f}\rho_{0}(\varepsilon_{f})S_{0}(\sigma_{i};\dot{\varepsilon_{i}},T|\varepsilon_{f})\right]^{\frac{V_{i}}{V_{0}}}.$
(11)

The quantity $\rho_{0}(\varepsilon_{f})$ corresponds to the Gumbel failure
strain probability density function

$\rho_{0}(\varepsilon_{f})=\frac{e^{\frac{E\varepsilon_{f}-\mu}{E\varepsilon_{0}}-e^{\frac{E\varepsilon_{f}-\mu}{E\varepsilon_{0}}}}}{\varepsilon_{0}},$
(12)

obtained by differentiating the expression 9
($\rho_{0}(\sigma)=-\frac{d\Sigma_{0}(\sigma)}{d\sigma}$) and substituting
$\sigma=E\varepsilon_{f}$ and $\beta=E\varepsilon_{0}$, with $E$ Young’s
modulus. The thermal factor $S_{0}$ can be derived from the Kramer’s theory
for the transition rate as [Sellerio  _et al._(2015)Sellerio, Taloni, and
Zapperi]

$S_{0}(\sigma_{i};\dot{\varepsilon_{i}},T|\varepsilon_{f})=e^{-\frac{\omega_{0}}{\sqrt{2\pi\dot{\varepsilon_{i}}}}\sqrt{\frac{EV_{0}}{k_{B}T}}\left[e^{-\frac{EV_{0}}{2k_{B}T}\left(\varepsilon_{f}-\sigma_{i}/E\right)^{2}}-e^{-\frac{EV_{0}}{2k_{B}T}\varepsilon_{f}^{2}}\right]},$
(13)

where $\omega_{0}$ is a characteristic frequency and $k_{B}$ the Boltzmann’s
constant.

The data reported in [Fontana and Palffy-Muhoray(2020)] correspond to the
fracture forces $F_{i}$ for fibers of lengths $L_{i}$. Each experiment was
performed at a different strain rate, such that
$\dot{\varepsilon_{i}}=\frac{\alpha}{L_{i}}$, where $\alpha=6\ 10^{-3}m/sec$
if $L_{i}<1m$ or $\alpha=71\ 10^{-3}m/sec$ if $L_{i}\geq 1m$. Under these
conditions, the proper way of fitting the experimental fracture data is the
Maximum Likelihood (ML) method [Taloni  _et al._(2018)Taloni, Vodret,
Costantini, and Zapperi], which consists in determining the values of the
parameters $\mu,\varepsilon_{0},V_{0},\omega_{0}$ that maximize the function

$\ln\mathcal{L}=\sum_{i=1}^{N}\ln\rho_{V_{i}}(\sigma_{i};\dot{\varepsilon_{i}};T).$
(14)

$\rho_{V_{i}}(\sigma_{i};\dot{\varepsilon_{i}};T)=-\frac{d\Sigma_{V_{i}}(\sigma_{i};\dot{\varepsilon_{i}};T)}{d\sigma_{i}}$
is the probability that a wire of volume $V_{i}$ fails at stress $\sigma_{i}$,
under the specific experimental conditions $(\dot{\varepsilon_{i}};T)$. Each
volume was calculated as $V_{i}=\pi L_{i}\left(\frac{d}{2}\right)^{2}$ with
$d=0.25\ 10^{-3}m$ for polyester fibers and $d=0.22\ 10^{-3}m$ for polyamide.
Finally, we set the experimental stresses as $\sigma_{i}=\frac{4F_{i}}{\pi
d^{2}}$ and $E=920MPa$ or $E=2.5GPa$, for polyester and polyamide wires
respectively. The parameters ML estimates were $\mu=0.47Gpa$,
$\varepsilon_{0}=0.0176$, $V0=5.28\ 10^{-13}m^{3}$, $\omega_{0}=0.13sec^{-1}$
in case of polyester, and $\mu=1.06Gpa$, $\varepsilon_{0}=0.011$, $V_{0}=1.68\
10^{-13}m^{3}$, $\omega_{0}=0.002sec^{-1}$ in case of polyamide fibers.

In Fig.1 the experimental fracture stresses (symbols) are plotted together
with the average values
$\langle\sigma_{V_{i}}(\dot{\varepsilon_{i}};T)\rangle$(red dashed lines),
defined as

$\langle\sigma_{V_{i}}(\dot{\varepsilon_{i}};T)\rangle=\int_{0}^{\infty}d\sigma\Sigma_{V_{i}}(\sigma;\dot{\varepsilon_{i}},T).$
(15)

Figure 1: Fracture stresses of polyester (A) and polyamide (B) fibers from the
experiments in [Fontana and Palffy-Muhoray(2020)] (Fig.2). The ML parameters
estimates were used to calculate $\langle\sigma\rangle$ (red dashed line).

The second set of experiments are reported in Fig.2 (symbols). We have plotted
the average stress as a function of the rupture time $t_{r}up$, calculated
from Eq.15 by adopting the parameters ML estimates found in the previous
analysis. The wires’ lengths were $L=6.2m$ in case of polyester and $L=3.8m$
for polyamides, while the strain rates could only be estimated indirectly by
$1/t_{rup}^{i}$, as also acknowledged by the authors [Fontana()]. Our
theoretical prediction would indicate a very small strain dependence of the
rupture forces in agreement with the experiments on polyester (Fig.2 (A)). On
the other side the experiments on polyamide indicate that the exhibited strain
dependence is not captured by the theoretical average stress. A possible
reason for this discrepancy is that the precise value of the strain-rate is
not known and its estimate $1/t_{rup}$ is not correct.

Figure 2: Breaking stresses as a function of the rupture times for polyester
(A) and polyamide (B) wires (Fig.3 in [Fontana and Palffy-Muhoray(2020)]).
Dashed blue lines are $\langle\sigma\rangle$, evaluated using the same
parameters values fitted in the Fig.1, and using $V=0.3\ 10^{-6}m^{3}$
(polyester), $V=0.14\ 10^{-6}m^{3}$ (polyamide) and
$\dot{\varepsilon}_{i}=1/t^{i}_{rup}$.

##  3 Conclusions

We have shown how classical fracture mechanics, coupled to extreme value
theory, provides an adequate theoretical framework to explain the results of
tensile experiments performed on polyamide as well as on polyester wires
reported in [Fontana and Palffy-Muhoray(2020)]. We have shown that the Gumbel
distribution reproduces the logarithmic decay of the average failure stress as
a function of the wire length. We have used the generalized EVT to fit the
experimental failure stresses, showing no significant strain-rate dependence.

## References

  * [Griffith(1921)] A. A. Griffith, Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character 221, 163 (1921).
  * [Taloni  _et al._(2018)Taloni, Vodret, Costantini, and Zapperi] A. Taloni, M. Vodret, G. Costantini, and S. Zapperi, Size effects on the fracture of microscale and nanoscale materials, Nature Reviews Materials 3, 211 (2018).
  * [Ball(2020)] P. Ball, Classic fail, Nature Materials 19, 829 (2020).
  * [Fontana and Palffy-Muhoray(2020)] J. Fontana and P. Palffy-Muhoray, St. petersburg paradox and failure probability, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 124,“ “bibinfo pages 245501 (“bibinfo year 2020).
  * [Duxbury  _et al._(1987)Duxbury, Leath, and Beale] P. Duxbury, P. Leath, and P. D. Beale, Breakdown properties of quenched random systems: the random-fuse network, Physical Review B 36, 367 (1987).
  * [Fisher and Tippett(1928)] R. A. Fisher and L. H. C. Tippett, Limiting forms of the frequency distribution of the largest or smallest member of a sample, in  _Mathematical Proceedings of the Cambridge Philosophical Society_ , Vol. 24 (Cambridge University Press, 1928) pp. 180–190.
  * [Gnedenko(1943)] B. Gnedenko, Sur la distribution limite du terme maximum d’une serie aleatoire, Annals of mathematics , 423 (1943).
  * [Leadbetter  _et al._(2012)Leadbetter, Lindgren, and Rootzén] M. R. Leadbetter, G. Lindgren, and H. Rootzén,  _Extremes and related properties of random sequences and processes_ (Springer Science & Business Media, 2012).
  * [Manzato  _et al._(2012)Manzato, Shekhawat, Nukala, Alava, Sethna, and Zapperi] C. Manzato, A. Shekhawat, P. K. Nukala, M. J. Alava, J. P. Sethna, and S. Zapperi, Fracture strength of disordered media: Universality, interactions, and tail asymptotics, Physical review letters 108, 065504 (2012).
  * [Kunz and Souillard(1978)] H. Kunz and B. Souillard, Essential singularity in the percolation model, Physical Review Letters 40, 133 (1978).
  * [Gumbel(2004)] E. J. Gumbel,  _Statistics of extremes_ (Courier Corporation, 2004).
  * [Sellerio  _et al._(2015)Sellerio, Taloni, and Zapperi] A. L. Sellerio, A. Taloni, and S. Zapperi, Fracture size effects in nanoscale materials: the case of graphene, Physical Review Applied 4, 024011 (2015).
  * [Fontana()] J. Fontana, personal communication.