Tunneling conductance due to discrete spectrum of Andreev states

We study tunneling spectroscopy of discrete subgap Andreev states in a superconducting system. If the tunneling coupling is weak, individual levels are resolved and the conductance $G(V)$ at small temperatures is composed of a set of resonant Lorentz peaks which cannot be described within perturbation theory over tunnelling strength. We establish a general formula for their widths and heights and show that the width of any peak scales as normal-state tunnel conductance, while its height is $\lesssim 2e^2/h$ and depends only on contact geometry and the spatial profile of the resonant Andreev level. We also establish an exact formula for the single-channel conductance that takes the whole Andreev spectrum into account. We use it to study the interference of Andreev reflection processes through different levels. The effect is most pronounced at low voltages, where an Andreev level at $E_j$ and its conjugate at $-E_j$ interfere destructively. This interference leads to the quantization of the zero-bias conductance: G(0) equals $2e^2/h$ (or 0) if there is (there is not) a Majorana fermion in the spectrum, in agreement with previous results from $S$-matrix theory. We also study $G(eV>0)$ for a system with a pair of almost decoupled Majorana fermions with splitting $E_0$ and show that at lowest $E_0$ there is a zero-bias Lorentz peak of width $W$ as expected for a single Majorana fermion (a topological NS-junction) with a narrow dip of width $E_0^2/W$ at zero bias, which ensures $G(0)=0$ (the NS-junction remains trivial on a very small energy scale). As the coupling $W$ gets stronger, the dip becomes narrower, which can be understood as enhanced decoupling of the remote Majorana fermion. Then the zero-bias dip requires extremely low temperatures $T\lesssim E_0^2/W$ to be observed.


