Topological supercurrents interaction and fluctuations in the multiterminal Josephson effect

TL;DR

This paper investigates the topological properties and supercurrent fluctuations in multiterminal topological superconductor Josephson junctions, revealing protected band crossings, Majorana-Weyl nodes, and quantized conductance linked to topological invariants.

Contribution

It introduces a symmetry-constrained scattering matrix approach to analyze sub-gap states and topological features in multiterminal Josephson junctions, including novel Majorana-Weyl nodes and conductance quantization.

Findings

Identification of topologically protected band crossings as monopoles of Berry curvature.

Discovery of finite energy Majorana-Weyl nodes in four-terminal junctions.

Quantized nonlocal conductance related to the Chern number of bands.

Abstract

We study the Josephson effect in the multiterminal junction of topological superconductors. We use the symmetry-constrained scattering matrix approach to derive band dispersions of emergent sub-gap Andreev bound states in a multidimensional parameter space of superconducting phase differences. We find distinct topologically protected band crossings that serve as monopoles of finite Berry curvature. Particularly, in a four-terminal junction the admixture of $2\pi$ and $4\pi$ periodic levels leads to the appearance of finite energy Majorana-Weyl nodes. This topological regime in the junction can be characterized by a quantized nonlocal conductance that measures the Chern number of the corresponding bands. In addition, we calculate current-phase relations, variance, and cross-correlations of topological supercurrents in multiterminal contacts and discuss the universality of these transport characteristics. At the technical level these results are obtained by integrating over the group of a circular ensemble that describes the scattering matrix of the junction. We briefly discuss our results in the context of observed fluctuations of the gate dependence of the critical current in topological planar Josephson junctions and comment on the possibility of parity measurements from the switching current distributions in multiterminal Majorana junctions.

Equations

$$\hat{S}(\epsilon )=\mathcal{P}\hat{S}(-\epsilon ){\mathcal{P}}^{-1},$$

$$\mathrm{Det}[{\mathbb{I}}_{2n}-{\hat{R}}_A(\epsilon ,\hat{\theta }){\hat{S}}_N(\epsilon )]=0.$$

$${\hat{S}}_N(\epsilon )=\left[\begin{array}{cc}\hat{s}(\epsilon )& 0\\ 0& {\hat{s}}^{\ast }(-\epsilon )\end{array}\right],$$

$${\hat{R}}_A(\epsilon ,\hat{\theta })=e^{-i\arccos \epsilon }\left[\begin{array}{cc}0& e^{i\hat{\theta }}\\ -{\mathcal{P}}^2{e}^{-i\hat{\theta }}& 0\end{array}\right],$$

$$P_n(\gamma ;\hat{\theta },\epsilon )\equiv \mathrm{Det}\left[{\mathbb{I}}_n+{\mathcal{P}}^2\gamma (\epsilon )e^{i\hat{\theta }}{s}^{\ast }(-\epsilon )e^{-i\hat{\theta }}s(\epsilon )\right].$$

$${\epsilon }_{{\mathcal{P}}^2=+1}^2(\hat{\theta })+{\epsilon }_{{\mathcal{P}}^2=-1}^2(\hat{\theta })=1.$$

$${\hat{S}}_0\equiv {\hat{R}}_A(0,\hat{\theta }){\hat{S}}_N=i\left[\begin{array}{cc}0& -e^{i\hat{\theta }}{\hat{s}}^{\ast }\\ e^{-i\hat{\theta }}\hat{s}& 0\end{array}\right],$$

$$\mathrm{Det}[\hat{m}(\hat{\theta })]=0,\hat{m}(\hat{\theta })=e^{-i\hat{\theta }/2}\hat{s}{e}^{i\hat{\theta }/2}-e^{i\hat{\theta }/2}{\hat{s}}^T{e}^{-i\hat{\theta }/2}.$$

$${\mathrm{Pf}}_{n\in \mathrm{even}}[\hat{m}(\hat{\theta })]=0.$$

$$s=\left[\begin{array}{cc}\sqrt{1-T}{e}^{i{\phi }_{00}}& \sqrt{T}{e}^{i{\phi }_{01}}\\ \sqrt{T}{e}^{i{\phi }_{10}}& -\sqrt{1-T}{e}^{i({\phi }_{01}+{\phi }_{10}-{\phi }_{00})}\end{array}\right],$$

