Pfaffian formula for fermion parity fluctuations in a superconductor and application to Majorana fusion detection

TL;DR

This paper generalizes Kitaev's Pfaffian formula to include quantum fluctuations in fermion parity for open systems, providing a statistical framework to identify Majorana zero-modes in superconductor experiments.

Contribution

It extends the Pfaffian formula to open systems using full counting statistics, enabling better detection of Majorana zero-modes amidst quantum fluctuations.

Findings

Generalized Pfaffian formula for open systems

Statistical criteria for Majorana zero-mode detection

Application to superconductor fusion experiments

Abstract

Kitaev's Pfaffian formula equates the ground-state fermion parity of a closed system to the sign of the Pfaffian of the Hamiltonian in the Majorana basis. Using Klich's theory of full counting statistics for paired fermions we generalize the Pfaffian formula to account for quantum fluctuations in the fermion parity of an open subsystem. A statistical description in the framework of random-matrix theory is used to answer the question when a vanishing fermion parity in a superconductor fusion experiment becomes a distinctive signature of an isolated Majorana zero-mode.

Equations

$$H=$$

$$+\frac{1}{2}{\sum }_{n,m=1}^N\left({\Delta }_{nm}{a}_n{a}_m+{\Delta }_{nm}^{\ast }{a}_m^{†}{a}_n^{†}\right).$$

$$H=\frac{1}{2}\sum _{n,m=1}^N{\Psi }_n^{†}\cdot {\mathcal{B}}_{nm}\cdot {\Psi }_m,$$

$${\Psi }_n=\left(\begin{array}{c}{a}_n\\ a_n^{†}\end{array}\right),{\mathcal{B}}_{nm}=\left(\begin{array}{cc}{V}_{nm}& -{\Delta }_{nm}^{\ast }\\ {\Delta }_{nm}& -V_{nm}^{\ast }\end{array}\right).$$

$${\mathcal{B}}_{nm}\mapsto U{\mathcal{B}}_{nm}{U}^{†}\equiv {\mathcal{A}}_{nm},\text{with}U=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ -i& i\end{array}\right),$$

$${\mathcal{A}}_{nm}=\left(\begin{array}{cc}i\mathrm{Im}(V_{nm}+{\Delta }_{nm})& i\mathrm{Re}({\Delta }_{nm}+V_{nm})\\ i\mathrm{Re}({\Delta }_{nm}-V_{nm})& i\mathrm{Im}(V_{nm}-{\Delta }_{nm})\end{array}\right)=-{\mathcal{A}}_{mn}^{\mathrm{T}}.$$

$$\gamma =({\gamma }_1,{\gamma }_2,\dots {\gamma }_{2N}),\text{with}\left(\begin{array}{c}{\gamma }_{2n-1}\\ {\gamma }_{2n}\end{array}\right)=U\left(\begin{array}{c}{a}_n\\ a_n^{†}\end{array}\right),$$

$${\gamma }_n{\gamma }_m+{\gamma }_m{\gamma }_n={\delta }_{nm},{\gamma }_n^2=1/2.$$

$$\mathcal{P}=(-1{)}^{{\sum }_{n=1}^N{a}_n^{†}{a}_n^{}}=(-2i{)}^N{\gamma }_1{\gamma }_2\cdots {\gamma }_{2N}$$

$${\mathcal{P}}_0=\text{sign}\text{Pf}(-i\mathcal{A}),\text{for}H=\frac{1}{2}\gamma \cdot \mathcal{A}\cdot \gamma .$$

$$⟨{\mathcal{P}}_{\mathrm{L}}⟩=⟨{\Psi }_0|(-1{)}^{{\sum }_{n=1}^{N_{\mathrm{L}}}{a}_n^{†}{a}_n^{}}|{\Psi }_0⟩.$$

$${\mathcal{P}}_0=\text{Pf}(-i\bar{\mathcal{A}}).$$

$$⟨{\mathcal{P}}_{\mathrm{L}}⟩=\text{Pf}[-i\bar{\mathcal{A}}{]}_{\mathrm{LL}}.$$

$$⟨{\gamma }_1{\gamma }_2\cdots {\gamma }_{2s}⟩=\underset{1\le k<l\le 2s}{\mathrm{Pf}}⟨{\gamma }_k{\gamma }_l⟩.$$

$$\mathcal{A}=\left(\begin{array}{cc}0& i\Gamma \\ -i{\Gamma }^{\mathrm{T}}& 0\end{array}\right).$$

$${\mathcal{P}}_0=\mathrm{sign}\mathrm{Pf}(-i\mathcal{A})=-\mathrm{sign}\mathrm{Det}\Gamma .$$

$$\bar{\mathcal{A}}=\left(\begin{array}{cc}0& i{O}_1{O}_2\\ -i{O}_2^{\mathrm{T}}{O}_1^{\mathrm{T}}& 0\end{array}\right)\Rightarrow [\bar{\mathcal{A}}{]}_{\mathrm{LL}}=0\Rightarrow ⟨{\mathcal{P}}_{\mathrm{L}}⟩=0,$$

$$\mathcal{A}=\left(\begin{array}{cc}i\Omega & i\Gamma \\ -i{\Gamma }^{\mathrm{T}}& i\Omega \end{array}\right),\Omega =\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& -1& 0\end{array}\right),\Gamma =\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right),$$

