# Singularities of Andreev spectrum in multi-terminal Josephson junction

###### Abstract

The energies of Andreev bound states (ABS) forming in a -terminal junction are affected by independent macroscopic phase differences between superconducting leads and can be regarded as energy bands in periodic solid owing to the periodicity in all phases. We investigate the singularities and peculiarities of the resulting ABS spectrum combining phenomenological and analytical methods and illustrating with the numerical results. We pay special attention on spin-orbit (SO) effects. We consider Weyl singularities with a conical spectrum that are situated at zero energy in the absence of SO interaction. We show that the SO interaction splits the spectrum in spin like a Zeeman field would do. The singularity is preserved while departed from zero energy. With SO interaction, points of zero-energy form an dimensional manifold in dimensional space of phases, while this dimension is in the absence of SO interaction. The singularities of other type are situated near the superconducting gap edge. In the absence (presence) of SO interaction, the ABS spectrum at the gap edge is mathematically analogues to that at zero energy in the presence (absence) of SO interaction. We demonstrate that the gap edge touching (GET) points of the spectrum in principle form () dimensional manifold when the SO interaction is absent (present). Certain symmetry lines in the Brillouin zone of the phases are exceptional from this rule, and GET there should be considered separately. We derive and study the effective Hamiltonians for all the singularities under consideration.

###### pacs:

74.45.+c, 85.25.Cp, 03.65.Vf, 71.70.Ej## I Introduction

Superconducting junctions give rise to many interesting and unique physical phenomena, this being a base of the numerous applications in the field of quantum devices. A conventional Josephson junction with two superconducting leads hosts the Andreev bound states (ABS), that carry the supercurrent determined by the difference of macroscopic phases of the leads NazarovBlanter (). The properties of ABS may be altered by connecting the superconductors with special materials. For example, the exchange field in a ferromagnetic junction splits the ABS energies in spin. This may result in the -state Oboznov (). Recent studies address topologically protected bound states with zero-energy, called Majorana bound states, that occur in the 1D semiconductor nanowire with spin-orbit (SO) interaction, Zeeman splitting, and proximity-induced superconducting gap Mourik (); Rokhinson (); Das (); Deng (). The presence of Majorana bound states in the junction may double the period of current-phase relation FuKane (); Lutchyn (); Beenakker1 (). The coexistence of SO interaction and Zeeman effect breaks the spin-rotation and time-reversal symmetries. With symmetry broken, the Josephson current is not an odd function of the phase difference. This is called the anomalous Josephson effect Buzdin (); Reynoso1 (); Reynoso2 (); YEN1 (); YEN2 ().

The Josephson junctions involving various materials have been mostly investigated in two-terminal setups. There is a recent interest in multi-terminal Josephson junctions Heck (); Riwar (); Padurariu (). Such junctions have been realized, for instance, with crossed InSb/As nanowires Plissard (), where SO interaction is strong. Multi-terminal Josephson junction with superconducting leads is affected by independent phase differences. The energies of ABS are periodic in all phase differences. The system of energy levels of ABS can be regarded as a band structure in a dimensions. The phase differences play the role of quasimomenta Riwar (). The multi-terminal junctions may exhibit topological properties even if the superconducting leads and the connecting region are not made from topological or other exotic materials. In the case of two-terminal junction, the Andreev levels touch zero energy only when the transmission coefficient of the normal region is unity and the phase difference is . For the multi-terminal junctions, the Andreev levels can reach zero energy at some isolated points in dimensional space of phase differences Heck (); Riwar (); Padurariu (). Such points are topologically protected being the Weyl singularities studied theoretically in 3D solids Wan ().

The energy gap closes at a Weyl point and satisfies a conical dispersion in the vicinity of the point. The Berry curvature field is divergent at the Weyl points. They can be regarded as Dirac monopoles of the Berry curvature field bearing the topological charge, . A band structure with Weyl points can be continuously transformed to that without the points if two Weyl points with opposite topological charge meet each other to annihilate. In a 3D solid, the SO interaction and the inversion-symmetry breaking are essential for occurrence stable Weyl point Wan (); Burkov (); Murakami (). Very recently, the angle-resolved photoemission spectroscopy experiments confirmed the existence of the Weyl points in the 3D solids, such as TaAs SuYXu1 (); Lv (); Yang (), TaP NXu (), and NbAs SuYXu2 ().

