{'fraction_non_alphanumeric': 0.03288609616028431, 'fraction_numerical': 0.010095027169014448, 'mean_word_length': 4.2933478735005455, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Majorana state on the surface of a disordered 3D topological insulator

We study low-lying electron levels in an"antidot"capturing a coreless vortex on the surface of a three-dimensional topological insulator in the presence of disorder. The surface is covered with a superconductor film with a hole of size R larger than coherence length, which induces superconductivity via proximity effect. Spectrum of electron states inside the hole is sensitive to disorder, however, topological properties of the system give rise to a robust Majorana bound state at zero energy. We calculate the subgap density of states with both energy and spatial resolution using the supersymmetric sigma model method. Tunneling into the hole region is sensitive to the Majorana level and exhibits resonant Andreev reflection at zero energy.

Topological insulators and superconductors are very peculiar materials with a gap in the bulk electron spectrum and a low-lying branch of subgap excitations on their surface (see [1] for a review). This surface metallic state appears due to topological reasons and is robust with respect to any (sufficiently small) perturbations. In particular, topological properties prevent these surface states from Anderson localization. One common example of a topological insulator is a two-dimensional system in the integer quantum Hall effect regime. The bulk of such a system has a spectral gap between successive Landau levels and is hence an insulator. At the same time quantized Hall conductance appears due to a fixed integer number of chiral propagating edge modes on the background of the bulk gap. This type of materials are referred to as Z topological insulators.
Another type of topological insulator is realized in three-dimensional (3D) semiconductor crystals with sufficiently strong spin-orbit interaction (BiSb, BiTe, BiSe, strained HgTe etc). The spin-orbit interaction leads to inversion of the spectral gap. As a result subgap surface excitations appear with a dispersion of the massless Dirac type. The topological invariant in these materials has a Z 2 nature. When the number of surface states is odd, one of them always remains gapless due to topological protection.
A general classification of topological insulators was developed in Refs. [2,3] based on the symmetry of the underlying Hamiltonian. The quantum Hall effect is an example of a topologically nontrivial state of the unitary symmetry (class A of the Altland-Zirnbauer classification [4]) in 2D. Strong spin-orbit interaction leads to a topological state in the symplectic class (AII) in 3D. Another important example is the 1D topological superconductor of the class BD symmetry (superconductor with both time-reversal and spin rotation symmetries broken). The topologically protected mode in this case is zero-dimensional and is known as the Majorana bound state (MBS). It appears, in particular, in the core of an Abrikosov vortex in the spinless p-wave superconductor [5]. A similar MBS appears [6] in a vortex in an ordinary s-wave superconductor brought in contact with the surface of a 3D Z 2 topological insulator, which corresponds to the symmetry class AII. The MBS appears as a descendant of the topologically protected massless Dirac state on the free surface of the topological insulator.
There are two general methods to observe the Majorana level. One of them relates to the anomalous Josephson effect [7,8]. Another, and a more direct, way involves tunneling into the region where the MBS is supposed to be localized [9]. The differential conductance in such a tunneling experiment yields the local density of states with spatial and energy resolution. The Majorana state in the vortex core manifests itself by resonant Andreev reflection at zero energy [10].
In this Letter we study local tunneling conductance for a setup depicted in Fig. 1. A superconducting (with a spectral gap ∆) film with a circular hole of the radius R is deposited on the surface of a 3D topological insulator, e.g., Bi 2 Te 3 crystal. Perpendicular magnetic field applied to the system produces a vortex pinned to the hole. The radius R is supposed to be relatively large, R ≥ (l, ξ 0 ), where l is the mean free path for topological insulator surface states and ξ 0 is the superconducting coherence length of the film. The condition R ≥ ξ 0 allows one to avoid [11] the abundance of low-lying Carolide Gennes-Matricon states whose presence in the core of the Abrikosov vortex complicates observation of the E = 0 Majorana state. Under these conditions disorder is relevant for the states localized in the hole region. Although the Majorana state is protected against disorder, i.e. its energy stays exactly zero, the spatial distribution of its wavefunction ψ 0 (r), and thus the local tunneling conductance are sensitive to disorder. The aim of this Letter is to calculate the tunneling conductance in the "dirty-limit" with spatial and energy resolution, and to identify the effects of the Majorana zero-energy state in the presence of strong disorder.
The problem is described by the following Bogolyubovde Gennes (BdG) Hamiltonian in polar coordinates r, ϕ: (1) Here the Hamiltonian H 0 describes the dynamics of surface excitations in a topological insulator without a superconducting layer. The Fermi velocity of surface electrons is denoted by v 0 , s is the spin operator, and V (r) is the random disorder potential. The Fermi level is shifted from the Dirac point by µ. The vector potential term in H 0 is neglected due to smallness of magnetic field. We assume the dirty limit with a disorder induced mean free path l ≪ ξ = D/∆ and a relatively large hole with R ≫ ξ; here D = v 2 0 τ /2. These conditions allow us to use a step-like radial dependence of the order parameter: Inside the hole, the order parameter is zero and electron dynamics is governed by H 0 . This Hamiltonian possesses time-reversal symmetry of the symplectic type, H 0 = s y H * 0 s y , and hence belongs to the symplectic symmetry class AII. Proximity to the superconductor induces a gap in the electronic spectrum outside the hole. This gap effectively confines the low-lying excitations and imposes boundary conditions for the Hamiltonian H 0 at r = R. These boundary conditions break time-reversal symmetry due to the spatially rotating phase of the order parameter. At the same time the total BdG Hamiltonian H acquires a specific particle-hole symmetry: Here τ y is the Pauli matrix acting in Nambu-Gor'kov (NG) space. Identity (3) implies a symmetry of the electron spectrum. For each eigenstate ψ with energy E there is a conjugate eigenstate s y τ y ψ * with energy −E.
The particle-hole symmetry (3) defines the class BD. On the level of random matrices there are two distinct versions of this class with even and odd number of eigenstates referred to as D and B class, respectively. In class B, one unpaired eigenstate has exactly zero energy and is self-conjugate: ψ = s y τ y ψ * . The BdG Hamiltonian counts every physical excitation twice due to the doubling in NG space. An unpaired level is thus "half" of a true excitation -a Majorana state.
We will calculate the density of states inside the hole with the help of the supersymmetric non-linear sigma model. Two-dimensional Dirac fermions with potential disorder are described by a very peculiar model of the class AII with Z 2 topological term [12]. This topological term appears as a consequence of the chiral anomaly of Dirac fermions. We will consider the minimal model operating with the 8 × 8 supermatrix Q. Apart from Fermi-Bose superspace, this matrix operates in the space of retarded-advanced (RA) fields and in a specific timereversal (TR) space introduced to take into account Cooperons and diffusons on equal footing. RA space is completely analogous to the NG space (see below); we denote Pauli matrices in this space by τ . Notation σ is used for Pauli matrices in TR space. A detailed derivation of the sigma model can be found in [13].
In class AII, the matrix Q obeys the non-linear constraint Q 2 = 1 and the linear constraint As a result, Q contains 8 commuting and 8 anticommuting (Grassmann) variables in total. Commuting degrees of freedom parameterize diagonal fermion-fermion (F) and boson-boson (B) blocks of Q. The general Q-matrix can be decomposed as where the central part contains only commuting variables while Grassmann parameters define the unitary supermatrix U . We will use the following explicit form of the central part of Q in terms of eight angle parameters: This representation fulfills all the constraints imposed by the symmetry class AII of the Hamiltonian H 0 . The F and B sectors of the sigma model are compact and noncompact, respectively. This is achieved by demanding that the angles θ b , k b are imaginary while all other angles in Eqs. (6), (7) are real. Below we will find that the saddle point describing the density of states in the hole occurs with real θ b . This implies a proper shift of the integration contour for this angle. The sigma model of class AII is designed to study transport properties of a disordered system. This implies that the Q matrix operates, in particular, in the RA space allowing for averaging the product of retarded and advanced Green functions. We are interested in the density of states and hence it suffices to average just the single retarded function. At the same time, superconducting boundary conditions implemented in the BdG Hamiltonian (1) require to introduce an additional doubling of fields in the Nambu-Gor'kov space. This can be achieved within the standard AII class sigma model with the 8 × 8 supermatrix Q while the role of NG space is taken over by the RA structure of Q (for detailed discussion of the transformation from RA to NG representation see [13,14]).
Thus we can incorporate the superconducting order parameter directly into the action of the sigma model, Here ǫ = E + iG t δ(r − r 0 )/4πν is the sum of the energy E and the local dwell term describing the coupling to a tunnel tip with dimensionless conductance G t ≪ 1 [13]. The action (8) involves the Z 2 topological term S θ [Q] which appears due to massless Dirac nature of underlying electrons as was discussed above. The topological term involves only the compact part of the sigma-model manifold, i.e., its F sector M F . In the general version of the sigma model, that is capable of averaging several retarded and advanced Green functions, the homotopy group is π 2 (M F ) = Z 2 . However, in the minimal model we are considering, M F has the structure of the product of two spheres S 2 × S 2 as seen from Eq. (6). In this case the homotopy group is richer, Z×Z. Two integer topological invariants can be introduced counting the degree of covering of the two spheres by the mapping Q from real space to M F . This allows us to write the topological term explicitly although in a non-invariant form. With the parameterization (6), it reads (cf. Ref. [12]): Let us now analyze the minima of the action (8) for the setup Fig. 1. Circular symmetry of the problem fixes the phase equal to the polar angle, φ f,b = ϕ. The other parameters depend only on the radial coordinate r. Both angles θ in F and B sectors obey the Usadel equation Inside the hole ∆ = 0. At low energies we can also neglect the E term and the equation becomes independent of k.
The step-like dependence of the order parameter, Eq. (2), imposes the boundary condition θ(R) = π/2. There are two possible solutions to the Usadel equation with this boundary condition: Saddle point equations also require k f = 0 and hence the angle χ f drops from the matrix Q and from the action. The two remaining angles k b and χ b are free and can take any constant values. The spatial profile of Q is thus fixed by the Usadel equation. The solutions θ 1,2 represent two disconnected saddle points in the F sector while in the B sector only the saddle point θ b = θ 1 is reachable. If the integration contour for θ b , which runs along the imaginary axis, is shifted to the point θ 2 a divergence occurs in the k b integral. Thus the B sector is reduced to a hyperboloid parameterized by ik b > 0 and 0 < χ b < 2π, see Ref. [13]. This is exactly the structure of the sigma model of class BD as we anticipated from symmetry analysis. The distinction between even (D) and odd (B) versions of this class is related to the disjoint character of the manifold due to the discrete degree of freedom in the F sector. Namely, in class D (B) the two parts of the manifold contribute to the partition function with the same (opposite) sign. In our problem the odd symmetry class B occurs. To demonstrate it, we compare the value of the action (8) at the two minima in the F sector. These two solutions indeed contribute with the opposite sign since the corresponding values of the topological term (9) differ by exactly iπ.
The density of states is given by the integral At low energies, this integral is to be calculated over the saddle manifold. Apart from the two variables k b and χ b and two disconnected points in the F sector this manifold involves two Grassmann variables in the matrix U , Eq. (5). In order not to spoil the saddle point, these variables must be constant in space and fulfill the condition [U,∆] = 0. We introduce them according to Within this parameterization, we can rewrite the integral (13) in terms of k b , χ b , η, ζ. Thus the problem is reduced to the 0D sigma model of class B (the two disjoint points in the F sector contribute with opposite signs). Explicit calculation of the integral [13,15] yields the local DOS in the factorized form Here γ = G t n(r 0 )/2π ≪ 1 and the low-energy level spacing is given by The spatial profile of DOS, n(r), is fixed by the solution of the Usadel equation, while the energy dependence is characteristic for the B class and contains a narrow lorentzian peak at zero energy. This peak is the zero-energy MBS broadened due to the finite tunneling time. Note that the width γ is position-dependent.
Integrating ρ(r, E) over space yields the global DoS In the extreme tunneling limit γ → 0, this function acquires the contribution δ(E)/2 and coincides with the result [15] up to a factor 2 due to BdG double counting.
The Majorana state appears as a half of a fermionic level. Spatial and energy dependence of ρ(r, E) factorize at energies much less than the Thouless energy E Th = D/R 2 . At higher energies fluctuations of Q are not important and DOS is given just by a single saddle point. Approximate solution of the Usadel equation (10) in the limit E ≫ E Th yields [13] N (E ≫ E Th ) = πνR 2 Local DOS is close to the normal value ν everywhere except for a narrow vicinity of the hole boundary.
The local density is measured in the tunneling experiment. The tunneling current is determined by with f (E) being the equilibrium Fermi distribution function. At low temperatures and voltages, (T, eV ) ≪ ω 0 , we keep only the first lorentzian term in the energy dependence of the local DOS (17). The differential conductance exhibits a peak At very low temperatures the height of this peak is universal and equals e 2 /π . This signifies resonant Andreev reflection at the Majorana state. This effect in the absence of disorder was studied in Ref. [10]. At larger (and more realistic) temperatures, γω 0 ≪ T ≤ ω 0 , the height of the peak is parametrically small, dI/dV ∼ γω 0 /T . Noise power of the tunneling current in the same The noise produced by the Majorana level is T -and Vindependent as long as γω 0 ≪ T . The results (22) and (23) apply at low temperature and voltage. When temperature and/or voltage are higher than the level spacing ω 0 , non-zero-energy states also contribute to the tunneling current. Positions and widths of these states depend on the realization of disorder. For low temperature but high voltage, T < ω 0 < eV , narrow resonances similar to Eq. (22) will occur due to non-zero levels [16]. Positions of these resonances strongly depend on disorder realization; their heights are smaller but close to e 2 /π and widths are of order γω 0 . When temperature exceeds ω 0 , all the resonances get smeared and the normal average DOS ν is recovered.
To conclude, we have studied the local density of states in the superconducting vortex on the surface of a topological insulator in the superconducting antidot setup depicted in Fig. 1. The spatial and energy dependence of the density of states factorize at low energies and the latter is given by the 0D sigma model of symmetry class B. We have identified the zero-energy Majorana state occurring due to the topological properties of the system. This Majorana state exhibits itself via the resonant Andreev reflection at zero energy yielding the peak in differential conductance with the universal amplitude e 2 /π and width proportional to the normal conductance G t .
We are grateful to B. Sacepe for stimulating discussions. This work was supported by the RFBR grant 10-02-00554-a and by the German Ministry of Education and Research (BMBF).

Supplementary Information
Derivation of the sigma-model

Sigma model for hybrid structure
In this section we derive the sigma model describing the dynamics of the topological insulator surface excitations in the presence of a random potential. We assume the hole geometry depicted in Fig. 1 of the main text. The derivation starts with the microscopic Bogolyubov -de Gennes Hamiltonian (1). The sigma model is aimed at averaging the single retarded Green function determining the density of states. We will also demonstrate the equivalence of this sigma model to the model of the symplectic symmetry class AII describing the transport properties of massless Dirac fermions. The latter model is based on the normal, rather than superconducting, Hamiltonian and is capable of calculating averaged products of retarded and advanced Green functions. It is known that the sigma model for Dirac fermions possesses the specific Z 2 topological term. The equivalence of the two models will prove the appearance of the topological term in the superconducting case considered in the present paper.
Let us start with the retarded Green function in the system governed by the Hamiltonian (1). We will use the supersymmetric integral representation Here the vector fields Φ and Φ † contain 4 commuting and 4 Grassmann parameters each. Apart from the spin and Bogolyubov -de Gennes structure of the Hamiltonian, we also introduce the superstructure in order to get rid of the normalizing denominator and facilitate further disorder averaging. The pre-exponential factor contains the supermatrix k = {1, −1} FB and supertrace is defined as in Ref. [17]: str A = A FF − A BB . The superconducting symmetry of the Hamiltonian H, Eq. (3), gives rise to specific soft modes -Cooperons -of the type G R E G R −E . These modes are relevant for the average density of states since they are built out of retarded functions only. In order to include them into our effective theory, we transform the action by writing half of it in the form (24) and another half in the time-reversed (transposed) form: Here we have introduced the matrix σ z operating in the space of time-reversed blocks -TR space. The matrix k appeared due to anticommutation of Grassmann variables. We will denote the doubled vectors as They are no longer independent, as Φ and Φ * were, but rather obey the linear constraint At this stage, we are ready to perform the disorder averaging. We adopt the standard Gaussian white noise disorder model characterized by the correlator The single-particle density of states is defined as ν = |µ|/2πv 2 0 . Note that the right-hand side is twice bigger compared to the standard definition applied to normal metals. The reason for this extra factor 2 is the Dirac nature of electron spectrum and hence doubling of the phase space due to the two types of excitations (electrons and holes).
Disorder averaging produces the quartic term in the action, Here we use the notation Λ = σ z τ z . The interaction term is further decoupled with the help of the Hubbard-Stratonovich transformation by introducing an auxiliary matrix field Q. The construction of the nonlinear sigma model implies that the field Q contains all relevant slow modes of the disordered system -diffusons and Cooperons. There are in total three ways to decouple the four-fermion interaction [17]. One of them involves a scalar field ∼ Ψ Ψ analogous to the random potential V . This leads only to an irrelevant renormalization of the chemical potential. Two other decoupling schemes introduce a matrix field with the structure ΨΨ and ΨΨ T corresponding to diffusons and Cooperons, respectively. In order to include all relevant slow modes in the theory, we have to perform decoupling in both ways, which can be achieved with a single matrix Q. Details of this calculation can be found in Ref. [17]. The resulting action has the form: The vectors Ψ andΨ are related by Eq. (27). This allows us to limit the matrix Q by the linear constraint Q =Q ≡ CQ T C T . This relation keeps only those parameters in Q that couple to the product ΨΨ. Integrating out the field Ψ yields the action for the matrix Q: Derivation of the sigma model proceeds with the saddle-point analysis of the above action. In the dirty limit, the energy and∆ terms in the argument of the logarithm are relatively small compared to Q/τ . With these small terms neglected, the action possesses the uniform saddle point Q = Λ. This saddle point corresponds to the self energy of the average electron Green function in the self-consistent Born approximation. Other saddle points can be achieved by rotations Q = T −1 ΛT , if the matrix T commutes with the spin operator s (E and ∆ terms are still neglected). Rotations T define the target manifold of the non-linear sigma model. The effective action of the sigma model within this manifold is derived with the help of the gradient expansion of Eq. (31) allowing for slow spatial variation of Q and perturbative expansion to the linear order in E and∆. Since the Q matrix is trivial in the spin space, we can safely reduce its size to 8 × 8 keeping only Nambu, TR, and FB structure. Then the self-conjugacy relation acquires the form of Eq. (4).
The gradient expansion of Eq. (31) is a highly nontrivial procedure in view of the chiral anomaly of the Dirac fermions. The momentum integrals arising after the expansion of the logarithm are divergent in the ultraviolet limit and require a proper regularization. The result of the expansion is independent of a particular regularization scheme provided the gauge invariance is preserved. The anomaly affects only the imaginary part of the action and leads to the appearance of the topological term [12]. At the same time, the real part can be obtained in a straightforward way since all the arising momentum integrals are convergent. The result reads Here the diffusion coefficient is D = v 2 0 τ /2 and the matrix Q is reduced to the 8 × 8 size. The action (8), used in the main text, differs from Eq. (32) only by a π/2 rotation of the superconducting phase and by an imaginary contribution to the energy term. The latter corresponds to a finite dwell time of the electron due to the coupling to the tunneling microscope tip as we elaborate below.
In order to explain the emergence of the topological term in the imaginary part of the action, we will prove the equivalence of the sigma model, derived here for the hybrid structure of Fig. 1, to the symplectic sigma model for Dirac fermions derived in Ref. [12].

Sigma model for Dirac fermions
The sigma model obtained by the gradient expansion of Eq. (31) belongs to the symplectic symmetry class AII. This is quite natural since the Hamiltonian possesses only time-reversal symmetry in the absence of the∆ term. Let us demonstrate the equivalence between our sigma model (8) describing the density of states for the Hamiltonian H, Eq. (1), and the standard symplectic class sigma model describing transport via the surface states of the topological insulator with Hamiltonian H 0 . In the latter case, the sigma model yields averaged products of retarded and advanced Green functions with the energy difference ω. This requires introducing retarded-advanced (RA) structure of the action, Here the matrix τ z operates in the RA space. Subscript 0 is used to distinguish this model from the sigma model derived in the previous section and used in the main text. The time-reversal symmetry H 0 = s y H T 0 s y is taken into account by further doubling of variables in the TR space, Note that the matrix Λ 0 = τ z multiplying the frequency term is different from its counterpart Λ from the previous section. The vectors Ψ andΨ as well as the charge conjugation matrix C 0 are also defined differently, Already at this stage we can prove the equivalence of the two models by selecting a proper unitary rotation of the field vector Ψ. This unitary rotation should obey the following properties: These relations bring the charge conjugation matrix C 0 and the matrix Λ 0 to the representation used in the previous section for the sigma model of the hybrid system. There are many matrices U that fulfill identities (36). One possible choice is In the case of usual rather than Dirac fermions, an analogous equivalence between the sigma model for the density of states in a superconducting hybrid structure and the orthogonal class sigma model is discussed in Ref. [14].

Topological term
So far we have discussed the derivation of the real part of the sigma-model action (8) and also proved the equivalence of this model to the symplectic class sigma model for Dirac fermions (up to the boundary conditions involving the term∆). The latter model is known to possess the Z 2 topological term as an imaginary part of its action, see Ref. [12]. This allows us to include the topological term in the action (8) and thus complete its derivation.
The explicit form of the topological term is written in the main text, Eq. (9), in a noninvariant form, using the parameterization (5) - (7). In this section we will discuss an indirect but explicitly invariant representation of this topological term.
The target manifold of the symplectic class sigma model is fixed by the representation Q = T −1 ΛT with the constraint Q =Q. If we neglect the Grassmann parameters and consider the central part of Q, see Eq. (5), these conditions yield Q F ∈ O(4)/O(2)× O(2) and Q B ∈ Sp(2, 2)/Sp(2)× Sp (2). The topological term arises in the compact F sector of the model. An explicit expression for the topological term can be written within a construction similar to the Wess-Zumino-Witten term, cf. Ref. [18].
We first extend the target manifold by relaxing the condition Q 2 = 1. Let us introduce the matrix Q =T ΛT . The only restriction on this matrix is Q =Q. The F sector of the extended manifold is Q F ∈ O(4). It has trivial second homotopy group, π 2 [O(4)] = 0. This allows us to introduce the third, auxiliary, coordinate 0 ≤ t ≤ 1, such that the real 2D space corresponds to t = 1, and continuously extend the matrix field according to Such an extension assumes that the physical space has no boundary and can be viewed as a surface of a 3D ball with t being the radial coordinate. With this definition of Q, the Wess-Zumino-Witten term has the form Here the indices a, b, and c take three values corresponding to the three coordinates t, x, and y, and ǫ abc is the unit antisymmetric tensor. The integrand in this expression (39) is a total derivative therefore the result of the integration depends only on the value of Q at the boundary of integration domain, i.e., on the physical field Q = Q| t=1 . This can be checked by an explicit calculation of the variation of S θ : (40) A topological term does not change under continuous variations of the field Q but takes different values in different (disconnected) topological sectors of the model. The general Wess-Zumino-Witten term has a nonzero variation and hence does not obey this property. However, the term (39) is constrained by Q 2 = 1. Under this condition, the variation (40) is identically zero. Hence the Wess-Zumino-Witten term (39) indeed plays the role of a topological term.
The value of the topological term is quantized, yielding the topological charge of the field configuration, only in a system without boundary. This is not the case for the geometry considered in the paper. The matrix Q is defined in the finite region inside the hole (see Fig. 1 of the main text) with the boundary conditions fixed by the∆ term in the Hamiltonian (1). The value of the topological term is not completely fixed in this case. In particular, the expression (39) depends not only on physical values of Q but also on the way it is extended in the third dimension, Eq. (38). However, this uncertainty leads only to an uncontrolled imaginary constant in the action. This constant is the same in both topological sectors of the model and hence does not alter any observable quantities.
An alternative derivation of the sigma-model action, including the topological term in the form Eq. (39) is possible within the non-abelian bosonization formalism. Let us for a moment assume, that the Hamiltonian has an extra chiral symmetry, H = −s z Hs z . This situation corresponds to the symmetry class DIII and can be realized, e.g., at the Dirac point in the presence of a random velocity disorder. The sigma model of the class DIII has an extended target space, corresponding to the manifold of Q introduced above. This sigma model possesses the Wess-Zumino-Witten term (39), see Ref. [19], and is the result of the non-abelian bosonization of the initial fermionic problem. Non-zero chemical potential drives the system away from the Dirac point and breaks the chiral symmetry. This reduces the model to the symplectic symmetry class and restricts the field by the condition Q 2 = 1. The resulting action has the topological term in the form of the restricted Wess-Zumino-Witten term (39) as discussed above.

Tunneling coupling term
In this section we discuss the appearance of the imaginary contribution to the energy due to the presence of the tunneling tip. In order to include the coupling to an external metallic probe, we extend the Hamiltonian as Diagonal blocks of this matrix are the Hamiltonian H of the topological insulator surface (1) and the Hamiltonian H M describing electron states in the metallic tip. Off-diagonal elements t and t † are tunneling amplitudes between the two subsystems. We assume these tunneling amplitudes are local in space and couple the states at point r 0 on the sample surface to the states at the tip of the metallic probe.
Deriving the sigma model as described in previous sections we arrive at the expression similar to Eq. (31) but with an extended matrix structure: The quantities ν M , τ M , and Q M refer to the metallic tip. In the absence of tunneling, t = 0, the two subsystems decouple and we obtain independent sigma models describing topological insulator and the tip. We will assume the density of states ν M to be large, that will fix Q M = Λ. Tunneling coupling is assumed to be weak. This allows us to expand the logarithm in powers of t and t † . The lowest non-vanishing contribution appears in the second order of the perturbation theory and yields the tunneling term in the sigma model, In the last expression we have used the saddle-point relation between Q and the Green function (EΛ−τ z H +iQ/2τ ) −1 in coincident points. We introduce normal dimensionless (in units e 2 /h) tunnel conductance of the junction G t = π 2 νν M tr(tt † ) (note that ν and ν M are total rather than spinless densities of states at the Fermi energy). Assuming Q M = Λ, we rewrite the action as Thus the tunneling term has the same structure as the energy term and produces an imaginary contribution to E, see Eq. (8) of the main text.

Calculation of the density of states
Low energies E ≪ E Th The local density of states on the surface of the topological insulator is given by Eq. (13) of the main text. The integral over matrix Q is to be calculated over the saddle manifold fixed by the solution of the Usadel equation (10). There are two distinct solutions θ 1,2 , [see Eqs. (11), (12) of the main text] for the angle θ F . In the F sector of the model, all other parameters are fixed and thus θ 1,2 represent two disjoint parts of the integration manifold. At the same time, only the solution θ b = θ 1 is allowed in the B sector, and two other angles, k B and χ B , are free forming a hyperboloid (the parameter k B must be imaginary in order to ensure convergence of the integral).
To calculate the density of states from Eq. (13), we have to express Str(kΛQ(r)), S[Q] and DQ in terms of k b , χ b , η, and ζ for the two disjoint parts of the manifold. However, a direct use of our parameterization leads to an uncertainty in the integral over the part corresponding to θ f = θ 2 . The integral over dη dζ is zero, while the integral over k b diverges. We resolve this uncertainty by a trick proposed in Ref. [15]. The two disjoint parts of the integration manifold can be mapped onto each other by the similarity transformation We use the parametrization (6), (7) in the domain with θ f = θ 1 and then apply the transformation (45) to cover the second topological sector with θ f = θ 2 . Owing to the relation (45), integration measure is the same in both sectors. It is given by the superdeterminant of the metric tensor in the space of Q matrices, see Ref. [17]. The integration measure reads with κ = ik b (k b is imaginary to ensure the convergence of the integral). The integration runs over the hyperboloid κ > 0, 0 < χ b < 2π. The action in the two topological sectors θ f = θ 1,2 has the form Here we have introduced the complex dimensionless energy parameterx = d 2 r ǫ cos θ 1 = x + iγ, with the real part x = E/ω 0 [see Eq. (18) in the main text] and the imaginary part γ = G t n(r 0 )/2π. The latter appears due to the coupling to the tunneling tip. The terms ±iπ/2 in Eqs. (47), (48) appear due to the topological term S θ [Q], see Eq. (9). With the above expressions for DQ and the action we can calculate the total partition function of the theory, This is exactly what one expects for the partition function of a supersymmetric theory. Note that the topologically nontrivial sector does not contribute to this result since the action S 2 contains no Grassmann variables. In order to calculate the density of states, we need the pre-exponential factor from Eq. (13). In our parameterization, it has the following form in the two sectors: (50) In the last expression we have used the condition γ ≪ 1. The result (51) coincides with Eqs. (15) -(17) of the main text. In the tunneling limit γ → 0 the lorentzian term yields a delta function. This implies that the Majorana state is robust in a closed system and disorder does not smear it. Note, that for an alternative model of the symmetry class D, when the topological term is absent, the density of states does not contain the Majorana delta-peak. Two topological sectors contribute with the same sign yielding [15] ρ D = ν cos θ 1 1 + sin(2πx) 2πx . (52) High energies E Th ≪ E ≪ ∆ In the previous section, we have calculated the density of states at low energies. The result (51) exhibits a deltapeak at zero energy and oscillations which decay at the characteristic energy scale ω 0 . These oscillations appear due to repulsion between low-lying levels with energies E and −E. The scale ω 0 is the global mean level spacing inside the hole. At higher energies E ≫ ω 0 , the effect of level repulsion can be neglected and the density of states is given by the mean-field expression ρ(r, E) = ν Re cos θ(r, E), where θ is a solution of the Usadel equation (10) with k = 0. The result (53) can be obtained from the sigma model identity (13) by integrating over Q around the unique global minimum of the action given by the Usadel equation. Supersymmetry ensures that the integration with respect to small fluctuations around the minimum yields unity and hence Eq. (53) follows. The Usadel equation is a complicated non-linear equation that cannot be solved analytically for arbitrary E. At relatively low energies ω 0 ≪ E ≪ E Th , the solution can be found perturbatively in E. While linear correction in E to θ is purely imaginary and does not change the density of states, the second order calculation up to terms ∼ E 2 is  (16) of the main text is reproduced, while at energies well above E Th the density of states is depleted only close to the hole edge within the strip ∼ R E T h /E in agreement with Eq. (55). Remarkably, the density of states is first increased above its normal value before dropping to zero at r = R.
necessary to obtain the observable result. The spatial profile of this small energy correction does not factorize as in the main term, see Eq. (15) in the main text. Second order perturbation theory amounts to solving linear differential equations, obtained by linearizing Eq. (10), and leads to very complicated functions. We will not discuss any further details of this approach here.
As the energy is increased further and exceeds the Thouless energy, E Th ≪ E ≪ ∆, the Usadel equation admits another approximate solution. The angle θ is small everywhere except a narrow vicinity of the hole boundary. In this vicinity we can neglect the curvature of the superconductor edge and write an approximate one-dimensional Usadel equation D ∂ 2 θ ∂r 2 + 2iE sin θ = 0.
With the boundary conditions θ(R) = π/2 and θ(r ≪ R) → 0, solution to this equation reads The density of states is given by Eq. (53) with the above result for θ. We see that θ decreases from π/2 down to exponentially small values in a narrow strip of the width ∼ R E Th /E near the hole boundary. The density of states takes its normal value ν everywhere inside the hole except for this narrow strip where it is depleted down to zero. Spatial integration yields global density of states This is the result (20) of the main text. We have also solved the Usadel equation numerically for the whole range of energies below ∆. The spatial profile of the density of states is shown in Fig. 3 and the energy dependence of global DOS is depicted in Fig. 2 in the main text. It perfectly matches both the low-energy limit (up to delta peak and oscillatory term) of Eq. (15) and high-energy asymptotics (56).

Tunneling current and noise
Tunneling current and noise can be found using the Landauer formalism instead of the sigma model. Electrons incident from the metallic probe are either reflected normally or experience Andreev reflection. Amplitudes of these processes can be found from the microscopic Hamiltonian of the system and used to calculate the current and higher