Pfaffian formula for fermion parity fluctuations in a superconductor and application to Majorana fusion detection
A. Grabsch Y. Cheipesh C. W. J. Beenakker

Engineering passive swimmers by shaking liquids

    M Laumann^1, A Förtsch^1, E Kanso^2 and W Zimmermann^1
    July 9, 2019
==========================================================





§ 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 <cit.>. While the connection between fermion-parity switches and level crossings was noticed already in 1970 by Sakurai <cit.>, these only became a topic of intense research activity after Kitaev <cit.> 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 P=± 1, the fermion parity expectation value ⟨ P⟩ 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 ⟨ P⟩ 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 <cit.>, and identify the origin of a statistical peak at ⟨ P⟩=0 discovered recently in computer simulations <cit.>. These findings have implications for proposed experiments <cit.> 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. <ref>).



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 <cit.> can be seen either as an application of the Wick theorem for Majorana operators <cit.> (cf. a similar application in Ref. Bau18), or as an application of Klich's theory of counting statistics for paired fermions <cit.>. In Sec. <ref> 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. <ref> with a statistical description of the double quantum dot geometry of Fig. <ref>, 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. <ref> with the case of weak coupling in Sec. <ref>. In Sec. <ref> 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. <ref> we discuss the implications of our analysis for the detection of Majorana zero-modes by means of a fusion experiment.



§ PFAFFIAN FERMION-PARITY FORMULA





 §.§ Kitaev's formula for an isolated system



To set the stage we recall some basic facts <cit.> needed to present Kitaev's formula <cit.> 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^†, a,

    H=   ∑_n,m=1^NV_nm(a^†_n a^_m-12δ_nm)
       +12∑_n,m=1^N(Δ_nma_na_m+Δ_nm^∗ a_m^† a_n^†).

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

More compactly, Eq. (<ref>) can be written in the matrix form



    H=12∑_n,m=1^NΨ_n^†· B_nm·Ψ_m,
       Ψ_n=[   a_n; a_n^† ],   B_nm=[    V_nm -Δ^∗_nm;    Δ_nm -V_nm^∗ ].


The 2N× 2N Hermitian matrix B is called the Bogoliubov-De Gennes (BdG) Hamiltonian <cit.>. Its eigenvalues come in pairs ± E_1, ± E_2,…± 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 

    B_nm↦ U B_nmU^†≡ A_nm,  with  U=1/√(2)[  1  1; -i  i ],

maps B onto the 2N× 2N imaginary antisymmetric matrix A with elements

    A_nm=[ i  Im (V_nm+Δ_nm) i  Re (Δ_nm+V_nm); i  Re (Δ_nm-V_nm) i  Im (V_nm-Δ_nm) ]=- A_mn^ T.

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

The transformed state

    γ=(γ_1,γ_2,…γ_2N),  with  [ γ_2n-1;   γ_2n ]=U[   a_n; a^†_n ],

contains 2N Hermitian operators γ_n=γ_n^†, with anticommutator

    γ_nγ_m+γ_mγ_n=δ_nm,  γ_n^2=1/2.

This is the Clifford algebra of Majorana operators.

The global fermion parity operator 

    P=(-1)^∑_n=1^N a^†_n a^_n=(-2i)^Nγ_1γ_2⋯γ_2N

commutes with H, so energy eigenstates have a definite fermion parity ± 1. Kitaev's formula <cit.> equates the fermion parity P_0 of the ground state to the Pfaffian[Wikipedia has a helpful collection of Pfaffian formulas.] (Pf) of the Hamiltonian in the Majorana basis,

    P_0=sign Pf (-i A),  for  H=12γ· A·γ.
















 §.§ Pfaffian formula for a subsystem



Our objective is to calculate the ground-state expectation value of the fermion parity P_ L of an open subsystem, say the left quantum dot with N_L fermionic modes in Fig. <ref> .