Riwar et al. Riwar () demonstrate the presence of Weyl points in the multi-terminal short Josephson junctions. In contrast to solids, the Weyl point does not require SO interaction. The Andreev levels are Kramers-degenerate for the whole ()-dimensional -space. They discuss the transconductance to detect the Chern number due to the Weyl points by using one phase as a control parameter to switch the topological state. To reveal experimental signatures of the topology associated with Weyl points, the authors of Ref. Riwar, propose the following scheme. They consider a 2D band structure that depends on the two phases, and . The property of this band structure can be tuned by the remaining phase, . The 2D band structure is characterized by a Chern number that is proportional to the flux of Berry curvature field through a () 2D plane. The Chern number is changed by one anytime the plane crosses the position of Weyl singularity. The Chern number is observed as a quantized transconductance between the leads one and two, in similarity with the quantum Hall effect.

Hech et al. Heck () and Padurariu et al. Padurariu () also study the Andreev levels at energies close to zero in three-terminal Josephson junction. The authors claim that the zero energy states in such junctions may open opportunity for a single fermion manipulation. In the three-terminal junctions, the Weyl singularity is generally absent although the energy of the ABS can pass zero. The authors obtain a condition for zero energy ABS and study the density of states in detail. When SO interaction is present, the energy levels of the ABS are split in spin.

In this study, we investigate theoretically the singularities of the ABS spectrum in four-terminal Josephson junction taking SO interaction into account. We thus attempt to formulate the full picture of such singularities combining phenomenological and analytical methods and illustrating it with numerical results. The ABS energies are found from the Beenakker’s determinant equation Beenakker2 () using the scattering matrix of the junction. Mostly we concentrate on the case of short junction where we can disregard the energy dependence of scattering matrix. Sometimes the absence of energy dependence leads to extra degeneracy in the spectrum. To lift those we take into account the energy dependence by perturbation theory.

Firstly, we concentrate on Weyl singularities that occur at in the absence of SO interaction. We consider the behavior of the singularities upon gradual increase of the strength of SO interaction. We have found that SO interaction splits the spectrum of ABS in spin. Very much like Zeeman effect would do. The conical points are departed from to mirror symmetric positive and negative energies. The Weyl singularities thus remain topologically protected. A small modification of scattering matrix by a parameter changes the position of the mirror symmetric conical points, not eliminating them. As we show in numerical illustration, an annihilation of Weyl points of opposite charge can take place upon the tuning of the scattering matrix. We derive the effective Hamiltonian describing the vicinity of the Weyl points. The cones in the vicinity intersect at a 2D surface in a 3D space of the phases. We prove that this is the general property of ABS spectrum in the presence of SO interaction.

Singularities of a different type arise when ABS energies approaches the gap edge, . Owing to mirror symmetry of Andreev spectrum, there is a level with at the same position in the space of the phases. We formulate a mathematical analogy that permits to map Weyl singularities at in the presence (absence) of the SO interaction to singularities at in the absence (presence) of the SO interaction. Since the ABS energies in the presence of SO interaction reach zero energy at 2D surface, we expect the gap touching point to form 2D surfaces in the absence of SO interaction. Indeed this can be seen in a concrete numerical calculation. Employing the same analogy from the fact that without SO interaction the ABS energies reach zero at isolated point only, we derive that the gap edge touching (GET) in the presence of SO interaction generally occurs only at the isolated point. This implies that even a weak SO interaction removes the GET. We construct the effective Hamiltonian to describe this situation. In this case, a weak energy dependence of scattering matrix can also become important. An important peculiarity of the ABS spectrum concerns the vicinity of symmetry lines in 3D elementary cell of the space of the phases. Three of the four superconducting phases are the same at a symmetry line. Therefore the four-terminal junction at a symmetry line can be regarded as a two-terminal junction with unequal number of conduction channels in the two leads. As mentioned in Ref. Beri, , for two-terminal short junction, the SO interaction is irrelevant to causing spin splitting. Thus the vicinity of symmetry lines requires a separate consideration. We derive an effective Hamiltonian to incorporate the details of the GET in the vicinity of the symmetry lines.

These two types of singular point would reveal the topological nature of multi-terminal Josephson junctions. It brings a goal to propose nanostructures as artificial exotic materials.

The structure of this article is as follows. In Sec. II, we explain the model for the multi-terminal Josephson junction and the equation to determine the ABS energies. In Sec. III, we describe the spin splitting of Weyl singularities. Section IV is devoted to the GET point in general. Here we formulate and employ the mapping between and . We also concentrate separately at the vicinity of the symmetry lines in this Section. We conclude in the last Section.

## Ii Model and Formulation

In this Section, we explain the model in use.

