A zero-voltage conductance peak from weak antilocalization in a Majorana nanowire

We construct the scattering matrix of the NS junction at the Fermi level by assuming a spatial separation of normal scattering in N and Andreev reflection in S. Within the excitation gap there is no transmission through the superconductor. The matrix r_A of Andreev reflection amplitudes from the superconductor is then a 2N × 2N unitary matrix. Mode mixing at the NS interface can be incorporated into the scattering matrix S_N of the normal region, so we need not include it in r_A. It has the block form [27, 28]

$$\begin{eqnarray} r_{\mathrm{A}}&=&\left(\begin{array}{@{}ll@{}} \Gamma & \Lambda\\ \Lambda^{\ast} & \Gamma \end{array}\right),\quad \Gamma=\bigoplus_{m=1}^{M}\left(\begin{array}{@{}ll@{}} \cos\alpha_{m} & 0\\ 0 & \cos\alpha_{m} \end{array}\right)\oplus \emptyset_{Q}\oplus\mathbbm{1}_{\zeta},\nonumber\\ \\ \Lambda&=&\bigoplus_{m=1}^{M}\left(\begin{array}{@{}ll@{}} 0 & -\mathrm{i}\sin\alpha_{m}\\ \mathrm{i}\sin\alpha_{m} & 0 \end{array}\right)\oplus \mathbbm{1}_{Q}\oplus\emptyset_{\zeta}.\nonumber \end{eqnarray} \tag{ 1 }$$

We have defined ζ = 0 if the difference N − Q is even and ζ = 1 if N − Q is odd, so that N − Q − ζ ≡ 2M is an even integer. The Andreev reflection eigenvalues ρ_m = sin²α_m that are not pinned at 0 or 1 are twofold degenerate [29]⁶.

The symbols $\mathbbm{1}_{n},\emptyset_{n}$ denote, respectively, an n × n unit matrix or null matrix for n ⩾ 1. The empty set is intended for n = 0. To make the notation more explicit, we give some examples of the direct sums,

$$\begin{eqnarray} & &\mathbbm{1}_{1}\oplus\emptyset_{1}=\left(\begin{array}{@{}ll@{}} 1 & 0\\ 0 & 0 \end{array}\right),\quad\mathbbm{1}_{2}\oplus\emptyset_{1}=\left(\begin{array}{@{}lll@{}} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0 \end{array}\right),\nonumber\\ \\ & &\mathbbm{1}_{2}\oplus\emptyset_{0}=\left(\begin{array}{@{}ll@{}} 1 & 0\\ 0 & 1 \end{array}\right),\quad \mathbbm{1}_{1}\oplus\emptyset_{0}=1,\quad \mathbbm{1}_{0}\oplus\emptyset_{1}=0.\nonumber\ \end{eqnarray} \tag{ 2 }$$

The normal region has a scattering matrix

$$\begin{equation} S_{\mathrm{N}}=\left(\begin{array}{@{}ll@{}} s_{0} & 0\\ 0 & s_{0}^{\ast} \end{array}\right),\quad s_{0}=\left(\begin{array}{@{}ll@{}} r' & t'\\ t & r \end{array}\right). \end{equation} \tag{ 3 }$$

The electron and hole blocks (with N × N reflection and transmission matrices r,r',t,t') are each other's complex conjugate at the Fermi level. The off-diagonal blocks of S_N vanish, because the normal metal cannot mix electrons and holes. The matrix s₀ is unitary, s₀s^†₀ = 1, without further restrictions in class D. In class BDI, chiral symmetry requires that s₀ = s^T₀ is also a symmetric matrix.

To separate the mixing of modes from backscattering, we make use of the polar decomposition

$$\begin{equation} s_{0}=\left(\begin{array}{@{}ll@{}} U & 0\\ 0 & V \end{array}\right)\left(\begin{array}{@{}ll@{}} -\sqrt{1-{\cal T}} & \sqrt{\cal T}\\ \sqrt{\cal T} & \sqrt{1-{\cal T}} \end{array}\right)\left(\begin{array}{@{}ll@{}} U' & 0\\ 0 & V' \end{array}\right). \end{equation} \tag{ 4 }$$

The matrices U,V,U' and V ' are N × N unitary matrices and ${\cal T}=\mathrm {diag}\,(T_{1},T_{2},\ldots , T_{N})$ is a diagonal matrix of transmission eigenvalues of the normal region. In class BDI, chiral symmetry relates U' = U^T, V ' = V^T.

We combine S_N and r_A to obtain the matrix r_he of Andreev reflection amplitudes (from electron e to hole h) of the entire system. This calculation is much simplified in the case of ζ = 0, ρ_m = 1 (m = 1,2,...,M) that all modes at the NS interface are Andreev reflected with unit probability. For this case Γ = 0, N − Q = 2M, we obtain

$$\begin{equation} r_{he}=t'^{\ast}\Lambda^{\ast}(1-r\Lambda r^{\ast}\Lambda^{\ast})^{-1}t,\quad \Lambda=\sigma_{y}^{\oplus M}\oplus \mathbbm{1}_{Q}. \end{equation} \noindent \tag{ 5 }$$

The notation σ^⊕M_y signifies the 2M × 2M matrix constructed as the direct sum of M Pauli matrices.

The Andreev reflection matrix determines the conductance

$$\begin{equation} G=G_{0}\,\mathrm{Tr}\,r_{he}^{\vphantom{\dagger}}r_{he}^{\dagger},\quad G_{0}=2e^{2}/h. \end{equation} \noindent \tag{ 6 }$$