A direct way to proceed, used for example in Ref. Cla17, is to calculate the many-particle ground state |Ψ_0⟩ in the basis of occupation numbers and evaluate

    ⟨ P_ L⟩=⟨Ψ_0|(-1)^∑_n=1^N_ La^†_n a^_n|Ψ_0⟩.

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 <cit.> has developed an efficient method, with a polynomial scaling in N, to calculate squares of expectation values of operators exp(iχ∑_na^†_n a^_n). This gives ⟨ P_ L⟩^2 if one sets χ=π and restricts the sum to indices n in L. In App. <ref> we show how the Klich method can be adapted to give also the sign of ⟨ P⟩_ L. That calculation is technically rather involved, but the final result can be easily understood as follows.

We make the flat-band transformation A↦A̅, which consists in replacing each of the 2N eigenvalues ± E_n of 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 (<ref>) invariant. And because the Pfaffian of -iA̅ can only equal ± 1 we no longer need to take the sign in Eq. (<ref>), hence the global fermion parity is 
    P_0=Pf (-iA̅).


At this point one may surmise that the desired subsystem generalization of Eq. (<ref>) simply amounts to taking the Pfaffian of the 2N_ L× 2N_ L submatrix [A̅]_ LL restricted to the subspace of modes in the left quantum dot,

    ⟨ P_ L⟩=Pf [-iA̅]_ LL.

This is indeed the correct expression, as one can see by application of the Wick theorem for Majorana operators <cit.>,

    ⟨γ_1γ_2⋯γ_2s⟩=1≤ k<l≤ 2s Pf ⟨γ_kγ_l⟩.

Substitution of P_ L=(-2i)^ N_ Lγ_1γ_2⋯γ_2N_ L on the left-hand-side and -2i⟨γ_kγ_l⟩=-iA̅_kl on the right-hand-side results in Eq. (<ref>). This is how an equivalent formula was derived recently for a different problem <cit.>.

Eq. (<ref>) is computationally efficient because the Pfaffian of an N× N matrix can be calculated in a time that scales polynomially with N <cit.>: It has the same O(N^3) complexity as the eigenvalue decomposition one needs for the flat-band transformation A↦A̅. Note that the flat-band transformation needs to be performed before the subblock restriction A̅↦[A̅]_LL — the two operations do not commute.



§ CONNECTION WITH THE MAJORANA FUSION RULE



As a fundamental application of Eq. (<ref>), consider the case that each quantum dot in Fig. <ref> has a single electronic mode (N_ L=N_  R=1), each consisting of two Majorana modes with inter-dot coupling matrix Γ 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

    A=[      0     iΓ; -iΓ^ T      0 ].

The global fermion parity is

    P_0= sign  Pf (-i A)=- sign  Det Γ.


To obtain the average local fermion parity we use that the real 2× 2 coupling matrix Γ has the singular value decomposition Γ=O_1  diag (κ_1,κ_2)O_2, with O_1,O_2 real orthogonal matrices and κ_1,κ_2>0. The eigenvalues of A are ±κ_1,±κ_2. In the flat-band transformation {κ_1,κ_2}↦{1,1}, which gives

    A̅=[              0        iO_1O_2; -iO_2^ TO_1^ T              0 ]⇒
    [A̅]_ LL=0⇒⟨ P_ L⟩=0,

so the average fermion parity in a single quantum dot vanishes. This is a manifestation of the Majorana fusion rule <cit.>: The fusion of the two Majorana zero-modes γ_1 and γ_2 produces an equal-weight superposition of a state of even and odd fermion parity.[The converse is not excluded:  ⟨ P_ L⟩= 0 without an isolated Majorana zero-mode is possible, for example for

A=i[    0  λ_1    0  λ_2; -λ_1    0 -λ_2    0;    0  λ_2    0  λ_1; -λ_2    0 -λ_1    0 ] with λ_1<λ_2.]