We consider a junction connected to superconducting leads. An example of a physical system of this sort is given in a Fig. 1(a). All microscopic detail of the junction can be incorporated into the scattering matrix of the electrons and holes (Fig. 1(b)). We assume certain numbers of spin-degenerate transport channels in each lead, so are matrices, . The symmetry of Bogoliubov de-Gennes equation implies the relation with being a matrix realizing a time-inversion in the spin space. In addition to this, the requirement of time reversibility implies . The superconducting leads do not provide extra potential scattering. They are described by the Andreev reflection amplitudes for converting electron to hole and hole to electron that do not change transport channel index and are presented by diagonal matrix . being a superconducting phase of a lead corresponding to the channel . is an energy-dependent phase of Andreev reflection. We assume the same material for all of the superconducting leads, so that their order parameters are the same . In this case, the energy-dependent phase is the same for the all leads and is given by . The eigenvectors of the ABS satisfy . Therefore the energies of ABS are determined from Beenakker2 ()

(1) |

with

(2) | |||||

(3) |

The hat symbol on and is omitted for simplification. Note that is a unitary matrix.

For numerical calculations, the scattering matrix is taken as a random matrix. The SO interaction is taken into account in . In the presence (absence) of SO interaction, random matrix is a member of symplectic (orthogonal) ensemble. We introduce a parameter tuning the strength of SO interaction YEN1 (). The parameter varies from to , providing a continuous transition between the orthogonal () and the symplectic () ensembles.

Let us discuss the ABS energies in the vicinities of zero energy, , and superconducting gap edge, , where the first term in the determinant in Eq. (1) becomes and , respectively. In the absence of SO interaction, the scattering matrix commutes with . Thus,

(4) |

When the SO interaction is present, we can introduce a unitary matrix . Since is a real matrix,

(5) |

Let us disregard energy dependence of scattering matrices. Suppose that for a sufficiently general unitary matrix , there is eigenvector . Such eigenvector would correspond to an ABS at in the presence of SO interaction and to an ABS at in the absence of SO interaction. Correspondingly, an eigenvector would be an ABS at in the absence of SO interaction and to an ABS at in the presence of SO interaction. We thus establish a mapping between the situation at () with SO interaction and at () without SO interaction.

## Iii Weyl Singularities at

In this Section, we concentrate on the Weyl singularities. In our numerical illustrations, we concentrate on four-terminal short junctions with a single channel in each lead, and . In this case, the ABS energies are periodic function of three independent phases so we can restrict ourselves to the Brillouin zone for . The time-reversibility manifests itself as the inversion symmetry in this Brillouin zone, . To start with, we demonstrate that in the absence of SO interaction the energy levels of ABS exhibit the band gap closing points at zero energy Riwar (). Next, we continuously change the scattering matrix increasing the parameter , that is strength of SO interaction. We demonstrate that the SO interaction splits the conical spectrum in spin. The conical point departs from and the energy levels cross zero energy at 2D surface rather than isolated point. This proves topological protection of the Weyl singularity. To prove the generality of our conclusions, we derive an effective Hamiltonian from Eq. (1) that is valid in the vicinity of the singularity and at weak SO interaction.

### iii.1 Energies of the Andreev bound state

We obtain the ABS energies from Eq. (1). The scattering matrix is random. In the absence of SO interaction, we chose the random matrix from the circular orthogonal ensemble and disregard its energy-dependence.

We examine the spectrum for many random matrices. About 6% of them show a gap closing of the ABS energies indicating the Weyl singularities. The singularities always come in groups of four. Figure 2(a) gives the positions of singularities for a random matrix of choice. The time-reversal invariance guarantees that the singularity at is accompanied by the singularity of the same topological charge at . Owing to this, the positions of four singularities and the center of Brillouin zone are in the same 2D plane. Figure 2(b) shows the lowest positive energy of the ABS in this plane. The spectrum is symmetric with respect to phase inversion . In the origin, this energy reaches maximum . It drops down to two valleys close to the edges of the Brillouin zone. A zoom into a valley (Fig. 2(c)) shows that the energy actually reaches zero in two isolated points.

In Figs. 2(d) and (e), we show the ABS energies versus at fixed . In Fig. 2(d), the choice of is such that a singularity is reached at some . In Fig. 2(e), the line spanned by changing passes close to a singularity. We see the spectrum is conical near the singularity. For in all four leads, they are four positive and four negative ABS energies. Since the SO interaction is absent, the energies of ABS are doubly degenerate. In Figs. 2(d) and (e), the second ABS band is close to .

At the same choice of random scattering matrix, we increase the parameter thereby together continuously increasing the strength of SO interaction. As we see in Fig. 3(a), the SO interaction splits the ABS energies in spin. The absence of time-reversibility required for such splitting comes about non zero . As in Fig. 2, we plot in Figs. 3(b) and (c) the lowest positive energy of the Andreev states in the plane that passes through the Weyl singularities. We see in Fig. 3(b) that overall energy landscape have not changed significantly in comparison with Fig. 2(b). However, as seen in Fig. 3(c), the landscape has changed drastically in the vicinity of the singularities. The gap is closed at the closed contour encircling the singularities. The singularities are shifted to non zero energy while the spectrum remains conical in the vicinity of a singularity (Fig. 3(a)). In Fig. 3(d), we plot energy difference between the second and the lowest positive energy levels and observed that it goes to zero at the position of the singularities.

