Andreev reflection in Euler materials

Many previous studies of Andreev reflection have demonstrated that unusual effects can occur in media which have a nontrivial bulk topology. Following this line of investigation, we study Andreev reflection by analysing a simple model of a bulk node with a generic winding number n > 0, where the even cases directly relate to topological Euler materials. We find that the magnitudes of the resultant reflection coefficients depend strongly on whether the winding is even or odd. Moreover this parity dependence is reflected in the differential conductance curves, which are highly suppressed for n even but not n odd. This gives a possible route through which the recently discovered Euler topology could be probed experimentally.


I. INTRODUCTION
The study of topological materials is currently a very active field of research.Investigations in this area range from broad theoretical analyses to direct experimental investigation of the properties of particular materials.Recent advancements in the characterization of topological band structures using symmetry eigenvalue methods have been significant [1][2][3][4][5][6][7][8][9][10].However, there is growing interest in multi-gap topological phases [11] since they generically cannot be explained within this paradigm.
A prominent example entails Euler topology [7,[12][13][14][15], which depends on the conditions between multiple gaps in the band structure of a C 2 T or PT symmetric material.Under these conditions, band nodes in two-dimensional materials carry non-Abelian 'frame charges' in momentum space [12,16,17], akin to π disclination defects in bi-axial nematics [6,18,19]; phases with a non-trivial Euler topology can be formed by braiding such degeneracies between successive bands [12,13,15,16].The ability of such nodes lying in a patch of the Brillouin zone, D ⊆ BZ, to annihilate is encoded in the Z-valued Euler class invariant χ, which is the real counterpart of the Chern number of a complex vector bundle [11,12].Since a region containing no nodes has χ = 0, a nonzero value of χ indicates a topological obstruction to the possibility of the nodes gapping out in this region.Around each node, a non-zero Euler class χ manifests itself as a winding w = 2χ in the two band subspace carrying the node.The braiding of nodes in reciprocal space therefore provides a means through which nodes carrying higher winding numbers could be realised in real materials.Importantly, such multi-gap phases are increasingly being related to novel physical effects.Examples include monopole-anti-monopole generation [14] that have been seen in trapped-ion experiments [20] or novel anomalous phases [21], and are increasingly gaining attention in different contexts that range from phononic systems and cold atom simulators to acoustic and photonic metamaterials [17,20,[22][23][24][25][26][27][28][29][30][31][32][33][34][35].
One way in which the unusual properties of topological materials can become manifest is in their Andreev reflection characteristics.Andreev reflection is the pro-cess through which an electronic excitation in a normal metal can scatter into a hole and/or Cooper pair when incident on the boundary between the normal metal and a superconductor.In the historically well-understood case of a quadratic band dispersion, an incident electron which scatters into a hole is retro-reflected, and travels away from the boundary along the same direction from which it arrived [36].In contrast, in graphene an electron incident on a normal-superconductor boundary can in addition undergo specular Andreev reflection.Such unusual scattering characteristics of the single-particle excitations in graphene can ultimately be attributed to the properties of the nodes (Dirac cones) within the band structure of monolayer graphene.The nearest-neighbour tight-binding model of this material famously exhibits degeneracies at two points K and K in the Brillouin zone; low-energy excitations about these points behave like relativistic particles with an approximately linear dispersion relation and, due to the non-zero winding numbers carried by the nodes, these excitations are chiral [37,38].These properties enable a condensed matter realisation of the Klein paradox, which in some cases can lead to perfect Andreev reflection [39,40].Studies have also shown that Rashba spin-orbit interactions can play an important role in Andreev reflection in graphene [41,42].
In some cases it has been found that the Andreev scattering matrix can be directly related to the topological invariants in the bulk of a material.For example, in topological superconductors with chiral symmetry, the quantum number Q ∈ Z of the BDI class is equal to the trace of the Andreev S-matrix r he [43][44][45].Moreover, in certain topological superconductors the interaction of Majorana bound states with electronic excitations can have a significant impact on the Andreev reflection characteristics [46][47][48][49][50].In such materials the differential conductance curves depend on the number of vortices present, and differ in the cases that this number is even and odd [51][52][53].The investigation of Andreev reflection in Weyl semimentals has also shown that the (pseudo-)spin structure of the Hamiltonian in the vicinity of a band node can strongly influence the directional dependence of Andreev reflection [54][55][56].
Motivated by the unusual Andreev reflection proper-arXiv:2302.07094v1[cond-mat.mes-hall]14 Feb 2023 ties of topological materials as described above, in this work we explore low-energy Andreev reflection in Euler materials.To do so, we make use of a simple model of a node which carries a generic integer winding number n ≥ 0 since, as previously mentioned, this is precisely the property exhibited by Euler materials around a band node.(Note that graphene has band nodes with winding numbers w = ±1; our results agree with previous findings in this limit.)We find that, when n is even, the ability of an electron to scatter into a hole is strongly suppressed at all angles, while for n odd the probability of Andreev reflection at normal incidence is always exactly unity.Moreover, for n even the degree of suppression increases with the strength of an externally applied gate potential U 0 (though this effect is less pronounced for large n).The differential conductance curves resulting from this behaviour are distinctive and could be measured experimentally and used as an indicator of the presence of Euler topology.

II. MODEL
We now describe and make use of a model to investigate the properties of Andreev reflection in topological Euler materials.Although the braiding of band nodes carrying non-Abelian frame charges requires the Bloch Hamiltonian to have three or more bands, we can obtain an effective description of the low energy physics around any given node by projecting the Hamiltonian to the twoband subspace containing the degeneracy.In this space, a non-zero Euler class is manifested as a winding of the Bloch Hamiltonian around the node.More precisely, the Euler class of a generic real two-band Bloch Hamiltonian in two dimensions where k = (k x , k y ) ∈ B.Z., is given by χ = −w/2, where the total winding number and the region D i ⊂ B.Z. contains the i th node of H.More generally, χ = D Eu − ∂D a, where the Euler form Eu and the connection a may be computed from the Berry-Wilczek-Zee connection 11,13,16,29].
With regard to the above Hamiltonian we note previous insighful reports that considered the above Hamiltonian in the context of stable winding and band degeneracies [57] as well as its use as the simplest ansatz supporting a finite Euler class [13].The point from a multi-gap perspective [11,12,14,29] is however that the nodes formed by a two-band subspace on an isolated patch of the Brillouin zone are stable as long as these bands remain disconnected from the other bands of the manyband context.Indeed, to induce a finite Euler class in a lattice model one necessarily needs a multi-gap structure, as can be directly analyzed [12] also in agreement with general homotopical arguments [11,15], and perform a braid between the nodes of adjacent energy gaps.We stress that this general multi-gap perspective sheds light on addressing physical features as recently discovered [14,21,24,29,35].
A minimal model of a single node with winding number w = n, which for simplicity we choose to be at the origin k = 0, is given by setting θ(k) = n arg(z) = arg(z n ) and r(k) = |z| n = |k| n , where z = k x + ik y .This makes the components of H homogeneous polynomials in k x , k y , since it then follows that The first few polynomials are x k y − k 3 y .In particular, when n = 1 the Hamiltonian is H 1 = k x σ x + k y σ z , which is related to the low-energy graphene Hamiltonians H ±K = k x σ x ±k y σ y via a unitary transformation.
Although the Euler class of the entire Brillouin zone must be quantised to an integer, it is possible to have isolated nodes with arbitrary winding numbers.For this reason, in the following we will allow n to be any integer rather than restricting it to be even.
We now suppose that we have a material which contains a node with a generic winding number n within its bulk band structure; from now on, we will refer to this as an Euler material.If, in addition, a position dependent superconducting pairing potential ∆(x) is induced in the material, for example via the proximity effect, then the quasiparticle excitations in the system can be described by a Bogoliubov-de Gennes (BdG) Hamiltonian of the form where H n (k) = P + n (k)σ x + P − n (k)σ z is the two-band Euler Hamiltonian with winding number n described above.In Eq. ( 4) we have also allowed for the possibility of an externally applied electrostatic potential U (x).When the pairing and electrostatic potentials vary as a function of position, the momentum k should be interpreted as a derivative in real space, and the excitations of the system may be determined by solving the PDE for positive eigenvalues ε ≥ 0, subject to appropriate boundary conditions.Suppose now that a normal Euler material fills the semi-infinite plane x < 0, while in the region x > 0 both the pairing and electrostatic potentials are non-zero and the material is superconducting.In particular, suppose that ∆ and U are uniform in each of these regions, Euler Superconductor Normal Euler Figure 3.In addition to the propagating wave solutions that describe electron-and hole-like quasiparticle excitations inside the bulk of the normal and superconducting Euler materials, when the winding number is n there exist 2(n − 1) modes localised at the boundary.
When an electron propagating in the normal region is incident on the boundary x = 0, it may scatter into an electron or a hole in the normal region, or a Cooper pair in the superconducting region, each with a certain probability.To determine the amplitudes for these various processes, it is necessary to solve the real-space BdG equation in the normal and superconducting regions.By finding the eigenstates of this equation in the normal and superconducting regions, and then matching these solutions at the boundary x = 0, we determine the reflection and transmission coefficients.The results of this calculation are shown in the Appendix.
In Fig. 2 the probability R = |r| 2 for an electron to reflect off the boundary is shown as a function of the incident angle α, the winding number n, and the strength U 0 of the applied electrostatic potential.The energy ε of the incoming electron is here less than the superconducting gap, ε < ∆ 0 , so the sum of the reflection and Andreev reflection coefficients R and R A is exactly equal to 1.This is because the excitation energy ε is less than the energy ∆ 0 required to create a Cooper pair out of the vacuum, so the electron can scatter only into states on the normal side.Above the critical angle there are no hole states available for the electron to scatter into, and the probability of Andreev drops to zero (equivalently, R = 1).We note that, for fixed ε, E F the angle α c increases monotonically with n, so that the range of angles over which Andreev reflection can take place is larger for greater values of n.
Interestingly, the form of the angular dependence of the reflection probability is qualitatively different when the winding number n is even or odd.For n odd, the probability for Andreev reflection is always exactly 1 at normal incidence (α = 0), independent of the strength of the potential U 0 .As U 0 is varied, the value of R increases uniformly across α.However, even for large changes in  U 0 , the change effected in R is small.On the other hand, for n even the value of R is not fixed to zero at α = 0, i.e. there is a non-zero probability for the electron to reflect as an electron even at normal incidence.Moreover, R shows a significant variation with the applied potential: as U 0 is increased, the value of R increases and tends towards the value R(α) = 1 as U 0 → ∞.The speed at which R increases depends on the magnitude of n: for larger n, the suppression of R A is smaller and it requires a stronger potential to send R → 1.These unusual features are qualitatively reproduced in the super-gap region ε > ∆ 0 , though in this latter case it no longer holds that Another interesting feature of the problem is that, in addition to the propagating wave states that, far into the bulk, represent the usual electron-and hole-excitations, for n > 1 there are 2(n − 1) evanescent solutions that are localised at the boundary x = 0. Though they do not influence the properties of the propagating wave states away from the boundary, these modes nonetheless make a significant contribution to the dynamics of scattering.Moreover, since these boundary modes exist only for n > 1, they are a signature of the higher bulk winding number of the Euler material.
The differential conductance of the interface may be calculated using the Blonder-Tinkham-Klapwijk formalism as described in [36]: The results of this computation for n = 1, . . ., 6 are shown for a range of Fermi and excitation energies in Fig. 4. Note that, as it must, the case n = 1 agrees with Fig. 4 of [37].Most notable however is the result for n = 2, which shows highly suppressed differential conductance over a wide range of excitation energies; indeed, it is only notably different from zero when the Fermi energy is very large compared to the superconducting gap and the excitation energy is slightly above this same energy ∆ 0 .It is qualitatively clear from Fig. 4 that the parity dependence displayed in the Andreev reflection in an Euler material (Fig. 2) is also manifest in the physically measurable quantities ∂I/∂V .For example, on the n = 3 graph in Fig. 4 it can be seen that the magnitude of the differential conductance tends towards approximately the same value in the limits ε → 0 and ε ∆ 0 .In contrast, the differential conductance in the case n = 6 tends towards to quite different values on either side of the line ε = ∆ 0 .
The plots shown in Fig. 4 are for a very large applied gate potential U 0 = 10 8 .For U 0 10 4 the corresponding plots are qualitatively similar, but the suppression in the case of n even is less severe.

III. CONCLUSION
We have demonstrated that the parity of the winding number of a band node has a strong influence on the scattering properties of the quasiparticle excitations at the boundary between a normal and a superconducting region.The suppression of the Andreev reflection coefficients in Euler materials with nodes carrying even winding numbers leads to a clear signature in the differential conductance curves that could in principle be probed experimentally.

A. S. M. is funded by an EPSRC PhD studentship (Project reference 2606546). A. B. was funded by a
Marie-Curie fellowship, grant no.101025315.R. J. S acknowledges funding from a New Investigator Award, EPSRC grant EP/W00187X/1, as well as Trinity college, Cambridge.The authors would like to thank X. Feng for helpful discussions on Andreev reflection, and for reading and providing feedback on a draft of the manuscript.definition of the polynomials P ± n to or, equivalently, Note that P ) is equal to the real (imaginary) part of the complex number (z 1 + iz 2 ) n only when both z 1 and z 2 are real.For instance if a, b, c ∈ R then