Several recent experimental proposals <cit.> are based on the connection between the Majorana fusion rule and vanishing average fermion parity. The implication “isolated Majorana zero-modes ⇒ ⟨ P_ L⟩=0” holds if there are only two pairs of Majorana zero-modes. For N_ L or N_ R greater than 1 the implication breaks down, as is demonstrated by the following counterexample for N_ L=N_ R=2: 



    A=[     iΩ     iΓ; -iΓ^ T     iΩ ],  Ω=[  0  0  0  0;  0  0  0  0;  0  0  0  1;  0  0 -1  0 ],  Γ=[ 0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0 ],
       ⇒A̅=1/√(5)[     iΩ'     iΓ'; -iΓ'^ T     iΩ' ],  Γ'=2Γ,  Ω'=[  0 -1  0  0;  1  0  0  0;  0  0  0  1;  0  0 -1  0 ],
       ⇒[A̅]_ LL=i/√(5)Ω'⇒ Pf [-iA̅]_ LL=-1/ 5,


and hence ⟨ P_ L⟩=-1/ 5 does not vanish even though each quantum dot has a pair of Majorana zero-modes without intra-dot coupling (γ_1 and γ_2 in the left dot, γ_5 and γ_6 in the right dot).

Since Pf (-i A)=+1 the global fermion parity is even, hence the negative sign for ⟨ P_ L⟩ 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 γ_1 and γ_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 |Ψ_0⟩ of the Hamiltonian H=12γ· A·γ. We find

    |Ψ_0⟩=√(5)/10[2i(a_1^† a_2^†+a_3^† a_4^†)   -(1+√(5))a_1^† a_3^†
       -(1-√(5))a_2^† a_4^†]|0⟩,

which indeed gives ⟨ P_ L⟩=-1/ 5 upon calculation of the expectation value (<ref>). 

In this case with N=N_ L+N_ 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. (<ref>) or on Eq. (<ref>) are equally efficient. For larger N the single-particle approach based on the Pfaffian formula has the more favorable scaling (polynomial instead of exponential).



§ RANDOM-MATRIX THEORY



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



 §.§ Skew Circular Real Ensemble



The matrix [-iA̅]_ LL which in view of Eq. (<ref>) determines the local fermion parity is a 2N_ L× 2N_ L submatrix of a matrix S=-iA̅ that is an antisymmetric (skew-symmetric) element of the real orthogonal group O(2N), with N=N_ L+N_ R. The corresponding ensemble from RMT is the class-DIII circular ensemble, which differs from the class-D circular ensemble by the antisymmetry restriction <cit.>.  The latter is called the Circular Real Ensemble (CRE) and we will refer to the former as the skew-Circular Real Ensemble (skew-CRE).[The qualifier “real” for the O(N) ensemble is used instead of “orthogonal” because the name Circular Orthogonal Ensemble (COE) was already used by Dyson <cit.> for the coset U(N)/ O(N).] The switch from symmetry class D to DIII is remarkable, because class DIII was originally introduced <cit.> 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 A with independent Gaussian elements on the upper diagonal (zero mean and unit variance), and perform the flat-band transformation to obtain S=-iA̅.

  * Draw a random element O from O(2N), uniformly with the invariant Haar measure, and construct

    S=O[  0_N× N  1_N× N; -1_N× N  0_N× N ]O^ T.


The two methods are equivalent because the distribution P( A)∝exp(14 Tr  A^2) as well as the flat-band transformation A↦A̅ are invariant under orthogonal transformations A↦ O AO^ T, so the matrix O in the decomposition (<ref>) is distributed according to the invariant Haar measure.

The matrix S has the block decomposition

    S=[ S_ LL S_ LR; S_ RL S_ RR ],   S_ LL=[-iA̅]_ LL,

with S_XY a matrix of dimension N_X× N_Y. In the context of scattering problems, where the skew-CRE ensemble was studied previously <cit.>, this is analogous to a decomposition of the scattering matrix into reflection and transmission matrices. In that context the eigenvalues ± iλ_n of the upper-left submatrix S_ LL correspond to reflection amplitudes.[Eq. (<ref>) follows from equation 5 of Ref. Dah10 upon change of variables from transmission probabilities T_n to reflection amplitudes λ_n=√(1-T_n).] Their joint probability distribution in the skew-CRE is known <cit.>,

    P(λ_1,λ_2,…λ_N_ min)∝∏_n(1-λ_n^2)^|N_ L-N_ R|∏_j<k(λ_k^2-λ_j^2)^2,
            N_ min=min(N_ L,N_ R),  0≤λ_n≤ 1.