The SO interaction changes the position of the Weyl singularities while preserving their topological charge and the conical dispersion. When increases gradually, the four Weyl singularities move in the 3D space of the phases. Figure 4(a) shows trajectories of their positions for arranging from to . Solid and dashed curves give the position of the singularities with positive and negative charge. For this particular choice of the scattering matrix, the singularities of the opposite charge get close to each other upon increasing the SO strength and eventually annihilate at , so the junction is not topological any more (Fig. 4(a)).

We compute the energy of the singularity and the plots result in Fig. 4(b). We see that this energy is zero in the absence of SO interaction. The energies of the singularities of different topological charge increase and become different with increasing . At , the energies come close together merging at the annihilation point .

### iii.2 Effective Hamiltonian

To prove the generality of the above results, we derive an effective Hamiltonian for ABS that is valid in the vicinity of a singularity. When the SO interaction is absent, the singularity is at . At the position of singularity, , the scattering matrix should have an eigenvalue . However, is a rather special matrix: at it can be presented as (cf. Eq. (4)). Here we introduce . We assume nothing about regarding it as an arbitrary unitary matrix. This implies that the eigenvalues of come in complex-conjugated pairs. If is an eigenvector of with eigenvalue , is an eigenvector with eigenvalue . This implies that the two orthogonal eigenvectors corresponding to eigenvalue ,

(6) |

To obtain the effective Hamiltonian, we project the matrix in Eq. (1) on to the space spanned by these two eigenvectors. For small deviations of the phases from , . We expand the scattering matrix as , where is a Hermitian matrix proportional to . Up to the first order in and , we find

(7) | |||||

(8) |

With this, the determinant equation (1) can be presented as an eigenvalue equation for a effective Hamiltonian Riwar (),

(9) |

where are the Pauli matrices in the basis of and . Using the relation , we prove that and real parameters are given simply by and . The Hamiltonian in Eq. (9) is a case of Weyl Hamiltonian.

It has opposite eigenvalues, . Expanding in , , we obtain

(10) |

For 3D space of the phases, the matrix is positively defined. We reproduce a conical spectrum of ABS in the vicinity of singularity. For dimensional space, the matrix has zero eigenvalues. We stay at zero energy if we depart from in this directions. So that, the singularities are concentrated at dimensional manifold.

Let us take into account weak SO interaction, expanding

(11) | |||||

(12) | |||||

with being Pauli matrices in spin space, being associated Hermitian matrices in channel space. We project on four dimensional space of spins and vectors . We observe that the structure of matrices is quit difference from . Using the relations between and , we prove and . With this, the effective Hamiltonian becomes:

(13) |

The BdG symmetry thus guarantees a special structure of this Hamiltonian where spin and orbital degree of freedom are totally separated. plays a role of an effective Zeeman field that splits the original conical spectrum.

(14) |

with being the spin projection on the axis of the effective Zeeman field. The SO interaction does not remove the conical point but rather shift it in energy by . In a 3D space, the zero energy is achieved at 2D surface of ellipsoid defined by . The singularity is enclosed by the surface. This consideration shows generality of our numerical results. The size of the ellipsoid enclosing a singularity increases with increasing SO interaction. In Fig. 3(c), we see that the ellipsoids enclosing each singularity have already merged together at .

## Iv Gap Edge Touching

In this Section, we consider the ABS in the vicinity of the superconducting gap edge . We show that the ABS energies reach the GET in the absence of SO interaction at a 2D surface in the 3D space of the phases. The SO interaction lifts the GET almost everywhere except particular manifolds: symmetry lines and isolated points. We investigate these cases separately and establish effective Hamiltonians.

### iv.1 GET at symmetry lines

To understand the peculiarities of the GET for multi-terminal junction, let us first consider two-terminal one with unequal number channels in the left and right lead, . The estimation of a number of localized Andreev states is somehow ambiguous. From one hand, one may argue that there are only such state because only that many states are sensitive to the superconducting phase difference between two leads. From the other hand, the full number of Andreev states is given by provided the total number of channels is even. The remaining states are pinned to the gap edge with no regard for superconducting phase difference. In a two-terminal junction, these states are indistinguishable from the states of the continuous spectrum. This is not a case of multi-terminal junction. We note that a multi-terminal junction in fact becomes two-terminal one if the leads are separated in two groups with superconducting phases are the same within each group. For instance, in our four-terminal set up, one can choose and . Such setting defines a symmetry line in multi-dimensional space of the phases. For our example with one channel in each lead, we find extra state pinned at the gap edge along the symmetry line. In distinction for a two-terminal case, this extra state can not be attributed to the continuous spectrum since it departs from the gap edge if we go off the symmetry line.