I. INTRODUCTION
Scanning Tunneling Microscopy (STM) is a powerful method of studying electronic surface properties of condensed matter systems. Electrons from the metallic tip of the microscope tunnel into the system, creating a stationary current at fixed voltage. At low temperatures the conductance G of the tunnelling contact is proportional to the surface density of states ρ beneath the tip, and the coupling strength between the tip and the system G(V ) ∝ |t 2 |ρ(r, eV ). (1) This relation allows to probe local density of states with energy and spatial resolution. In recent works 1-14 transport measurements have been suggested to probe the Majorana fermion in superconducting structures with non-trivial topological properties, such as nanowires combining strong spin-orbit, induced superconductivity and magnetic field 2,4,5 , and vortex cores on topological insulator surfaces with induced s-wave superconductivity 7,10,15,16 . There are two major transport phenomena suggested, the first one is the fractional Josephson effect 1,7,14,17 : the current-phase relation in an SNS-junction hosting a pair of Majorana fermions is 4π-periodic with phase if the fermionic parity in the junction is conserved 18 . This effect also implies the doubling of Shapiro steps in an ac measurement 11 , which has recently been reported 19 . A second, and more direct effect of the Majorana state on transport is the peak in zero-bias conductance. Transport experiments have been done on nanowire-based setups [20][21][22] , showing a zero-bias peak, that is robust to the change of magnetic field (in both magnitude and direction) and gate voltages. Un-der the conditions of these experiments, the Majorana state is localized and thus is expected to produce a conductance peak proportional to the peak in the density of states, in agreement with Eq. (1). For the vortex core, direct STM measurements have been proposed 10,16 with a similar prediction of a zero-bias conductance peak. The idea of probing the Majorana bound state with tunneling transport measurements motivated us to revisit the conventional relation (1) for the case of a Majorana level and, more generally, for an arbitrary discrete spectrum of Andreev bound states. In this work we study the conductance through a tunneling contact between a normal metal and a superconductor hosting a subgap spectrum of Andreev bound states E j . While the G(V ) characteristic is expected to reflect the structure of the spectrum and exhibit peaks at voltages close to E j /e, it is clear, that the classical expression (1) cannot accurately describe the current in our system of interest. Since we consider a discrete spectrum, the density of states ρ(E) is a set of δ-functions. The conductance however, cannot be a set of δ-functions as it is limited by the number of channels in the contact. Moreover, expression (1) is a Fermi Golden rule result, that describes electrons slowly leaking the probed system. In our superconducting system there are no extended states at subgap energies and single electrons cannot leak into it. The only possibility for a stationary current is due to Andreev reflection: the process of electrons being reflected from the NS-junction as holes. The discrete Andreev levels in the superconducting system act as quasi-stationary levels for this process: an incident electron may occupy an energetically favorable level and is then emitted back into the metal tip either as an electron or as a hole. At voltages close to one of the levels E j0 this process should exhibit resonant behavior, and cannot be described perturbatively in terms of the coupling strength. Our aim is to describe these resonances for a generic tunneling contact.
Another interesting problem of the discrete spectrum is interference of Andreev reflection processes involving different levels. There is no interference involved in the simple formula (1), where the conductance is directly proportional to the number of available energy levels on the probed surface. The process of Andreev reflection is different from conventional tunneling as it involves tunneling into the probed system and then back into the metal. Reflection of the electron into a certain hole channel can happen in multiple ways, that differ in the Andreev level E j they tunnel into. The amplitudes of these different ways have to be summed together, which allows for interference between different levels. The interference effects should be important for a few-channel contact, especially a single-channel one, where all amplitudes are added together. Most interesting is the case of low voltages, where the BdG-symmetry of the spectrum becomes important and defines the interference between conjugate levels ±E j .
In section II we formally introduce the system we study and the formalism we will use. This includes the Landauer formula for superconductors, a Kubo-type formula for Andreev reflection, and theŤ -matrix expression of exact Green's functions in terms of the tunneling Hamiltonian and unperturbed Green's functions of the connected systems. In section III we study the resonance peaks that occur at voltages close to a single Andreev level E j0 . We derive general formulae describing the peaks and study the particular cases of a peak on a Majorana level, formulae for a point-contact and the conductance at finite temperatures. In section IV we establish an exact formula for conductance in the single-channel contact case that takes the whole spectrum E j into account. We then study interference of different levels and observe quantization of zero-bias conductance. Next, we study conductance at finite voltages. In particular, we consider a system with a pair of weakly coupled Majorana levels. We find that a zero-bias peak in this system coexists with an extremely narrow dip at E = 0 that shows that the system is still topologically trivial. We discuss the experimental requirements to observe the dip structure of this peak and show that, surprisingly, they become more severe when the tunnelling coupling to the probe increases. In section V we discuss some experimentally realized systems where level spacing of Andreev states is large enough for the individual levels to be observed.

II. SETUP
The system we are studying consists of a metallic tip described by the Hamiltonian H M , a superconducting system H S and a tunneling Hamiltonian H T that allows electrons to tunnel from the tip onto the system.
We are interested in voltages and temperatures much lower than the superconducting gap and thus only take the Andreev levels into account: The contact is described by the tunnelling Hamiltonian where |θ m , |σ s describe some electronic wave functions in the metallic tip and superconductor, correspondingly.
In the Bogolyubov-de Gennes formalism we will use, H T has to be extended into Nambu spacě whereΘ is the time-reversal operator, and |θ * m is the particle-hole-conjugated partner of |θ m , i.e. |θ * m = C|θ m with C denoting particle-hole conjugation. Throughout the paper the check sign indicates a structure in Nambu space.
At subgap energies E < ∆ electrons incident from the metallic tip are reflected back either normally or as holes. The latter process is Andreev reflection and produces a current through the NS-junction, which can be described in the Landauer formalism. The formula describing the stationary current at some voltage bias reads (we derive it from the general Lesovik-Levitov formula 23 for the BdG-Hamiltonian) Note that the factor of 2e is due to the quantum of transferred charge in the Andreev process. f 0 (E) denotes the equilibrium Fermi distribution function, and T A (E) is the sum of Andreev reflection eigenvalues: where S eh (E) is the electron-hole block of the S-matrix of the N-S junction. With the help of the relations T A (E) = T A (−E) and f 0 (E) + f 0 (−E) = 1 formula (6) can be rewritten in a form similar to the one for the standard N-N contact current: Below we mainly study the zero-temperature conductance G(eV ) = T A (E = eV ) · 2e 2 /h which is just the sum of Andreev transmission eigenvalues times the relevant conductance quantum 2e 2 /h. Conductance at T > 0 can then be easily obtained from the zero-temperature result via the convolution (9). The trace (7) can be written as which can be understood as an Andreev-reflection analog of the Kubo formula for the linear response.v is the velocity operator for electrons andv * =ΘvΘ −1 is its time-reversal. The explicit matrices inside the trace (10) act in Nambu space, andǦ A,R E are exact advanced and retarded Green's functions at energy E. The exact Green's functions in the metal are exressed through unperturbed Green' here and in what follows Green's functions without index are assumed retarded): The Green's functions here act in the Nambu space, sǒ We use g to denote normal metal Green's functions without Nambu structure and g * =ΘgΘ −1 for its timereversal image. The superconducting Green's function isǦ The sum over j > 0 reflects the symmetry of the levels -each positive level |j has a particle-hole conjugated partner |j * with opposite energy. G S has poles at E = ±E j and theŤ -matrix cannot be treated perturbatively in their vicinity. All terms of the geometrical sequence have to be included in the expression (12). This can be done in a compact way in two cases we consider below.