Substitution of the polar decomposition (4) gives the compact expression

$$\begin{equation} \begin{array}{l} G/G_{0}=\mathrm{Tr}\,{\cal T}{\cal M}{\cal T}{\cal M}^{\dagger},\\ \hspace{1.4pc}{\cal M}=(1-\Omega^{\ast}\sqrt{1-{\cal T}}\Omega\sqrt{1-{\cal T}})^{-1}\Omega^{\ast},\quad \Omega=V'\Lambda V^{\ast}. \end{array} \end{equation} \noindent \tag{ 7 }$$

This is the zero-temperature conductance at the Fermi level, in the limit of zero bias voltage. Away from the Fermi level particle–hole symmetry is broken, so the electron and hole blocks in S_N are distinct unitary matrices s_e and s_h. If the bias voltage V₀ remains small compared to the excitation gap, we can keep the same r_A. The finite-voltage differential conductance $\tilde {G}=\mathrm {d}I/\mathrm {d}V_{0}$ is then given by

$$\begin{eqnarray} &&\begin{array}{l} \tilde{G}/G_{0}=\mathrm{Tr}\,{\cal T}_{\mathrm{h}}\tilde{\cal M}{\cal T}_{\mathrm{e}}\tilde{\cal M}^{\dagger},\\ \tilde{\cal M}=(1-\Omega_{\mathrm{h}}^{\ast}\sqrt{1-{\cal T}_{\mathrm{e}}}\Omega_{\mathrm{e}}\sqrt{1-{\cal T}_{\mathrm{h}}})^{-1}\Omega_{\mathrm{h}}^{\ast},\\ \Omega_{\mathrm{e}}=V'_{\mathrm{e}}\Lambda V_{\mathrm{h}}^{\ast},\;\;\Omega_{\mathrm{h}}=V'_{\mathrm{h}}\Lambda V_{\mathrm{e}}^{\ast}. \end{array} \end{eqnarray} \noindent \tag{ 8 }$$

The electron matrices are evaluated at energy eV₀ above the Fermi level and the hole matrices at energy −eV₀ below the Fermi level. Chiral symmetry remains operative away from the Fermi level; hence V '_e = V^T_e, V '_h = V^T_h ⇒Ω_h = Ω^†_e in class BDI. We will apply equation (8) to voltages large compared to the Thouless energy, when the electron and hole matrices may be considered to be statistically independent.

Isotropic mixing of the modes by scattering in the normal region means that the unitary matrices in the polar decomposition (4) are uniformly distributed in the unitary group ${\cal U}(N)$. We can calculate the average conductance for a given set of transmission eigenvalues by integration over ${\cal U}(N)$ with the uniform (Haar) measure. A full average would then still require an average over the T_n's, but if these are dominated by a tunnel barrier they will fluctuate a little and the partial average over the unitary matrices is already informative.

The calculation is easiest if all T_n's have the same value 0 ⩽ T ⩽ 1. The average zero-voltage conductance 〈G〉 is then given by the integral

$$\begin{equation} \langle G\rangle=T^{2}G_{0}\int_{0}^{2\pi}\mathrm{d}\phi\,\rho(\phi)\left|1-(1-T)\mathrm{e}^{\mathrm{i}\phi}\right|^{-2}, \end{equation} \noindent \tag{ 9 }$$

with $\rho (\phi )=\langle \sum _{n}\delta (\phi -\phi _{n})\rangle$ the density on the unit circle of the eigenvalues e^iϕ_n of ΩΩ*. The corresponding finite-voltage expression has a uniform ρ = N/2π, leading to

$$\begin{equation} \langle\tilde G\rangle=NG_{0}T/(2-T), \end{equation} \noindent \tag{ 10 }$$

irrespective of the symmetry class and independent of the topological quantum number Q.

The zero-voltage average (9) does depend on Q and is different for class D and BDI. The calculations are given in the appendix. Explicit expressions in class D are