### iv.2 ABS energies near the gap edge: general

Let us consider GET at general position in 3D space of phases. Let us note mathematical analogy between the spectrum at and . The spectrum is determined by properties of the scattering matrix in Eq. (1). Zero energies () correspond to the eigenvalue of of the matrix while GET correspond to the eigenvalue of . From the other hand, as we have seen in Section II, in the absence and presence of SO interaction, we can represent as and , respectively, being a general unitary matrix. This establishes a rather unexpected mapping of the spectrum at in the presence (absence) of the SO interaction to the spectrum at in the absence (presence) of the SO interaction. We have derived in Section III.1 that the ABS reach at 0D or 2D manifolds in the 3D spaces in the absence or presence of SO interaction, respectively. This implies that the GET occurs at 2D and 0D manifolds in the absence and presence of SO interaction, respectively.

Let us see this in numerical results. In the absence of SO interaction, scattering matrix is chosen randomly from the circular orthogonal ensemble. Its energy-dependence is disregarded. We examine the spectrum for many random scattering matrices. The GET is observed for all the matrices. For a 3D space of the phases, the ABS energy reaches at a 2D surface. This is illustrated in Fig. 5. We fix and find the GET points forming a 1D curve. Interestingly, the curve passes symmetry lines where we indeed expect GET. For instance, in Fig. 5(a), the curve passes the symmetry line at and passes the symmetry line at . As is tuned from to (Fig. 5(b)), the positions of the first passing is gradually sifted to . At , the curve approaches two symmetry lines in this plane, and . For positive , the spectrum is obtained from the inversion symmetry. We see that the GET occurs at a single connected 2D surface that includes all four symmetry lines. In principle, there is nothing to for bit more sophisticated topology of the surface. However, in a dozen of samples we have explored, we have found no complex topology.

Figures 6(a) and (a’) show the ABS energies with the GET in the absence of SO interaction. The energies are plotted versus while and are fixed. The GET is found at . In the vicinity of the touching point, the spectrum of ABS energy is parabolic rather than conical. A mirror symmetry of the ABS energies guarantees the GET at at the same position in the 3D space of the phases.The SO interaction (Figs. (b) 6and (b’)) lifts a degeneracy of the energies in spin. In Fig. 6(b’), one spin-resolved ABS energy comes very close to the gap edge. However, as seen in the inset, the touching does not take place in the presence of SO interaction conform to our expectation.

In Fig. 7, we illustrate the situation in the vicinity of symmetry lines. We plot the ABS energies along the line . This line crosses the symmetry lines. The crossing points are indicated by arrows in the Figure. We see that at these points, , the GET survives in the presence of SO interaction. Also, the spin splitting of the lowest ABS vanishes at this point. Then, this proves that the SO interaction is irrelevant at the symmetry lines.

### iv.3 effective Hamiltonian in general case

To prove the generality our numerical results concerning the GET, we derive an effective Hamiltonian for ABS near taking into account weak SO interaction and energy-dependence of the scattering matrix. We assume no vicinity of symmetry lines. This specific situation will be considered separately in the next Section.

Let us first neglect SO interaction and energy-dependence of the scattering matrix. For a GET point at , the matrix has to have an eigenvalue . We put . This eigenvalue is double degenerate; if is an eigenvector of with the eigenvalue , is also an eigenvector of ,

(15) |

We project the matrix in Eq. (1) on to together the space spanned by these two eigenvectors. Next we expand full scattering matrix in Eq. (1) with respect to small phase deviation from , weak SO interaction, and weak energy-dependence. The first order expansion in (deviation) and (SO interaction) is the same as given by in Eqs. (7) and (12), respectively. The energy-dependence is expanded with a Hermitian matrix ,

(16) | |||||

(17) | |||||

The causality of scattering matrix implies that all eigenvalue of are positive.

The first order correction terms are similar to those for Weyl singularity at . However, since the sign of the eigenvalue is opposite, the selection rules for matrix elements representing the deviation and SO interaction are interchanged. For small phase deviations, ; . The expansion of the weak SO interaction gives three independent vectors in spin space, and . For weak energy-dependence, . No restriction applies to the off-diagonal matrix elements. So we define . The positivity of implies .

To deviate an effective Hamiltonian, we introduce the energy deviation from the edge . Up to the first order in , we obtain . Summarizing all these, we obtain that the spectrum is defined by the following eigenvalue equation

(18) |

