Topological phase diagram of a three-terminal Josephson junction: From the conventional to the Majorana regime

TL;DR

This paper investigates how the transconductance in three-terminal Josephson junctions changes during a topological phase transition from s-wave to p-wave (Majorana) induced by an in-plane magnetic field, revealing fractional quantization due to Majorana quasiparticles.

Contribution

It provides a comprehensive analysis of the topological transition in Josephson junctions, highlighting the fractionalization of transconductance quantization caused by Majorana states.

Findings

Transconductance quantization becomes fractional in the Majorana regime.

The transition is characterized as a function of magnetic field and flux.

Majorana quasiparticles lead to effective half-Cooper pair charge pumping.

Abstract

We study the evolution of averaged transconductances in three-terminal Josephson junctions when the superconducting leads are led throughout a topological phase transition from an $s$-wave to a $p$-wave (Majorana) phase by an in-plane magnetic field $B_x$. We provide a complete description of this transition as a function of $B_x$ and a magnetic flux $\Phi$ threading the junction. For that we use a spinful model within a formalism that allows to treat on equal footing the contribution to the transconductance from both the Andreev subgap levels and the continuum spectrum. We unveil a fractionalization in the quantization of the transconductance due to the presence of Majorana quasiparticles, reflecting the effective pumping of half a Cooper pair charge in the $p$-wave regime.

Figures

Equations

$$\hat{H}_{\nu}\!=\!\frac{1}{2}\!\sum_{j=0}^{N-1}\!\hat{\Psi}^{\nu\dagger}_{j}\mathcal{H}_{\nu}\hat{\Psi}^{\nu}_{j}+\!\frac{1}{2}\!\sum_{j=0}^{N-2}\!\Big{[}\hat{\Psi}^{\nu\dagger}_{j}T_{\nu}\hat{\Psi}^{\nu}_{j+1}+\hat{\Psi}^{\nu\dagger}_{j+1}T_{\nu}^{\dagger}\hat{\Psi}^{\nu}_{j}\Big{]},$$

$${\mathcal{H}}_{\nu }$$

$$T_{\nu }$$

$${\hat{H}}_{\mathcal{T}}$$

$${\mathcal{T}}_{\nu {\nu }^{\prime }}$$

$$G_{\nu {\nu }^{\prime }}=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }d{\phi }_1{\int }_0^{2\pi }d{\phi }_2\frac{\partial ⟨{\mathcal{J}}_{\nu }⟩}{\partial V_{{\nu }^{\prime }}}=-\frac{2e^2}{h}{W}_{\nu {\nu }^{\prime }}$$

$$\displaystyle W_{\nu\nu^{\prime}}\!\!=\!\!\frac{1}{8\pi^{2}}\!\!\!\!\int\!\!\!\!d^{2}\varphi\!\!\int\!\!\!d\omega\!\mathrm{Tr}\Big{[}\mathcal{G}^{-1}\!\frac{\partial\mathcal{G}}{\partial{\varphi_{\nu}}}\!\cdot\mathcal{G}^{-1}\!\frac{\partial\mathcal{G}}{\partial{\varphi_{\nu^{\prime}}}}\cdot\!\mathcal{G}^{-1}\!\frac{\partial\mathcal{G}}{\partial{\omega}}-\!\nu\leftrightarrow\nu^{\prime}\!\Big{]},$$

$$H_T^{\mathrm{eff}}({\phi }_1,{\phi }_2)=i{t}_{12}{\hat{\gamma }}_1{\hat{\gamma }}_2+i{t}_{23}{\hat{\gamma }}_2{\hat{\gamma }}_3+i{t}_{31}{\hat{\gamma }}_3{\hat{\gamma }}_1,$$

$$C_{-}\!=\!\frac{1}{8\pi}\!\int_{0}^{4\pi}\!d\varphi_{1}\!\int_{0}^{4\pi}\!\!d\varphi_{2}\Big{[}i\langle\partial_{\varphi_{1}}\psi_{-}|\partial_{\varphi_{2}}\psi_{-}\rangle-(\varphi_{1}\!\leftrightarrow\!\varphi_{2})\Big{]}$$

$$[i{\partial }_t-H(t)]\hat{\mathcal{G}}(t,t^{\prime })=\delta (t-t^{\prime }),$$

$$\hat{\mathcal{G}}(t,t^{\prime })$$

$$\tilde{\mathcal{G}}(\omega ,t_{\text{av}})$$

$$[\omega +\frac{i}{2}{\partial }_{t_{\text{av}}}-e^{\frac{1}{2i}{\partial }_{\omega }^{\mathcal{G}}{\partial }_{t_{\text{av}}}^H}H(t_{\text{av}})]\tilde{\mathcal{G}}(\omega ,t_{\text{av}})=\mathcal{I},$$