$$\Rightarrow \bar{\mathcal{A}}=\frac{1}{\sqrt{5}}\left(\begin{array}{cc}i{\Omega }^{\prime }& i{\Gamma }^{\prime }\\ -i{\Gamma }^{\prime \mathrm{T}}& i{\Omega }^{\prime }\end{array}\right),{\Gamma }^{\prime }=2\Gamma ,{\Omega }^{\prime }=\left(\begin{array}{cccc}0& -1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& -1& 0\end{array}\right),$$

$$\Rightarrow [\bar{\mathcal{A}}{]}_{\mathrm{LL}}=\frac{i}{\sqrt{5}}{\Omega }^{\prime }\Rightarrow \mathrm{Pf}[-i\bar{\mathcal{A}}{]}_{\mathrm{LL}}=-\frac{1}{5},$$

$$\displaystyle|\Psi_{0}\rangle=\frac{\sqrt{5}}{10}\bigl{[}2i(a_{1}^{\dagger}a_{2}^{\dagger}+a_{3}^{\dagger}a_{4}^{\dagger})$$

$$\displaystyle-(1-\sqrt{5})a_{2}^{\dagger}a_{4}^{\dagger}\bigr{]}|0\rangle,$$

$$\mathcal{S}=O\left(\begin{array}{cc}{0}_{N\times N}& 1_{N\times N}\\ -1_{N\times N}& 0_{N\times N}\end{array}\right)O^{\mathrm{T}}.$$

$$\mathcal{S}=\left(\begin{array}{cc}{\mathcal{S}}_{\mathrm{LL}}& {\mathcal{S}}_{\mathrm{LR}}\\ {\mathcal{S}}_{\mathrm{RL}}& {\mathcal{S}}_{\mathrm{RR}}\end{array}\right),{\mathcal{S}}_{\mathrm{LL}}=[-i\bar{\mathcal{A}}{]}_{\mathrm{LL}},$$

$$P({\lambda }_1,{\lambda }_2,\dots {\lambda }_{N_{\min }})\propto \prod _n(1-{\lambda }_n^2{)}^{|N_{\mathrm{L}}-N_{\mathrm{R}}|}\prod _{j<k}({\lambda }_k^2-{\lambda }_j^2{)}^2,$$

$$N_{\min }=\min (N_{\mathrm{L}},N_{\mathrm{R}}),0\le {\lambda }_n\le 1.$$

$${\lambda }_c=(2/N)(N_{\mathrm{L}}{N}_{\mathrm{R}}{)}^{1/2},$$

$$\rho (\lambda )=\frac{1}{{\delta }_{\mathrm{eff}}}+\frac{\sin (2\pi \lambda /{\delta }_{\mathrm{eff}})}{2\pi \lambda },\lambda \lesssim 1/{\delta }_{\mathrm{eff}}.$$

$$|⟨{\mathcal{P}}_{\mathrm{L}}⟩|=\prod _{n=1}^{N_{\min }}{\lambda }_n=\sqrt{\mathrm{Det}{\mathcal{S}}_{\mathrm{LL}}}.$$

$$P(⟨{\mathcal{P}}_{\mathrm{L}}⟩=p)=\frac{45}{32}\left(1-p^4+4p^2\ln |p|\right),N_{\mathrm{dot}}=2,|p|\le 1.$$

$$\mathrm{Var}⟨{\mathcal{P}}_{\mathrm{L}}⟩=\frac{(2N_{\mathrm{dot}}){!}^3}{(N_{\mathrm{dot}}){!}^2(4N_{\mathrm{dot}})!}=\frac{\sqrt{2}}{4^{N_{\mathrm{dot}}}}[1+\mathcal{O}(1/N_{\mathrm{dot}})],$$

$$P_{\mathrm{uniform}}(⟨{\mathcal{P}}_{\mathrm{L}}⟩=p)=\frac{(-\ln |p|{)}^{N_{\mathrm{dot}}-1}}{2(N_{\mathrm{dot}}-1)!},|p|\le 1,$$

$$\mathcal{A}=\left(\begin{array}{cc}i{\Omega }_{\mathrm{L}}& i\Gamma \\ -i{\Gamma }^{\mathrm{T}}& i{\Omega }_{\mathrm{R}}\end{array}\right).$$

Full text

\shortabstract

1 Introduction

The pairing interaction in a superconductor favors a ground state with an even number of electrons, but when both time-reversal and spin-rotation symmetry are broken the ground state may have odd parity — for example when a magnetic impurity binds an unpaired electron [1]. While the connection between fermion-parity switches and level crossings was noticed already in 1970 by Sakurai [2], these only became a topic of intense research activity after Kitaev [3] made the connection with topological phase transitions and Majorana fermions: The absence of level repulsion at a fermion-parity switch indicates a change in a topological quantum number, which Kitaev identified as the sign of the Pfaffian of the Hamiltonian in the basis of Majorana fermions.