Here . The full analysis of this Hamiltonian is involved. Here, we discuss several simple cases. Let us first neglect both SO interaction and energy-dependence. The eigenvalues of the Hamiltonian are in this case simply . The negative eigenvalue does not lead to any localized state. The positive eigenvalue provides since is linear in phase deviation this defines a parabolic spectrum touching the gap edge at . We see that the GET requires one condition to be fulfilled in distinction from three condition required for Weyl singularity. This is why the manifold of GET point in -dimensional space of phases generally has dimension in the absence of SO interaction and energy-dependence of scattering matrix.

If we take energy-dependence into account, the eigenvalues are given by . The plot of eigenvalues and spectrum is given in Fig. 8(a). We see that the energy-dependence modifies the GET in a rather complex way. The spectrum in the vicinity of the GET line forms two bands. While one of the band never touches the edge, the other band exist only in the vicinity of the line and merges with the continue upon increasing .

Let us neglect now the energy-dependence and take into account SO interaction. When terms are taken into account, the four eigenvalues of are

(19) |

with . Two positive eigenvalues define two spin-split bands that never touch the edge (Fig. 8(b)) unless a special condition discussed below is fulfilled. The minimal energy distance to the gap edge is not achieved at but rather is shifted depending on the parameters of the SO interaction.

If the energy-dependent terms and the spin-splitting are the same order of magnitude, the picture of the GET can become complex. For instance for an example given in Fig. 8(c), there are only three bands. One existing only in the vicinity of the gap edge line.

Let us now discuss a special condition of the GET in the presence of sufficiently strong SO interaction. To analyze this, let us derive an effective Hamiltonian with strong SO interaction assuming a GET point to be present at . Strong SO interaction guarantees that only two spin-dependent eigenvectors are important instead of the four as in the previous consideration. The two eigenvectors and satisfy

(20) |

with . The problem is thus mathematically equivalent to our consideration of Weyl singularity in Sec. III.2. The first order expansion in (deviation) and (energy-dependence) results in correction terms and , respectively. We prove the selection rules for elements of matrices and ; . The effective Hamiltonian thus reads

(21) |

with , , and . We keep here energy-dependent term that is absent in the Hamiltonian (9). We have just found a Weyl singularity with a conical point near the gap edge. The conical point requires three conditions to be fulfilled . Therefore we expect the points to form dimensional manifold in the dimensional spaces. The energy-dependent term shifts the energy of the conical point from the gap edge similar to the effect of the SO interaction in Eq. (13) (Fig. 8(d)). We have tried to find this singularities in our numerical simulation. So far, we have found none. The reason of this is not completely clear for us. We hypothesize that probability to find a random scattering matrix with such singularity is low because the GET point in the presence of SO interaction tend to stick to the symmetry line.

### iv.4 The vicinity of a symmetry line

Let us consider an effective Hamiltonian in the vicinity of symmetry lines. The SO interaction is assumed to be strong. We concentrate on a four-terminal junction. At the symmetry lines, three of the four superconducting phases are equal. The electron scattering matrix can be presented in a block structure,

(22) |

Assuming there are channels in three leads having the same phase and channels in the other lead. is a reflection matrix for three leads having the same phase. is that of the other lead, and are the transmission matrices between the three leads and the other lead, their dimensions are , , correspondingly. At the symmetry line, we have independent vectors satisfying . These states are disconnected from the other lead. Their energies are precisely at the gap edge, therefore

(23) |

is satisfied. Similar to previous considerations, these eigenvectors can be arranged into conjugated pairs: if is an eigenvector, is also an orthogonal eigenvector. Note that is a unitary matrix.

To derive the effective Hamiltonian, we project Eq. (1) onto the subspace of these eigenvectors. The consideration can be done for arbitrary dimension, but for the sake of comprehensibility, we concentrate on this situation with a single channel in each lead. In this case, and and we project on four orthogonal states . For a small phase deviation from the symmetry line, . The matrix element can be represented as

(24) |

Here, is the diagonal matrix of the phase deviations. Note that we use to arrive at this. A matrix representation of this correction term for the four bases gives a Hamiltonian in the vicinity of the symmetry line

(25) |

Here, we use . The parameters in block diagonal components are defined as , , and the same way for . Those in off-diagonal components are , and . The eigenvalues of are

(26) |

with and . This is similar to Eq. (19), however, in the case under consideration, all elements of the effective Hamiltonian are proportional to the phase deviations. This results in a linear splitting of four-hold degenerate eigenvalue and two bands touching the gap edge at the symmetry line (Fig. 9(a)).

The energy-dependence is taken into account the same way as Sec. IV.3. It gives a correction term in the Hamiltonian that importantly does not vanish at ,

(27) |

with and the same for , , and . Here, due to the causality. This guarantees the positive eigenvalues at the symmetry lines, , given by . The eigenvalues are doubly degenerate. The deviation from the symmetry line gives rise to a linear splitting of the eigenvalues at further increase . The eigenvalues approach a linear asymptotics given by Eq. (26). We illustrate the spectrum in Fig. 9(b).