$$\displaystyle\int_{-\infty}^{\infty}e^{i\omega t_{\text{rel}}}H(t_{\text{av}}+t_{\text{rel}}/2)\hat{\mathcal{G}}\big{(}t_{\text{av}}+\frac{t_{\text{rel}}}{2},t_{\text{av}}-\frac{t_{\text{rel}}}{2}\big{)}dt_{\text{rel}}=$$

$$\displaystyle\int_{-\infty}^{\infty}e^{i\omega t_{\text{rel}}}e^{\frac{t_{\text{rel}}}{2}\partial_{t_{\text{av}}}^{H}}H(t_{\text{av}})\hat{\mathcal{G}}\big{(}t_{\text{av}}+\frac{t_{\text{rel}}}{2},t_{\text{av}}-\frac{t_{\text{rel}}}{2}\big{)}dt_{\text{rel}}=$$

$$\displaystyle\int_{-\infty}^{\infty}e^{\frac{1}{2i}\partial_{\omega}\partial_{t_{\text{av}}}^{H}}e^{i\omega t_{\text{rel}}}H(t_{\text{av}})\hat{\mathcal{G}}\big{(}t_{\text{av}}+\frac{t_{\text{rel}}}{2},t_{\text{av}}-\frac{t_{\text{rel}}}{2}\big{)}dt_{\text{rel}}=$$

$$=e^{\frac{1}{2i}{\partial }_{\omega }^{\mathcal{G}}{\partial }_{t_{\text{av}}}^H}H(t_{\text{av}})\tilde{\mathcal{G}}(\omega ,t_{\text{av}}).$$

$$\Big{(}\omega+\frac{i}{2}\partial_{t_{\text{av}}}-H(t_{\text{av}})-\frac{1}{2i}\frac{\partial H(t_{\text{av}})}{\partial t_{\text{av}}}\partial_{\omega}\Big{)}\widetilde{\mathcal{G}}(\omega,t_{\text{av}})=\mathcal{I}.$$

$$\displaystyle\Big{(}\omega-H\Big{)}\widetilde{\mathcal{G}}^{(0)}$$

$$\displaystyle\Big{(}\omega-H\Big{)}\widetilde{\mathcal{G}}^{(1)}$$

$${\tilde{\mathcal{G}}}^{(0)}$$

$${\tilde{\mathcal{G}}}^{(1)}$$

$${\partial }_{\omega }{\tilde{\mathcal{G}}}^{(0)}$$

$${\partial }_{t_{\text{av}}}{\tilde{\mathcal{G}}}^{(0)}$$

$$⟨{\varphi }_0|\hat{\mathcal{O}}(t)|{\varphi }_0⟩$$

$$⟨\mathcal{O}(t)⟩$$

$$⟨{\mathcal{J}}_{\nu }^{(1)}(t)⟩$$

$$G_{\nu {\nu }^{\prime }}=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }d{\phi }_1{\int }_0^{2\pi }d{\phi }_2\frac{\partial ⟨{\mathcal{J}}_{\nu }^{(1)}⟩}{\partial V_{{\nu }^{\prime }}},$$

$$G_{\nu {\nu }^{\prime }}$$

$$\hat{H}\!=\!\frac{1}{2}\!\sum_{j=0}^{N-1}\!\hat{\Psi}^{\dagger}_{j}\mathcal{H}\hat{\Psi}_{j}+\!\frac{1}{2}\!\sum_{j=0}^{N-2}\!\Big{[}\hat{\Psi}^{\dagger}_{j}T\hat{\Psi}_{j+1}+\hat{\Psi}^{\dagger}_{j+1}T^{\dagger}\hat{\Psi}_{j}\Big{]},$$

$$\mathcal{H}$$

$$T$$

$$g_n^L(\omega )$$

$$g_n^B(\omega )$$

$$g_n^R(\omega )$$

$$T_n$$

Full text

I Introduction.

The quest for finding novel topological properties in condensed matter systems has drawn increasing attention both in theoretical and experimental research. Much of the progress in the field has relied on the idea of the emergence of topological features by engineering heterostructures of conventional materials Alicea (2012); Aguado (2017); Sato and Ando (2017); Ozawa et al. (2019), driving them out of equilibrium Oka and Aoki (2009); Kitagawa et al. (2010); Lindner et al. (2011) or by the direct design of topological phases with cold atoms systems Bloch et al. (2012); Jotzu et al. (2014); Aidelsburger et al. (2014); Cooper et al. (2019). In this way, it is now possible to artificially mimic topological matter without the need of working with exotic materials per se.