$$\begin{equation} \frac{\langle G\rangle_{\mathrm{D}}}{G_{0}}=\left\{ \begin{array}{@{}l@{\quad}llll@{}} 2T, & \mathrm{for}\;\;N=2,\;\;Q=0,\\ 1+2T^{2}, & \mathrm{for}\;\;N=3,\;\;Q=1,\\ 2T(2-T+T^2), & \mathrm{for}\;\;N=4,\;\;Q=0,\\ 1+2T^{2}(3-2T+T^{2}), & \mathrm{for}\;\;N=5,\;\;Q=1. \end{array}\right. \end{equation} \noindent \tag{ 11 }$$

The Q dependence appears to second order in the reflection probability R = 1 − T, while the general first-order result

$$\begin{equation} \langle G/G_{0}\rangle_{\mathrm{D}}=N(1-2R)+2R+{\cal O}(R^{2}) \end{equation} \noindent \tag{ 12 }$$

is Q independent. The corresponding expressions in class BDI are more lengthy, and we only record the small-R result

$$\begin{equation} \langle G/G_{0}\rangle_{\mathrm{BDI}}=N(1-2R)+2R\frac{Q^{2}+N}{N+1}+{\cal O}(R^{2}), \end{equation} \noindent \tag{ 13 }$$

to show that it is Q-dependent already to first order in R. These are all finite-N results. In the large-N limit the Q dependence is lost,

$$\begin{equation} \langle G/G_{0}\rangle=\frac{NT}{2-T}+\frac{2(1-T)}{(2-T)^{2}}+{\cal O}(N^{-1}), \end{equation} \noindent \tag{ 14 }$$

irrespective of the symmetry class.

As illustrated in figure 2, for the case when all T_n's have the same value T the difference $\delta G=\langle G\rangle -\langle \tilde {G}\rangle$ is positive, corresponding to weak antilocalization and a conductance peak. The sign of the effect may change if the T_n's are very different, in particular, in class BDI—which has δG < 0 in a quantum dot geometry (circular ensemble) [28]. It is a special feature of quantum interference in a magnetic field that the distinction between weak localization and antilocalization is not uniquely determined by the symmetry class [21, 31, 32].

Zoom In Zoom Out Reset image size

Following  [6, 7], we consider a conducting channel parallel to the x-axis on a substrate in the x–y-plane (width W and Fermi energy E_F), in a magnetic field B (orientation $\hat {n}$ and Zeeman energy $E_{\mathrm {Z}}=\frac {1}{2}g_{\mathrm {eff}}\mu _{\mathrm {B}} B$), with Rashba spin–orbit coupling (characteristic energy E_so = m_effα²_so/ℏ² and length l_so = ℏ²/m_effα_so) and induced s-wave superconductivity (excitation gap Δ₀). The Hamiltonian is

$$\begin{equation} \begin{array}{l} \displaystyle {\cal H}=\left(\begin{array}{@{}ll@{}} H_{0}-E_{\mathrm{F}} & \Delta\sigma_{y}\\ \Delta^{\ast}\sigma_{y} & E_{\mathrm{F}}-H_{0}^{\ast} \end{array}\right),\\ \displaystyle H_{0}=\frac{p_{x}^{2}+p_{y}^{2}}{2m_{\mathrm{eff}}}+U(x,y)+\frac{\alpha_{\mathrm{so}}}{\hbar}(\sigma_{x}p_{y}-\sigma_{y}p_{x})+E_{\mathrm{Z}}\hat{n}\cdot\boldsymbol{\sigma}. \end{array} \end{equation} \noindent \tag{ 15 }$$

The electrostatic potential U = U_gate + δU contains the gate potential U_gate that creates the tunnel barrier and the impurity potential δU that varies randomly from site to site on a square lattice (lattice constant a), distributed uniformly in the interval (−U_disorder,U_disorder). The disordered region is −L_N < x < L_S, an NS interface is constructed by increasing the pair potential Δ from 0 to Δ₀ at x = 0 and a rectangular barrier of height U_barrier, thickness δL_barrier is placed at x = −x_barrier. The conductance of the normal region (x < 0) contains a contribution G_disorder from disorder and G_barrier from the barrier.

The orientation of the magnetic field plays an important role [6, 7]: it lies in the x–y-plane to eliminate orbital effects on the superconductor and we will only include its effect on the electron spin (through the Zeeman energy). A topologically nontrivial phase needs a nonzero excitation gap for E_Z > Δ₀, which requires a parallel magnetic field B_∥ ($\hat {n}=\hat {x}$). We will consider that case in the next subsection and then discuss the case of a perpendicular magnetic field B_⊥ ($\hat {n}=\hat {y})$ in section 3.3.

To avoid the complications from chiral symmetry, we first focus on a relatively wide junction, W = 3 l_so, when symmetry class D (rather than BDI) applies [28]. (We turn to class BDI in the next subsection.) The normal region has N = 8 propagating modes (including spin) in zero magnetic field, for E_F = 12 E_so. The topological quantum number Q was determined both from the determinant of the reflection matrix [33]⁷ and independently by counting the gap closings and reopenings upon increasing the magnetic field. A transition from Q = 0 to 1 is realized by increasing B_∥ at fixed Δ₀ = 8 E_so.

The results are shown in figure 3 (solid curves) for two geometries, one with the tunnel barrier far from the NS and another with the barrier close to the interface⁸.

The disorder-averaged conductance shows a zero-voltage peak in a magnetic field, regardless of whether the nanowire is topologically trivial (Q = 0) or nontrivial (Q = 1). The peak disappears in zero magnetic field and instead a conductance minimum develops, indicative of an induced superconducting minigap in the normal region. The two geometries in panels 3(a) and (b) show comparable results, the main difference being a broadening of the zero-bias peak when the tunnel barrier is brought closer to the NS interface—as expected from the increase in Thouless energy⁹. The shallow maximum which develops around zero voltage in the B = 0 curve of panel 3(b) is a precursor of the reflectionless tunneling peak, which appears in full strength when the barrier is placed at the NS interface [5].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All this applies to the average conductance in an ensemble of disordered nanowires. Individual members of the ensemble show mesoscopic, sample-specific conductance fluctuations, in addition to the systematic weak antilocalization effect. For some disorder realizations the zero-voltage conductance peak remains clearly visible, see figure 4. The peak sticks to zero bias voltage over a relatively wide magnetic field range, even though the superconductor is topologically trivial (Q = 0). The appearance and disappearance of the peak are not associated with the closing and reopening of an excitation gap, so it cannot produce Majorana fermions [34]¹⁰.

So far we have considered a class D nanowire with a magnetic field B_∥ parallel to the wire axis. In a perpendicular magnetic field B_⊥ (perpendicular to the wire in the plane of the substrate) the symmetry class remains D (broken time-reversal and spin–rotation symmetry), although the topologically nontrivial phase disappears [6, 7]. We therefore expect the class D zero-bias peak to persist in a perpendicular field as a result of the weak antilocalization effect.

This expectation is borne out by the computer simulations, see the dashed curves in figure 3. A zero-bias peak exists for both B_⊥ and B_∥. If the nanowire is topologically trivial, there is not much difference in peak height for the two magnetic field directions (compare the blue solid and dashed curves). In contrast, if the nanowire is topologically nontrivial for a parallel field, then the peak is much reduced in a perpendicular field (red solid versus dashed curves). The disappearance of the Majorana zero mode and the collapse of the zero-bias peak may be accompanied by the appearance of propagating modes in the superconducting part of the nanowire. This explains the increased background conductance in the red dashed curve of figure 3(a).

The effect of a magnetic field rotation is entirely different when W ≲ l_so and the symmetry class is BDI rather than D [26, 28]. The term σ_xp_y in the Hamiltonian (15) can then be neglected, so that ${\cal H}$ commutes with σ_y in a perpendicular magnetic field ($\hat {n}=\hat {y}$). The two spin components along $\pm \hat {y}$ decouple and for each spin component separately the particle–hole symmetry is broken. We therefore expect both the Majorana resonance and the weak antilocalization peak to disappear in a perpendicular magnetic field for sufficiently narrow wires.

This is demonstrated by the computer simulations shown in figure 5, for the average conductance in a topologically trivial wire of width W = 0.3 l_so. The main difference with the data in figure 3 is that the symmetry class is now BDI rather than D, because of the narrower wire. This change of symmetry class does not significantly affect the weak antilocalization peak in a parallel magnetic field. But if the magnetic field is rotated to a perpendicular direction, the peak disappears—as expected for a class BDI nanowire.

Zoom In Zoom Out Reset image size

All the results presented so far are in the zero-temperature limit. We calculate the temperature dependence of the differential conductance from the finite-T₀ and finite-V₀ generalization of equation (6),

$$\begin{equation} G=\frac{2e}{h}\,\int_{-\infty}^{\infty}\mathrm{d}\varepsilon\,\frac{\mathrm{d}f(\varepsilon-eV_{0})}{\mathrm{d}V_{0}}\,\mathrm{Tr}\,r_{he}^{\vphantom{\dagger}}(\varepsilon)r_{he}^{\dagger}(\varepsilon), \end{equation} \tag{ 16 }$$

$$\begin{equation} f(\varepsilon)=\frac{1}{1+\exp(\varepsilon/k_{\mathrm{B}}T_{0})}. \end{equation} \noindent \tag{ 17 }$$

Thermal averaging at a nonzero temperature T₀ broadens the conductance peak around V₀ = 0 and reduces its height, at constant area $\int G\,\mathrm {d}V_{0}$ under the peak.

This effect of thermal averaging applies to both the weak antilocalization peak and the Majorana resonance, but the characteristic temperature scale is different, as shown in figure 6. The Majorana zero mode is more sensitive to thermal averaging because it is more tightly bound to the NS interface, with a smaller Thouless energy and therefore a smaller characteristic temperature.

Zoom In Zoom Out Reset image size

We employ Dyson's Brownian motion approach [37], which sets up a stochastic process for the unitary matrix U whose stationary distribution coincides with the Haar measure on ${\cal U}(N)$. In each infinitesimal step of the process, $U\to U\exp (\mathrm {i} H)$, where H is a Hermitian matrix from the Gaussian unitary ensemble, with identically normal distributed complex numbers H_lm = H*_ml (l ⩽ m), $\overline {H_{lm}} =0$, $\overline { H_{kl}H_{mn}} =\delta _{kn}\delta _{lm}\tau$, $\overline { H_{kl}H_{mn}^*} =\delta _{km}\delta _{ln}\tau$; the limit τ → 0 is implied to generate infinitesimal increments.

The corresponding increments δμ_n can be calculated in perturbation theory. The drift coefficients $c_l=\lim _{\tau \to 0}\tau ^{-1} \overline {\delta \mu _l}$ and the diffusion coefficients $c_{lm}=\lim _{\tau \to 0}\tau ^{-1} \overline {\delta \mu _l\delta \mu _m}$ follow by averaging over the random variables in H. As we will see, the symmetries in the classes D and BDI are restrictive enough so that these coefficients can be expressed in terms of the quantities μ_n alone, without requiring data from the eigenvectors of X. Thus, the stochastic process for these quantities closes.

Introducing a fictitious time t, the evolution of the joint probability distribution is governed by a Fokker–Planck equation,

$$\begin{equation} \frac{\partial P}{\partial t}=\left[-\sum_l\frac{\partial}{\partial \mu_l} c_l+\frac{1}{2}\sum_{l,m}\frac{\partial}{\partial \mu_l}\frac{\partial}{\partial \mu_m}c_{lm}\right]P\left(\{\mu_n\},t\right). \end{equation} \noindent \tag{ A.1 }$$

The stationary solution P({μ_n}), for which the right-hand side of the Fokker–Planck equation vanishes, is the required eigenvalue distribution.

In class D, we have X = UU* with U uniformly distributed in ${\cal U}(N)$. Note that the operation of complex conjugation is basis dependent; if B = A* in one basis, then this relation is only preserved under orthogonal transformations, but not under general unitary transformations. Thus, we work in a fixed basis |r〉 (at most permitting orthogonal basis changes), and define for any $|\psi \rangle =\sum _r\psi _r|r\rangle$ a complex-conjugated vector $|\psi ^*\rangle \equiv \sum _r\psi _r^*|r\rangle$. As usual, $\langle \psi |=\sum _r\psi _r^*\langle r|$; thus $\langle \psi ^*|=\sum _r\psi _r\langle r|$.

The matrices X and U are unitary and obey Det X = |Det U|² = 1. Moreover, the matrix X* has the same eigenvalues x₁,x₂,...,x_N as the matrix X. For even N, it follows that all eigenvalues appear in complex-conjugated pairs; every eigenvalue x_k has a partner $x_{\bar k}= x^{*}_{k}=x^{-1}_{k}$. For odd N, in addition to such pairs there is a single unpaired eigenvalue, denoted as x_N, which (because of the constraint on the determinant) is pinned at x_N = 1. The paired eigenvectors are related according to

$$\begin{equation} |\bar k\rangle=\xi_{k}U|k^*\rangle. \end{equation} \noindent \tag{ A.2 }$$

Here we have to set ξ_k such that ξ²_k = λ_k; this guarantees that the relation between both eigenvectors in a pair is reciprocal, $|\bar {\bar k}\rangle =|k\rangle$. Observing that the eigenvectors form an orthogonal basis, we find the matrix elements

$$\begin{equation} \langle k |U | l^*\rangle=\xi_k \delta_{k \bar l} =(\langle k^* | U^* | l\rangle)^*=\langle l | U^{\mathrm{T}} | k^*\rangle. \end{equation} \noindent \tag{ A.3 }$$

With the help of these matrix elements we can now evaluate the drift and diffusion coefficients. In second-order perturbation theory,

$$\begin{equation} \delta x_l= \langle l |\delta X| l\rangle+{\sum_k}'\frac{\langle l | \delta X | k\rangle\langle k|\delta X | l \rangle}{x_l-x_k}, \end{equation} \noindent \tag{ A.4 }$$

where the prime restricts the sum to k ≠ l, while

$$\begin{equation} \delta X= \mathrm{i}UHU^*-\mathrm{i}XH^*+UHU^*H^* -\frac{1}{2}UH^2U^*-\frac{1}{2}X{H^*}^2 \end{equation} \noindent \tag{ A.5 }$$

is the increment of X to leading order in τ. The Gaussian averages are now carried out according to the rules

$$\begin{equation} \langle k | A H B |l \rangle\langle m | C H D |n \rangle=\tau \langle k | A D |n \rangle\langle m | C B |l \rangle, \end{equation} \tag{ A.6 }$$

$$\begin{equation} \langle k | A H B |l \rangle\langle m | C H^* D |n \rangle=\tau \langle k | A C^{\mathrm{T}} | m^* \rangle\langle n^* | D^{\mathrm{T}} B |l \rangle. \end{equation} \tag{ A.7 }$$

In particular, $\overline {H^2}= N\tau$, $\overline {UHU^*H^*}=\tau UU^\dagger =\tau$, and

$$\begin{eqnarray} &&\fl \overline{\langle l | UHU^*-XH^* | k\rangle\langle k| UHU^*-XH^* | l \rangle}\nonumber\\ &&=2\tau\langle l |X| l \rangle \langle k |X| k \rangle -\tau\langle l | U^{\mathrm{T}} | k^*\rangle \langle l^* |U^*| k\rangle-\tau\langle k | U^{\mathrm{T}} | l^*\rangle \langle k^* |U^*| l\rangle\nonumber \\ &&=2\tau x_l x_k-\tau \delta_{l\bar k}(x_l+ x_{\bar l}), \end{eqnarray} \noindent \tag{ A.8 }$$

where we invoked equation (A.3). We thus obtain

$$\begin{equation} \overline{\delta x_l}= \tau -N\tau x_l -\tau {\sum_k}'\frac{2 x_lx_k-\delta_{l\bar k}(x_l+x_{\bar l}) }{x_l-x_k} . \end{equation} \tag{ A.9 }$$

Analogously, we find that

$$\begin{equation} \overline{\delta x_l\delta x_m} = \overline{\langle l |\delta X| l\rangle\langle m |\delta X| m\rangle} =-2\tau\delta_{lm}x_l^2+2\tau\delta_{l\bar m}. \end{equation} \tag{ A.10 }$$

Note that these expressions only depend on the eigenvalues. We remark that for the pinned unpaired eigenvalue x_N = 1, occurring if N is odd, these relations deliver $\overline {\delta x_N} =\overline {(\delta x_N)^2}=0$.

We now pass over to the quantities $\mu _l=(x_l+x_{\bar l})/2$, and restrict the index l such that it enumerates the pairs of eigenvalues. For even N we then find that

$$\begin{equation} \overline{\delta \mu_l} =\tau -2\tau \mu_l -2\tau(\mu_l^2-1) {\sum_k}'\frac{1}{\mu_l-\mu_k}, \end{equation} \tag{ A.11 }$$

while for odd N we have

$$\begin{equation} \overline{\delta \mu_l} =-3\tau \mu_l -2\tau(\mu_l^2-1) {\sum_k}''\frac{1 }{\mu_l-\mu_k}, \end{equation} \tag{ A.12 }$$

where the double prime excludes the pinned eigenvalue. Furthermore,

$$\begin{equation} \overline{\delta \mu_l\delta \mu_m} =2\tau(1-\mu_l^2)\delta_{lm}. \end{equation} \tag{ A.13 }$$

The stationarity condition of the associated Fokker–Planck equation (A.1) can be expressed as

$$\begin{equation}\frac{\partial}{\partial \mu_l}\overline{\delta \mu_l}P=\frac{1}{2}\frac{\partial^2}{\partial \mu_l}\overline{(\delta \mu_l)^2}P. \end{equation} \tag{ A.14 }$$

For even N = 2M, this is solved by

$$\begin{equation} P(\mu_{1},\mu_{2},\ldots,\mu_{M})\propto\prod_{k=1}^{M}\frac{1+\mu_k}{\sqrt{1-\mu_k^2}}\,\prod_{l<m=1}^{M}(\mu_l-\mu_m)^2, \end{equation} \tag{ A.15a }$$

up to a normalization constant. Each of the μ_n's (n = 1,2,...,M) is twofold degenerate. For odd N = 2M + 1 one eigenvalue is pinned at +1, and the remaining ones are twofold degenerate with the distribution

$$\begin{equation} P(\mu_{1},\mu_{2},\ldots,\mu_{M})\propto\prod_{k=1}^{M}\sqrt{1-\mu_k^2}\,\prod_{l<m=1}^{M}(\mu_l-\mu_m)^2. \end{equation} \tag{ A.15b }$$

This concludes our derivation of the eigenvalue distribution of UU* with U uniform in ${\cal U}(N)$. We have not found the result (A.15) in the literature, but there is a curious correspondence with the known [36, 38] eigenvalue distribution of orthogonal matrices (uniformly distributed according to the Haar measure). An (N + 1) × (N + 1) orthogonal matrix O with determinant −1 has one eigenvalue pinned at −1. If we exclude that eigenvalue, the remaining N eigenvalues of O have the same probability distribution as the N eigenvalues of UU*.

As an independent demonstration of this correspondence between the eigenvalue distributions of UU* and O, we have investigated the Brownian motion of orthogonal matrices. Let O be a random (N + 1) × (N + 1)-dimensional matrix in the orthogonal group, constrained to the sector Det O = −1.

The Brownian motion is induced by O(1 + A + A²/2), where (in the fixed basis) A = −A^T is a real antisymmetric matrix, with $\overline {A_{lm}^2}=\tau$. Due to the condition on the determinant, there is always one eigenvalue pinned at x_N+1 = −1, while an additional eigenvalue is pinned at x_N = 1 if N is odd. All other eigenvalues appear in pairs x_l, $x_{\bar l}$, with $|\bar l\rangle =| l^*\rangle$ (no additional factors are required).

We calculate the increments and average:

$$\begin{eqnarray} \delta x_l&=&\frac{1}{2}\langle l |O A^2| l \rangle +\sum_{k\neq l}\frac{\langle l |O A| k \rangle\langle k |O A| l \rangle}{x_l-x_k}\nonumber\\ \\ &\Rightarrow& \overline{\delta x_l}=-\frac{1}{2}\tau Nx_l+\tau\sum_{k\neq l}\frac{x_l x_k(\delta_{k\bar l}-1)}{x_l-x_k},\nonumber \end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray} \delta x_l\delta x_k&=&\langle l |O A| l \rangle\langle k |O A| k \rangle\nonumber\\ \\ &\Rightarrow &\overline{\delta x_l\delta x_k}=\tau x_l x_k (\delta_{l\bar k}-\delta_{lk})=\tau (\delta_{l\bar k}-x_l^2\delta_{lk}).\nonumber \end{eqnarray} \tag{ A.17 }$$