## V Conclusions and Discussion

We have studied the singularities and peculiarities in the ABS spectrum of a Josephson junction connected to superconducting leads. The ABS energies in such junctions depend on independent superconducting phase differences, being a periodic function of all phases. Therefore, they can be regarded as energy bands in the dimensional periodic solid, if one associates with quasimomenta. We have concentrated on the singularities related to topological properties and use numerical illustrations of the spectrum and as well as derive effective Hamiltonians to describe the vicinity of the singularities. The illustrations are made for a four-terminal short junction. In this case, the energies of ABS correspond to the bands in a 3D solid. The ABS energies are calculated from Beenakker’s determinant equation using scattering matrix.

We reveal the singularities in the vicinity of zero energy and near the gap edge. We establish a mathematical analogy between the spectrum at and .

First, we have considered Weyl singularities near zero energy. When the SO interaction is absent, the singularities are found at accompanying conical spectrum. The Weyl singularities occur at isolated 0D points in the 3D space of the phases. The SO interaction splits the singular points to mirror symmetric positive and negative energies in spin. A small modification of scattering matrix only shifts the position of the Weyl singular points but does not eliminate those since they are topologically protected. In the presence of SO interaction, zero energy points in the vicinity of Weyl singularities form a 2D manifold in the 3D space of the phases. This 2D manifold encircles the singular point. To prove our numerical results, we have derived an effective Hamiltonian that is valid in the vicinity of the singularity and at weak SO interaction. Eigenvalues of the Hamiltonian reproduce a conical singularity in the spectrum. They reach zero at an ellipsoid enclosing the singular point.

Exploiting the mentioned analogy between the ABS at and , we have investigated the spectrum in the vicinity of the gap edge. The analogy implies the GET occurs at a 2D surface in the 3D space in the absence of SO interaction and at isolated points in the presence of SO interaction. Thus, the SO interaction generally lifts the GET except specific situations. We have established effective Hamiltonians for two specific cases: symmetry lines and isolated points. We have also taken into consideration a weak energy dependence of the scattering matrix relevant for fine structure of GET. The effective Hamiltonians derived prove the generality of our numerical results. For the GET in a strong SO interaction, the Hamiltonian indicates a Weyl singularity with a conical spectrum. The singularity is shifted from the edge by the energy-dependent term similar to the effect of SO interaction at . However, we do not find the isolated GET point in our numerical simulations.

At a symmetry line, since three of the four phases are equal, the four-terminal junction can be regarded as a two-terminal junction with unequal channel numbers in the leads. Then, only a single ABS is sensitive to the phase differences while the other stick to the gap edge. We have also derived an effective Hamiltonian in the vicinity of the symmetry lines. The energy-dependence of the scattering matrix lifts the GET and adds another Andreev state, which is localized in the vicinity of the symmetry line.

Our study provides a numerical and analytical evidence for the Weyl points in multi-terminal Josephson junctions. These points are different from Weyl singular points from that predicted and recently found in 3D solids, such as TaAs SuYXu1 (); Lv (); Yang (), TaP NXu (), and NbAs SuYXu2 (). In the 3D solids, the absence of inversion symmetry in the material and removing the spin degeneracy are essential for a stable Weyl point. However, multi-terminal Josephson junction shows stable Weyl points even though the ABS energies are not spin-split. Our study thus facilitates an alternative method of realization of Weyl fermions in condensed matter.

The Weyl points in multi-terminal Josephson junction should provide positive and negative topological charges in the space of the phases. The charges are sources of the Berry curvature fields. Riwar et al. Riwar () have proposed an experimental setup to detect the Chern number by measuring the quantized transconductance, in similarity with the quantum Hall effect. In their scheme, the Chern number is defined as the integral of the Berry curvature field on a 2D plane in the space. Physically, the integration is achieved by sweeping the phases by applying finite bias voltages to the superconducting leads. Riwar et al. considered a situation without SO interaction where the Weyl singular points appear always at E=0 and are doubly degenerate with respect to spin. Then the lowest positive (and the highest negative) levels of the ABS are relevant to the Chern number. Near the Weyl point, the quantization of transconductance may be sensitive to the temperature of the junction owing to the undesired thermal activation of a quasiparticle in the lowest subband. In this paper, we take into account the SO interaction and consider Weyl singularities of two types. Although the SO interaction shifts the conical point from E=0, it does not influence the Chern number arising from the singularity. Therefore, the topological charge may be still detected by measuring the quantized transconductance. The same applies to the singularity at the gap edge. The SO shift of the singular point from the Fermi energy implies that the probability of undesired thermal activation of a quasiparticle state near the point of the singularity is small at sufficiently low temperatures, this facilitates the observation of the topological effect.