Following this line of thought, multi-terminal Josephson junctions of at least four conventional superconductors coupled through a normal scattering region have been put forward as an exciting platform to engineer Weyl-topological phases (Riwar et al., 2016). The existence of Weyl singularities in the spectrum of multiply connected networks can be traced back to the seminal work by Avron et al Avron et al. (1988), who envisioned the possibility of non trivial quantization of averaged transconductances (ATCs) in these devices. Different Josephson junction (JJ) circuits were also suggested as practical realizations to experimentally probe topologically protected adiabatic transport in these artificial setups Avron and Cross (1989).

More recently, it was shown that three terminal structures can also develop topological properties when piercing the central scattering region with an external magnetic flux Meyer and Houzet (2017); Xie et al. (2017) or when driven with an external microwave field Peralta Gavensky et al. (2018). There have also been proposals to measure the Fubini-Study geometric tensor Klees et al. (2018) and the Berry curvature Peralta Gavensky et al. (2018); Klees et al. (2018) of a multi-terminal JJ using microwaves. Remarkably, all this theoretical activity has been fuelled by several stimulating experiments on multi-terminal superconducting devices Strambini et al. (2016); Vischi et al. (2017); Draelos et al. (2019); Pankratova et al. (2018).

In a multi-terminal JJ, ATCs can be experimentally probed by voltage biasing one of the superconducting reservoirs as a function of a control parameter and measuring the phase-averaged current in a different lead. For conventional $s$-wave pairing, this would lead to a transconductance quantization in units of $4e^{2}/h$ Eriksson et al. (2017). This outcome can be simply thought of as a non-trivial adiabatic charge pumping between leads, say $\nu$ and $\nu^{\prime}$, a product of the existence of Weyl sources in the spectrum Riwar et al. (2016); Avron et al. (1988). For vanishing voltage $V$, the ATC results in $G_{\nu\nu^{\prime}}=Q_{P}/(V\Delta T)=n_{\nu\nu^{\prime}}2eq^{\star}/h$, where the pumping period $\Delta T$ is merely the Josephson period $h/2eV$ and the phase-averaged pumped charge $Q_{P}$ over this period is an integer multiple ($n_{\nu\nu^{\prime}}\,\epsilon\,\mathbb{Z}$) of the effective charge of the Cooper pair $q^{\star}=2e$.

A natural follow up question is whether multi-terminal JJs with $p$-wave pairing that host Majorana fermions may also develop a quantized ATC. Interestingly enough, this would provide a robust measurement of the effective tunneling charge in the topological superconductors, which is expected to be $q^{\star}=e$ as in the case of the $4\pi$-periodic fractional Josephson effect  Oreg et al. (2010); Lutchyn et al. (2010); Wiedenmann et al. (2016); Laroche et al. (2019). A partial answer using the Kitaev spinless model was recently found in Ref. [Houzet and Meyer, 2018] for the case of a four-terminal device (see also [Xie and Levchenko, 2019]). Yet, a complete description of the transition from the $s$-wave to the $p$-wave regime within a spinful approach in the entire parameter space is still missing. This is a necessary condition to obtain a unified topological phase diagram that could be experimentally tested. Even more, the calculation of the topological invariant associated with the ATC has so far been mainly based on the Chern numbers of the Andreev bound states (ABS) Riwar et al. (2016); Xie et al. (2017); Xie and Levchenko (2019); Houzet and Meyer (2018); Deb et al. (2018) despite the fact that in many regions of the parameter space they cannot be disentangled from the continuum spectrum above the superconducting gap (which might also contribute to the ATC Meyer and Houzet (2017)).

In this work, we present the complete topological phase diagram of a three-terminal JJ obtained by a direct calculation of the adiabatic corrections to the Josephson current. We use the general approach developed in Ref. [Peralta Gavensky et al., 2018] which includes the contributions from both the ABS and the continuum spectrum. To describe the s-wave to p-wave transition, we consider a junction formed by spinful one-dimensional ($1D$) wires with Rashba spin-orbit coupling, proximity induced superconductivity and an in plane magnetic field $B_{x}$ (see Fig. 1). This model, introduced in Refs. [Oreg et al., 2010] and [Lutchyn et al., 2010], is known to have a topological phase transition that allows Majorana fermions to develop at the end of the wires for $B_{x}>B_{c}$, with $B_{c}$ a critical field. An external magnetic flux $\Phi$ threading the central region of the device allows a fine tuning of the topological properties. For $B_{x}<B_{c}$ ($s$-wave phase) we recover the usual quantization ($4e^{2}/h$) whenever the system is gapped, while for $B_{x}>B_{c}$ ($p$-wave phase) we find a quantization in $2e^{2}/h$ that we identify with the development of effective spinless-like degrees of freedom, which make the measurement of a fractionalized pumped charge ($q^{\star}=e$) possible. All these results are summarized in Fig. 2.