If N_ L>N_ R there are additionally 2(N_ L-N_ R) trivial eigenvalues pinned at ± 1, not included in the distribution (<ref>).

Symmetry class DIII has the ℤ_2 invariant Pf  S=± 1, which in view of Kitaev's formula (<ref>) is the global fermion parity P_0. This does not enter in Eq. (<ref>) because in the skew-CRE the distribution of the λ_n's is independent of the ℤ_2 invariant <cit.>.

The density ρ(λ) of the nontrivial eigenvalues has ±λ symmetry with a three-peak structure: There are two peaks at  the band edges ±λ_c, with <cit.>

    λ_c=(2/N)(N_ LN_ R)^1/2,

and a peak at the band center[For the density profile near λ=0 we can approximate the distribution (<ref>) by P({λ})∝∏_j<k(λ_k^2-λ_j^2)^2 and ignore the restriction |λ_n|≤ 1. The distribution of the λ_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 (<ref>). 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 ρ(λ).] described by <cit.>

    ρ(λ)=1/δ_ eff+sin(2πλ/δ_ eff)/2πλ,  λ≲ 1/δ_ eff.

The parameter δ_ eff=π/2N_ min is the mean eigenvalue spacing in the center of the band. The peak at λ=0 is a weak antilocalization effect in the scattering context <cit.>. 



Fig. <ref> shows the eigenvalue density for N_ L=N_ R=N_ dot ranging from 1 to 6. The three-peaked structure is evident except for N_ dot=1, when the density profile is flat.



 §.§ Distribution of the local fermion parity in the skew-CRE





The peak at λ=0 in the eigenvalue density ρ(λ) increases the probability for vanishing local fermion parity, since

    |⟨ P_ L⟩|=∏_n=1^N_ minλ_n=√( Det  S_ LL).

Indeed, as shown in Fig. <ref>, while the distribution of ⟨ P_ L⟩ in the skew-CRE is broad for a single electronic level N_ dot=1 in each quantum dot, it quickly narrows to a sharp peak at ⟨ P_ L⟩=0 with just a few levels — in accord with numerical calculations reported by Clarke, Sau, and Das Sarma <cit.>.

The peak at zero ⟨ P_ L⟩≡ p appears as a sharp cusp in Fig. <ref>, it has a logarithmic singularity ∝ (p^2ln|p|)^N_ dot-1, for example

    P(⟨ P_ L⟩=p)=4532(1-p^4+4 p^2 ln|p|),   N_ dot=2,  |p|≤ 1.

For large-N_ dot the width of the distribution becomes exponentially small, as follows from the variance

    Var ⟨ P_ L⟩=(2N_ dot)!^3/(N_ dot)!^2(4N_ dot)!= √(2)/4^N_ dot[1+ O(1/N_ dot)],

see App. <ref>.

We may quantify the effect of the spectral peak in ρ(λ) on the distribution of the local fermion parity by comparing with a set of independent λ_n's with uniform density. In that uniform case one would have the fermion parity distribution

    P_ uniform(⟨ P_ L⟩=p)=(-ln|p|)^N_ dot-1/2(N_ dot-1)!,  |p|≤ 1,

with a variance 3^-N_ dot that decays less rapidly than Eq. (<ref>).



 §.§ 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=12γ· A·γ of the two coupled quantum dots of Fig. <ref> has the block structure

    A=[  iΩ_ L     iΓ; -iΓ^ T  iΩ_ R ].

The real antisymmetric matrices Ω_X of size 2N_X× 2N_X, with X∈{R,L}, describe the left and right quantum dot in isolation, while the 2N_ L× 2N_ R real matrix Γ describes the coupling via a quantum point contact (QPC) with N_ QPC propagating fermionic modes. In what follows we take N_ L=N_ R=N_ dot.

The number N_ dot counts the number of electronic modes in each quantum dot. One electronic mode a_n corresponds to two Majorana modes γ_2n-1 and γ_2n, according to

    a_n=(γ_2n-1+iγ_2n)/√(2),