(Note that 〈l|A|l〉 does not vanish if |l〉 is complex, as is generally the case for the unpinned eigenvalues.)

As before, in passing to μ_l we restrict indices to enumerate different pairs. For N even, we find that (considering that the restricted sum has (N − 2)/2 terms)

$$\begin{eqnarray} \overline{\delta \mu_l}&=&\frac{1}{2}\tau-\frac{1}{2}\tau N\mu_l- \tau\sum_{k\neq l,N+1}\frac{\mu_l\mu_k-1}{\mu_l-\mu_k}\nonumber\\ \\ &=&\frac{1}{2}\tau-\tau\mu_l- \tau(\mu_l^2-1)\sum_{k\neq l,N+1}\frac{1}{\mu_l-\mu_k},\nonumber\ \end{eqnarray} \tag{ A.18 }$$

while if N is odd (where the restricted sum has (N − 3)/2 terms),

$$\begin{eqnarray} \overline{\delta \mu_l}&=&-\frac{1}{2}\tau N\mu_l- \tau\sum_{k\neq l,N,N+1}\frac{\mu_l\mu_k-1}{\mu_l-\mu_k}\nonumber\\ \\ &=&-\frac{3}{2}\tau\mu_l- \tau(\mu_l^2-1)\sum_{k\neq l,N,N+1}\frac{1}{\mu_l-\mu_k}.\nonumber \end{eqnarray} \tag{ A.19 }$$