II Model Hamiltonian.

The Hamiltonian of the Josephson tri-junction is given by $\hat{H}=\sum_{\nu=1}^{3}\hat{H}_{\nu}+\hat{H}_{\mathcal{T}}$, where $\hat{H}_{\nu}$ is the Hamiltonian of the lead $\nu$, and $\hat{H}_{\mathcal{T}}$ describes the Josephson tunneling. We use a spatially discretized version of the continuum model Oreg et al. (2010); Lutchyn et al. (2010) with $N\rightarrow\infty$ lattice sites to describe the leads:

[TABLE]

with the four-component spinor in site $j$ defined as ${\hat{\Psi}^{\nu}_{j}}{}=(c_{j\uparrow}^{\nu},c_{j\downarrow}^{\nu},c_{j\downarrow}^{\nu\dagger},-c_{j\uparrow}^{\nu\dagger})^{\mathrm{T}}$ in the Nambu representation and

[TABLE]

where $t$ is the hopping matrix element, $\mu$ is the chemical potential, $\Delta$ is the superconducting order parameter, $\alpha$ is the spin-orbit coupling, and $B_{x}$ is the Zeeman field not . The Pauli matrices $\tau_{a}$ ($\sigma_{a}$) and the identity $\tau_{0}$ ($\sigma_{0}$) act in particle-hole (spin) space. The junction’s tunneling Hamiltonian $\hat{H}_{\mathcal{T}}$ couples the leads only through their end site. To simplify the notation, we redefine $\hat{\Psi}_{N-1}^{1}\equiv\hat{\Psi}^{1}$, $\hat{\Psi}_{0}^{2}\equiv\hat{\Psi}^{2}$ and $\hat{\Psi}_{N-1}^{3}\equiv\hat{\Psi}^{3}$ and then write

[TABLE]

with $\mathcal{T}_{\nu\nu^{\prime}}^{\dagger}=\mathcal{T}_{\nu^{\prime}\nu}$ and $\epsilon_{\nu\nu^{\prime}}$ the totally antisymmetric tensor such that $\epsilon_{12}=\epsilon_{23}=\epsilon_{31}=1$. The dimensionless parameter $\Phi$ describes a magnetic flux threading the junction (in units of the flux quantum $\Phi_{0}=h/e$) and $U$ is an on-site surface energy not .

III Topological averaged transconductance.

We are interested in computing the current $\langle\mathcal{J}_{\nu}\rangle$ flowing in lead $\nu$ when driving adiabatically the phase of a different lead $\nu^{\prime}$ with a time dependence given by Josephson’s relation $\dot{\varphi}_{\nu^{\prime}}=2eV_{\nu^{\prime}}/\hbar$. The adiabatic correction to the current expectation value can be calculated using Wigner’s representation of the two-time out-of-equilibrium Green’s functions and performing a gradient expansion in the slow time scale (see Appendix A for details). The ATC results in

[TABLE]

where

[TABLE]

is the winding number associated to the Green’s function $\mathcal{G}$, which is topologically quantized whenever there is a gap at the Fermi energy. The trace is performed on both spin and particle-hole spaces. In the above expressions, $\mathcal{G}(\omega,\{\varphi_{\nu}\})$ is the Feynman Green’s function of the entire Hamiltonian in the total adiabatic limit—it follows the perturbation, $\{\varphi_{\nu}(t)\}$, instantaneously at time $t$—. This expression can be manipulated so as to be written in terms of the projected junction’s Green function involving only the end sites of the nanowires ($\mathcal{G}_{J}$), see Appendix A. We compute $\mathcal{G}_{J}$ numerically by a simple decimation procedure of the Hamiltonian defined in Eq. (1) (see Appendix B). Eqs. (4) and (5), originally derived in Ref. [Peralta Gavensky et al., 2018], generalize the ones derived in Ref. [Riwar et al., 2016] since they provide, on equal footing, the contribution to the ATC of both the ABS and the continuum states above the superconducting gap. In a more recent article, a similar expression in terms of the scattering matrix was obtained (Repin et al., 2019). There are two points worth mentioning: (i) If both spin projections are degenerate, then the trace over these variables would simply result in a factor $2$, which means a transconductance quantization in units of $4e^{2}/h$, and (ii) the Green’s function formalism implicitly assumes a thermodynamic average (taken here at zero temperature), which means that if there is a band crossing of the ABS at zero energy, the transconductance will not necessarily be quantized.

III.1 Topological phase diagram.