cf. Eq. (<ref>). Because of this double-counting, the mean level spacing δ_0 of eigenstates of Ω_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 Ω_X. Each upper-diagonal matrix element has zero mean and variance  2N_ dotδ_0^2/π^2, corresponding to superconductors in symmetry class D (broken time-reversal and broken spin-rotation symmetry) <cit.>. 

Following Refs. Bar00,Die06, the quantum dots are coupled by a Gaussian random matrix Γ of rank N_ QPC, with elements[The coupling matrix (<ref>) describes a ballistic point contact. For tunnel coupling, rather than ballistic coupling, the coupling strength δ_0/π is to be multiplied by T_n^-1(2-T_n-2√(1-T_n)), with T_n the tunnel probability of the n-th mode in the QPC, see Ref. Bee15.]

    Γ_kl=2N_ dotδ/π∑_n=1^2N_ QPCv^(n)_k w^(n)_l,

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

For the weak coupling regime we focus on the case of a single propagating electronic mode in the point contact, N_ QPC=1, corresponding to two propagating Majorana modes. We do not have an analytical solution, so we show numerical results in Fig. <ref> for the probability distribution of ⟨ P_ L⟩= Pf (-iA̅) in the ensemble of random matrices Ω_ L, Ω_ R, and Γ. The variance of the distribution is compared with that in the skew-CRE in Fig. <ref>. 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.



§ EFFECT OF AN ISOLATED MAJORANA ZERO-MODE



The random Hamiltonians of the previous section do not contain isolated Majorana zero-modes: the 2N_ 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 Ω_ L or Ω_ R in the Hamiltonian (<ref>). (The row and column number should be the same to preserve the antisymmetry of Ω_X.) The effect on the distribution of the local fermion parity is shown in Fig. <ref>. 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.



§ CONCLUSION


	
In summary, we have studied the fusion of Majorana zero-modes using a generalization of Kitaev's Pfaffian formula <cit.> 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. (<ref>), and an equivalent formulation from Ref. Bau18, 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. <ref> 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 <cit.>. 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 <cit.>, 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 ⟨ P_ L⟩≈ 0 of the local fermion parity. The Majorana fusion rule provides a fundamentally different, topological mechanism <cit.>: 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 ⟨ P_ L⟩ because the two fusion channels, with or without an unpaired quasiparticle, have equal weight.

To carry out such a fusion experiment it is proposed <cit.> that one would rapidly decouple the subsystems, on a time scale τ_c sufficiently short that quasiparticles from the environment cannot leak in. The complication <cit.> is that even if there are isolated Majorana zero-modes, the presence of even a small number N_ dot of higher levels at energies below ħ/τ_c may hide the presence of the zero-modes by favoring ⟨ P_ L⟩≈ 0 (see Fig. <ref>).

Fig. <ref> 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.


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






§ 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 <cit.>, to reproduce his result for ⟨ P_ L⟩^2. Then we will resolve the sign ambiguity to arrive at Eq. (<ref>) for ⟨ P_ L⟩. An equivalent formula is obtained by a different method in Ref. Bau18, Appendix B.

The superconductor in Fig. <ref> 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 (<ref>) 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→ 0 limit at the end of the calculation in order to describe an isolated system.

At inverse temperature β=1/k_BT the average fermion parity ⟨ P_L⟩ of subsystem L (the left quantum dot in Fig. <ref>) is given by the trace of the equilibrium density matrix 

    ρ_eq=1/Z e^-β H,  Z=Tr ρ_eq,

acting on the fermion parity operator

    P_L=exp(iπ∑_n∈L a^†_na^_n).

Because H=12∑_n,m A_nmγ_nγ_m in the basis of Majorana operators γ_n, and a^†_na^_n=iγ_2n-1γ_2n+12, this can be written as

    ⟨ P_L⟩=   e^iπ N_L/2/Z Tr [exp(-12β∑_n,m A_nmγ_nγ_m).
       ×.exp(-12iπ∑_n,m(σ_y⊗ P_L)_nmγ_nγ_m)].

The matrix σ_y is a Pauli matrix and the operator P_L projects onto N_L fermionic modes in subsystem L.