Furthermore,

$$\begin{equation} \overline{\delta \mu_l\delta \mu_k}=\tau(1-\mu_l^2)\delta_{lk}. \end{equation} \tag{ A.20 }$$

Comparison with equations (A.11)–(A.13) shows that these are the same average increments, if we rescale τ by a factor of 2. The eigenvalues of UU* and O therefore execute the same Brownian motion process, with the same stationary solution (A.15).

In class BDI, we have X = UλU^†U*λU^T, with U uniform in ${\cal U}(N)$ and λ a fixed diagonal matrix with entries ±1 that sum up to Q. Since here the matrix X is symmetric, X = X^T, it is now diagonalized by an orthogonal transformation; thus, the eigenvectors |k〉 = |k*〉 are real. As in class D, eigenvalues appear in complex-conjugate pairs, apart from eigenvalues pinned at 1. We observe that Ω mediates between the associated eigenvector, $|\bar k\rangle =\xi _k\Omega |k\rangle =\xi _k^*\Omega ^*|k\rangle$. In order to treat the partners symmetrically we have to require that $|\bar k\rangle$ is also real, so ξ_k compensates for any complex overall factor. It then follows that 〈k|ΩΩ*|k〉 = ξ²_k = λ_k, and thus the coefficients ξ_k are related to the eigenvalues as in class D.