Figure 2 summarizes the main result of this work: (a) a color map of the thermal transconductance as a function of the external flux $\Phi$ and the in-plane Zeeman field $B_{x}$; (b) representative cuts of this color map for fixed $B_{x}$. For $B_{x}=0$ both spin projections are degenerate and four zero energy Weyl points develop as a function of the external flux at $\Phi_{c}^{(i)}$ for $i=1,2,3,4$, marked with arrows in Fig. 2(b), making the ATC jump at these singular points between quantized values of $\pm 4e^{2}/h$ and zero. As $B_{x}$ increases, the quantization does not persist at all flux values, which can only be due to band overlaps developing in the spectrum at zero energy, as we shall show below. At this point, we would like to emphasize that in this entire region ($B_{x}<B_{c}$) the independent contributions to the ATC of both the ABS and the continuum spectrum are non trivial and that our method adds them up to give the correct value. In some cases, it is possible to separate them by resorting to the so-called topological Hamiltonian, defined as $H_{\mathrm{top}}^{J}(\varphi_{1},\varphi_{2})\equiv-\mathcal{G}_{J}^{-1}(\omega=0)$, which has been proven (Wang and Zhang, 2012; Wang and Yan, 2013) to have the same topological properties as the system described by the entire $\mathcal{G}_{J}(\omega)$. Indeed, the Chern numbers calculated from the spectrum of $H_{\mathrm{top}}^{J}$ do agree (see Appendix C) with the results presented in Fig. 2 and, for the particular case of $B_{x}=0$, with those obtained in Ref. [Meyer and Houzet, 2017].

The most striking result shown in Fig. 2 is manifested for $B_{x}>B_{c}$, where the ATC only takes quantized values in $2e^{2}/h$. For these fields, Majorana degrees of freedom are expected to emerge in the system, with an effectively spinless superconductivity that causes the effective charge $q^{\star}$ to be that of the single electron.

III.2 Junction’s spectral density.

In order to understand these results, we show in Fig. 3 the spectral density of the tri-junction $\mathcal{A}(\omega)=-1/\pi\,\mathrm{Im}[\mathrm{Tr}(\mathcal{G}^{r}_{J}(\omega))]$ along the path $\varphi_{2}=-\varphi_{1}$ for different values of $B_{x}$, where $\mathcal{G}^{r}_{J}(\omega)$ stands for the retarded junction Green’s function. The external magnetic flux was taken to be $\Phi_{c}^{(1)}=0.53$, which corresponds to the first critical flux where a Weyl point appears in the ABS spectrum at zero magnetic field. As $B_{x}$ increases, the spin degeneracy is lifted and the ABSs overlap–as discussed in Ref. [Yokoyama and Nazarov, 2015]–causing the ATC to be not quantized. At the critical field $B_{c}=\sqrt{\Delta^{2}+\mu^{2}}$ the superconducting gap closes, and for $B_{x}>B_{c}$ a zero energy flat band emerges. The existence of this dispersionless state can be traced, as discussed below, to the fact that in the $p$-wave regime, Majorana modes $\hat{\gamma}_{\nu}$ develop at the end of each lead $\nu$. Only two effectively spinless dispersive ABSs coexist inside the gap with the Majorana band. Notice that in this three-terminal setup, no further requirement beyond thermal occupation is needed in order to observe the $2e^{2}/h$ quantization since for any flux aside from the critical ones the two ABSs never touch in the Majorana phase.

In Fig. 4 we show the spectral density along a different path ($\varphi_{1}=2\pi$) as a function of the Zeeman field for no external flux ($\Phi=0$). Here, it becomes clear that the jump in the ATC from $-2e^{2}/h$ to $2e^{2}/h$ for $B_{x}>B_{c}$ can be traced back to the emergence of a pseudo-spin $1$ Dirac-Weyl point, as it is apparent in the figure for $B_{x}=1.5\,B_{c}$. In this condition, when the external magnetic flux is turned on, it breaks the degeneracy between the two ABS, making them detach from the zero energy (Majorana) solution.

We notice that there are situations in which the ABSs are not detached from the continuum but they merge into it in the form of resonances (see Figs. 3 and 4)—this is so even in the Majorana regime (not shown). These mergings do not change the Green’s function winding number since when two states become degenerate the Chern number of each one may change but their sum remains constant Avron et al. (1983). Hence, all changes in the quantized ATC presented in the topological phase diagram of Fig. 2 are due to the occurrence of crossings or Weyl points in the Andreev spectrum at zero energy.

In the $p$-wave phase ($B_{x}>B_{c}$), the topological Hamiltonian $H_{\mathrm{top}}$ is not well defined due to the presence of the Majorana bound state. However, for $B_{x}\gg B_{c}$, i.e. deep into the Kitaev regime, a simple low-energy model helps to describe and analyze the results. In this case, the tunneling Hamiltonian takes the form