Application of the identity <cit.> 

    [ Tr ∏_ke^γ· O_k·γ]^2=e^∑_k Tr O_k  Det (1+∏_ke^O_k-O_k^ T),

results in

    ⟨ P_L⟩^2   =e^iπ N_L  Det [1+exp(-β A)exp(-iπσ_y⊗ P_L)]/ Det [1+exp(-β A)]
       =(-1)^N_L  Det [1-2/1+exp(β A) (σ_0⊗ P_L)].

In the second equality we made use of the identity

    e^iχ σ_y⊗ P_L=1+σ_0⊗ P_L(cosχ-1)+iσ_y⊗ P_Lsinχ,

with χ=π. (The matrix σ_0=σ_y^2 is the 2× 2 unit matrix.) Note that, in a basis of energy eigenstates of the BdG Hamiltonian, the operator (1+e^β A)^-1 is the Fermi function f(E)=(1+e^β E)^-1. 

Eq. (<ref>) is Klich's result for the square of the average fermion parity (equation 84 in Ref. Kli14). Klich shows how the sign of ⟨ P_ L⟩ 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× 2N imaginary anti-symmetric matrix A can be decomposed as

    A=iO(J⊗ E)O^T,  J=[  0  1; -1  0 ],

where O is a 2N× 2N real orthogonal matrix and
 E=diag (E_1,E_2,… E_N) is an N× N real diagonal matrix. Substitution into Eq. (<ref>) gives

    ⟨ P_L⟩^2=(-1)^N_L  Det [1-O2/1+exp(iβ J⊗ E) O^T(σ_0⊗ P_L)]
         =(-1)^N_L  Det [1-O[1-iJ⊗tanh(12β E)] O^T(σ_0⊗ P_L)].


This may be written in a more compact form by defining the restriction [M]_LL of a 2N× 2N matrix M to the 2N_L× 2N_L submatrix of modes in region L,

    ⟨ P_L⟩^2   =(-1)^N_L  Det [O[iJ⊗tanh(12β E)]O^T]_LL
       =Det [O[J⊗tanh(12β E)]O^T]_LL.

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

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

    ⟨ P_L⟩=Pf [O[J⊗tanh(12β E)]O^T]_LL.

In the zero-temperature, β→∞ limit this reduces to

    ⟨ P_L⟩=Pf [O[J⊗( sign  E)]O^T]_LL,

which is Eq. (<ref>) with -iA̅=O[J⊗(sign  E)]O^T. Kitaev's formula (<ref>) for the global ground-state fermion parity is recovered when L is the entire isolated system. This correspondence also identifies √()Det with +Pf rather than with -Pf.



§ MOMENTS OF DETERMINANTS OF ANTISYMMETRIC RANDOM MATRICES



In Sec. <ref> 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 <cit.>. We therefore give the derivation in this appendix, and for completeness and reference also derive the corresponding result for antisymmetric Hermitian matrices.



 §.§ Principal minor of antisymmetric orthogonal matrix


Consider a 2N× 2N antisymmetric real orthogonal matrix S, with a uniform distribution in O(2N) subject to the antisymmetry constraint. This is the class-DIII circular ensemble of RMT <cit.>, referred to as the skew-Circular Real Ensemble (skew-CRE) in the main text.[The antisymmetric orthogonal matrices form a disconnected set in O(2N), distinguished by the sign of the Pfaffian. For the probability distribution (<ref>) it does not matter whether or not we restrict the ensemble to Pf  S=± 1.]

The 2N_ L× 2N_ L upper-left submatrix S_ LL has eigenvalues ± i λ_n, 0≤λ_n≤ 1. Denoting N_ R = N-N_ L and N_ min=min(N_ L,N_ R), we have that N-N_ min of the λ_n's are pinned to +1. The set {λ_n}={λ_1,λ_2,…λ_N_ min} can vary freely in the interval [0,1], with joint probability distribution <cit.>

    P({λ_n }) ∝∏_n (1-λ_n^2)^|N_ L - N_ R|∏_i<j (λ_i^2 - λ_j^2)^2
      .