To identify the pinned eigenvalues note that Ω = Ω^† = Ω⁻¹ is both Hermitian and unitary, and thus has eigenvalues ±1. Let Ω_± be the eigenspace for each set of eigenvalues, and Ω*_± the analogous eigenspace for Ω*, which is spanned by the complex-conjugated vectors. We denote ξ = sign Q. The space [span(Ω_−ξ,Ω*_−ξ)]^⊥ is then of dimension q = |Q| (barring accidental degeneracies), and all of the vectors in this space obey X|k〉 = |k〉. Thus X has q = |Q| eigenvalues pinned at 1. For each pinned eigenvalue, insisting that $|\bar k\rangle =|k\rangle$ implies Ω|k〉 = Ω*|k〉 = ξ|k〉, ξ = sign Q = ± 1 (consistent with the property that these states lie in the joint subspace of Ω_ξ and Ω*_ξ).

With these additional properties in hand, the evaluation of drift and diffusion coefficients can proceed along the same steps as before. With the specified form of X, the incremental step of U carries over to an increment

$$\begin{eqnarray} \delta X&=& \mathrm{i}U[H,\lambda]U^\dagger \Omega^* -\mathrm{i}\Omega U^*[H^*,\lambda]U^{\mathrm{T}}\nonumber\\ \\ & &+\tau Q (\Omega^*+\Omega)+2\tau(1-X)-2N\tau X,\nonumber\ \end{eqnarray} \tag{ A.21 }$$