An open subsystem need not be in a state of definite fermion parity ${\cal P}=\pm 1$, the fermion parity expectation value $\langle{\cal P}\rangle$ may take on any value in the interval $[-1,1]$. Here we generalize Kitaev’s Pfaffian formula so that it can describe both closed and open systems. This generalization has a computational as well as a conceptual merit. Computationally, it reduces the complexity of a calculation of $\langle{\cal P}\rangle$ for $N$ levels from order $2^{N}$, when all possible occupation numbers are enumerated, down to order $N^{3}$. Conceptually, it allows us to make contact with the random-matrix theory of topological superconductivity [4, 5], and identify the origin of a statistical peak at $\langle{\cal P}\rangle=0$ discovered recently in computer simulations [6]. These findings have implications for proposed experiments [7] to search for signatures of isolated Majorana zero-modes in the fermion parity of two superconductors that have first been fused and then decoupled (see Fig. 1).

The outline of the paper is as follows. In the next section we derive the Pfaffian formula for the average subsystem fermion parity. This generalization of Kitaev’s formula [3] can be seen either as an application of the Wick theorem for Majorana operators [8, 9, 10] (cf. a similar application in Ref. 11), or as an application of Klich’s theory of counting statistics for paired fermions [12]. In Sec. 3 we use the fermion parity formula to establish the connection between vanishing average fermion parity and the presence of isolated Majorana zero-modes in the decoupled quantum dot. We continue in Sec. 4 with a statistical description of the double quantum dot geometry of Fig. 1, by identifying the random-matrix ensemble in symmetry class DIII that describes the fermion parity fluctuations. We contrast the case of strongly coupled quantum dots in Sec. 4.2 with the case of weak coupling in Sec. 4.3. In Sec. 5 we show how weak coupling by a single-mode quantum point contact can distinguish quantum dots with or without isolated Majorana zero-modes. In the concluding Sec. 6 we discuss the implications of our analysis for the detection of Majorana zero-modes by means of a fusion experiment.

2 Pfaffian fermion-parity formula

2.1 Kitaev’s formula for an isolated system

To set the stage we recall some basic facts [13] needed to present Kitaev’s formula [3] for the ground-state fermion parity of an isolated superconductor.

At the mean-field level the Hamiltonian of a superconductor is a Hermitian quadratic form in the fermion creation and annihilation operators $a^{\dagger}$, $a$,

[TABLE]

The indices $n,m$ label spin and orbital degrees of freedom of $N$ fermionic modes. The $N\times N$ Hermitian matrix $V$ represents the kinetic and potential energy and the antisymmetric matrix $\Delta$ is the pair potential.

More compactly, Eq. (1) can be written in the matrix form

[TABLE]

The $2N\times 2N$ Hermitian matrix ${\cal B}$ is called the Bogoliubov-De Gennes (BdG) Hamiltonian [14]. Its eigenvalues come in pairs $\pm E_{1},\pm E_{2},\ldots\pm E_{N}$ of opposite sign, with the positive entries equal to the single-particle excitation energies of the many-particle Hamiltonian $H$.

The unitary transformation

[TABLE]

maps ${\cal B}$ onto the $2N\times 2N$ imaginary antisymmetric matrix ${\cal A}$ with elements

[TABLE]

The superscript T denotes the transpose. An antisymmetric matrix is also referred to as “skew-symmetric”.

The transformed state

[TABLE]

contains $2N$ Hermitian operators $\gamma_{n}=\gamma_{n}^{\dagger}$, with anticommutator

[TABLE]

This is the Clifford algebra of Majorana operators.

The global fermion parity operator

[TABLE]

commutes with $H$, so energy eigenstates have a definite fermion parity $\pm 1$. Kitaev’s formula [3] equates the fermion parity ${\cal P}_{0}$ of the ground state to the Pfaffian111Wikipedia has a helpful collection of Pfaffian formulas. (Pf) of the Hamiltonian in the Majorana basis,

[TABLE]

2.2 Pfaffian formula for a subsystem

Our objective is to calculate the ground-state expectation value of the fermion parity ${\cal P}_{\rm L}$ of an open subsystem, say the left quantum dot with $N_{\text{L}}$ fermionic modes in Fig. 1 .

A direct way to proceed, used for example in Ref. 6, is to calculate the many-particle ground state $|\Psi_{0}\rangle$ in the basis of occupation numbers and evaluate

[TABLE]

Since the Fock space of occupation numbers has dimension $2^{N}$, this direct approach scales exponentially with system size and is therefore prohibitively expensive for large systems.

Klich [12] has developed an efficient method, with a polynomial scaling in $N$, to calculate squares of expectation values of operators $\exp(i\chi\sum_{n}a^{\dagger}_{n}a^{\vphantom{\dagger}}_{n})$. This gives $\langle{\cal P}_{\rm L}\rangle^{2}$ if one sets $\chi=\pi$ and restricts the sum to indices $n$ in L. In App. A we show how the Klich method can be adapted to give also the sign of $\langle{\cal P}\rangle_{\rm L}$. That calculation is technically rather involved, but the final result can be easily understood as follows.

We make the flat-band transformation ${\cal A}\mapsto\bar{\cal A}$, which consists in replacing each of the $2N$ eigenvalues $\pm E_{n}$ of ${\cal A}$ by their sign. (We assume that no eigenvalue is identically zero, meaning that we are not precisely at a fermion-parity switch.) Since no eigenvalue crosses zero when it is replaced by its sign, the flat-band transformation leaves the sign of the Pfaffian (8) invariant. And because the Pfaffian of $-i\bar{\cal A}$ can only equal $\pm 1$ we no longer need to take the sign in Eq. (8), hence the global fermion parity is

[TABLE]