III. SINGLE-LEVEL RESONANCE FOR ARBITRARY HT
TheŤ -matrix and consequently T A can be found explicitly in a closed form in the resonance situation |E − E j0 | ≪ |E − E j =j0 | when the energy of the incident electron is closer to the Andreev level E j0 than to any other level E j . We only keep the (E − E j0 ) −1 term inǦ S in this case, and writě Here |τ is a function of the electron coordinates in the metallic tip, which is determined by the contact Hamiltonian H T and the resonant Andreev state |j 0 . It has both electron and hole components which we denote as |τ e and |τ h for later use. The simple projector form of the operator (15) leads tǒ The real part of τ |Ǧ 0 |τ is zero, and we denote its imaginary part as W . Substituting theŤ -matrix (16) into formula (11) and the result into the Kubo formula (10), we find (here |τ * h =Θ|τ h is just the time-reversed wave function) In the case of a clean metallic tip (and assuming that Fermi energy E F is the largest energy scale), the gvg-terms can be simplified in terms of the n e,h densities: Formula (19) constitutes our main result for the singlelevel resonance, and we will now discuss it in more detail. The V -dependence of (19) has a Lorentz shape, which is generic: at resonance, the S-matrix is governed by the pole at E = E j0 + iW . The width W describes the broadening of the Andreev level at the energy E j0 due to coupling to the tunneling probe. T * A ≤ 1 is the probability of Andreev reflection for an incident electron with E = E j0 . It is defined by the relation between the electron and hole components of |j 0 and the geometry of the contact. If |j 0 was a purely electronic state T * A ∼ n h = 0 and there would be no Andreev reflection from |j 0 and hence no conductance. For deep subgap states (E 0 ≪ ∆) electron and hole components are very similar (e.g. connected by phase rotation for a uniform ∆(r)) so we expect 1 − T * A ≪ 1. However, such a feature may be absent when measuring tunneling conductance at points where the gap ∆ is locally suppressed, ∆(r) ≤ E 0 , e.g. near the origin of a vortex core, or in the middle of an SNS-junction. In such regions electron and hole components may differ considerably, leading to a small T * A . The presence of the T * A factor in (19) is an important difference from the classical relation (1): Andreev reflection amplitude depends on the e−h structure of the wave function, and not just on the overall probability density.

A. Majorana peak
An important special case of the single-level resonance described by (19) is the resonance on a Majorana bound state. In this case |j 0 is self-conjugate, so that |τ e = |τ * h . Hence n e = n h and T * A = 1, This agrees with earlier derivations of the Majoranainduced zero-bias tunneling conductance peak 3,6 . The width W of the peak still depends on contact and wave function geometry, while the height is fixed to 2e 2 /h by the symmetry of the Majorana state. Other energy levels |j = j 0 generally lead to a negative correction to the conductance of the order of O(W 2 /ω 2 0 ) where ω 0 is the typical level spacing of the Andreev spectrum in Eq.(3).