The determinant of S_ LL is a principal minor given by

    Det  S_ LL
      = ∏_n=1^N_ L (i λ_n) (-i λ_n)
      = ∏_n=1^N_ minλ_n^2
      .

We seek the moments μ_q=𝔼[( Det  S_ LL)^q] of this determinant in the skew-CRE.

For that purpose we make a change of variables from λ_n to R_n = λ_n^2∈[0,1],
with distribution

    P({ R_n }) ∝∏_n R_n^-1/2 (1-R_n)^|N_ L-N_ R|∏_i<j (R_i-R_j)^2
      .

We can then compute the moments
of the determinant from

    μ_q
      = ∫_0^1 d{R_n}∏_i<j (R_i-R_j)^2
        ∏_n (1-R_n)^|N_L-N_R| R_n^q-1/2/∫_0^1 d {R_n}∏_i<j (R_i-R_j)^2
        ∏_n (1-R_n)^|N_L-N_R| R_n^-1/2
      ,

where we abbreviated ∫_0^1 d{R_n}= ∫_0^1 d R_1 ⋯∫_0^1 d R_N_ min.

These socalled Selberg integrals have a closed-form expression <cit.>,

    μ_q=
        ∏_j=0^N_ min-1Γ( max(N_ L,N_ R)+j+1/2)
          Γ( q+j+1/2)
        /Γ( max(N_ L,N_ R)+q+j+1/2)
          Γ( j+1/2)
        
        .

For the first few moments, Eq. (<ref>) reduces to

    μ_1= 
      (2 N_ L)! (2N_ R)! N!/N_ L! N_ R! (2N)!
      ,
       μ_2
      = (2N_ L+1)(2N_ R+1)/2N+1μ_1^2
      .

Eq. (<ref>) in the main text is Eq. (<ref>) for N_ L=N_ R=N_ dot=N/2.



 §.§ Antisymmetric Hermitian matrix


A similar calculation can be carried out for moments of the determinant of a 2N× 2N antisymmetric Hermitian matrix 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 ±λ_n. The N eigenvalues λ_n≥ 0 have the joint distribution <cit.>

    P({λ_n }) ∝∏_n  e^- λ_n^2/2∏_i<j (λ_i^2 - λ_j^2)^2
    .

The determinant is

    Det  A =(-1)^N ∏_n=1^N λ_n^2
      .


Let us introduce the variables x_n =λ_n^2/2≥ 0, with distribution

    P({ x_n }) ∝∏_n
      x_n^-1/2 e^- x_n∏_i<j (x_i - x_j)^2.

The q-th moment μ_q of the determinant of A is given by

    μ_q= (-2)^Nq∫_0^∞ d{x_n}∏_i<j (x_i-x_j)^2 ∏_n x_n^q - 1/2 e^-x_n/∫_0^∞ d{x_n}∏_i<j (x_i-x_j)^2 ∏_n x_n^- 1/2 e^-x_n
      ,

with ∫_0^∞ d{x_n}=∫_0^∞ dx_1⋯∫_0^∞ dx_N.
This is the ratio of normalisation constants of Laguerre
distributions, which is known <cit.>.
We thus obtain

    μ_q
        = (-2)^Nq∏_j=0^N-1Γ(q+N-j-1/2)/Γ(N-j-1/2)
        .


For q=1,2 this reduces to

    μ_1= (-1)^N
      (2N)!/2^N N!,  μ_2= 
      (2N+1)!(2N)!/2^2N (N!)^2,
       ⇒ Var ( Det  A) = 2 N [𝔼( Det  A)]^2.

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

    μ_1=√(2)(-2/e)^N  N^N[1+ O(1/N)],

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