At this point one may surmise that the desired subsystem generalization of Eq. (8) simply amounts to taking the Pfaffian of the $2N_{\rm L}\times 2N_{\rm L}$ submatrix $[\bar{\cal A}]_{\rm LL}$ restricted to the subspace of modes in the left quantum dot,

[TABLE]

This is indeed the correct expression, as one can see by application of the Wick theorem for Majorana operators [8, 9, 10],

[TABLE]

Substitution of ${\cal P}_{\rm L}=(-2i)^{N_{\rm L}}\gamma_{1}\gamma_{2}\cdots\gamma_{2N_{\rm L}}$ on the left-hand-side and $-2i\langle\gamma_{k}\gamma_{l}\rangle=-i\bar{\cal A}_{kl}$ on the right-hand-side results in Eq. (11). This is how an equivalent formula was derived recently for a different problem [11].

Eq. (11) is computationally efficient because the Pfaffian of an $N\times N$ matrix can be calculated in a time that scales polynomially with $N$ [15, 16]: It has the same ${\cal O}(N^{3})$ complexity as the eigenvalue decomposition one needs for the flat-band transformation ${\cal A}\mapsto\bar{\cal A}$. Note that the flat-band transformation needs to be performed before the subblock restriction $\bar{\cal A}\mapsto[\bar{\cal A}]_{\rm{LL}}$ — the two operations do not commute.

3 Connection with the Majorana fusion rule

As a fundamental application of Eq. (11), consider the case that each quantum dot in Fig. 1 has a single electronic mode ($N_{\rm L}=N_{\rm R}=1$), each consisting of two Majorana modes with inter-dot coupling matrix $\Gamma$ but vanishing intra-dot coupling — so these become fully isolated zero-modes when the quantum dots are decoupled. The Hamiltonian in the Majorana basis is

[TABLE]

The global fermion parity is

[TABLE]

To obtain the average local fermion parity we use that the real $2\times 2$ coupling matrix $\Gamma$ has the singular value decomposition $\Gamma=O_{1}\,{\rm diag}\,(\kappa_{1},\kappa_{2})O_{2}$, with $O_{1},O_{2}$ real orthogonal matrices and $\kappa_{1},\kappa_{2}>0$. The eigenvalues of ${\cal A}$ are $\pm\kappa_{1},\pm\kappa_{2}$. In the flat-band transformation $\{\kappa_{1},\kappa_{2}\}\mapsto\{1,1\}$, which gives

[TABLE]

so the average fermion parity in a single quantum dot vanishes. This is a manifestation of the Majorana fusion rule [17]: The fusion of the two Majorana zero-modes $\gamma_{1}$ and $\gamma_{2}$ produces an equal-weight superposition of a state of even and odd fermion parity.222The converse is not excluded: $\langle{\cal P}_{\rm L}\rangle=0$ without an isolated Majorana zero-mode is possible, for example for

${\cal A}=i\tiny{\begin{pmatrix}0&\lambda_{1}&0&\lambda_{2}\\ -\lambda_{1}&0&-\lambda_{2}&0\\ 0&\lambda_{2}&0&\lambda_{1}\\ -\lambda_{2}&0&-\lambda_{1}&0\end{pmatrix}}$ with $\lambda_{1}<\lambda_{2}$.

Several recent experimental proposals [6, 7, 18] are based on the connection between the Majorana fusion rule and vanishing average fermion parity. The implication “isolated Majorana zero-modes $\Rightarrow$ $\langle{\cal P}_{\rm L}\rangle=0$” holds if there are only two pairs of Majorana zero-modes. For $N_{\rm L}$ or $N_{\rm R}$ greater than 1 the implication breaks down, as is demonstrated by the following counterexample for $N_{\rm L}=N_{\rm R}=2$:

[TABLE]

and hence $\langle{\cal P}_{\rm L}\rangle=-1/5$ does not vanish even though each quantum dot has a pair of Majorana zero-modes without intra-dot coupling ($\gamma_{1}$ and $\gamma_{2}$ in the left dot, $\gamma_{5}$ and $\gamma_{6}$ in the right dot).

Since ${\rm Pf}\,(-i{\cal A})=+1$ the global fermion parity is even, hence the negative sign for $\langle{\cal P}_{\rm L}\rangle$ means that the states with odd-odd occupation numbers in the left and right quantum dot have a greater weight in the ground state than the states with even-even occupation numbers — even though the fusion of the Majorana modes $\gamma_{1}$ and $\gamma_{2}$ would favor equal weight of even and odd fermion parity.

As a check on the formalism, we have also calculated the average fermion parity directly from the many-particle ground state wave function $|\Psi_{0}\rangle$ of the Hamiltonian $H=\tfrac{1}{2}\gamma\cdot{\cal A}\cdot\gamma$. We find

[TABLE]

which indeed gives $\langle{\cal P}_{\rm L}\rangle=-1/5$ upon calculation of the expectation value (9).

In this case with $N=N_{\rm L}+N_{\rm R}=4$ electronic levels the size $2^{N-1}=8$ of the basis of many-particle states in the even-parity sector is the same as the size $2N=8$ of the basis of single-particle states, so the two calculations based on Eq. (9) or on Eq. (11) are equally efficient. For larger $N$ the single-particle approach based on the Pfaffian formula has the more favorable scaling (polynomial instead of exponential).

4 Random-matrix theory