Appendix B: Derivation of scattering coefficients
In the normal and superconducting regions x < 0 and x > 0 we search for wave-like scattering solutions ψ = exp(ik x x + ik y y)(u, v) of Eq. 5 .In each case this leads to a particular eigenvalue equation for the excitation energies ε, and for the vectors u, v.

Normal side
On the normal side the pairing potential ∆ vanishes and the gate potential U = 0.The BdG equation to be solved is therefore where the excitation energy ε > 0. The characteristic equation of this block diagonal matrix factorises into the product of the characteristic equations of each block, so the eigenstates are divided into electronic states (with u = 0, v = 0), and hole states (with u = 0, v = 0).
Explicit expressions for all of these states are as follows: where in Eqs.B2c and B2f the index q = 1, 2, . . ., n − 1, and factors of e ikyy common to all terms have been omitted since they cancel in the final equations.In the above expressions the letters 'e' and 'h' indicate whether the state corresponds to an electron or a hole, and the label 'N' indicates that the wavevectors are valid in the normal region.The states ψ e,+ n and ψ h,+ n describe excitations propagating in the positive x-direction, while the states ψ e,− n and ψ h,− n propagate in the negative xdirection.The electron and hole wavevectors for these states are real and given by k where S = sgn(ε − E F ). On the other hand, the wavevectors k (N, e) x;n,q = Sn,q (ε + E F ) 1/n r n,q (α) e iφn,q(α)/2 (B3c) α )e iφn,q(α )/2 (B3d) for the states ψ e n,q and ψ h n,q are complex, so these states are evanescent and describe modes localised at the boundary.The signs ensure that these wavevectors all have positive imaginary parts so that the states decay to zero in the region x < 0, rather than blowing up as x → −∞.The existence of these n − 1 additional states can be traced back to the fact that the BdG Hamiltonian depends on |k| n , which in real space translates to an n th order differential equation -precisely n states for each energy.The speeds of propagation of the electron and hole states in the x-direction are respectively given by the states in Eqs.B2a, B2e, B2d, and B2b have been normalised by their respective speeds to ensure that they carry a single particle.(Since the states B2c and B2f are not propagating their normalisation is unimportant so has been omitted).
The angle α at which the incident electron propagates with respect to the x-axis may be expressed in terms of the conserved quantities k y and ε as and the corresponding angle for the propagating hole state is The remaining terms are c (N, e) n,q = P + n Sn,q r n,q (α) e iφn,q(α)/2 , sin(α) = P − n Sn,q r n,q (α) e iφn,q(α)/2 , sin(α) (B7b) n,q = P + n Sn,q r n,q (α )e iφn,q(α )/2 , S sin(α ) (B7c) n,q = P − n Sn,q r n,q (α )e iφn,q(α )/2 , S sin(α ) where the polynomials P ± (z 1 , z 2 ) are defined in Eq. 3, and The function arctan(x, y) is here defined to give the angle subtended between the point (x, y) and the positive xaxis, and to have a single branch cut along the line y = 0, x < 0.