[TABLE]

with $t_{12}=\tilde{\Gamma}_{12}\cos((\varphi_{1}-\varphi_{2})/2-\Phi/3)$, $t_{23}=-\tilde{\Gamma}_{23}\cos(\varphi_{2}/2-\Phi/3)$, $t_{31}=-\tilde{\Gamma}_{13}\sin(\varphi_{1}/2+\Phi/3)$ and $\tilde{\Gamma}_{ij}=\Gamma_{ij}/2$. It possesses two particle-hole symmetric eigenvectors $|\psi_{\pm}\rangle$ with energies $E_{\pm}=\pm\sqrt{t_{12}^{2}+t_{13}^{2}+t_{23}^{2}}$, which correspond to a finite energy fermion (two Majorana degrees of freedom are paired up), and one zero-energy solution (the flat band) that is merely an “unpaired” Majorana mode living in the junction. Within this approximation it is much easier to show that all the winding comes from the finite energy fermionic band whose Chern number

[TABLE]

changes from $-1$ to $1$ at $\Phi=0$ and $\Phi=\pi$. In this regime, the continuum does not contribute. To further check these results, we have also made use of an analytic expression for the boundary Green’s function for a Kitaev chain in the wide band limit Zazunov et al. (2016) to calculate the winding number $W_{\nu\nu^{\prime}}$. Interestingly, we do recover the same results as in the spinful model, even though in this case, the continuum has a smoother dependence in energy than in the complete model.

IV Conclusions.

By expressing the average transconductance of a multi-terminal JJ in terms of the Green’s function of the junction, we have calculated the topological phase diagram of a three-terminal junction that completely describes the transition from the $s$-wave to the $p$-wave regime within a simple spinful model. In the former case we found that transconductance is quantized in units of $4e^{2}/h$, whenever the ABS spectrum is gapped, while in the latter it is in units of $2e^{2}/h$. This is consistent with a change of the pumped charge over a Josephson period from $2e$ to $e$, a distinctive signature of the fractional Josephson effect. Notably, in our odd-terminal setup this occurs under thermal occupation conditions—this is not the case in $4$-terminal JJ, see [Houzet and Meyer, 2018].

As mentioned above, the Green’s function formalism used to calculate the winding number turns out to be essential to properly account for the contribution of both the Andreev subgap levels and the continuum states, which are, in general, unavoidably mixed. In addition, this method possesses several important advantages. On the one hand, it can in principle incorporate the effect of interactions, provided a good approximation for the single-particle Green’s function is obtained. It can also properly handle the presence of a uncoupled Majorana state always pinned at zero energy. On the the other hand, the method could easily incorporate microscopic details of the junction through more realistic calculations of the junction’s Green’s functions that might be relevant for a proper comparison with experimental data.

Acknowledgements.

We acknowledge financial support from ANPCyT (grants PICTs 2013-1045 and 2016-0791), from CONICET (grant PIP 11220150100506) and from SeCyT-UNCuyo (grant 06/C526).

Appendix A Adiabatic expansion of the Josephson’s current.

We present here, for the sake of completeness, the expression derived in Ref. [Peralta Gavensky et al., 2018] to calculate the topological transconductance of a multi-terminal Josephson junction. In general, the adiabatic expansion of the Josephson current can be done directly with the Green’s function of the junction’s Hamiltonian. We begin with Dyson’s equation of motion of the Green’s function of the entire (non-interacting) Hamiltonian $H(t)$ describing the junction and the leads written in the Bogoliubov-de Gennes (BdG) representation of particle-hole and spin spaces used in the main text

[TABLE]

where the time dependence appears only on the tunneling Hamiltonian throughout the superconducting phases. The approach is based on the separation of slow and fast timescales in Dyson’s equation by taking advantage of Wigner’s representation. Under this formalism, two time variables, the relative time $t_{\text{rel}}=t-t^{\prime}$ and the average one $t_{\text{av}}=\frac{t+t^{\prime}}{2}$, are handled. Time derivatives of the former are used as a small parameter, making feasible a perturbation scheme. Wigner’s representation is given by

[TABLE]

Taking into account that $\partial_{t}=\frac{1}{2}\partial_{t_{\text{av}}}+\partial_{t_{\text{rel}}}$, Eq. (8) can be written to all orders as [Kershaw and Kosov, 2017]

[TABLE]

where we applied the following manipulation

[TABLE]

We are interested in a first order approximation, where the Green’s function is perturbed up to first derivatives of the average time. In that case, Eq. (10) is given by

[TABLE]

Introducing a series expansion for the Green’s function $\widetilde{\mathcal{G}}(\omega,t_{\text{av}})=\sum_{n}\widetilde{\mathcal{G}}^{(n)}$, we can obtain the next system of equations