B. Point-contact
The case most relevant to a scanning tunneling microscope setup is the point contact H T = |θ t 0 σ| + h.c. where wave functions |θ and |σ are localized at some points r M and r S in the metallic tip and superconductor correspondingly. For simplicity we assume here that H T commutes with spin, but our results are easily generalized to a spin-sensitive setup.
Here ν M is the density of states in the metallic tip per spin projection and |ψ e (h) 2 | = |ψ 2 e(h)↑ | + |ψ 2 e(h)↓ | is the probability density of the electron (hole) component of the resonant Andreev level. We can rewrite (21) for the point-contact as which has the form of a standard expression for normal conductance times the resonant reflection probability T * A . ν in Eq. (24) is a broadened density of states of the Andreev level |j 0 : The point-contact expression (24) can be generalized to an arbitrary contact: whereǦ S is the broadened Green's function of the probed system:Ǧ If instead ofǦ S in (26) we substitute the Green's function corresponding to the normal state of our probed system, and put T * A = 1, the expression (24) will yield the normal-state conductance. In the superconducting regime in presence of the Andreev spectrum, the Green's function pole has to be broadened according to Eq. (27), so that δ-peaks in G R − G A are smeared into Lorentz peaks of the width W . We remind that the parameter W is different for different peaks according to Eq. (18). Thus, if W is calculated for |j 0 , then expression (27) with that W only works for the peak at eV = E j0 .

C. Temperature dependence
Given the zero-temperature G(T = 0, eV ) from (19), the finite-temperature conductance is readily obtained via formula (9) by convoluting it with the thermal distribution function. The shape of G(T ) starts to differ from the G(T = 0) case as soon as T becomes comparable to W . At this point the peaks become broader, while retaining spectral weight. At ω 0 ≫ T ≫ W (we denote the level spacing by ω 0 ) the temperature-broadened resonant peak looks like Here W enters as the spectral weight of the original Lorentz peak (19), while the peak shape is defined by the temperature T only. If instead of expression (19) we use expression (24) or (26) in the convolution (9), we find that the broadening of the δ-peaks by W is much less than thermal broadening and thus drops out. For the point-contact IV. SINGLE-CHANNEL CONTACT,

MULTI-LEVEL SYSTEM
Another case which can be solved exactly is presented by a single-channel contact H T coupled to an arbitrary system of Andreev levels H S . In this section we account for all Andreev levels simultaneously in an exact way, which allows us to study important interference effects due to Andreev reflection via different levels of the discrete spectrum.
An arbitrary single-channel tunneling Hamiltonian can be written as h T = |θ t σ|+h.c. with a real t and normalized electronic wave functions |θ and |σ that are defined in the metal and superconductor respectively. To include h T into formula (12), we extend it into Nambu space: Thus, in the BdG formalismȞ T appears as a two-channel contact. The calculation of theŤ -matrix is similar to that of the expression (16). Using (30) and (14) we writě Here N 0 = Im θ|g R E |θ is a normalization factor (for a point-contact N 0 = πν). Σ A has a simple physical meaning: 2Σ A is the amplitude of Andreev reflection (up to some phase factor) in the lowest order in H T . Similarly, Σ e,h describe the corrections to electron-electron and hole-hole reflection due to the presence of H T . Using Eq. (31) we now construct theŤ -matrix: where the 2 × 2 matrix in Eq.(36) reads aš Finally we substituteǦ =Ǧ 0 +ŤǦ 0 (andǦ A =Ǧ R † ) withŤ -matrix (36) into the Kubo formula (10). Only off-diagonal (in Nambu space) terms ofǦ R,A and consequently only off-diagonal terms of (37) enter the trace, yielding Voltage eV enters the formula throughΣ where E is set equal to eV . At resonance, |eV − E j0 | ≪ W , formula (38) reproduces our single-level result (19) with n e = | σ|j | 2 t 2 N 0 /π and n h = | σ|j * | 2 t 2 N 0 /π. Away from any resonance,Σ ≪ 1, so that G(eV ) ≃ 2e 2 /h × 4|Σ A | 2 , which is the lowest order perturbative result in powers of H T . Σ A is a sum of Andreev amplitudes over different levels j, and at some voltages several levels contribute to Σ A terms of comparable magnitudes, which may interfere. These interference effects are most drastic at lowest voltages, where interference pushes G(0) to be equal either to 0 or to 2e 2 /h, as we discuss below.