where we have already averaged terms of second order in H; in particular, terms such as $\overline {UH\lambda HU^\dagger U^*\lambda U^{\mathrm {T}}}=\tau Q \Omega ^*$ produce the topological invariant Q. The associated eigenvalue increment averages to

$$\begin{eqnarray} &&\fl \overline{\delta x_l}= -2N\tau x_l +\tau Q \langle l |\Omega^*+\Omega|l\rangle +2\tau(1-x_l)\nonumber \\ &&\fl -{\sum_k}'(x_l-x_k)^{-1}\overline{\langle l| U[H,\lambda]U^\dagger \Omega^* -\Omega U^*[H^*,\lambda]U^{\mathrm{T}}|k\rangle\langle k| U[H,\lambda]U^\dagger \Omega^* -\Omega U^*[H^*,\lambda]U^{\mathrm{T}}|l\rangle}\nonumber\\ &&\fl = -2N\tau x_l +2\tau q \delta_{l\bar l} +2\tau(1-x_l) -4\tau{\sum_k}'\frac{ x_lx_k-\delta_{l\bar l}\delta_{k\bar k} -\delta_{k\bar l}(x_l+x_{\bar l})/2 +x_lx_k\delta_{lk} }{x_l-x_k}, \end{eqnarray} \noindent \tag{ A.22 }$$

where the δ_lk term can be dropped because of the constraint k ≠ l on the sum. Note how Q changes to q = |Q| because of the sign of the matrix element involving pinned eigenvalues.

Again we find that eigenvalues at unity remain pinned. For the other eigenvalues, we separate out from the sum the q eigenvalues that are pinned, and sum over the M = (N − q)/2 pairs of unpinned eigenvalues,

$$\begin{eqnarray} \overline{\delta x_l} &=& -2N\tau x_l +2\tau(1-x_l) -2\tau\frac{2-(x_l+x_{\bar l})}{ x_l-x_{\bar l}}\nonumber\\ \\ & &-4\tau q\frac{x_l }{x_l-1} -4\tau x_l {\sum_k}'' \frac{ x_k+x_{\bar k}-2x_{\bar l} }{x_l+x_{\bar l}-x_k-x_{\bar k}},\nonumber\ \end{eqnarray} \tag{ A.23 }$$

where the double prime again indicates the exclusion of the pinned eigenvalues. Furthermore,

$$\begin{eqnarray} \overline{\delta x_l\delta x_m}&=& 8\tau(\delta_{l\bar m}-\delta_{l m}x_l^2). \end{eqnarray} \noindent \tag{ A.24 }$$

For the quantities $\mu _l=(x_l+x_{\bar l})/2$, this gives

$$\begin{equation} \overline{\delta\mu_l} = -2q\tau (\mu_l+1) +2\tau(1-3\mu_l) -4\tau {\sum_k}'' \frac{ \mu_l^2-1 }{\mu_l-\mu_k}, \end{equation} \tag{ A.25 }$$

$$\begin{equation} \overline{\delta\mu_l\delta\mu_m} =8\tau(1-\mu_l^2)\delta_{lm}. \end{equation} \tag{ A.26 }$$

The stationarity condition (A.14) is now fulfilled for

$$\begin{equation} P(\mu_{1},\mu_{2},\ldots,\,\mu_{M})\propto\prod_{k=1}^{M}(1-\mu_k)^{(q-1)/2}\,\prod_{l<m=1}^{M}|\mu_l-\mu_m|, \end{equation} \noindent \tag{ A.27 }$$

which gives the joint probability distribution of the twofold degenerate, unpinned eigenvalues μ_n (n = 1,2,...,M).

The probability distributions (A.15) and (A.27) are both of the form

$$\begin{equation} P(\mu_1,\mu_2,\ldots, \mu_M)\propto\prod_{k=1}^{M}(1+\mu_k)^{a}(1-\mu_k)^{b}\,\prod_{l<m=1}^{M}|\mu_l-\mu_m|^{\beta}, \end{equation} \tag{ A.28 }$$