For a statistical description of the fermion parity fluctuations we apply the methods of random-matrix theory (RMT). In Sec. 4.2 we assume a strong mixing of the states in the two quantum dots of Fig. 1, and then in Sec. 4.3 we consider the opposite regime of weakly coupled quantum dots. We will need results [19] from the RMT in symmetry class DIII, which we summarize in Sec. 4.1.

4.1 Skew Circular Real Ensemble

The matrix $[-i\bar{\cal A}]_{\rm LL}$ which in view of Eq. (11) determines the local fermion parity is a $2N_{\rm L}\times 2N_{\rm L}$ submatrix of a matrix ${\cal S}=-i\bar{\cal A}$ that is an antisymmetric (skew-symmetric) element of the real orthogonal group ${\rm O}(2N)$, with $N=N_{\rm L}+N_{\rm R}$. The corresponding ensemble from RMT is the class-DIII circular ensemble, which differs from the class-D circular ensemble by the antisymmetry restriction [5]. The latter is called the Circular Real Ensemble (CRE) and we will refer to the former as the skew-Circular Real Ensemble (skew-CRE).333The qualifier “real” for the ${\rm O}(N)$ ensemble is used instead of “orthogonal” because the name Circular Orthogonal Ensemble (COE) was already used by Dyson [20] for the coset ${\rm U}(N)/{\rm O}(N)$. The switch from symmetry class D to DIII is remarkable, because class DIII was originally introduced [4] in superconductors with preserved time-reversal symmetry — which is broken in our physical system.

Two equivalent methods to randomly choose a matrix from the skew-CRE are:

Generate a real antisymmetric matrix $-i{\cal A}$ with independent Gaussian elements on the upper diagonal (zero mean and unit variance), and perform the flat-band transformation to obtain ${\cal S}=-i\bar{\cal A}$. 2. 2.

Draw a random element $O$ from ${\rm O}(2N)$, uniformly with the invariant Haar measure, and construct

[TABLE]

The two methods are equivalent because the distribution $P({\cal A})\propto\exp(\tfrac{1}{4}{\rm Tr}\,{\cal A}^{2})$ as well as the flat-band transformation ${\cal A}\mapsto\bar{\cal A}$ are invariant under orthogonal transformations ${\cal A}\mapsto O{\cal A}O^{\rm T}$, so the matrix $O$ in the decomposition (18) is distributed according to the invariant Haar measure.

The matrix ${\cal S}$ has the block decomposition

[TABLE]

with ${\cal S}_{XY}$ a matrix of dimension $N_{X}\times N_{Y}$. In the context of scattering problems, where the skew-CRE ensemble was studied previously [5], this is analogous to a decomposition of the scattering matrix into reflection and transmission matrices. In that context the eigenvalues $\pm i\lambda_{n}$ of the upper-left submatrix ${\cal S}_{\rm LL}$ correspond to reflection amplitudes.444Eq. (20) follows from equation 5 of Ref. 19 upon change of variables from transmission probabilities $T_{n}$ to reflection amplitudes $\lambda_{n}=\sqrt{1-T_{n}}$. Their joint probability distribution in the skew-CRE is known [19],

[TABLE]

If $N_{\rm L}>N_{\rm R}$ there are additionally $2(N_{\rm L}-N_{\rm R})$ trivial eigenvalues pinned at $\pm 1$, not included in the distribution (20).

Symmetry class DIII has the $\mathbb{Z}_{2}$ invariant ${\rm Pf}\,{\cal S}=\pm 1$, which in view of Kitaev’s formula (10) is the global fermion parity ${\cal P}_{0}$. This does not enter in Eq. (20) because in the skew-CRE the distribution of the $\lambda_{n}$’s is independent of the $\mathbb{Z}_{2}$ invariant [19].

The density $\rho(\lambda)$ of the nontrivial eigenvalues has $\pm\lambda$ symmetry with a three-peak structure: There are two peaks at the band edges $\pm\lambda_{c}$, with [19]

[TABLE]

and a peak at the band center555For the density profile near $\lambda=0$ we can approximate the distribution (20) by $P(\{\lambda\})\propto\prod_{j<k}(\lambda_{k}^{2}-\lambda_{j}^{2})^{2}$ and ignore the restriction $|\lambda_{n}|\leq 1$. The distribution of the $\lambda_{n}$’s is then identical to the distribution of the energy levels of a Hermitian matrix in symmetry class D, which has the spectral peak (22). A Hermitian matrix in class DIII, rather than class D, has a vanishing density of states at zero energy, but this is not relevant for $\rho(\lambda)$. described by [4, 21, 22]

[TABLE]

The parameter $\delta_{\rm eff}=\pi/2N_{\rm min}$ is the mean eigenvalue spacing in the center of the band. The peak at $\lambda=0$ is a weak antilocalization effect in the scattering context [23].

Fig. 2 shows the eigenvalue density for $N_{\rm L}=N_{\rm R}=N_{\rm dot}$ ranging from 1 to 6. The three-peaked structure is evident except for $N_{\rm dot}=1$, when the density profile is flat.

4.2 Distribution of the local fermion parity in the skew-CRE

The peak at $\lambda=0$ in the eigenvalue density $\rho(\lambda)$ increases the probability for vanishing local fermion parity, since

[TABLE]