99
Bal06 A. V. Balatsky, I. Vekhter, and Jian-Xin Zhu, Impurity-induced states in conventional and unconventional superconductors, Rev. Mod. Phys. 78, 373 (2006).
Sak70 A. Sakurai, Comments on superconductors with magnetic impurities, Prog. Theor. Phys. 44, 1472 (1970).
Kit01 A. Yu. Kitaev, Unpaired Majorana fermions in quantum wires, Phys. Usp. 44 (suppl.), 131 (2001).
Alt97 A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).
Bee15 C. W. J. Beenakker, Random-matrix theory of Majorana fermions and topological superconductors, Rev. Mod. Phys. 87, 1037 (2015).
Cla17 D. J. Clarke, J. D. Sau, and S. Das Sarma, Probability and braiding statistics in Majorana nanowires, Phys. Rev. B 95, 155451 (2017).
Aas16 D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon, M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flensberg, and J. Alicea, Milestones toward Majorana-based quantum computing, Phys. Rev. X 6, 031016 (2016).
Per77 J. H. H. Perk and H. W. Capel, Time-dependent xx-correlation functions in the one-dimensional xy-model, Physica 89A, 264 (1977).
Per84 J. H. H. Perk, H. W. Capel, G. R. W. Quispel, and F. W. Nijhoff, Finite-temperature correlations for the Ising chain in a transverse field, Physica 123A, 1 (1984).
Bra05 S. Bravyi, Lagrangian representation for fermionic linear optics, Quantum Inf. Comp. 5, 216 (2005).
Bau18 B. Bauer, T. Karzig, R. V. Mishmash, A. E. Antipov, and J. Alicea, Dynamics of Majorana-based qubits operated with an array of tunable gates, SciPost Phys. 5, 004 (2018).
Kli14 I. Klich, A note on the full counting statistics of paired fermions, J. Stat. Mech. P11006 (2014). When comparing formulas, note that Klich has a factor of two in the anticommutator of Majorana operators.
Bud13 J.-C. Budich and E. Ardonne, Equivalent topological invariants for one-dimensional Majorana wires in symmetry class D, Phys. Rev. B 88, 075419 (2013).
DeG66 P. G. De Gennes, Superconductivity of Metals and Alloys (Benjamin, New York, 1966).
Rub11 J. Rubow and U. Wolff, A factorization algorithm to compute Pfaffians, Comp. Phys. Comm. 182, 2530 (2011).
Wim12 M. Wimmer, Efficient numerical computation of the Pfaffian for dense and banded skew-symmetric matrices, ACM Trans. Math. Software 38, 30 (2012).
Nay08 C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
Bee19 C. W. J. Beenakker, A. Grabsch, and Y. Herasymenko, Electrical detection of the Majorana fusion rule for chiral edge vortices in a topological superconductor, SciPost Phys. 6, 022 (2019).
Dah10 J. P. Dahlhaus, B. Béri, and C. W. J. Beenakker, Random-matrix theory of thermal conduction in superconducting quantum dots, Phys. Rev. B 82, 014536 (2010).
Dys62 F. J. Dyson, The threefold way: Algebraic structure of symmetry groups and ensembles in quantum mechanics, J. Math. Phys. 3, 1199 (1962).
Mehta M. L. Mehta, Random Matrices (Elsevier, Amsterdam, 2004).
Iva02 D. A. Ivanov, The supersymmetric technique for random-matrix ensembles with zero eigenvalues, J. Math. Phys. 43, 126 (2002). 
Pik12 D. I. Pikulin, J. P. Dahlhaus, M. Wimmer, H. Schomerus, and C. W. J. Beenakker, A zero-voltage conductance peak from weak antilocalization in a Majorana nanowire, New J. Phys.14, 125011 (2012).
Bar00 C. I. Barbosa, T. Guhr, and H. L. Harney, Impact of isospin breaking on the distribution of transition probabilities, Phys. Rev. E 62, 1936 (2000).
Die06 B. Dietz, T. Guhr, H. L. Harney, and A. Richter, Strength distributions and symmetry breaking in coupled microwave billiards, Phys. Rev. Lett. 96, 254101 (2006).
Forrester P. J. Forrester, Log-Gases and Random Matrices (Princeton, 2010).
Girko V. L. Girko, Theory of Random Determinants (Springer, Berlin, 1990). 
For09 P. J. Forrester and E. Nordenstam, The anti-symmetric GUE minor process, Mosc. Math. J. 9, 749 (2009).
Ede15 A. Edelman and M. La Croix, The singular values of the GUE (Less is More), Random Matrices: Theory and Applications 04, 1550021 (2015).