[TABLE]

where the arguments $(\omega,t_{\text{av}})$ are omitted for the sake of brevity. Solving Eq. (13) and replacing into Eq. (14) we get

[TABLE]

where \Big{[}\cdot,\cdot\Big{]} is the commutator. We’ve also used the properties of the derivatives of inverse functions $\partial_{\lambda}A^{-1}=-A^{-1}\partial_{\lambda}AA^{-1}$, meaning in this case that

[TABLE]

The mean value of a single particle operator written in the BdG basis $\hat{\mathcal{O}}(t)=\frac{1}{2}\sum_{\alpha\beta}\hat{\Psi}^{\dagger}_{\beta}(t)\mathcal{O}_{\beta\alpha}\hat{\Psi}_{\alpha}(t)$ is given by

[TABLE]

By writing Eq. (19) in Wigner’s represention and expanding the Green’s function up to first order we obtain

[TABLE]

By specifying $\mathcal{O}$ as the current operator $\mathcal{J}_{\nu}=2e\partial_{\varphi_{\nu}}H$ and using that $\lim_{t^{\prime}-t\to\epsilon^{+}}\partial_{t_{\text{av}}}H=\sum_{\nu^{\prime}}\partial_{\varphi_{\nu^{\prime}}}H\dot{\varphi}_{\nu^{\prime}}(t)$, we get the first order correction to the current expectation value given by

[TABLE]

The averaged transconductance then simply results in

[TABLE]

which leads to Eq. (4) used to obtain our main results. We then recover an expression of this adiabatic transport coefficient in terms of the Green’s function winding number. By noticing that $\widetilde{\mathcal{G}}^{(0)}(\omega,t_{\mathrm{av}})$ depends on the average time $t_{\mathrm{av}}$ only throughout the superconducting phases we can define $\mathcal{G}(\omega,\{\varphi_{\nu}\})=\widetilde{\mathcal{G}}^{(0)}(\omega,\{\varphi_{\nu}(t_{\mathrm{av}})\})$. Furthermore, since the tunneling Hamiltonian acts locally on the end sites of each lead, this expression can be simplified so as to be written only in terms of the projected junction Green’s function $\mathcal{G}_{J}$. The averaged transconductance can then be written as

[TABLE]

where the phase derivatives have been explicitly written in terms of the tunneling Hamiltonian $\mathcal{T}$ as defined in Eq. (3), and the trace $\text{Tr}_{J}$ acts locally on the junction’s states. We emphasize that the result does not depend explicitly on the slow time-scale any more, since the instantaneous superconducting phases are integrated out.

Appendix B Boundary Green’s functions of the uncoupled leads by decimation.

The spinful boundary Green’s function of each one of the leads is obtained by a numerical decimation procedure applied to the Hamiltonian

[TABLE]

where the four component spinor in site $j$ is defined as ${\hat{\Psi}_{j}}{}=(c_{j\uparrow},c_{j\downarrow},c_{j\downarrow}^{\dagger},-c_{j\uparrow}^{\dagger})^{\mathrm{T}}$ in the Nambu representation and

[TABLE]

where $t$ is the hopping matrix element, $\alpha$ the spin orbit-coupling and $B_{x}$ the Zeeman field. The Pauli matrices $\tau_{a}$ ($\sigma_{a}$) and the identity $\tau_{0}$ ($\sigma_{0}$) act in particle-hole (spin) space.

The decimation procedure iteratively eliminates at each step of the Green’s functions equations of motion of even sites of the chain, eventually encapsulating numerous sites through a recursion relation using renormalized parameters. Assuming that the initial undressed Green’s functions of the left boundary $g^{L}_{0}(\omega)$, the bulk $g^{B}_{0}(\omega)$, and the right boundary $g^{R}_{0}(\omega)$ are all the same and equal to $(\omega-\mathcal{H})^{-1}$, at step $n\geq 1$ of the decimation they take the form

[TABLE]

where the renormalized hoppings are given by

[TABLE]

and $\Sigma_{n-1}(\omega)=T_{n-1}g_{n-1}^{B}(\omega)T_{n-1}^{\dagger}+T_{n-1}^{\dagger}g_{n-1}^{B}(\omega)T_{n-1}$. The initial tunneling is just the undressed one $T_{n=0}=T$ [see Eq. (25)]. For sufficiently large $n$, each of these Green’s functions accurately represent the left and right boundary of a semi-infinite chain (generally not equal in topological systems) and the bulk of an infinite one.