Indeed, as shown in Fig. 3, while the distribution of $\langle{\cal P}_{\rm L}\rangle$ in the skew-CRE is broad for a single electronic level $N_{\rm dot}=1$ in each quantum dot, it quickly narrows to a sharp peak at $\langle{\cal P}_{\rm L}\rangle=0$ with just a few levels — in accord with numerical calculations reported by Clarke, Sau, and Das Sarma [6].

The peak at zero $\langle{\cal P}_{\rm L}\rangle\equiv p$ appears as a sharp cusp in Fig. 3, it has a logarithmic singularity $\propto(p^{2}\ln|p|)^{N_{\rm dot}-1}$, for example

[TABLE]

For large-$N_{\rm dot}$ the width of the distribution becomes exponentially small, as follows from the variance

[TABLE]

see App. B.

We may quantify the effect of the spectral peak in $\rho(\lambda)$ on the distribution of the local fermion parity by comparing with a set of independent $\lambda_{n}$’s with uniform density. In that uniform case one would have the fermion parity distribution

[TABLE]

with a variance $3^{-N_{\rm dot}}$ that decays less rapidly than Eq. (25).

4.3 RMT model of weakly coupled quantum dots

The RMT description in terms of the skew-CRE from the previous subsection assumes a strong (chaotic) mixing in the entire phase space, appropriate for strongly coupled quantum dots. To describe also the weakly coupled regime, we consider an alternative approach where the RMT ensemble is applied to the two quantum dots individually, rather than to the system as a whole.

In the Majorana representation, the Hamiltonian $H=\tfrac{1}{2}\gamma\cdot{\cal A}\cdot\gamma$ of the two coupled quantum dots of Fig. 1 has the block structure

[TABLE]

The real antisymmetric matrices $\Omega_{X}$ of size $2N_{X}\times 2N_{X}$, with $X\in\{\text{R,L}\}$, describe the left and right quantum dot in isolation, while the $2N_{\rm L}\times 2N_{\rm R}$ real matrix $\Gamma$ describes the coupling via a quantum point contact (QPC) with $N_{\rm QPC}$ propagating fermionic modes. In what follows we take $N_{\rm L}=N_{\rm R}=N_{\rm dot}$.

The number $N_{\rm dot}$ counts the number of electronic modes in each quantum dot. One electronic mode $a_{n}$ corresponds to two Majorana modes $\gamma_{2n-1}$ and $\gamma_{2n}$, according to

[TABLE]

cf. Eq. (5). Because of this double-counting, the mean level spacing $\delta_{0}$ of eigenstates of $\Omega_{X}$ is one half the electronic mean level spacing of a quantum dot (taken the same in each dot, for simplicity).

For a statistical description we take independent Gaussian distributions for the two matrices $\Omega_{X}$. Each upper-diagonal matrix element has zero mean and variance $2N_{\rm dot}\delta_{0}^{2}/\pi^{2}$, corresponding to superconductors in symmetry class D (broken time-reversal and broken spin-rotation symmetry) [4, 5].

Following Refs. 24, 25, the quantum dots are coupled by a Gaussian random matrix $\Gamma$ of rank $N_{\rm QPC}$, with elements666The coupling matrix (29) describes a ballistic point contact. For tunnel coupling, rather than ballistic coupling, the coupling strength $\delta_{0}/\pi$ is to be multiplied by $T_{n}^{-1}(2-T_{n}-2\sqrt{1-T_{n}})$, with $T_{n}$ the tunnel probability of the $n$-th mode in the QPC, see Ref. 5.

[TABLE]

in terms of $2N_{\rm QPC}$ real Gaussian vectors $v^{(n)}$ and $w^{(n)}$ of unit average length (each element independently distributed with zero mean and variance $1/2N_{\rm dot}$).

For the weak coupling regime we focus on the case of a single propagating electronic mode in the point contact, $N_{\rm QPC}=1$, corresponding to two propagating Majorana modes. We do not have an analytical solution, so we show numerical results in Fig. 4 for the probability distribution of $\langle{\cal P}_{\rm L}\rangle={\rm Pf}\,(-i\bar{\cal A})$ in the ensemble of random matrices $\Omega_{\rm L}$, $\Omega_{\rm R}$, and $\Gamma$. The variance of the distribution is compared with that in the skew-CRE in Fig. 5. The two figures show that the distribution of the local fermion parity is much broader when the coupling is via a single-mode point contact.

5 Effect of an isolated Majorana zero-mode

The random Hamiltonians of the previous section do not contain isolated Majorana zero-modes: the $2N_{\rm dot}$ Majorana modes in each quantum dot have intradot coupling as well as interdot coupling. We may introduce a pair of isolated Majorana zero-modes in a quantum dot by setting to zero one row and one column of the submatrix $\Omega_{\rm L}$ or $\Omega_{\rm R}$ in the Hamiltonian (27). (The row and column number should be the same to preserve the antisymmetry of $\Omega_{X}$.) The effect on the distribution of the local fermion parity is shown in Fig. 6. The distribution of the local fermion parity is strongly peaked at zero if and only if there is a pair of isolated Majorana zero-modes in each of the two quantum dots.

6 Conclusion