## Acknowledgment

We appriciate the fruitful discussion with Roman-Pascal Riwar, Manuel Houzet, Julia S. Meyer in Université Grenoble Alpes. This work has been partially supported by JSPS Postdoctoral Fellowships for Research Abroad and the Nanosciences Foundation in Grenoble, in the frame of its Chair of Excellence program grand in Grenoble.

## References

- (1) Y. V. Nazarov and Y. M. Blanter, Quantum Transport: introduction to nanoscience, (Cambridge University Press, Cambridge, 2009).
- (2) V. A. Oboznov, V. V. Bol’ginov, A. K. Feofanov, V. V. Ryazanov, and A. I. Buzdin, Phys. Rev. Lett. 96, 197003 (2006).
- (3) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
- (4) L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nature Phys. 8, 795 (2012).
- (5) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nature Phys. 8, 887 (2012).
- (6) M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Nano Lett. 12, 6414 (2012).
- (7) L. Fu and C. L. Kane, Phys. Rev. B 79, 161408(R) (2009).
- (8) R. M. Lutchyn, J. D. Sau, and S. Das Sarma Phys. Rev. Lett. 105, 077001 (2010).
- (9) C. W. J. Beenakker, D. I. Pikulin, T. Hyart, H. Schomerus, and J. P. Dahlhaus Phys. Rev. Lett. 110, 017003 (2013).
- (10) A. Buzdin, Phys. Rev. Lett. 101, 107005 (2008).
- (11) A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Phys. Rev. Lett. 101, 107001 (2008).
- (12) A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Phys. Rev. B 86, 214519 (2012).
- (13) T. Yokoyama, M. Eto, and Yu. V. Nazarov, J. Phys. Soc. Jpn. 82, 054703 (2013).
- (14) T. Yokoyama, M. Eto, and Yu. V. Nazarov, Phys. Rev. B 89, 195407 (2014).
- (15) B. van Heck, S. Mi, and A. R. Akhmerov, Phys. Rev. B 90, 155450 (2014).
- (16) R.-P. Riwar, M. Houzet, J. S. Meyer, and Yu. V. Nazarov, arXiv: 1503.06862v2.
- (17) C. Padurariu, T. Jonckheere, J. Rech, R. Mélin, D. Feinberg, T. Martin, and Yu. V. Nazarov, arXiv: 1508.03289v1.
- (18) S. R. Plissard, I. van Weperen, D. Car, M. A. Verheijen, G. W. G. Immink, J. Kammhuber, L. J. Cornelissen, D. B. Szombati, A. Geresdi, S. M. Frolov, L. P. Kouwenhoven, and E. P. A. M. Bakkers, Nature Nanotech. 8, 859 (2014).
- (19) X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
- (20) A. A. Burkov and L. Balents, Phys. Rev. Lett. 107, 127205 (2011).
- (21) S. Murakami, New J. Phys. 9, 356 (2007).
- (22) S.-Y. Xu, I. Belopolski. N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Science DOI: 10.1126/science.aaa9297 ¡http://www.sciencemag.org/content/early/2015/07/15/science.aaa9297.full.pdf¿ (2015).
- (23) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Phys. Rev. X 5, 031013 (2015).
- (24) L. Yang, Z. Liu, Y. Sun, H. Peng, H. Yang, T. Zhang, B. Zhou, Y. Zhang, Y. Guo, M. Rahn, D. Prabhakaran, Z. Hussain, S.-K. Mo, C. Felser, B. Yan, Y. Chen, Nature Phys. 11, 728 (2015).
- (25) N. Xu, H. M. Weng, B. Q. Lv, C. Matt, J. Park, F. Bisti, V. N. Strocov, D. Gawryluk, E. Pomjakushina, K. Conder, N. C. Plumb, M. Radovic, G. Autès, O. V. Yazyev, Z. Fang, X. Dai, G. Aeppli, T. Qian, J. Mesot, H. Ding, M. Shi, arXiv:1507.03983.
- (26) S.-Y. Xu, N. Alidoust, I. Belopolski, C. Zhang, G. Bian, T.-R. Chang, H. Zheng, V. Strokov, D. S. Sanchez, G. Chang, Z. Yuan, D. Mou, Y. Wu, L. Huang, C.-C. Lee, S.-M. Huang, B. Wang, A. Bansil, H.-T. Jeng, T. Neupert, A. Kaminski, H. Lin, S. Jia, and M. Z. Hasan, arXiv:1504.01350v2.
- (27) C. W. J. Beenakker, Phys. Rev. Lett. 67, 3836 (1991); Phys. Rev. Lett. 68, 1442(E) (1992).
- (28) B. Bèri, J. H. Bardarson, and C. W. J. Beenakker, Phys. Rev. B 77, 045311 (2008).