These boundary Green’s functions are used to construct the three terminal junction Green’s function (see Fig. 1) as $\mathcal{G}_{J}^{r/a}(\omega)=(g^{-1}(\omega\pm i\eta)-{\mathcal{T}})^{-1}$, where the supraindexes $r/a$ stand for retarded and advanced respectively, $g(\omega)$ is a $12\times 12$ block diagonal matrix of the form

[TABLE]

and $\mathcal{T}$ is the tunneling between the leads as described in Eq. (3) in the main text. The Feynman Green’s function used to calculate the winding number is then obtained as

[TABLE]

with $f(\omega)$ the Fermi-Dirac distribution function.

Appendix C Topological Hamiltonian for $B<B_{c}$: contributions from the ABS and the continuum states.

It was recently stated that, under quite general assumptions, all the topological information of the total $\mathcal{G}(\omega)$ of a system is encoded in the the so called “topological Hamiltonian” defined as $H_{\mathrm{top}}\equiv-\mathcal{G}^{-1}(\omega=0)$ (Wang and Zhang, 2012; Wang and Yan, 2013). Notice this is a well defined hermitian operator, since the self energy $\Sigma(\omega=0)$ has no imaginary part. The key idea is that the entire $\mathcal{G}(\omega)$ can be smoothly deformed to $\mathcal{G}_{\text{eff}}(\omega)=[\omega+\mathcal{G}^{-1}(\omega=0)]^{-1}$ preserving its topology. In this case, it is easy to prove that the winding number of the Green’s function in a given parameter space $(\varphi_{1},\varphi_{2})$ can be calculated with the Chern numbers associated to each of the eigenstates $|\psi_{n}\rangle$ of $H_{\text{top}}$ as

[TABLE]

where the sum runs over all occupied states and

[TABLE]

This theoretical tool provides a much simpler way of obtaining the topological invariants without the need of an integral over frequency space (see Eq. (5) in the main text).

In the case of Josephson junctions, the eigenstates resulting from diagonalizing $H_{\mathrm{top}}^{J}\equiv-\mathcal{G}_{J}^{-1}(\omega=0)$ should then reflect the topological nature of both the Andreev subgap structure and the continuum states above the superconducting gap, as we shall show below. There is only one caveat: the method is only applicable in our system for in plane magnetic fields below the critical value $B_{x}<B_{c}$, since for larger fields –after the topological transition– an unpaired Majorana emerges at zero energy (the “flat band”). In this case, the Green’s function is singular at $\omega=0$ and hence, not invertible.

In the upper panels of Fig. 5 we present the spectral density of the trijunction $\mathcal{A}(\omega)=-1/\pi\,\mathrm{Im}[\mathrm{Tr}(\mathcal{G}_{J}^{r}(\omega))]$ along the path $\varphi_{2}=-\varphi_{1}$ for different values of the fluxes threading the central scattering region between the superconducting leads. The in plane magnetic field was chosen to be zero ($B_{x}=0$) so that the eigenenergies are spin degenerate. The lower panels show the outcome of diagonalizing the topological Hamiltonian for the corresponding fluxes. Six particle-hole symmetric and spin-degenerate bands are present. The eigenvalues in solid lines, which are closer to the Fermi level, should possess the topological information of the Andreev subgap states, while the dashed eigenvalues should represent the continuum states. Note that even though the energy scale is very different from the original one, the touching of the lowest ABS with the continuum at $\varphi_{1}=\varphi_{2}=0$ and zero external flux $\Phi=0$ [see Figs. (a) and (d)] is accurately described, as well as the closing of the gap between the ABS at the first critical flux $\Phi=\Phi_{c}^{(1)}=0.53$ [see Figs. (c) and (f)].

In Fig. 6 we show the Chern numbers of these bands as a function of the external flux $\Phi$. The Chern associated to the ABS (upper panel in Fig. 6) is obtained from the band closer to the Fermi level, while the one related to the continuum (middle panel in Fig. 6) is just the sum of the Cherns of the two lower bands in Fig. 5. All of them are spin degenerate and, as a result, they lead to a simple factor 2 in their respective topological invariants. The sum of this contributions (lower panel in Fig. 6) is precisely the winding number of the full frequency-dependent Green’s function obtained in the main text, revealing the validity of the method. Notice also that the contribution from the ABS agree well with the calculation of Ref. [Meyer and Houzet, 2017].

This same procedure can be done for every $0<B_{x}<B_{c}$, where the spin degeneracy is lifted. The topological Hamiltonian $H_{\text{top}}^{J}$ possesses in this case 12 different particle-hole symmetric eigenvalues. In this regime, the contribution from the continuum states and the Andreev levels cannot be separated, since the gap closing makes one of the ABSs merge with the continuum. Nonetheless, the sum of all contributions up to the Fermi level correctly recovers the $s$-wave part of the phase diagram presented in Fig. 2(a) of the main text.