A. Zero-bias conduction quantization
The fact, that the zero-bias conductance of a singlechannel NS-junction is quantized follows from general symmetry consideration 24,25 . The S-matrix describing an NS-junction has BdG-symmetry which relates S * E and S −E . In conjunction with unitarity S † S = 1 this fixes det S 0 = ±1, dividing possible S-matrices into two disconnected classes: trivial (det S 0 = 1) and non-trivial (det S 0 = −1) . The non-trivial case is represented by topological superconductor-metal junctions supporting a Majorana level on the junction 2 , (this Majorana level is a well-defined quasistationary Andreev level if there is a barrier in the metal in the vicinity of the contact to superconductor, otherwise the Majorana level is leaked into the normal metal). In the single-channel case S is a 2 × 2 matrix and G(0) = (1 − det S 0 )e 2 /h. Thus G(0) is 0 (2e 2 /h) in the trivial (non-trivial) case.
Let us now show that expression (38) is indeed quantized at E = 0. Suppose first there are no zero levels among E j . Then, Σ A = 0 since terms with |j and |j * in (35) cancel each other out at E = 0, and consequently, The two above cases (no Majorana states, or a single one) describe all physical situations, as any higher number of zero modes would split. Let us show how it happens in somethat more details. One possible case corresponds to a non-Majorana zero level with E 0 = 0 and |j 0 = |j * 0 , which can also be understood as a pair of uncoupled Majorana fermions |j 0 = e iϕ (|γ 1 + i|γ 2 ). Depending on the choice of ϕ, the state |j 0 can be decomposed into different pairs of Majorana fermions γ 1 , γ 2 . If | σ|j 0 | = | σ|j * 0 | we choose e 2iϕ = σ|j 0 / σ|j * 0 and find that the correspoding γ 2 is not coupled to the contact, σ|γ 2 = 0. Thus, γ 2 does not enter the hamiltonian H = H T + H M + H S and we end up with the case of a single Majorana fermion γ 1 . In the second case, | σ|j 0 | = | σ|j * 0 | both γ 1 and γ 2 enter H T for any choice of ϕ. The Majorana fermions are effectively coupled by the contact into a finite-energy fermionic level, which drives our NS-junction into the trivial regime. Indeed, if | σ|j 0 | = | σ|j * 0 | the 1/E 4 singularities in the denominator of expression (38) do not cancel out so that G(0) = 0, as expected for a trivial junction. Finally, one may ask what happens if there are multiple zero levels j 0 , j ′ 0 . . . . This problem reduces to those considered above by changing the basis of the zero levels so that only one level couples to the contact, which is always possible for a single-channel contact.
Note that exact quantization only happens in the single-channel case. For instance the exact cancellation of |j and |j * Andreev amplitudes is only possible when the incoming electron and outgoing hole are in the same channel. This suggests that there may be some traces of the destructive interference for few-channel junctions, which become negligible as the number of channels becomes large.