$$\epsilon (\vartheta )=\pm \{\begin{array}{ll}\sqrt{T}\cos (\vartheta /2),& {\mathcal{P}}^2=+1,\\ \sqrt{1-T{\sin }^2(\vartheta /2)},& {\mathcal{P}}^2=-1,\end{array}$$

$$J(\vartheta )=\pm \frac{e\Delta }{\hslash }\times \{\begin{array}{ll}\sqrt{T}\sin (\vartheta /2),& {\mathcal{P}}^2=+1,\\ T\sin \vartheta /4\epsilon (\vartheta ),& {\mathcal{P}}^2=-1.\end{array}$$

$${\epsilon }_{\pm }(\mathbf{\theta })=\pm \sqrt{\frac{1-B_3(\mathbf{\theta })}{2}},{\epsilon }_0(\mathbf{\theta })=0.$$

$$B_3=\frac{1}{2}[2a^2-(1+a^2)(b^2+c^2-2b^2{c}^2)$$

$$-4abc\sqrt{(1-b^2)(1-c^2)}\cos \phi ]$$

$$+bc(1-a^2)\cos {\vartheta }_1+(1-a^2)\sqrt{(1-b^2)(1-c^2)}\cos {\vartheta }_2$$

$$+\left[bc(1+a^2)\sqrt{(1-b^2)(1-c^2)}+a(b^2+c^2-2b^2{c}^2)\cos \phi \right]$$

$$\times \cos ({\vartheta }_1-{\vartheta }_2)+a(b^2-c^2)\sin \phi \sin ({\vartheta }_1-{\vartheta }_2).$$

$$\epsilon (\mathbf{\theta })=\pm \sqrt{\frac{4-A_4\pm \sqrt{A_4^2-4B_4+8}}{8}},$$

$$A_4=$$

$$B_4=$$

$$+2\sum _{jkl\in P_{123}}\Re [{\mathsf{B}}_{jkl}{e}^{-i{\theta }_{jkl}}].$$

$${\mathsf{A}}_0=\sum _{j=0}^3|s_{jj}{|}^2,{\mathsf{A}}_j=s_{0j}^{\ast }{s}_{j0},{\mathsf{A}}_{jk}=s_{kj}^{\ast }{s}_{jk},$$

$${\mathsf{B}}_0=\sum _{j<k}|{\mathsf{s}}_{jj;kk}{|}^2,{\mathsf{B}}_j=\sum _{k\ne 0,j}{\mathsf{s}}_{0k;kj}^{\ast }{\mathsf{s}}_{jk;k0},$$

$${\mathsf{B}}_{jk}=\sum _{l\ne j,k}{\mathsf{s}}_{kl;lj}^{\ast }{\mathsf{s}}_{jl;lk},{\mathsf{B}}_{jkl}={\mathsf{s}}_{0k;lj}^{\ast }{\mathsf{s}}_{k0;jl},$$

$$\mathrm{Pf}_{4}(\bm{\theta})=\sum_{jkl\in P_{123}}\big{[}\mathsf{s}_{0k;lj}e^{i\theta_{jkl}/2}+\mathsf{s}_{jl;k0}e^{-i\theta_{jkl}/2}\big{]}=0.$$

$$\rho (\vec{x})\equiv \sqrt{|\mathrm{Det}\hat{M}(\vec{x})|},M_{\mu \nu }\equiv \sum _{i,j}\frac{\partial s_{ij}}{\partial x_{\mu }}\frac{\partial s_{ij}^{\ast }}{\partial x_{\nu }}.$$

$$p(\lambda )=⟨e^{i\lambda Q(\vec{x})}⟩,P(Q)=\frac{1}{2\pi }\int d\lambda e^{-i\lambda Q}p(\lambda ),$$

$$⟨J^m⟩=$$

$$\times {}_2{F}_1\left[\frac{m}{2},m+1;m+2;{\sin }^2\left(\frac{\vartheta }{2}\right)\right],{\mathcal{P}}^2=-1,$$

$$\mathrm{Var}J=$$

$$P(J;\vartheta )=\frac{2J}{J_{\mathrm{c}}^{+}|J_{\mathrm{c}}^{+}|}\Theta \left(|J_{\mathrm{c}}^{+}|-\left|2J-J_{\mathrm{c}}^{+}\right|\right),$$

$$P(J;\vartheta )=\frac{8\Theta \left(|J_{\mathrm{c}}^{-}|-\left|2J-J_{\mathrm{c}}^{-}\right|\right)}{\left|\sin \vartheta \right|K\left(J\tan \frac{\vartheta }{2}\right)\left[1+K^2\left(J\tan \frac{\vartheta }{2}\right)\right]},$$

$$\rho (\vec{x})=\mathcal{N}ab\sqrt{\frac{(2-a^2)[(1-a^2)(2-b^2)+a^2{b}^2{\sin }^2\phi ]}{(1-a^2)(1-b^2)(1-c^2)}},$$

Full text

I Introduction

The universality of conductance fluctuations (UCF) is the hallmark of mesoscopic physics BLA-JETPLett85 ; BLA-DEK-JETPLett86 ; Lee-Stone-PRL85 ; Lee-Stone-Fukuyama-PRB87 . This phenomenon emerges from the quantum coherence of electron trajectories and is sensitive to changes in external magnetic field or gate voltage. At temperatures below the Thouless energy, $T<E_{\mathrm{Th}}$, which is related to the inverse dwell time for an electron to diffuse across the sample $E_{\mathrm{Th}}=D/L^{2}$, the root-mean-square (rms) value of conductance fluctuations saturates to the universal value of the order conductance quantum $\sim e^{2}/h$ as long as the characteristic sample size $L$ is smaller than the dephasing length $L<L_{\phi}$. Interaction effects in normal metals barely change the magnitude and universality of conductance fluctuations, although they are crucially important in determining the temperature dependence of dephasing effects and, in particular, $L_{\phi}$ Aleiner-Blanter-PRB02 . The robustness of UCF can be rooted to the random matrix theory description of Wigner-Dyson statistics of electron energy levels in disordered conductors Dyson . Indeed, in the Landauer picture of transport across a mesoscopic sample, conductance is given by $e^{2}/h$ times the number of single-particle levels within the energy strip of the width of Thouless energy. While the average number of such levels depends on the dimensionality, random matrix theory predicts that their fluctuation is universally of the order of one Altshuler-Shklovskii-JETP86 ; Mehta .

When superconductivity is induced at the boundary of the mesoscopic sample via the proximity effect, the universality of fluctuations remains intact Exp-UCF-SN-1 ; Exp-UCF-SN-2 . Indeed, the magnitude of sample-to-sample conductance fluctuations changes only by a numerical factor of the order of unity whose value depends on the underlying symmetry Brouwer-PhD ; Beenakker-RMP . Interestingly, the universality of fluctuations extends beyond conductance as it also manifests in the Josephson current of a superconducting-normal-superconducting (SNS) bridge. Indeed, extending the original ideas of Altshuler and Spivak BLA-BZS-JETP87 , who argued that random shifts of sub-gap energy levels with superconducting phase difference would alter the current, Beenakker showed Beenakker-PRL91 that in short junctions, $L\ll\xi$, where $\xi$ is the superconducting coherence length, the rms value of critical current fluctuations saturates to a universal bound $\sim e\Delta/h$ determined only by the superconducting energy gap $\Delta$ in the leads. Further a complete characterization of the supercurrent variance as a function of phase across the point contact Josephson junction was computed by Chalker and Macêdo Chalker . In long junctions, $L\gg\xi$, supercurrent fluctuations cease to be universal and scale with $\sim eE_{\mathrm{Th}}/h$. However, a remarkable property of these fluctuations is that there is a regime where the entire critical current through the junction can be determined by the mesoscopic contribution when the average current is suppressed.

In recent years the interest in Josephson physics has shifted towards junctions whose elements include topological materials Sacepe ; Veldhorst ; Bestwick ; Mason ; Finck-PRX ; Kurter-PRB ; Kurter-NC ; Sochnikov ; Stehno or where topological properties are enabled by a specific design of the hybrid junction with otherwise conventional materials Frolov ; Shtrikman ; Marcus ; Moler ; Kouwenhoven . These possibilities and advances motivate our work to investigate how universal mesoscopic effects manifest in topological Josephson junctions that, in particular, host Majorana states (see the review in Beenakker-RMP-MF and references therein). We carry out this analysis in the context of multiterminal devices that were brought into the spotlight of recent theoretical attention with the observation that they can emulate topological matter Riwar ; Eriksson ; Houzet ; HYX1 ; HYX2 ; Qi ; Nazarov ; Deb ; Meyer , which triggered experimental efforts in realizing these systems in various proximitized circuits Giazotto ; Heiblum ; Finkelstein ; Manucharyan ; Pribiag .

The rest of the paper is organized as follows. In Sec. II we briefly review symmetry-constrained scattering matrix transport formalism in application to the Josephson effect in multiterminal circuits. In Sec. III we apply these methods to two-terminal junctions as a benchmark and then extend our analysis to three- and four-terminal devices, for which we compute the emergent band structure of sub-gap states, investigate their topology, and derive transport characteristics such as transconductance and supercurrent. In Sec. IV we focus our attention on the statistical properties of topological supercurrents and obtain analytical results for variance that takes a universal form and also inherits the $4\pi$ periodicity of the Majorana Josephson effect.

II Scattering matrix formalism

Consider a Josephson junction (JJ) where $n$ superconducting (S) terminals are connected through the common normal (N) region, thus forming a multiterminal SNS contact. To keep the presentation simple, we assume that each superconducting lead is coupled by only a single conducting channel in the normal region and both time-reversal and chiral symmetries are broken (unconventional classes D and C Beenakker-RMP-MF ). Formation of the sub-gap bound states in the JJs is the result of coherent Andreev reflections that describe electron-to-hole conversion at the superconductor-normal interface. In $n$-terminal junctions an elastic scattering event at energy ${\varepsilon}$ is characterized by a scattering matrix $\hat{S}({\varepsilon})\in\mathrm{U}(2n)$, where “$2$” denotes the particle-hole degrees of freedom. In what follows we assume that all leads have the same superconducting gap $\Delta$ and normalize all energies in units of $\Delta$. The particle-hole (PH) symmetry is represented by

[TABLE]

where the antiunitary PH transform $\mathcal{P}$ falls into two categories, $\mathcal{P}^{2}=\pm 1$. For example, for $s$-wave paring, $\mathcal{P}=\hat{\tau}_{1}\mathcal{K}$ ($\mathcal{P}^{2}=+1$) in the spin-nondegenerate case, and $\mathcal{P}=i\hat{\tau}_{2}\mathcal{K}$ ($\mathcal{P}^{2}=-1$) in the spin-degenerate case, where $\hat{\tau}_{1,2,3}$ are the Pauli matrices acting in particle-hole space and $\mathcal{K}$ denotes the complex conjugation. The Andreev bound state energies are determined by the determinant equation Beenakker-PRL91

[TABLE]

Here $\hat{S}_{N}({\varepsilon})$ is the scattering matrix of the normal region, and $\hat{R}_{A}({\varepsilon},\hat{\theta})$ is the scattering matrix describing Andreev reflections, where $\theta_{\alpha}\in\{\theta_{0},\theta_{1},\cdots,\theta_{n-1}\}$ is the diagonal matrix of superconducting phases. We set $\theta_{0}=0$ owing to global gauge invariance. Due to the PH symmetry Eq. (1) these scattering matrices take the block-diagonal forms

[TABLE]

where $\hat{s}({\varepsilon})\in\mathrm{U}(n)$. The determinant in Eq. (2) simplifies further to a degree-$n$ characteristic polynomial of $\gamma({\varepsilon})\equiv e^{-2i\arccos{{\varepsilon}}}$,

[TABLE]

which is (anti)palindromic $\mathcal{P}^{2n}\gamma^{n}P_{n}(\gamma^{-1})=P_{n}(\gamma)$. Importantly, from Eq. (4) we observe that, for a fixed normal-region scattering matrix $\hat{s}$, the Andreev bands of $\mathcal{P}^{2}=\pm 1$ symmetry classes are dual via the relation

[TABLE]

Previously, we extensively discussed the $\mathcal{P}^{2}=-1$ scenario in application to three- and four-terminal junctions HYX1 ; HYX2 . In this work we primarily focus on the Andreev spectrum of $\mathcal{P}^{2}=+1$ junctions that can support zero-energy Majorana modes. In what follows, we also assume energy-independent scattering matrices $\hat{s}$ that correspond to, for example, weak links where the length of the junction is small compared to the superconducting coherence length, $L\ll\xi$, so that retardation effects of traveling quasiparticles can be neglected. We note that the existence of Majorana zero modes does not depend on this assumption.

To study the energy spectra of emergent states in junctions with $\mathcal{P}^{2}=+1$ terminals, we introduce the scattering matrix at ${\varepsilon}=0$,

[TABLE]

which belongs to the circular real ensemble since $\mathrm{Det}\hat{S}_{0}=(-1)^{n}$. Via Eq. (2) the zero-energy Majorana modes are determined by the determinant equation of an antisymmetric matrix $\hat{m}(\theta)$,

[TABLE]

From here we draw important properties. (i) For $n\in\mathrm{odd}$, Eq. (7) is generally satisfied for any scattering matrices $\hat{s}$ and phases $\hat{\theta}$. This implies that Andreev-Majorana zero modes are present at any phases and robust to elastic scattering and superconducting order parameter nonuniformity. These nondispersive flat bands do not contribute to Josephson currents. (ii) For $n\in\mathrm{even}$, the Andreev bands cross at zero energy at phases determined by the Pfaffian equation

[TABLE]

Based on our study of two- and four-terminal junctions, we conjecture that there always exists a pair of Majorana zero-modes modes on an $(n-2)$-dimensional hypersurface in the $\bm{\theta}=(\theta_{1},\cdots,\theta_{n-1})$ space described by Eq. (8). Next we reveal the energy spectrum of the junction for several concrete forms of the scattering matrix.

III Multiterminal Josephson effect

III.1 Two-terminal junctions

We first study two-terminal junctions as a benchmark. We parametrize the $2\times 2$ unitary matrix $\hat{s}$ by four independent parameters,

[TABLE]

where $T\in[0,1]$, representing the normal-region transmission and scattering phases $\varphi_{00,01,10}\in[0,2\pi]$. The sub-gap spectrum of excitations is determined by the $n=2$ characteristic polynomial (4) via the equation $P_{2}(\gamma)=\gamma^{2}+2B_{2}\gamma+1=0$, where the $B_{2}$-function takes the form $B_{2}(\theta)=1-2T\sin^{2}(\vartheta/2)$, with $\vartheta\equiv\theta-\phi+\pi$ and $\phi\equiv\varphi_{10}-\varphi_{01}$. The two branches of dispersive solutions are given by

[TABLE]

where for comparison we recall the results for the conventional $\mathcal{P}^{2}=-1$ junctions. For finite transmission $T\neq 0$, the $n=2$ Pfaffian equation (8) reduces to $\mathrm{Pf}_{2}(\theta)=m_{01}\sim\cos(\vartheta/2)=0$, so that a Majorana crossing occurs at $\vartheta=(2k+1)\pi$, with $k\in\mathbb{Z}$. The zero-temperature Josephson current $J(\theta)\equiv(2e\Delta/\hbar)\partial_{\theta}{\varepsilon}$ takes the form

[TABLE]

The typical energy dispersion and supercurrent-phase relation are shown in Fig. 1. We note that for $\phi=\pi$ doubly degenerate Majorana states emerge at $\theta=\pi$ and the energy and supercurrent exhibit $4\pi$ periodicity in $\theta$. In addition, the bound states are detached from the continuum with a minimal gap $1-\sqrt{T}$ at $\theta=0$ and $2\pi$. Equations (10) and (11) are consistent with the prior results (e.g., Ref. Kwon ).

III.2 Three-terminal junctions

For $n=3$ the spectrum of localized states is determined by the palindromic polynomial $P_{3}(\gamma)=(\gamma+1)(\gamma^{2}+2B_{3}\gamma+1)=0$ and composed of three bands,

[TABLE]

Adopting the same parametrization of the scattering matrix as in Ref. HYX1 , the $B_{3}$-function can be found in the closed analytical form

[TABLE]

Consequently, there are only six independent parameters of the scattering matrix $\{a,b,c,\varphi,\phi_{1,2}\}$ that enter the spectrum of Andreev bound states (ABS). Furthermore, scattering phases $\phi_{1,2}$ only shift the phases of the leads $\vartheta_{1,2}=\theta_{1,2}-\phi_{1,2}$.

Depending on the choice of scattering matrix parameters, we find rich behavior of the energy bands. For the special case $c=\sqrt{1-b^{2}}$ and $\phi=\pi$ the spectrum exhibits nontrivial topology, as shown in Fig. 2. Zero-energy Weyl points appear at $\vartheta_{1,2}=0$ for $b=b_{0}=1/\sqrt{2}$ [Fig. 2(b)-(d)]. As shown in Fig. 2(a), the Chern number of the corresponding band structure exhibits a sign jump $C_{12}=\mathrm{sgn}(b_{0}-b)$ for $b\to b_{0}$. We also note that the other topological phase transitions for $b\approx 0.28$ and $0.96$ are related to the gap closing/reopening at the Andreev band edge ${\varepsilon}=1$. Figure 2(e) displays Josephson currents $J_{1,2}$ in two terminals when the system is tuned to the nodal gapless states. In Fig. 2(f) the series of one-dimensional cuts in either the $\theta_{1}$ phase or $\theta_{2}$ shows how Josephson currents change as one tunes phases to the vicinity of nodal points. We observe that moving across the node currents exhibit discontinuous jumps. Note that the Majorana flat band ${\varepsilon}_{0}=0$ does not contribute to the Josephson current.

III.3 Four-terminal junctions

The energy spectrum of four-terminal junctions can host Majorana zero modes and Weyl nodes simultaneously. The four Andreev bands determined by the palindromic equation $P_{4}(\gamma)=\gamma^{4}+A_{4}\gamma^{3}+B_{4}\gamma^{2}+A_{4}\gamma+1=0$ are given explicitly by the following expressions:

[TABLE]

where the $A_{4}$- and $B_{4}$-functions are defined by

[TABLE]

Here we have used short-hand notations for phases $\theta_{jk}\equiv\theta_{j}-\theta_{k}$, and $\theta_{jkl}\equiv\theta_{j}+\theta_{k}-\theta_{l}$, permutations $P_{123}\in\{123,312,231\}$, and $\Re[\cdot]$ denotes the real part of a complex number. Additionally, parameters $\mathsf{A}$ and $\mathsf{B}$ are functions of the scattering matrix elements $\{s_{jk}\}$. Specifically,

[TABLE]

and

[TABLE]

with $\mathsf{s}_{jk;lq}\equiv s_{jk}s_{lq}-s_{jq}s_{lk}$. The scattering matrix is parametrized by sixteen real parameters as in Ref. Dita , where $\{a,b,c,d,h,f\}\in[0,1]$ and $\varphi_{00,11,01,10,02,20,03,30,12,21}\in[0,2\pi]$. An inspection of these expressions reveals that despite the fact that we need ten independent phases to parametrize the scattering matrix, only six effective angles, $\phi_{1}\equiv\varphi_{12}-\varphi_{21}$, $\phi_{2}\equiv\varphi_{13}-\varphi_{31}$, $\phi_{3}\equiv\varphi_{14}-\varphi_{41}$, $\phi_{4}\equiv\varphi_{22}-\varphi_{23}$, $\phi_{5}\equiv\varphi_{22}-\varphi_{32}$, and $\phi_{6}\equiv\varphi_{11}+\varphi_{22}-\varphi_{12}-\varphi_{21}$, affect the Andreev spectrum in Eq. (14). The zero-energy states are determined by the $n=4$ Pfaffian equation (8),

[TABLE]

Via the unitary condition of $\hat{s}$, Eq. (18) implies that $\sum_{jkl\in P_{123}}C_{jkl}\cos[(\theta_{jkl}-\zeta_{jkl})/2]=0$, where $C_{ijk}$ and $\zeta_{ijk}$ are real functions of $\{s_{jk}\}$. Most importantly, this determines a Majorana-crossing surface in $\bm{\theta}$ space given by $\theta_{jkl}=\zeta_{jkl}\pm\frac{\pi}{2}(\mathrm{mod}\,2\pi)$.

As a practical example, we study the energy bands of this model for the choice of incommensurate parameters: $a=1/4$, $b=1/\sqrt{3}$, $c=1/5$, $d=1/2$, $f=1/3$, $h=4/5$, $\phi_{1}=\phi_{2}=\pi$, $\phi_{3}=0$, $\phi_{4}=\phi_{5}=-\pi/3$, and $\phi_{6}=\pi/6$. The energy spectrum, corresponding Chern number, and Josephson currents are shown in Fig. 3. We observe that the lower bands exhibit $4\pi$-periodicity due to the Majorana crossings described by Eq. (18). Moreover, at $\theta_{3}=\theta_{3}^{\ast}\approx 0.6990$ [Fig. 3(b)] and $5.7883$ [Fig. 3(c)] finite-energy Weyl nodes form at $(\theta_{1}^{\ast},\theta_{2}^{\ast})\approx(3.7785,3.2853)$ and $(3.2046,3.0051)$, respectively, between one of the higher and lower bands. The appearance of these nodal points is signaled by a change in the Chern number [Fig. 3(a)]. At this point we comment that it was recently shown that such Majorana-Weyl crossings occur in a different model of a four-terminal junction formed between the end-states of one-dimensional topological superconductors (TS) of class D Meyer . It has been pointed out that a finite Chern number $C_{ij}$ in this regime is associated with a quantized transconductance $G_{ij}=(2e^{2}/h)C_{ij}$. We confirm this result in our scattering matrix model and remark that the extra phase transitions in Fig. 3(a) are related to gap closing/reopening at the band edge ${\varepsilon}=1$ that may not be stable since the higher bands can strongly hybridize with the continuum ${\varepsilon}>1$.

The one-dimensional cut of the spectrum in Fig. 3(b) along $\theta_{1}$ at $\theta_{2}=\theta_{2}^{\ast}\approx 3.2853$ is shown in Fig. 3(d). The Josephson currents $J_{1,2}$ as functions of $\theta_{1,2}$ corresponding to the spectrum in Fig. 3(b) are shown in Fig. 3(e), where the hedgehog-like singularities are present at the Weyl nodes $(\theta_{1}^{\ast}+2\pi,\theta_{2}^{\ast})$ and $(\theta_{1}^{\ast},\theta_{2}^{\ast}+2\pi)$. We note that the other two nodal points at $(\theta_{1}^{\ast},\theta_{2}^{\ast})$ and $(\theta_{1}^{\ast}+2\pi,\theta_{2}^{\ast}+2\pi)$ do not induce current singularities since the higher- and lower-band contributions cancel each other. This can be observed in Fig. 3(f), which displays $J_{1}(\theta_{1})$ along the cut at $\theta_{2}=\theta_{2}^{\ast}$.

IV Fluctuations in topological junctions

In the previous section we studied Josephson current for the given realization of the scattering matrix. As alluded to in the Introduction, this current is expected to display reproducible sample-to-sample fluctuations and it is thus of interest to study its statistical properties. We primarily focus on its variance and also on the cross-correlation function that can be experimentally accessed in the multiterminal devices. As is known from quantum transport theoretical approaches, statistical transport properties of phase-coherent mesoscopic systems can be conveniently computed by means of averaging over a random-matrix that describes the system. In open systems, the averaging is done over the scattering matrix and one typically considers two models of junctions: chaotic cavities or disordered contacts. The former case is more suitable for the model considered in this work. We thus follow classical works by Baranger and Mello Baranger , and Jalabert, Pichard, and Beenakker Jalabert who studied conduction through a chaotic cavity on the assumption that the scattering matrix is uniformly distributed in the unitary group, restricted only by symmetry. This is the circular ensemble, introduced by Dyson, and shown to apply to a chaotic cavity by Blumel and Smilansky Smilansky . In other words we consider a SNS junction where the normal region is a chaotic quantum dot Whisler .

The probability density $\rho(\vec{x})$, an invariant Haar measure, of the $\hat{s}$-matrix parameters $\vec{x}$ is given by

[TABLE]

The distribution function of an observable $Q(\vec{x})$, defined as $P(Q)\equiv\int\,d\vec{x}\rho(\vec{x})\delta(Q-Q(\vec{x}))$, is in practice calculated by the characteristic function

[TABLE]

where $\langle\cdots\rangle\equiv\int\!d\vec{x}\rho(\vec{x})(\cdots)$, denoting the circular unitary ensemble (CUE) average.

For a benchmark we first study the statistics of the Josephson current $J(\vartheta)$ in the two-terminal junctions and take $2\Delta/\hbar$ as the units of $J$ in the following discussion. From Eqs. (9) and  (19) we obtain a constant invariant measure $\rho(T)=1$. We recall that this simplicity is specific to the unitary case; for instance, in orthogonal symmetry the probability density is not flat in $T$ even for a single-channel limit. All the moments as well as the distribution function of the Josephson current can be obtained analytically. The $m$-moment is given by the expression

[TABLE]

where $m\in\mathbb{N}$ and ${}_{2}F_{1}(a,b;c;z)$ is the hypergeometric function. Therefore, the variance for $\mathcal{P}^{2}=+1$ reads

[TABLE]

as depicted in Fig. 4(a). In the topological regime the variance inherits $4\pi$ periodicity and has a remarkably simple form. In the non-topological regime, our result is similar to that of Chalker and Macêdo Chalker , albeit with different numerical coefficients as they considered the multi-mode disordered junction model where averaging is done over the Dorokhov distribution of transmission eigenvalues. Finally, the Josephson-current distribution function takes the form, for $\mathcal{P}^{2}=+1$,

[TABLE]

where $J_{\mathrm{c}}^{+}(\vartheta)=\sin(\vartheta/2)$ and $J_{\mathrm{c}}^{-}(\vartheta)=\sin(\vartheta/2)\,\mathrm{sgn}[\cos(\vartheta/2)]/2$ are the critical currents, $\Theta(x)$ is the Heaviside step function, and the function $K(x)=|x|+\sqrt{1+x^{2}}$. As shown in Fig. 4(b), in both $\mathcal{P}^{2}=\pm 1$ classes the relation $P(J;-\vartheta)=P(-J;\vartheta)$ is satisfied. (i) For $\mathcal{P}^{2}=+1$, $P(J)$ is a linear function of $J$ for which $P(0)=0$ and the slope is defined by the phase $\vartheta$. In particular $P(J)=\delta(J)$ for $\vartheta=2k\pi$ with $k\in\mathbb{Z}$. (ii) For $\mathcal{P}^{2}=-1$, $P(J)$ is smaller for a larger current amplitude and $P(J)=\delta(J)$ for $\vartheta=k\pi$ with $k\in\mathbb{Z}$.

We proceed to study the Josephson-current statistics of the three-terminal junctions. From Eq. (19) we obtain the probability density of the effective parameters $\vec{x}=(a,b,c,\varphi)$,

[TABLE]

where $\mathcal{N}\approx 1.7671\times 10^{-2}$ is the normalization constant. The numerical results of the expectation value, variance, and covariance of $J_{1,2}(\vartheta_{1,2})$ are shown in Fig. 5, where the covariance is defined as $\mathrm{Cov}J_{1,2}\equiv\langle J_{1}J_{2}\rangle-\langle J_{1}\rangle\langle J_{2}\rangle$. The general relations $\langle J_{1}^{m}(\vartheta_{1},\vartheta_{2})\rangle=\langle J_{2}^{m}(\vartheta_{2},\vartheta_{1})\rangle$ and $\mathrm{Cov}J_{1,2}(\vartheta_{1},\vartheta_{2})=\mathrm{Cov}J_{1,2}(\vartheta_{2},\vartheta_{1})$ are satisfied. We observe that the variances and covariances distinguish the $\mathcal{P}^{2}=\pm 1$ junctions.

For $\mathcal{P}^{2}=+1$ three-terminal junctions, by integrating one and two terminals we can construct an effective two-terminal S-TS junction Ioselevich which supports a single channel on one lead (topological), with phase $\theta_{0}=0$, and two channels on the other (conventional), with phase $\theta=\theta_{1}=\theta_{2}$. Defining $\vartheta=\vartheta_{1}$ and $\phi=\phi_{2}-\phi_{1}$ in Eqs. (12) and (13), we obtain the Josephson current through the two leads,

[TABLE]

which depends on five independent parameters $\vec{x}^{\prime}=(a,b,c,\varphi,\phi)$. We present the numerical results of the statistical properties of $J(\vartheta)$ for such a configuration in Fig. 6. The expectation and variance of $J$ as functions of $\vartheta$ are shown in Fig. 6(a). The characteristic function $p(\lambda;\vartheta)$ for various $\vartheta$ is shown in Fig. 6(b). We observe that, for $\lambda\gg 1$, $p(\lambda)\propto\lambda^{-\alpha}$ with the exponent $\alpha$ being $\vartheta$-independent. We calculate the $\vartheta$-averaged characteristic function $\bar{p}(\lambda)\equiv\frac{1}{2\pi}\int_{0}^{2\pi}\!d\vartheta p(\lambda;\vartheta)$ and fix the exponent $\alpha\sim 0.55$. In particular, this analysis enables us to extract the asymptotic behavior of the full distribution function $P(J)$ which is found to exhibit a universal power-law scaling $P(J)\propto J^{-(1-\alpha)}$ in the limit $J\ll 1$.

V Discussion and outlook

In this work we applied methods of scattering matrix theory to study the transport properties of multiterminal Josephson junctions of topological superconductors. We have examined the spectrum of sub-gap states in two-, three-, and four-terminal configurations and determined that the texture of resulting Andreev bands in the multidimensional parameter space of superconducting phases can produce nonvanishing fluxes of the Berry curvature. These properties translate into the quantized nonlocal conductances of these devices. We have also studied the current-phase relationships and interaction of supercurrents, as well as their mesoscopic statistical properties. In particular, we discussed the universal regime of current fluctuations and computed supercurrent variance as well as the current cross-correlation function in the topological regime. We close this work with a few comments in relation to existing and possible future experiments in which the fundamental physics of multiterminal Josephsonic devices could be further explored.

Recently, Josephson supercurrent and conductance were measured as a function of geometry, temperature, and gate voltage in proximitized planar junction devices composed of superconductors and surface states of a topological insulator (S-TI-S junctions) in order to determine the nature of the electronic transport in these systems. The supercurrent was found to exhibit a sharp drop as a function of gate voltage, see Figs. (2) and (6) of Ref. Kurter-PRB , superimposed with reproducible noise whose magnitude was a fraction of the critical current. The systematic trend in the critical current dependence was explained by a mechanism related to the relocation of the topological surface state with respect to trivial conducting two-dimensional states formed by band-banding near the surface. In real space, a negative gating potential pushes the trivial state below the topological surface states, exposing the topological state to the disordered surface of the TI. As a result, the magnitude of the supercurrent changes sharply. The noise was attributed to the percolation effects as near the voltage threshold it is likely that local charge fluctuations cause the path of the supercurrent to be highly meandering. We wish to point out that there is possibly an alternative picture as this noise could be of mesoscopic origin. This evidence is further supported by observed similar reproducible noise features in Fraunhofer magneto-oscillations of the critical current. While our model is not directly applicable to S-TI-S junctions we make the observation that the magnitude of current fluctuations is consistent with the expectations that disorder scattering causes observed mesoscopic effects.

In addition, we wish to nore that related statistical properties of supercurrents can also be studied by measuring switching current distributions. In particular, for topological Josephson devices, the critical current measurements can potentially enable determining the parity state of a Majorana fermion (pair) in a junction since the supercurrent acquires an anomalous fractional component due to Majorana modes, $\pm\sin(\vartheta/2)$, where the sign encodes the parity. The typical switching measurement is performed by ramping the bias current through the junction to detect the current value at which the junction jumps to the finite voltage state. By repeating this protocol many times and accumulating statistics of random supercurrent switching events (as previously successfully implemented in various mesoscopic proximity circuits, e.g. nanowires and graphene layers Coskun ; Murphy ), one expects to reveal a bimodal distribution indicating the two parity states. If the separation of the two peaks in this distribution is wide enough, one can detect (with some fidelity) the parity state. Mutiterminal devices considered in this work can provide an actual hardware platform to conduct such experiments and our transport theory will be useful in modeling future measurements. In particular, knowledge of the current-phase relationship is needed for determining the energy barrier of a phase slip that triggers the switching. Furthermore, these developments are also inspired by the potential application of multiterminal devices in the design of protected superconducting qubits.

Acknowledgment

We would like to thank Dale Van Harlingen for the discussions, especially on the topic of the gate dependence of the critical current fluctuations in topological Josephson junctions Kurter-PRB and the feasibility of switching current experiments in the tri-junction configuration. The work of H.-Y. X was supported by National Science Foundation Grants No. DMR-1653661 and DMR-1743986. The work of A.L. was supported in part by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DESC0017888, and the Army Research Office, Laboratories for Physical Sciences Grant No. W911NF-18-1-0115. H.-Y. X. acknowledges the hospitality of the Kavli Institute of Theoretical Sciences (KITS) where this work was partially done.