with β = 2, a = 1/2, b = |Q| − 1/2 in class D and β = 1, a = 0, b = |Q|/2 − 1/2 in class BDI. These are called Jacobi distributions, because the eigenvalue density ρ(μ) can be written in terms of Jacobi polynomials [36].

For small N it is easier to calculate the eigenvalue density by integrating out all μ_n's except a single one. Keep in mind that |Q| of the μ_n's are pinned at unity and that the N − |Q| = 2M unpinned μ_n's are twofold degenerate. (The products in equation (A.28) run only over these M unpinned pairs.) The eigenvalue density $\rho (\mu )=\langle \sum _{n=1}^{N}\delta (\mu -\mu _{n})\rangle$ is then given by

$$\begin{eqnarray} \rho(\mu)&=&|Q|\delta(\mu-1)+2Mp(\mu),\nonumber\\ p(\mu)&=&\int_{-1}^{1}\mathrm{d}\mu_{1}\int_{-1}^{1}\mathrm{d}\mu_{2}\\ & &\cdots\int_{-1}^{1}\mathrm{d}\mu_{M}\,\delta(\mu-\mu_{1})P(\mu_1,\mu_2,\ldots, \mu_M).\nonumber \end{eqnarray} \tag{ A.29 }$$

The delta functions satisfy $\int _{-1}^{1}\delta (\mu \pm 1)\mathrm {d}\mu =1$. The average conductance follows from the eigenvalue density according to equation (9),

$$\begin{equation} \langle G\rangle=T^{2}G_{0}\int_{-1}^{1}\mathrm{d}\mu\,\rho(\mu)[1+(1-T)^2-2(1-T)\mu]^{-1}. \end{equation} \tag{ A.30 }$$

This gives the small-N results in equation (11) and figure 2.

The large-N limit (14) is obtained from an integral equation for the eigenvalue density in the Jacobi ensemble [5, 39],

$$\begin{eqnarray} & &M\int_{-1}^{1}\mathrm{d}\mu\,p(\mu')\ln|\mu-\mu'|=-\frac{1}{2}(1-2/\beta)\ln p(\mu)\nonumber\\ &&\quad-\frac{a}{\beta}\ln(1+\mu)-\frac{b}{\beta}\ln(1-\mu)+C+{\cal O}(1/M). \end{eqnarray} \tag{ A.31 }$$

The constant C is determined by the normalization

$$\begin{equation} \int_{-1}^{1}\mathrm{d}\mu\,p(\mu)=1. \end{equation} \tag{ A.32 }$$

The solution is

$$\begin{eqnarray} Mp(\mu)&=&\frac{\tilde{M}}{\pi\sqrt{1-\mu^{2}}}-\frac{a}{\beta}\delta(\mu+1)-\frac{b}{\beta}\delta(\mu-1)\nonumber\\ & &+\frac{1}{4}(1-2/\beta)[\delta(\mu+1)+\delta(\mu-1)]+{\cal O}(1/M), \end{eqnarray} \tag{ A.33 }$$

$$\begin{equation} \tilde{M}=M+(a+b)/\beta-\frac{1}{2}(1-2/\beta). \end{equation} \tag{ A.34 }$$

Upon substitution of the values for a,b,β into the two symmetry classes and transforming back from p to ρ, we find that

$$\begin{equation} \rho(\mu)=\frac{N}{\pi}\frac{1}{\sqrt{1-\mu^{2}}}+\frac{1}{2}\delta(\mu-1)-\frac{1}{2}\delta(\mu+1)+{\cal O}(1/M), \end{equation} \tag{ A.35 }$$

independent of Q and for both symmetry classes D and BDI. The corresponding result for the conductance is equation (14), to order 1/N if the limit N → ∞ is taken at a fixed Q.

For completeness we also give the derivation of the large-voltage limit (10) of the average conductance. We need to evaluate

$$\begin{equation} \langle \tilde{G}\rangle=T^{2}G_{0}\int_{0}^{2\pi}\mathrm{d}\phi\,\tilde{\rho}(\phi)\left|1-(1-T)\mathrm{e}^{\mathrm{i}\phi}\right|^{-2}, \end{equation} \tag{ A.36 }$$

with $\tilde {\rho }(\phi )=\langle \sum _{n}\delta (\phi -\phi _{n})\rangle$ the density on the unit circle of the eigenvalues e^iϕ_n of a unitary matrix $\tilde {\Omega }$.

In class D, the matrix $\tilde {\Omega }\equiv U$ is uniformly distributed in ${\cal U}(N)$. This is the circular unitary ensemble (CUE, β = 2). In class BDI, the chiral symmetry enforces that $\tilde {\Omega }$ is unitary symmetric, $\tilde {\Omega }=UU^{\mathrm {T}}$ with U uniform in ${\cal U}(N)$. This is the circular orthogonal ensemble (COE, β = 1). Unlike the probability distributions we needed for the zero-voltage limit, these two distributions are in the literature [36],

$$\begin{equation} P(\phi_1,\phi_2,\ldots,\phi_N)\propto\prod_{k<l=1}^{N}|\mathrm{e}^{\mathrm{i}\phi_{k}}-\mathrm{e}^{\mathrm{i}\phi_{l}}|^\beta. \end{equation} \tag{ A.37 }$$

The corresponding density

$$\begin{equation} \tilde{\rho}(\phi)=N/2\pi,\;\;0<\phi\leqslant 2\pi \end{equation} \tag{ A.38 }$$

is uniform irrespective of the value of β and without any finite-N corrections. Substitution into equation (A.37) gives the result (10).