In summary, we have studied the fusion of Majorana zero-modes using a generalization of Kitaev’s Pfaffian formula [3] for the global fermion parity of the superconducting ground state, to include local fermion parity fluctuations in an open subsystem. The Pfaffian formula in Eq. (11), and an equivalent formulation from Ref. 11, is computationally efficient since it works with the single-particle (Bogoliubov-De Gennes) Hamiltonian rather than with the many-particle Hamiltonian. One limitation of the single-particle formulation is that it is limited to a mean-field description of the superconductor — in particular we are assuming that the quantum dots in the geometry of Fig. 1 have a sufficiently large capacitance that Coulomb charging energies can be neglected.

The Pfaffian fermion parity formula is particularly suited to an analysis in terms of random-matrix theory, in an ensemble of antisymmetric matrices [5]. For strongly coupled quantum dots the circular ensemble in symmetry class DIII is the appropriate ensemble, which allows for analytical results for the statistical distribution of the local fermion parity. There is no eigenvalue repulsion at the particle-hole symmetry point in such an ensemble [4], and the resulting accumulation of near-zero eigenvalues enforces a nearly equal-weight superposition of even and odd fermion parity in a quantum dot.

This is a nontopological mechanism for vanishing expectation value $\langle{\cal P}_{\rm L}\rangle\approx 0$ of the local fermion parity. The Majorana fusion rule provides a fundamentally different, topological mechanism [17]: The merging or “fusion” of two isolated Majorana zero-modes (“isolated” in the sense of zero intradot coupling, while allowing for interdot coupling) also favors a vanishing $\langle{\cal P}_{\rm L}\rangle$ because the two fusion channels, with or without an unpaired quasiparticle, have equal weight.

To carry out such a fusion experiment it is proposed [7] that one would rapidly decouple the subsystems, on a time scale $\tau_{c}$ sufficiently short that quasiparticles from the environment cannot leak in. The complication [6] is that even if there are isolated Majorana zero-modes, the presence of even a small number $N_{\rm dot}$ of higher levels at energies below $\hbar/\tau_{c}$ may hide the presence of the zero-modes by favoring $\langle{\cal P}_{\rm L}\rangle\approx 0$ (see Fig. 3).

Fig. 6 illustrates our proposal to distinguish the two mechanisms for vanishing local fermion parity: A low-rank coupling between the quantum dots, via a single-mode quantum point contact, suppresses the nontopological effect from levels at nonzero energy, without affecting the topological effect from the fusion of isolated Majorana zero-modes.

Acknowledgements.

We have benefited from discussions with A. R. Akhmerov and Y. Herasymenko. J. H. H. Perk alerted us to the connection between Eqs. (11) and (12). This research was supported by the Netherlands Organization for Scientific Research (NWO/OCW) and by the European Research Council.

Appendix A Derivation of the Pfaffian formula from Klich’s counting statistics theory

We follow the steps of Klich’s theory of counting statistics of paired fermions [12], to reproduce his result for $\langle{\cal P}_{\rm L}\rangle^{2}$. Then we will resolve the sign ambiguity to arrive at Eq. (11) for $\langle{\cal P}_{\rm L}\rangle$. An equivalent formula is obtained by a different method in Ref. 11, Appendix B.

The superconductor in Fig. 1 is assumed to be an isolated system, so that the global fermion parity does not fluctuate. For the derivation of the subsystem fermion parity formula (11) it is convenient to start from the more general case that the superconductor is in contact with a reservoir in thermal equilibrium at temperature $T$. We will then take the $T\rightarrow 0$ limit at the end of the calculation in order to describe an isolated system.

At inverse temperature $\beta=1/k_{{\text{B}}}T$ the average fermion parity $\langle{\cal P}_{\text{L}}\rangle$ of subsystem L (the left quantum dot in Fig. 1) is given by the trace of the equilibrium density matrix

[TABLE]

acting on the fermion parity operator

[TABLE]

Because $H=\tfrac{1}{2}\sum_{n,m}{\cal A}_{nm}\gamma_{n}\gamma_{m}$ in the basis of Majorana operators $\gamma_{n}$, and $a^{\dagger}_{n}a^{\vphantom{\dagger}}_{n}=i\gamma_{2n-1}\gamma_{2n}+\tfrac{1}{2}$, this can be written as

[TABLE]

The matrix $\sigma_{y}$ is a Pauli matrix and the operator $P_{\text{L}}$ projects onto $N_{\text{L}}$ fermionic modes in subsystem L.

Application of the identity [12]

[TABLE]

results in

[TABLE]

In the second equality we made use of the identity

[TABLE]

with $\chi=\pi$. (The matrix $\sigma_{0}=\sigma_{y}^{2}$ is the $2\times 2$ unit matrix.) Note that, in a basis of energy eigenstates of the BdG Hamiltonian, the operator $(1+e^{\beta{\cal A}})^{-1}$ is the Fermi function $f(E)=(1+e^{\beta E})^{-1}$.

Eq. (34) is Klich’s result for the square of the average fermion parity (equation 84 in Ref. 12). Klich shows how the sign of $\langle{\cal P}_{\rm L}\rangle$ can be recovered if the determinant is known analytically as a function of the matrix elements. Here we take a different route, more suitable for numerical calculations, which gives the sign directly upon evaluation of a Pfaffian instead of a determinant.

Any $2N\times 2N$ imaginary anti-symmetric matrix ${\cal A}$ can be decomposed as

[TABLE]

where $O$ is a $2N\times 2N$ real orthogonal matrix and