Superconducting side
On the superconducting side the pairing potential ∆ = 0 and the gate potential U = −U 0 is finite, so the BdG equation becomes (B9) The characteristic equation now no longer factorises, indicating that the basic quasiparticle excitations of the system consist of coherent superpositions of electrons and holes together.We now list all of the scattering states on the superconducting side: e −iφ (s where again q = 1, 2, . . ., n − 1, and the letter 'S' indicates that we are now dealing with the superconducting side.The letters 'e' and 'h' now indicate whether the excitations are 'electron-like' or 'hole-like', that is, whether they propagate in the same or the opposite directions to their group velocities.These names are also motivated by the fact that ϕ e → ψ e and ϕ h → ψ h in the limit that ∆ 0 → 0. Again, there are two real wavevectors k (S, e) x;n,0 corresponding to the forwards-and backwardspropagating electron-and hole-like solutions ϕ e± n and ϕ h± n .The remaining wavevectors k (S, e) x;n,q = − Sn,q E 0 + ε 2 − ∆ 2 0 1/n r n,q (γ) e iφn,q(γ)/2 (B11c) r n,q (γ ) e iφn,q(γ )/2 (B11d) are complex, implying that the states ϕ e n,q and ϕ h n,q are evanescent waves.The signs Sn,q now ensure that these states decay to zero as x → +∞.

Matching
To determine the probability that the incident electron will be scattered into an outgoing electron or hole on the normal side, or an electron-or hole-like quasiparticle in the superconducting side, we propose a trial ansatz on both sides and determine the amplitudes of each state in the solution.The incoming state on the normal side is a sum of a right-moving electron, a left-moving electron, a left-moving hole, and all evanescent states: x;n,0 + rψ e− n e −ik (N, e) x;n,0 + r A ψ h− n e −ik (N, h) x;n,0 a q ψ e n,q e ik (N, e) x;n,q + b q ψ h n,q e ik (N, h) x;n,q .(B14a) On the superconducting side, the ansatz is a sum of an electron-like and a hole-like quasiparticle propagating to the right, and all evanescent states: ϕ n = tϕ e+ n e ik (S, e) x;n,0 + t ϕ h+ n e ik (S, h) x;n,0 a q ϕ e n,q e ik (S, e) x;n,q + b q ϕ h n,q e ik (S, h) x;n,q .
(The factor e ikyy common to all terms has again been omitted in Eqs.B14.)To determine the coefficients, we apply the boundary/matching conditions a q ik (N, e) x,q p ψ e,q + b p ik (N, h) x,q p ψ h,q (B16) = t ik c q ik (S, e) x,q p φ e,q + d q ik (S, h) x,q p φ h,q for p = 0, 1, . . ., n − 1; this is a set of 4n equations which may be solved to yield the 4n coefficients r, r A , a 1 , . .(B18) and the 4 × (4n + 1)-dimensional matrix of states Ψ = ψ e,+ ψ e,− ψ h,− • • • φ h,p .The vector v T 0 is used to set the normalisation I 0 = 1.The matrix A is non-singular so this yields a unique solution for each of the coefficients r, r A , etc.The expressions resulting from the solution of this set of equations are rather involved, so we omit them; they may be obtained using a computer algebra program.Importantly, when the incident angle α is greater than a critical angle the hole states are no longer valid scattering states.After the final result has been obtained, all of the coefficients r A , b 1 , . . ., b n−1 must therefore be multiplied by θ(α−α c ), where θ(x) is the Heaviside step function.

Figure 1 .
Figure 1.Low energy excitation spectrum of the Hamiltonian Hn(k) of Eq. 1 for n = 1 (left) and n = 2 (right).The field in the kx-ky plane displays the vector (cos θ(k), sin θ(k)), which winds n times around the origin.

Figure 4 .
Figure 4. Variation of the differential conductance ∂I/∂V with the Fermi and excitation energies EF and ε with a large external gate potential U0 = 10 8 .