B. Conductance at finite voltage and interference
The quantization discussed above only concerns E = eV = 0 and cannot provide any information about finite energies. Our formula (38) allows to find how narrow the quantization region is, i.e. at which voltages (or temperatures) it should be visible. The study of (38) at finite E also resolves an apparent paradox similar to that discussed in 25 : consider a finite superconducting system designed to host Majorana fermions. e.g. 2,15 . It is topologically trivial, because there is an even number of Majorana fermions across the system and there is finite energy splitting between them. So, det S 0 = 1 and G(0) = 0. Let us now extend the system to move one Majorana level to infinity. It is now uncoupled from the system and we are left with a topological junction which should show a zero-bias conductance peak. However, from continuity S 0 should remain in the trivial class and G(0) should stay at zero. Below we show that both statements are true: there is a peak at low energies, and, at the same time, The low-voltage conductance is governed by the lowest energy level |j and |j * so we omit other levels in what follows -they only lead to minor corrections in G. We find The shift in the resonance energy in (40) is negligible, as it is smaller than E 0 by a factor of W/ω 0 (except for the exotic case mentioned earlier when there are two uncoupled Majorana fermions in the superconducting system, which are then hybridized by the tunneling contact; in this caseẼ 0 is the effective energy splitting and cannot be neglected). The two-level formula (39) agrees with what has been derived earlier in a different way for the special case when H T couples the metal tip to a single Majorana operator 6 . If E 0 ≫ W (gray curves on Fig.  1,2), the function (39) shows two well-separated Lorentz peaks at ±E 0 with a width W each. If we decrease E 0 , these two peaks move towards each other, but G(0) stays zero. When E 0 becomes comparable to W (blue curves on Fig. 1,2) the shape does not look like a pair of Lorentz peak anymore. In the region |eV | < E 0 the conductance is suppressed due to destructive interference, while at eV > |E 0 | the conductance is enhanced: at eV ≫ E 0 one gets G(eV ) = 2e 2 /h · 4W 2 /(eV ) 2 which is twice the sum of two Lorentz tails. Finally, at E 0 ≪ W (red curves on Fig. 1,2), we see one accurate Lorentz peak of double width 2W with a parametrically narrow dip in the region E < ∼ E 2 0 /W all the way to 0. While quantization G(0) = 0 is fulfilled, it only affects a very narrow region of energies. Note that throughout the change of E 0 the current I(V ≫ E 0 , W ) = ∞ 0 G(V )dV (proportional to the area under the curves on Fig. 1) remains the same, and equals 2e/h·πW . Hence, at high temperatures T ≫ E 0 , W the interference effect upon the conductance is exponentially small.
The above consideration solves the apparent paradox: the region of voltages and temperatures where the trivial character of the system is observable is very narrow, eV, T < ∼ E 0 . Interestingly, increasing the coupling strength W to the probe makes the observation of the conductance dip even more difficult -at W ≫ E 0 one needs to resolve energies of the order of E 2 0 /W instead of just E 0 , as can be seen from Eq. (39) and Figs. 1,2,3. This can be understood in terms of the specific topological transition of the S-matrix described by Pikulin and Nazarov 25 . Our S-matrix has two poles (which are also poles of function (39)) at E ± = ± E 2 0 − W 2 + iW . For a weak contact (gray curve on Fig. 2) they are connected by BdG-symmetry, E + = −E * − and represent a quasistationary Andreev level at E + and its image at E − . As W is increased and the coupling gets stronger the poles move towards the imaginary axis. At W = E 0 (blue curve on Fig. 2) they meet at the imaginary axis and  At T = 0.1 (blue curve) which is of the order of the zerotemperature dip width E 2 0 /W , the dip is significantly smeared out, while the wings of the Lorentz peak remain unchanged. At T = 1 (gray curve) which is of the order of E0 there is only a very shallow dip left. Finally, at T = 4 (black dashed curve), which is comparable to W = 10, there is no visible trace of the dip and the Lorentz peak starts to broaden. the transition 25 happens. At higher W one of the poles moves up the axis, the other down and each of them is now BdG-self-conjugate: E ± = −E * ± . At very high W (red curve on Fig. 2) the lower pole is "buried" 25 -it sits parametrically close to the real axis, E − = iE 2 0 /2W , and is very weakly coupled to the probe. This means, that the better the contact, the more our NS-junction looks topological from the electronic transport point of view. This is an important observation -while the quantization of G(0) happens for any W , a coupling stronger than the splitting E 0 diminishes the temperature and voltage range in which the conductance dip (and hence the remote second Majorana fermion) is observable, as shown in Fig. 3.