${\cal E}=\text{diag}\,(E_{1},E_{2},\ldots E_{N})$ is an $N\times N$ real diagonal matrix. Substitution into Eq. (34) gives

[TABLE]

This may be written in a more compact form by defining the restriction $[M]_{\text{LL}}$ of a $2N\times 2N$ matrix $M$ to the $2N_{\text{L}}\times 2N_{\text{L}}$ submatrix of modes in region L,

[TABLE]

Note that, because of the submatrix restriction, the product rule ${\text{Det}}\,(AB)=({\text{Det}}\,A)({\text{Det}}\,B)$ cannot be applied to ${\text{Det}}[AB]_{\text{LL}}$, so the orthogonal matrix $O$ cannot be cancelled with the inverse $O^{\text{T}}$.

We have now arrived at the determinant of a real antisymmetric matrix, hence we can take the square root without introducing branch cuts,

[TABLE]

In the zero-temperature, $\beta\rightarrow\infty$ limit this reduces to

[TABLE]

which is Eq. (11) with $-i\bar{\cal A}=O[J\otimes(\text{sign}\,{\cal E})]O^{\text{T}}$. Kitaev’s formula (8) for the global ground-state fermion parity is recovered when L is the entire isolated system. This correspondence also identifies $\sqrt{\text{Det}}$ with $+\text{Pf}$ rather than with $-\text{Pf}$.

Appendix B Moments of determinants of antisymmetric random matrices

In Sec. 4.2 we used a formula for the average determinant of a submatrix (a principal minor) of an antisymmetric real orthogonal matrix. This would seem like a classic result in RMT, but we have not found it in the literature on such matrices [27, 28, 29]. We therefore give the derivation in this appendix, and for completeness and reference also derive the corresponding result for antisymmetric Hermitian matrices.

B.1 Principal minor of antisymmetric orthogonal matrix

Consider a $2N\times 2N$ antisymmetric real orthogonal matrix ${\cal S}$, with a uniform distribution in ${\rm O}(2N)$ subject to the antisymmetry constraint. This is the class-DIII circular ensemble of RMT [4, 5], referred to as the skew-Circular Real Ensemble (skew-CRE) in the main text.777The antisymmetric orthogonal matrices form a disconnected set in ${\rm O}(2N)$, distinguished by the sign of the Pfaffian. For the probability distribution (41) it does not matter whether or not we restrict the ensemble to ${\rm Pf}\,{\cal S}=\pm 1$.

The $2N_{\rm L}\times 2N_{\rm L}$ upper-left submatrix ${\cal S}_{\rm LL}$ has eigenvalues $\pm i\lambda_{n}$, $0\leq\lambda_{n}\leq 1$. Denoting $N_{\rm R}=N-N_{\rm L}$ and $N_{\rm min}=\min(N_{\rm L},N_{\rm R})$, we have that $N-N_{\rm min}$ of the $\lambda_{n}$’s are pinned to $+1$. The set $\{\lambda_{n}\}=\{\lambda_{1},\lambda_{2},\ldots\lambda_{N_{\rm min}}\}$ can vary freely in the interval $[0,1]$, with joint probability distribution [19]

[TABLE]

The determinant of ${\cal S}_{\rm LL}$ is a principal minor given by

[TABLE]

We seek the moments \mu_{q}=\mathbb{E}\bigl{[}({\rm Det}\,{\cal S}_{\rm LL})^{q}\bigr{]} of this determinant in the skew-CRE.

For that purpose we make a change of variables from $\lambda_{n}$ to $R_{n}=\lambda_{n}^{2}\in[0,1]$, with distribution

[TABLE]

We can then compute the moments of the determinant from

[TABLE]

where we abbreviated $\int_{0}^{1}d\{R_{n}\}=\int_{0}^{1}dR_{1}\cdots\int_{0}^{1}dR_{N_{\rm min}}$.

These socalled Selberg integrals have a closed-form expression [26],

[TABLE]

For the first few moments, Eq. (45) reduces to

[TABLE]

Eq. (25) in the main text is Eq. (46) for $N_{\rm L}=N_{\rm R}=N_{\rm dot}=N/2$.

B.2 Antisymmetric Hermitian matrix

A similar calculation can be carried out for moments of the determinant of a $2N\times 2N$ antisymmetric Hermitian matrix ${\cal A}$, in the Gaussian ensemble of independent upper-diagonal elements with a normal distribution (zero mean and unit variance).

The $2N$ eigenvalues come in pairs $\pm\lambda_{n}$. The $N$ eigenvalues $\lambda_{n}\geq 0$ have the joint distribution [21]

[TABLE]

The determinant is

[TABLE]

Let us introduce the variables $x_{n}=\lambda_{n}^{2}/2\geq 0$, with distribution

[TABLE]

The $q$-th moment $\mu_{q}$ of the determinant of ${\cal A}$ is given by

[TABLE]

with $\int_{0}^{\infty}d\{x_{n}\}=\int_{0}^{\infty}dx_{1}\cdots\int_{0}^{\infty}dx_{N}$. This is the ratio of normalisation constants of Laguerre distributions, which is known [26]. We thus obtain

[TABLE]

For $q=1,2$ this reduces to

[TABLE]

The average determinant of antisymmetric Hermitian matrices increases exponentially with $N$,

[TABLE]

in contrast to the exponential decay for antisymmetric orthogonal matrices, cf. Eq. (25).