V. DISCUSSION
We now discuss real systems where the effects studied can be observed. A natural object to look at is a core of an Abrikosov vortex in a quasi-2D superconductor. The above results are applicable if the Caroli-De Gennes-Matricon (CdGM) levels 30,31 accessible by the tunnel contact on the top plane of the superconductor are not spreaded into the continuous band dispersing along the vortex core. The CdGM spectrum of a clean vortex core is described by with ω 0 ∼ ∆ pF ξ and m * = m z (p F ⊥ ξ) 2 , where µ is the angular momentum, ξ is the coherence length and m z is the effective mass along the z-direction assumed to be larger than (x, y)-plane effective mass m in an anisotropic superconductor. In the presence of even a low concentration of impurities, CdGM states in a neighboring atomic planes of layered material will be displaced by the amount of energy ∼ ω 0 , whereas bandwidth of each level is about µω 0 m/m z . Therefore, in strongly anisotropic materials with m z ≫ m individual low-lying CdGM levels are well localized 32,33 and tunnelling from the STM tip into the surface of such a material may show the features discussed above in this paper. However, in order to discern these features, one needs to work at very low temperatures T ≤ ω 0 . Currently, it should be possible to reach this condition while dealing with a material like NbSe 2 , whereas older data 34,35 were taken at too high temperatures.
Another possibility is to use genuine two-dimensional electron systems like graphene or a surface of topolgical insulator. In these materials there are no transverse quantum numbers and the Andreev spectrum is discrete with level spacing ω 0 (in graphene the spectrum is fourfold degenerate due to valley and spin degrees of freedom). An additional advantage of these two-dimensional systems is a controllable E f . The charge neutrality point in graphene coincides with the K-point of its Dirac spectrum, which allows to have ω 0 comparable to E f . An important issue is the implementation of superconductivity in these two-dimensional systems. A possible solution is to cover them with a conventional superconductor to induce superconductivity by proximity effect 15,[26][27][28] . However, one needs to make sure that there are no parasite Andreev states in the covering superconductor. This can be done by making a hole in the superconductor of a radius of the order of the coherence length 10,16,29 physically removing the normal vortex core with all its low-lying Andreev states from the superconductor. In addition, such a hole should facilitate the application of a tunneling microscope to the graphene or topological insulator surface.

VI. CONCLUSION
In conclusion, we have studied resonant Andreev reflection in an arbitrary tunneling N-S contact between a metal and a superconductor hosting a discrete Andreev spectrum E j . If the level spacing is larger than the inverse dwell time W , the subgap conductance G(V ) consists of Lorentz peaks centered at eV = E j . We have derived formulae describing these peaks in the general case. The heights of the peaks do not scale with the coupling strength and are generally very close to 2e 2 /h, while widths scale as normal conductance of the contact. We also established an exact formula describing the Andreev conductance for a generic single-channel contact that takes into account all Andreev levels at once. This formula reveals the role of interference between Andreev reflection processes involving different levels E j . At eV = 0 interference between the level E j and its particle-hole image −E j is destructive, so that G(0) is exactly 2e 2 or 0 depending on whether there is a Majorana fermion in the spectrum, in agreement with Smatrix theory of topological NS-junctions. We have also studied G(V ) for finite energies for different W and lowest energy level E 0 and explained what happens when the system is gradually driven from the trivial to the topological phase by removing a Majorana fermion from the system and thus driving E 0 to zero. We have found that at E 0 ≪ W the conductance shows a zero-bias Lorentz peak of the usual width 2W , with a very narrow dip on top of it, which secures G(0) = 0. The parametrically small width E 2 0 /W of the dip can be understood as the "burying" of the remote Majorana fermion: the strong coupling of the tunneling tip to one of the Majorana states effectively decreases the coupling of the system to the other Majorana state. To observe the dip (i.e. the "buried Majorana fermion") experimentally, temperature and energy resolution has to be smaller than E 2 0 /W . The effects we discussed should be observable in two-dimensional electronic systems including graphene or the surface of three-dimensional topological insulators with proximity-induced superconductivity. Another possibility is to study the level structure of Abrikosov vortices in a layered strongly anisotropic superconductor like NbSe 2 at ultra-low temperatures.
We thank Alexey Ioselevich, Pavel Ostrovsky and Dmitry Pikulin for very helpful discussions. This work was supported by the RFBR grant # 10-02-00554-a and by the RAS program "Mesoscopic and disordered systems".