Size Constraints on Majorana Beamsplitter Interferometer: Majorana Coupling and Surface-Bulk Scattering
SSize Constraints on Majorana Beamsplitter Interferometer:Majorana Coupling and Surface-Bulk Scattering
Henrik Schou Røising ∗ and Steven H. Simon Rudolf Peierls Center for Theoretical Physics, Oxford OX1 3NP, United Kingdom (Dated: March 6, 2018)Topological insulator surfaces in proximity to superconductors have been proposed as a way toproduce Majorana fermions in condensed matter physics. One of the simplest proposed experimentswith such a system is Majorana interferometry. Here, we consider two possibly conflicting constraintson the size of such an interferometer. Coupling of a Majorana mode from the edge (the arms) of theinterferometer to vortices in the center of the device sets a lower bound on the size of the device.On the other hand, scattering to the usually imperfectly insulating bulk sets an upper bound. Fromestimates of experimental parameters, we find that typical samples may have no size window inwhich the Majorana interferometer can operate, implying that a new generation of more highlyinsulating samples must be explored.
PACS numbers: 71.10.Pm, 74.45.+c, 73.23.-b
I. INTRODUCTION
There has been an ongoing search for Majoranafermions in condensed matter systems which has been in-tensified over several years. Vortices in a spinless p -wavesuperconductor have long been known to bind zero en-ergy Majorana modes. With p -wave superconductorsbeing very rare in nature, no experiment has convincinglyobserved such Majoranas yet. More recently it was pre-dicted that Majorana bound states also exist in vorticesin a proximity-induced superconductor on the surface ofa topological insulator (TI). A number of recent exper-iments on TIs in proximity to superconductors andother similar experimental systems have increasedthe interest in this possibility.
FIG. 1. The Majorana beamsplitter interferometer proposedby Fu and Kane, and Akhmerov et al.
It consists of a 3Dstrong TI in proximity to a superconductor and magnets ofopposite polarization.
The surfaces of TIs support gapless excitations withthe dispersion of a Dirac cone (in principle, TIs canhave any odd number of Dirac cones in the Brillouinzone but we consider the simplest case of a single conefor simplicity). The spectrum can be gapped either by applying a magnetic field to give the Dirac fermion ei-ther a positive or negative mass, or by placing a super-conductor in proximity to the surface. At interfaces be-tween different gapped regions, gapless one-dimensionalfermionic channels can develop. For example, an inter-face between two magnetically gapped regions with op-posite mass signs will contain a gapless and chiral one-dimensional Dirac fermion mode. An interface betweena magnetically gapped region and a superconducting re-gion will contain a gapless chiral Majorana mode. It is the physics of these modes that we are exploring inthe current paper.A very elegant experiment, building an interferome-ter out of these gapless chiral modes, was proposed byFu and Kane and simultaneously by Akhmerov et al. The device is depicted schematically in Fig. 1. Incom-ing particles or holes, biased at a low voltage, flow into aDirac channel between two oppositely polarized magneticregions. The Dirac fermion is split into two Majoranafermions upon hitting the superconductor, one flowing ineach direction around the superconducting region (drawnas a disk in Fig. 1). At the other end of the supercon-ducting region the two Majorana modes are re-combinedinto a Dirac mode. The differential conductance of thisdevice was predicted to take the values 0 or 2 e /h de-pending on whether the number of Φ = h/ (2 e ) vorticesin the superconductor is even or odd, respectively. Withthe exception of quantum Hall systems, this was the firstproposed realistic Majorana interferometry experiment,and it remains a good candidate for the first experimentto successfully establish the existence of chiral Majoranamodes (although, we note that promising evidence of chi-ral Majorana modes was reported very recently ).The most experimentally explored TI materials arethe Bismuth based compounds, including (amongmany others) Bi Se , Bi Te , and Bi Te Se. Within thisclass of materials, several experiments have successfullyformed some sort of superconducting interfaces.
Inthis paper we mainly have this type of material in mind.However, we note that the material SmB , which is a r X i v : . [ c ond - m a t . s t r- e l ] M a r possibly a topological Kondo-insulator, will be discussedin the conclusion.In this paper we study two effects that restrict the sizeof the interferometer device, as summarised in Fig. 2. Onthe one hand, we consider coupling of the Majorana edgemode to the Majorana modes trapped in the core of thevortex in the center of the interferometer. When thiscoupling is sufficiently strong, i.e. when the interferom-eter is sufficiently small, the conductance signal will bedistorted, rendering the interpretation of the experimentdifficult. The coupling is generally strong on the scale ofthe coherence length ξ = (cid:126) v F / ∆ where v F is the Fermivelocity and ∆ is the proximity gap. This length scalecan be on the order of a micron (we discuss materialsparameters in section III A, IV B, and IV C).Next, we consider surface-bulk scattering due to theunintentionally doped and poorly insulating bulk ofTIs. We find that when current leaks from the Ma-jorana edge channel to the ground, the signal obtained atthe end of the interferometer drops exponentially with alength scale set by the surface-to-bulk scattering lengthfrom disorder. For typical samples this length scale canbe shorter than a micron. Thus we have two potentiallyconflicting constraints on the size of the proposed inter-ferometer. We will discuss the possible directions forwardin the conclusion.The outline of this paper is as follows. In section II wereview the formalism of Majorana interferometry and thecharacteristic conductance signal of the experiment. Insection III we consider the impact of coupling betweenthe chiral Majorana and a bound Majorana mode in avortex core. We study how this limits the source-drainvoltage and sets a lower bound on the size of the system.In section IV we consider the leakage of current from thedevice to the TI substrate and we establish an upper limiton the system size due to this leakage. We conclude withnumerical estimates and an outlook for Majorana inter-ferometry on the surface of TIs. Appendix A containsa derivation of the average currents. In Appendix B weprovide generalizations of the single-point Majorana cou-pling displayed in Fig. 2 (a). Vortex-bound Majoranafermions in a TI/SC hybrid structure and the energysplitting between zero modes is discussed in AppendixC. In Appendix D we present details of our estimate ofthe surface-bulk scattering rate based on Fermi’s Goldenrule. Finally, in Appendix E signal loss due to scatteringwith acoustic phonons is briefly discussed.
II. BACKGROUND ON INTERFEROMETRYWITH MAJORANA FERMIONS
In this section, we first review the formalismneeded to calculate the conductance and interferometrycurrent.
In Fig. 2 (a) we show a top view ofthe interferometer that was described in the introduction.For the moment we ignore the Majorana coupled to thetop arm (marked as an “X”) at position 0. Charge trans-
S Dλ ↑ M ↓ ξ L ξ R ξ (a) S Dξ L ξ R ξ η η (b) FIG. 2. Top view of the interferometry device with a super-conductor (SC) and magnetic domains (with magnetizationM ↑ / ↓ ) that probes charge transport from the source (S) tothe drain (D) through a pair of chiral Majorana modes. (a)Single-point coupling of strength λ at coordinate x = 0. Here, ξ R ≡ ξ ( x = 0 + ) and ξ L ≡ ξ ( x = 0 − ). (b) The interfer-ometer coupled to a single-mode conducting lead representingleakage to the bulk of the TI. port from the source to the drain is computed by usinga transfer matrix which describes transport of particlesand holes from the source on the left ( L ), via the perime-ter of the superconductor, to the drain on the right ( R ),[ ψ e , ψ h ] TR = T [ ψ e , ψ h ] TL where T here means transpose.The matrix T can be decomposed into three pieces cor-responding to the three key steps between the source andthe drain T = S † P S = (cid:18) T ee T eh T he T hh (cid:19) . (1)The unitary matrix S relates the Majorana states [ ξ , ξ ]running along the upper and lower edge of the supercon-ducting disk to the electron and hole states [ ψ e , ψ h ] thatenter via the leads, [ ξ , ξ ] T = S [ ψ e , ψ h ] T . The matrix P contains plane wave phases that the low energy chiralMajorana modes accumulate as they move along the edgeof the superconducting disk. Finally, the matrix S † re-assembles the two Majoranas into outgoing electron andhole states that enter the drain.The matrices S and T are functions of energy. Dueto the particle-hole symmetry, we must have S ( E ) = S ∗ ( − E ) τ x , where τ x is a Pauli matrix in particle-holespace. At E = 0 these constraints fix S (0) up to anoverall phase that observables do not depend on, S (0) = 1 √ (cid:18) i − i (cid:19) (cid:18) e iα e − iα (cid:19) . (2)We will apply S (0) with α = 0 for convenience. As dis-cussed in Ref. 17 this form is exact when the systemhas a left-right symmetry. Even in cases where the sys-tem breaks this symmetry, corrections are O ( E ) andcan thus be ignored at low temperature and low voltage(see Appendix C 1). We study the symmetric situationwhere the magnetization is M ≡ M ↑ = − M ↓ through-out the main text. The magnetization enters the modelHamiltonian given in Eq. (C1).The transfer matrix is used to calculate the averagecurrent in the drain as the difference between the electronand hole current, with the source biased at voltage V (seeAppendix A for a derivation), I D = eh (cid:90) ∞ d E δf ( E ) (cid:0) | T ee | − | T eh | (cid:1) . (3)Here, δf = f e − f h with f e/h ( E ) = f ( E ∓ eV ), and f ( E ) = (1 + e βE ) − is the Fermi-Dirac distribution with β = 1 / ( k B T ) and E measured relative to the Dirac point.The incoming current is I S = eh (cid:90) ∞ d E δf ( E ) . (4)By current conservation, a net current of I SC = I S − I D is absorbed by the grounded superconductor. As followsfrom Eq. (3), (4), and unitarity of T the differential con-ductance measured in the grounded superconductor atzero temperature is: G SC ( V ) = d I SC d V (cid:12)(cid:12)(cid:12) T =0 = 2 e h | T eh ( eV ) | . (5)In order to calculate G SC ( V ) we need only establish theproperties of the propagation matrix P . If the arms ofthe interferometer are of length l and l we have P ( E ) = (cid:18) e ik ( E ) l +2 iφ ( E ) e ik ( E ) l (cid:19) . (6)Above the wavevector is of the form k ( E ) = E/v m dueto the linear dispersion of the Majorana modes where v m is the Majorana velocity. We have included an ad-ditional phase shift φ which may come from a numberof sources. In Refs. 16 and 17 the possibility of n vor-tices being added in the center of the superconductingregion was considered. In this case the additional phaseis e iφ = ( − n . Inserting P into Eq. (1) and (5) at zerotemperature yields the conductance G SC , ( V ) = 2 e h sin (cid:18) nπ eV δL (cid:126) v m (cid:19) . (7)Here, δL = l − l is the difference between the lengthsof the two arms. At δL = 0, G SC , ( V = 0) = 2 e /h for n odd, and is zero for n even. This would be a rather clearexperimental signature. For the same phase matrix P ,the drain current in Eq. (3) evaluates to I D, ( V ) = ( − n πk B T eh sin ( eV δL (cid:126) v m )sinh ( πδLk B T (cid:126) v m ) , (8)which holds in the low temperature and low voltage limit(compared to the bulk gap). We emphasise that theseresults are derived with zero coupling to the central Ma-jorana and to any other degrees of freedom (e.g. phononsor conducting bulk states), hence the subscript 0. III. EFFECT OF MAJORANA COUPLING
We now consider the effect of coupling the Majoranastrapped in the cores of vortices in the superconductor tothe chiral edge states. Each vortex traps a single Majo-rana mode. If a vortex is close to the edge (roughly withinthe coherence length) there will be tunnelling coupling asshown in in Fig. 2 (a) where the Majorana zero mode ismarked as an “X” and the (tunnelling) coupling matrixelement is of magnitude λ . The magnitude of the cou-pling drops exponentially with the distance between thevortex and the edge, see Appendix C 2 and C 3.We very generally describe the chiral Majorana on theupper interferometer arm, ξ , by a Lagrangian density L = iξ ( ∂ t + v m ∂ x ) ξ where x is the spatial coordinatealong the upper edge. A similar description holds for thestate on the lower edge, ξ . The vortex bound Majorana, ξ , is described by L = iξ ∂ t ξ . We add the couplingterm between the central bound state and the chiral mode L bulk − edge = 2 iλξ ( x = 0) ξ . (9)The equations of motion, following from the full La-grangian L + L + L bulk − edge , are given by ∂ t ξ = λ [ ξ R + ξ L ] , (10) v m ξ R = v m ξ L + λξ . (11)Here, the notation ξ R = ξ ( x = 0 + ) and ξ L = ξ ( x =0 − ) was introduced. A Fourier transformation yields aphase shift across the coupling point, ξ R ( ω ) = ω + iνω − iν ξ L ( ω ) = e iφ ξ L ( ω ) , (12)with ω the frequency, ν ≡ λ / (2 (cid:126) v m ), and φ ( ω ) =arctan( ν/ω ). This energy-dependent phase shift is in-serted into Eq. (6) and we obtain the zero-temperatureresult G SC ( V ) = 2 e h sin (cid:18) nπ eV δL (cid:126) v m + arctan (cid:16) νeV (cid:17)(cid:19) . (13)Observe that the even-odd effect undergoes a crossoverwhen the coupling strength is of the order λ ≈ √ (cid:126) v m eV .At zero voltage, or at infinite coupling strength, the even-odd effect is reversed from the value at high voltage or lowcoupling strength. This crossover is equivalent to shifting n by one in G SC , , causing ξ to acquire a phase shiftof π at low energy, and it is assigned the interpretationthat a vortex Majorana effectively is absorbed by theedge. Similar results are known from Quantum Hallinterferometers at filling fraction 5 / The originalconductance is recovered at high voltage (see Fig. 3).The above result applies when there is no position de-pendence in the coupling to the edge. If we instead con-sider a continuous bulk-edge coupling L bulk − edge = 2 i (cid:90) d x λ ( x ) ξ ( x ) ξ , (14) [ μ eV ] h e G ( V ) Even nOdd n
FIG. 3. The differential conductance from Eq. (13) as a func-tion of the voltage. Here, ν = 1 µ eV and (cid:126) v m /δL = 2 µ eV.The dotted lines show the conductance with ν = 0. At lowvoltage the even-odd effect is reversed. then the total phase shift is again given by Eq. (12) with λ → (cid:90) d x λ ( x ) e ikx (cid:90) x d x (cid:48) λ ( x (cid:48) ) e − ikx (cid:48) (15)in the numerator and a similar replacement in the de-nominator, see Appendix B 1. We note that if x c (cid:28) /k ,where x c is defined such that λ ( | x | > x c ) ≈
0, then the ef-fective λ above becomes (cid:104) λ ( x ) (cid:105) with (cid:104) λ ( x ) (cid:105) = (cid:82) d xλ ( x )up to corrections of order kx c .The scheme above can also be generalized to includemore complicated couplings to multiple (vortex) Majo-rana modes. So long as these modes are coupled to onlya single edge, and not to each other, each coupling causesa phase shift of arctan( ν/ ( eV )) in the phase of propaga-tion along the edge, see Appendix B 2. In the case thatmultiple vortex Majoranas are coupled to each other, aneven number of Majoranas will gap out, whereas an oddnumber will leave a single effective Majorana zero mode.In the above calculation we assumed Majorana cou-pling to one edge only. Since coupling varies exponen-tially with distance to the edge it is not unreasonablethat this will effectively be the case. However, it is alsorealistic that a vortex will be roughly equal distance from,and hence equally coupled to, both edges. This case isdiscussed in detail in Appendix B 3. While the generalresult becomes complicated, at least in the case of equalcouplings to both edges and both edges of equal length,the physics of the even-odd crossover found in this sec-tion remains unchanged. Finally, coupling between a sin-gle vortex Majorana and the edge at finite temperatureis discussed in Appendix B 4. A. Lower Bound on Size and Voltage Constraints
The above derived crossover makes the observation ofthe even-odd conductance effect impossible at low volt-age. The proposed experiment is to add a single vortexand observe a change in conductance (say, from zero to2 e /h ). However, if the Majorana is then effectively ab-sorbed into the edge, the conductance remains zero even once the vortex is added, destroying the predicted effect.This will occur for interferometers of size comparable tothe coherence length, i.e. the length scale of an order-parameter deformation. Assuming that ∆ = 0 . Se , the naive lower size bound for the disk is ξ = (cid:126) v F / ∆ (cid:39) µ m. The tunnelling coupling is expected to decay like λ ∝ µ exp( − R/ξ ) when R (cid:29) ξ (see Eq. (C5)) where µ is the chemical relative to the Dirac point. For disksof size R (cid:39) ξ the energy splitting is comparable to theenergy gap and the notion of stable edge/vortex statesbreaks down.We note, however, that the Majorana coupling termvanishes identically at µ = 0 in the topological insulatorsuperconductor hybrid structure. This is related toappearance of an additional symmetry at the Dirac point,bringing the Hamiltonian (Eq. (C1) in the absence of amagnetic field) from symmetry class D to the BDI in theAltland-Zirnbauer classification. If there is disorder in-ducing local fluctuations in the chemical potential (seesection IV C), say on some scale δµ , a random couplingterm of the type in Eq. (14) will be present and causeenergy splitting. Still, the scenario of having the aver-age (cid:104) λ ( x ) (cid:105) ∼ (cid:104) µ ( x ) e − r ( x ) /ξ (cid:105) ≈
0, which would make theparasitic phase shift vanish, is possible but becomes ex-tremely geometry sensitive (e.g. sensitive to the vortexposition) if R (cid:39) ξ . Moreover, we should expect (cid:104) λ ( x ) (cid:105) to be on the scale of δµ (cid:112) l/d , which in general will notbe small. Here, d is the length scale associated with theenergy fluctuations. Assuming that we cannot controldisorder, the only way to assure suppression of unwantedphase shifts and energy splitting is to increase R . Thus,we will use R = ξ as a strict lower size bound on thedevice.As far as chiral transport on the interferometer armsis concerned, it seems at this stage that one can operateat high voltage, eV (cid:38) ν to avoid the coupling. Natu-rally, the voltage is constrained from above by the globalbulk gap, eV (cid:46) min { M , ∆ } , to avoid excitation of non-topological states. Thus, if R (cid:39) ξ the remaining voltagerange, ν (cid:46) eV (cid:46) min { M , ∆ } , might be too limited tohave a clear experimental signature. Increasing the diskradius to suppress the coupling (and therefore lower ν )induces the problem of signal leakage to conducting bulkstates, which we estimate below.Finally, we note that the upper voltage bound in prac-tice can be lower. This is due to the Caroli-de Gennes-Matricon excited states of the vortex, characterised by aminigap, ∼ ∆ /µ , which can be less than one mK.The excited vortex bound states can be activated ther-mally or by tunnelling from the edge states, in analogy tothe zero mode tunnelling in the beginning of this section.By pinning the vortex to a hole in the superconductor,the minigap can be increased to a substantial fraction of∆ . Although the details of such a tunnelling processis outside the scope of this paper, additional resonancesand conductance phase shifts are expected as eV hits thebound state energies. IV. SURFACE-BULK SCATTERING
Topological qubits are intrinsically protected fromdecoherence. In protocols based on braiding with topo-logical qubits, leakage of current is harmful since it gener-ically causes entanglement with the environment that po-tentially corrupt the qubits. Although the experimentwe study here does not probe a topological qubit, bothtypes of experiment are sensitive to bulk leakage. Sincemost TIs are poorly insulating, bulk leakage is a rele-vant problem to consider.In this section we model leakage of Majorana modesfrom the the surface of the TI to its poorly insulatingbulk. This is done by coupling the interferometer armfirst to a single metallic lead. Then, we take multipleweakly coupled metallic leads to represent a continuouslyleaking environment. We combine the result of this scat-tering process with an estimate of the surface-to-bulkscattering rate from disorder in doped TIs. Our resultssuggest an upper size bound that potentially coincideswith the lower bound discussed in section III for manyunintentionally doped TIs. A. Scattering on Conducting Leads
Let the upper interferometer arm be coupled to ametallic lead as depicted in Fig. 2 (b). Referring to theformalism of Ref. 27, the lead fermions are transformedto a Majorana basis, [ η ( ± )1 , η ( ± )2 ] T = S [ ψ ( ± ) e , ψ ( ± ) h ] T withthe superscript indicating incoming ( − ) or outgoing (+)states. The S matrix here may differ from the one inEq. (2) at low energy only by having a different phase α ,which is irrelevant for observables. The scattering pro-cess in the Majorana basis is denoted by A ( E ), (cid:16) ξ R , η (+)1 , η (+)2 (cid:17) T = A ( E ) (cid:16) ξ L , η ( − )1 , η ( − )2 (cid:17) T . (16)By rotating the lead particles and holes into the appro-priate Majorana basis, the chiral Majorana mode canbe shown to decouple from one of the (artificial) leadMajoranas. In the low energy limit, the scattering ma-trix is a real rotation matrix, A ( E ) ∈ SO(3). If wedenote the local reflection amplitude of η ( η ) by r ( r = 1), this means that the scattering matrix can beparametrized by r at the junction only (the transmissionamplitude is t = (cid:112) − r ), A ( E ) = (cid:18) r − t t r (cid:19) ⊕ r . (17)Including the state ξ on the lower arm and the planewave phases acquired across the interferometer, we ob-tain the matrix P acting on [ ξ L , ξ , η ( − )1 , η ( − )2 ] T , P = r e ikl + inπ − t e i kl e ikl t e i kl + inπ r ⊕ . (18) Here, r = 1 was used. The transfer matrix can be com-puted, and it has the 2 × T = (cid:0) S † ⊕ S † (cid:1) P (cid:0) S ⊕ S (cid:1) = (cid:18) T S → D T (cid:96) in → D T S → (cid:96) out T (cid:96) in → (cid:96) out (cid:19) . (19)Above, the subscripts indicate the result of transportfrom/to the source ( S ), the drain ( D ), or the incom-ing/outgoing lead ( (cid:96) in / out ). The blocks in this transfermatrix are used to find the average current contributionas measured in contact β , with | ( T α → β ) ee | − | ( T α → β ) eh | = r ( − n cos (cid:16) EδLv m (cid:17) ( α, β ) = ( S, D ) r ( α, β ) = ( (cid:96) in , (cid:96) out )0 ( α, β ) = ( S, (cid:96) out ) , ( (cid:96) in , D ) . (20)Here, the first line gives the contribution from combiningMajoranas ξ and ξ , the second line from combining η and η , and the third line from combining ξ and η .Combining Majoranas from different sources always givea vanishing contribution to the average current. The drain current is thus I D ( V ) = r I D, ( V ), where I D, is defined in Eq. (8); the visibility is coherentlyreduced by r . The current in the outgoing lead is I (cid:96) out = r e V (cid:48) /h at T = 0 when the lead is biased at V (cid:48) .Subtracting the incoming current I (cid:96) in = e V (cid:48) /h yields anet current of I (cid:96) out − I (cid:96) in = ( r − e V (cid:48) /h in the conduct-ing lead. When decoupling the lead completely, r = 1,no net current goes in the lead.The matrix P can be trivially extended to include scat-tering on many single-mode leads. Repeating the calcu-lation above with two scattering leads of the same re-flection amplitude r gives the same reduction of currentin both leads, independent of which arms the leads arecoupled to. The drain current is reduced by r . Gener-alising these statements by induction, with N identicalscatterers of reflection amplitude r , the total conduc-tance from the collection of leads (representing the bulkof the TI) is G leads = e N (1 − r ) /h . Furthermore, the N scatterers reduce the visibility of the drain currentmultiplicatively as I D ( V ) = r N I D, ( V ).If we let l S denote the average length a chiral Majoranatravels before it scatters into the bulk, we may by defi-nition express the reflection amplitude as r = 1 − N ll S ,which is equivalent to defining the leakage conductanceinto the collection of leads by G leads = e l/ ( hl S ). Takingthe continuum limit of infinitely many weak scatterersmeans that lim N →∞ r N = exp( − l/l S ), and the current isexponentially suppressed in the drain, I D ( V ) = I D, ( V ) e − l/l S . (21)As for the differential conductance measured in thegrounded superconductor, the amplitude of the oscilla-tions are suppressed G SC ( V ) = G SC , ( V ) e − l/l S + e h (1 − e − l/l S ) , (22)and distinguishing even from odd n is rendered difficultas l surpasses l S . Thus, l S acts as an upper bound on thecircumference of the length of the interferometer arms. B. The Surface-Bulk Scattering Rate
The remaining step of our argument is to estimate thescattering length l S . We start by calculating the surface-bulk (SB) scattering length for a TI surface electron ( l ( e ) S )in the absence of any surface superconductor or magnet.Then, we use a similar calculation to estimate the Ma-jorana scattering length ( l ( m ) S ). We restrict ourselves toelastic (zero temperature) surface-bulk scattering, wherethe electron-phonon coupling is irrelevant (see subsectionIV D and Appendix E). Building on the formalism devel-oped in Ref. 44 we consider scattering via screened chargeimpurities. We note that the consideration in Ref. 44 isrestricted to point scatterers where the overall scatteringstrength is a priori unknown, whereas in our approachthe effective potential strength is fixed by the screenedCoulomb potential and the dopant concentration.Specifically, we consider a surface state ( S ) with initialwave vector k = ( k (cid:107) ,
0) that scatters to a lowered con-duction band (a charge puddle) n (cid:48) in the bulk ( B ) withfinal wave vector k (cid:48) = ( k (cid:48)(cid:107) , k (cid:48) z ). The incoming surfacestate has energy (cid:15) F = v F k (cid:107) . For numerical purposes wework with a TI slab of thickness L = 40 nm. In Fig.4 (a) we show the form of | Ψ( z ) | at zero parallel mo-mentum for the lowest three eigenstates (derived in Ref.45). One of the TI surface states, penetrating tens of˚As into the bulk, is seen in purple. Assuming randomlydistributed screened (dopant) charges with average con-centration n we get the scattering rate from Fermi’sGolden rule (cf. Appendix D):Γ impSB ( (cid:15) F ) = (23) n (cid:32) e π(cid:15) (cid:15) r k (cid:107) (cid:33) (cid:88) n (cid:48) ∈ B k (cid:48)(cid:107) |∇ ξ B,n (cid:48) ,k (cid:48)(cid:107) | (cid:90) ∞ d k (cid:48) z d σ ( n (cid:48) )long ( (cid:15) F )d k (cid:48) z . Here, ξ B,n (cid:48) ,k (cid:48)(cid:107) is the dispersion of bulk band n (cid:48) definingthe outgoing wave vector k (cid:48)(cid:107) by ξ B,n (cid:48) ,k (cid:48)(cid:107) = (cid:15) F . The dif-ferential scattering rate d σ ( n (cid:48) )long ( (cid:15) F ) / d k (cid:48) z is the screenedCoulomb coupling convoluted with the overlap betweenthe bulk and surface states (cf. Eq (D5)), see Fig. 4 (c).In most topological insulators the Fermi level resides inthe bottom of the conduction band or the top of the va-lence band, and a significant doping of acceptors/donorsis needed to reach a bulk insulating state. For Bi Se films (typically being n -doped) we use the dopant con-centration n = 10 cm − . Such doping causesfluctuations in the average (doping) concentration, mak-ing the conduction band bend and induce electron andhole puddles, i.e. effectively filled pockets from theconduction band at the Fermi level. As a simple model ofa typical situation we imagine that chemical doping of the - -
100 100 200 z [ Å ] | ψ ( z ) | [ Å - ] | ψ S ( z ) | | ψ B , n = ( z ) | | ψ B , n = ( z ) | yMajorana profile Δ M (a) - - k ∥ [ Å - ] ξ n ( k ∥ ) [ eV ] (b) k z ' [ Å - ] σ / dk z ' ϵ F = ϵ F = ϵ F = (c) ϵ F [ eV ] Γ imp [ μ eV ] L =
20 nmL =
30 nmL =
40 nm (d)
FIG. 4. (a) The three lowest-lying eigenstates in a TI filmof thickness L = 40 nm. The surface state localized on theTI top is displayed in purple; the bottom surface state isnot shown. Inset: A sketch of the typical transverse surfaceMajorana wavefunction. Its effect on the scattering rate isdiscussed in Appendix D 1. The model parameters are ad-justed to those of Bi Se except for the gap, which is hereset to ∆ TI = 0 . (b) The dispersion of the states displayed in (a) butwith two additional bulk states. In black: The typical Ma-jorana surface dispersion. (c) The differential scattering rated σ ( n (cid:48) =1)long ( (cid:15) F ) / d k (cid:48) z defined in Eq. (D5) shown for three valuesof the Fermi energy (measured in eV). (d) The electronic scat-tering rate as given in Eq. (23) with n = 10 cm − . VanHove singularities are seen whenever the Fermi energy hitsthe bulk bands, but they are softened by the k (cid:48)(cid:107) factor thatgoes to zero at these points. Here, we used (cid:15) r = 100. TI bulk lowers the bulk Fermi level to lie in some region0 .
10 eV (cid:46) (cid:15) F (cid:46) .
18 eV relative to the Dirac point. Thiswould naively be bulk insulating since the Dirac point inBi Se is separated from the conduction band by a gap of0 .
28 eV. Puddles coming from stretching the bottom ofthe conduction band are simply modelled by setting theTI gap to ∆ TI = 0 . ). Let us thereforeassume that the surface Fermi energy is tuned near theDirac point. The electronic scattering rates as a functionof the bulk Fermi energy are shown in Fig. 4 (d). We notethat the rate stays similar in magnitude for a large TI filmthickness range. When the film is made thinner, there arefewer bulk states to scatter into, but this is compensatedby an increased surface-bulk overlap. Moreover, the ratesin Fig. 4 (d) have the same qualitative features as seenfor point source scattering. The (elastic) scattering lifetime of the surface electronsis τ ( e )SB = 1 / Γ impSB , and the scattering length is given by l ( e ) S = τ ( e )SB v F . As a typical value of the rates in Fig. 4(d) we use Γ impSB = 3 µ eV, and we arrive at the unusu-ally long scattering time τ ( e )SB ≈ . l ( e ) S ≈ . − longer than forbulk scattering lifetimes seen in experiment. Yet, ourresults could be consistent with what was attributed tobe unusually long surface lifetimes as observed after op-tically exciting bulk states.
In Ref. 51 one deducedthat τ ( e )BS >
10 ps after seeing a stable population of thesurface state induced by elastic bulk-surface scattering inBi Se (the samples were kept at T = 70 K) after subse-quent inelastic decays associated with much shorter timescales. It would be desirable to see similar experimentsengineered to measure the elastic surface-bulk scatteringrate in doped compounds conducted at lower tempera-tures. C. Velocity Suppression and ContradictingRequirements
We now use the same expression in Eq. (23) to extractinformation about the Majorana lifetime and scatteringrate. In doing so we ignore the (confined) transverse pro-file of the Majorana surface state. This is a valid approxi-mation because the reciprocal coherence length is far lessthan the maximal transverse scattering momentum in thestudied energy regime, see Appendix D 1. Moreover, weassume that the induced superconductivity does not gapthe TI bulk such that elastic surface-bulk scattering isnot precluded. The Majorana scattering rate will typi- cally be larger or equal to the electronic scattering ratedue to a suppression of the Majorana velocity relative tothe Fermi velocity.The Majorana velocity is calculated from the low-energy chiral solution of the Hamiltonian in Eq. (C1): v m /v F = (cid:112) − ( µ/M ) [1 + ( µ/ ∆ ) ] − . Here, µ is thesurface chemical potential relative to the Dirac point.From the k (cid:107) ∝ /v F dependency of Γ impSB in Eq. (23) wededuce that l ( m ) S = τ ( m )SB v m is suppressed approximativelyas l ( m ) S ∼ ( v m /v F ) α l ( e ) S with α ≈ Correspondingly, τ ( m )SB ∼ ( v m /v F ) α − τ ( e )SB .A non-zero size window for the interferometer size is re-quired by imposing l ( m ) S (cid:38) ξ , which in turn means that v m /v F (cid:38) (cid:16) ξ/l ( e ) S (cid:17) / . (24)Inserting the electronic rate from the end of the last sub-section and ξ (cid:39) µ m in this yields v m /v F (cid:38) .
4. Hence,we must require a very fine tuning of µ , | µ | (cid:46) . (cid:39) . M (cid:29) µ .When the TI bulk is increasingly doped with screenedCoulomb charges to make the Fermi energy approach theDirac point, fluctuations in the surface electrical poten-tial energy are enlarged. The spatial dependence of thelocal density of states broadens the Fermi energy intoa Gaussian distribution of finite width. The standarddeviation of this energy smearing has been estimated intheory and measured in experiment to be δµ (cid:39) − . δµ decays inversely proportional to the average chemicalpotential further away). For Fermi energies close to theaverage Dirac point the notion of a uniform local densityof states breaks down. Tuning of the surface chemicalpotential is consequently not globally possible within anenergy resolution set by the distribution width. Impor-tantly, the smearing greatly exceeds the tuning requiredabove, | µ | ∼ δµ (cid:29) ∆ , even if we relax the assumption ofa very small ∆ by increasing it one order of magnitude.Finally, the spatial scale of the chemical potential fluc-tuations is n − / ≈ n − / (cid:28) ξ . Thus, unless these spatialfluctuations in the local density of states can be broughtunder control, unintentionally doped TIs are left unsuitedfor Majorana interferometry. D. Other Limitations
Majorana interactions and coupling between Ma-jorana modes and other degrees of freedom, such asphonons, are other potential sources of decoherence. Lo-cal interactions between chiral Majorana fermions are ex-pected to be heavily suppressed at low temperatures andmomenta with the leading order term going like O ( k ). The chiral Majoranas on the interferometer arms can alsoexcite phonons if an electron-phonon coupling is present;see Appendix E for a short discussion. For spatial inver-sion symmetric materials, e.g. Bi Se and Bi Te , wherethis coupling is dictated by a deformation potential, thedecay rate of quasiparticles due to scattering on acous-tic phonons at T = 0 exhibits a Γ ph ∼ ( eV ) behaviour.This is a posteriori expected to be small compared to bulkleakage at low voltage. Without spatial inversion sym-metry, a piezoelectric interaction can cause the electron-phonon coupling to follow a reciprocal power law in q ,making scattering an increasing problem for small mo-menta. V. CONCLUSIONS
Two limiting effects in Majorana interferometers areconsidered. We include a previously neglected couplingbetween chiral Majoranas and vortex-pinned modes. Atlow voltage this coupling yields a crossover in the con-ductance even-odd effect, distorting the interpretation ofthe experiment. We also find that surface-to-bulk scat-tering in the TI sets an upper bound on the size of thedevice. With a proximity gap of ∆ = 0 . ξ = (cid:126) v F / ∆ (cid:39) µ m for Bi Se . Due to the doping needed to reach a bulkinsulating state in many TIs, conduction band puddleslead to large fluctuations in the surface Fermi energy. Inturn, this is in conflict with the surface chemical potentialfine tuning required to have a non-zero size window forthe interferometer. This leaves the possibility of probingthe experiment in many poorly bulk insulating Bismuthcompounds potentially extremely difficult.The most natural ways to overcome the restrictionsconsidered here are (i) to find superconductors with ex-cellent contact to the TI such that ∆ can be made (ide-ally orders of magnitude) larger, or (ii) to pursue a searchfor TIs with a highly insulating bulk. Most Bismuthbased TIs are poorly insulating and are unintentionallydoped, which results in them being unsuited for Ma-jorana interferometry. However, recently reported mixedBismuth compounds have shown evidence of an appre-ciably insulating bulk, which could potentially open awindow of opportunity for this material. We note thatthe recommendation (i) above is similar to the need forhigh-quality interfaces in superconductor-semiconductornanowire devices. Another possible material system to consider is theputative topological Kondo insulator SmB , which givesstrong sign of surface conduction. The bulk resistiv-ity of this material can reach several Ωcm at tempera-tures below a few Kelvins. Very recently, evidence ofthe superconducting proximity effect in Nb/SmB bi-layers was reported, hence providing a possible plat-form for this experiment. However, the physics of highlyinteracting topological Kondo insulators is still poorlyunderstood, and it is unclear how much of the de-tails of the simple non-interacting TI surface physics will carry through to this more complicated case. ACKNOWLEDGMENTS
We thank Mats Horsdal, Dmitry Kovrizhin, RamilAkzyanov for useful discussions and Kush Saha for feed-back on the numerics. H. S. R. acknowledges Aker Schol-arship. S. H. S. is supported by EPSRC grant numbersEP/I031014/1 and EP/N01930X/1.
Appendix A: Derivation of the Average Current
The derivation of Eq. (3) and (4) goes along the linesof Ref. 59. The scattering states of the interferometerare (chiral) plane waves in the one-dimensional channelconvoluted with a transverse wavefunction (we let k > ψ eL = e ik L x ϕ eL ( y ) , (A1) ψ hL = e − ik L x ϕ hL ( y ) , (A2) ψ eR = T ee e ik R x ϕ eR ( y ) + T eh e − ik R x ϕ hR ( y ) , (A3) ψ hR = T hh e − ik R x ϕ hR ( y ) + T he e ik R x ϕ eR ( y ) . (A4)We have assumed that the contacts α ∈ { L, R } containonly single mode states and that the incident state inthe L contact is either an electron or a hole. One can T Lψ eL ∼ e ikx ψ hL ∼ e − ikx Rψ eR ∼ T ee e ikx + T eh e − ikx ψ hR ∼ T hh e − ikx + T he e ikx FIG. 5. The scattering states at the two contact points. construct arbitrary states by expanding in the scatteringstates above. This is incorporated in the (second quan-tized) field operatorˆΨ α,σ ( r , t ) = 1 √ π (cid:90) d E α √ (cid:126) v α ψ σα ( E α , r )ˆ a α,σ ( E α ) e − iω α t , (A5)where we introduced ω α = E α / (cid:126) = v α k α and the annihi-lation operator ˆ a α,σ ( E α ) of type σ ∈ { e, h } satisfying (cid:8) ˆ a † α,σ ( E ) , ˆ a β,σ (cid:48) ( E (cid:48) ) (cid:9) = δ α,β δ σ,σ (cid:48) δ ( E − E (cid:48) ) . (A6)The current operator of type σ in contact α is defined byˆ I α,σ = (cid:126) e im (cid:90) d r ⊥ ,α (cid:104) ˆΨ † α,σ ∂ x ˆΨ α,σ − ( ∂ x ˆΨ † α,σ ) ˆΨ α,σ (cid:105) . (A7)Here, m is the effective mass, mv α = (cid:126) k α . Assumingfurther that the contacts act as thermal reservoirs kept atequal temperature, we can average the density operatorby (cid:10) ˆ a † α,σ ( E )ˆ a α (cid:48) ,σ (cid:48) ( E (cid:48) ) (cid:11) = δ α,α (cid:48) δ σ,σ (cid:48) δ ( E − E (cid:48) ) f σ ( E ) , (A8)where f e/h ( E ) = [exp( β [ E ∓ µ ])+1] − is the Fermi func-tion for particles and holes. Finally, we assume the trans-verse wavefunctions to be orthonormal, (cid:90) d y ( ϕ σα ( y )) ∗ ϕ σ (cid:48) α (cid:48) ( y ) = δ α,α (cid:48) δ σ,σ (cid:48) . (A9)The source and the drain current are defined by the aver-age particle minus hole current, I D ≡ (cid:10) ˆ I D,e − ˆ I D,h (cid:11) and I S ≡ (cid:10) ˆ I S,e − ˆ I S,h (cid:11) . Calculating these explicitly by usingEq. (A7), (A8), (A9), and unitarity of T leads exactly tothe expressions in Eq. (3) and (4). Appendix B: Charge Transport with MajoranaCoupling Disorder1. One Majorana with Smeared Coupling to theEdge
Consider the case where the edge Majorana is coupledcontinuously to the bound state as in Eq. (14). Usingthe ansatz ξ ( x ) = f ( x ) e i ( kx − ωt ) with f ( −∞ ) being theMajorana (fermion) field on the far left and ω = k/v m forthe linearly dispersing modes. The equations of motioncan be combined to give the continuum version of Eq.(12): − iv m ω∂ x f ( x ) = λ ( x ) e − ikx (cid:90) d x (cid:48) λ ( x (cid:48) ) f ( x (cid:48) ) e ikx (cid:48) . (B1)Here, the plane wave phase contribution that can be ne-glected in the discrete case where λ ( x ) = λδ ( x ) is in-cluded. The equation above can be solved as follows.Integrating both sides gives the solution implicitly by f ( x ) = f ( −∞ ) − ζiωv m (cid:90) x d x (cid:48) λ ( x (cid:48) ) e − ikx (cid:48) , (B2)with ζ = (cid:82) d xλ ( x ) f ( x ) e ikx . Inserting this back into (B1)yields ζ = f ( −∞ ) (cid:82) d x λ ( x ) e ikx iωv m (cid:82) d x λ ( x ) e ikx (cid:82) x d x (cid:48) λ ( x (cid:48) ) e − ikx (cid:48) , (B3)Finally, this expression is used in (B2), and we obtain f (+ ∞ ) = (B4) ω + i (cid:126) v m (cid:82) d x λ ( x ) e ikx (cid:82) x d x (cid:48) λ ( x (cid:48) ) e − ikx (cid:48) ω − i (cid:126) v m (cid:82) d x λ ( x ) e ikx (cid:82) x d x (cid:48) λ ( x (cid:48) ) e − ikx (cid:48) f ( −∞ ) . Comparing this to Eq. (12) proves the statement in Eq.(15).
2. Multiple Majoranas Coupled to Each Edge
In the main text, we obtained a phase contributionof φ = arctan (cid:0) νeV (cid:1) in the differential conductance whenthe edge was coupled to one vortex (Eq. (13)). This canbe generalized if the chiral Majoranas have single-pointcouplings to several vortices. In that case φ is replacedby (cid:80) i φ i − (cid:80) j φ j , where φ i = arctan (cid:0) ν i eV (cid:1) comes fromvortex-edge-coupling with the upper edge and φ j fromcoupling with the lower edge. Hence, multiple phasecrossovers will occur if several vortices are located closeto the edge.
3. One Majorana Coupled to Both Edges
Another generalization is to study point-couplings toboth the lower and the upper arm, in which case the cou-pling term in the Lagrangian is L bulk − edge = 2 iλ ξ ( x =0) ξ + 2 iλ ξ ( x = 0) ξ . The corresponding equations ofmotion are2 ∂ t ξ = λ [ ξ R + ξ L ] + λ [ ξ R + ξ L ] , (B5) v m ξ R = v m ξ L + λ ξ , (B6) v m ξ R = v m ξ L + λ ξ . (B7)Here, we again use the notation ξ R = ξ ( x = 0 + ), ξ R = ξ ( x = 0 + ), ξ L = ξ ( x = 0 − ), and ξ L = ξ ( x = 0 − ).In frequency space we find a relation between ξ and ξ across the coupling points, [ ξ , ξ ] TR = U ( ω )[ ξ , ξ ] TL ,where U ( ω ) is the unitary matrix U ( ω ) = 1 W ( ω ) (cid:32) − ν + ν + iω − √ ν ν − √ ν ν ν − ν + iω (cid:33) . (B8)Above, ν i ≡ λ i / (2 (cid:126) v m ) and W ( ω ) = ν + ν + iω . No-tice how U ( ω ) is off-diagonal in the low-energy limit forthe configuration ν = ν . The chiral Majoranas on theperimeter therefore switch place by tunnelling across thevortex Majorana in this case. The phase matrix is givenby P = (cid:32) e i kl e i kl (cid:33) U ( ω ) (cid:32) e i kl + inπ e i kl (cid:33) , (B9)which leads to the differential conductance G SC ( V ) = 2 e h eV ) + ( ν + ν ) (cid:104) ( ν − ν ) + (cid:2) ( eV ) − ( ν − ν ) (cid:3) sin (cid:16) nπ eV δL (cid:126) v m (cid:17) + eV ( ν − ν ) sin (cid:16) nπ + eV δL (cid:126) v m (cid:17) (B10)+ 2 ν ν (1 + ( − n ) (cid:105) . This reduces to Eq. (13) when ν = 0. For δL = 0the expression above is identical to Eq. (13) with thereplacement ν → ν + ν .0
4. One Majorana Coupled to One Edge at FiniteTemperature
We return to the case with a one-point Majorana cou-pling (Fig. 2 (a)). The drain current can generally bere-expressed as a residue sum over poles z + j in the uppercomplex half-plane. Rewriting Eq. (3), I D ( ν ) = ( − n e i sinh( βeV ) (cid:88) j Res { h ( z ) , z + j } , (B11)where h ( z ) is the function h ( z ) = 1 z + ν ( z − ν ) cos (cid:16) δLzv m (cid:17) − zν sin (cid:16) δLzv m (cid:17) cosh[ β ( z − eV )] cosh[ β ( z + eV )] . (B12)The set of simple poles of h ( z ) in the upper half-plane is z + j ∈ { iν } ∪ { iπ (2 m − /β + eV } ∞ m =1 . For ν = 0, thesum is obtainable in closed form and stated in Eq. (8).For weak coupling at finite temperature, ν (cid:28) k B T, eV ,we expand h ( z ) in powers of βν , with β − = k B T .For convenience we define the reduced current as R ≡ I D ( ν ) /I D, . To second order in βν with δL = 0 we find R = 1 − νπeV tanh ( βeV /
2) (B13)+ βν πeV Im (cid:110) ψ (1) (cid:16) − i βν π (cid:17)(cid:111) + O ( βν ) . Above, ψ (1) ( z ) = d d z log Γ( z ) is the trigamma func-tion. In the infinite coupling limit the sign of the draincurrent is flipped, R → −
1, which can be seen from Eq.(B12). The reduced current decreases monotonically, afeature of setting δL = 0, with coupling strength untilthe vortex Majorana is fully absorbed by the edge. At - FIG. 6. The reduced drain current R as function of βν forsymmetric arms δL = 0. The dotted lines represent the weakcoupling result from Eq. (B13) and the full lines are numericalresults for βeV being 5 (purple), 10 (orange), and 25 (green). T = 0 the current is found to be (with (cid:126) = 1 here):( − n πe I D = v m δL sin (cid:16) δLeVv m (cid:17) + 2 ν (cid:104) cosh (cid:16) δLeVv m (cid:17) − sinh (cid:16) δLeVv m (cid:17)(cid:105) (B14) × (cid:104) Im (cid:110) Ci (cid:16) δLv m ( eV + iν ) (cid:17)(cid:111) + Im (cid:110) Ci (cid:16) − iν δLv m (cid:17)(cid:111) − Re (cid:110) Si (cid:16) δLv m ( eV + iν ) (cid:17)(cid:111)(cid:105) . Here, Ci and Si are trigonometric integral functions. Forsymmetric arms the above result simplifies to give thereduced current R = 1 + 2 νeV [arctan ( νeV ) − π ]. Fromthis, the aforementioned result lim ν →∞ R = − Appendix C: Majorana Fermions in a TI/SC HybridStructure
Superconductivity induced on the surface of a strongTI support chiral Majorana fermions on domain wallsand Majorana bound states localized in vortices. The Fuand Kane Hamiltonian of a TI/SC hybrid structure witha single Dirac-like dispersion is given by H = ( v F σ · p − µ ) τ z + M ( r ) σ z +∆( r ) τ + +∆ ∗ ( r ) τ − , (C1)where τ j and σ j are Pauli matrices acting in particle-holeand spin space, respectively. Moreover, τ ± = ( τ x ± iτ y ) / M is the Zeeman energy associated with a magneticfield applied in the z direction. The standard procedureis to introduce u σ ( r ) and v σ ( r ) that define the quasipar-ticles, conveniently arranged in Ψ = [ u ↑ , u ↓ , v ↑ , − v ↓ ] T satisfying the BdG equations H Ψ = E Ψ .
1. Corrections to the S Matrix
One may solve Eq. (C1) on a magnetic domain wall inthe absence of a superconductor. There are then two chi-ral solutions: one in the particle sector, | τ z = +1 (cid:105) , andone in the hole sector, | τ z = − (cid:105) , both localized at theinterface. On a magnet-superconductor interface thereexists a single chiral Majorana solution with linear dis-persion. By the definition of the S matrix, (cid:32) ξ ξ (cid:33) = (cid:32) S ee S eh S he S hh (cid:33) (cid:32) ψ e ψ h (cid:33) , (C2)one can estimate the scattering elements at the tri-junction as overlaps between the chiral states, S ee = (cid:104) τ z = +1 | ξ (cid:105) , S eh = (cid:104) τ z = − | ξ (cid:105) , etc. Here, ξ ( ξ ) isthe chiral Majorana fermion on the upper (lower) inter-ferometer arm. The Majorana wavefunction decays ex-ponentially with the length scale v F / (cid:112) M − µ ( v F / ∆ )in the magnetic (superconducting) region (inset of Fig. 41(a)). Interestingly, if the magnetization on each side ofthe Dirac channel are of equal magnitude (with oppositepolarization), compactly expressed in terms of the sym-metry H ( − y ) = M − H ( y ) M where M = iσ y , thezero energy result in Eq. (2) is obtained at all energies.Experimentally, this is likely a weak assumption, and wecan therefore safely apply S ( E = 0) throughout. If thetwo magnetic fields in the Dirac region are of differentstrengths there will be corrections to S that are smallwhen E (cid:28) v m v F min { M , ∆ } . One can heuristically adda correction term to S ( E ) proportional to E/ Λ (whereΛ ∼ v m v F min { M , ∆ } is some phenomenological scale)and impose unitarity to order O ( E ) and particle-holesymmetry. This leads to order O ( E ) corrections to thedrain current.
2. Majorana Bound States and Energy Splitting
If a vortex is present in the system described by Eq.(C1), ∆( r ) = ∆( r ) e i(cid:96)θ , a single Majorana bound statewill be localized in the vortex core when the vorticity (cid:96) is odd. Following the procedure in Ref. 60, the zeroenergy BdG equations for the model in Eq. (C1) with acentral vortex of vorticity (cid:96) = 1 are expressed in termsof two real and coupled radial equations, (cid:32) − µ η ∆( r ) + v F (cid:0) ∂ r + r (cid:1) − η ∆( r ) − v F ∂ r − µ (cid:33) (cid:32) u ↑ ( r ) u ↓ ( r ) (cid:33) = 0 . (C3)Here, v σ ( r ) = ηu σ ( r ) with η = ± R ≈ ξ , ∆( r ) = ∆ Θ( r − R ), the resulting zeromode is expressed in terms of Bessel functions, (cid:16) u ↑ ( r ) , u ↓ ( r ) (cid:17) T = (C4) N (cid:32) J ( µrv F ) J ( µrv F ) (cid:33) (cid:104) Θ( R − r ) + e − ∆0 vF ( r − R ) Θ( r − R ) (cid:105) . Above, N is given by normalization of the two-component spinor. Two vortex bound Majoranas sep-arated by a distance R will generally lift from zero en-ergy by an energy splitting exponentially small in R/ξ when R (cid:29) ξ . If µ is tuned close to the Dirac point (typi-cally achieved by bulk doping of screened charges ), thissplitting has previously been estimated to asymptoticallyapproach ε + ∼ − µ ( R/ξ ) / exp( − R/ξ ) up to some pref-actor of order one. Similarly, a superconducting disk of radius R (with acentral vortex) deposited on a TI supports a Majoranaedge state located exponentially close to the edge. Forweak magnetic fields the splitting between the centralbound state and the edge state has been estimated todecay like ε + ∼ − µl B e − Rξ /ξ, (C5)with l B the magnetic length.
3. Toy Model Calculation: Energy Splitting in aSpinless p -wave Superconductor For completeness, and serving as an illuminating ex-ample not found elsewhere, we explain how the energysplitting between edge states in a Corbino geometry canbe calculated in the spinless p + ip superconductor. Seee.g. Ref. 1 for general aspects of this system and Ref.61 for a presentation of the similar intervortex splitting.Finally, the problem of the Majorana energy hybridiza-tion in superconductor-semiconductor hybrid structuresis addressed in Ref. 62.The spinless p -wave superconductor realises a topo-logical phase for µ > µ <
0. The model has the corresponding zero energy BdGequations (cid:32) − m ∇ − µ p F { ∆( r ) , ∂ z ∗ }− p F { ∆ ∗ ( r ) , ∂ z } m ∇ + µ (cid:33) (cid:32) u ( r ) v ( r ) (cid:33) = 0 , (C6)with the operator ∂ z ∗ = e iθ ( ∂ r + ir ∂ θ ). We let a ra-dial annulus geometry between R and R be super-conducting with a vortex located in the center hole,∆( r ) = ∆ e iθ Θ( r − R )Θ( R − r ), with chemical potential2 mv F µ > ∆ (causing the wavefunctions to oscillate) and µ < R ≡ R − R (cid:29) ξ , we treatthe two single-edged systems separately and constructthe ground state candidates φ ± = √ ( ψ ± i ψ ), where ψ ( ψ ) is the solution localized on the inner (outer)edge. The two solutions are exponentially damped awayfrom their respective edges and they oscillate with fre-quency k = (cid:112) mµ − (∆ /v F ) . Moreover, if both R , R (cid:29) /k , the Bessel functions that appear as radialsolutions can be expanded asymptotically. The splittingintegral can be calculated, and we send µ → −∞ outsidethe annulus in the end. Upon evaluating the splittingintegrals and then taking the limit, we obtain ε ± = (cid:10) φ ± (cid:12)(cid:12) H (cid:12)(cid:12) φ ± (cid:11) ≈ ∓ µv F k sin ( kR ) e − R/ξ . (C7)The splitting is zero for particular values of the domainseparation. Whenever R = πn/k for n ∈ N , there aretwo degenerate ground states. The zero mode conditiongives associations of an interference phenomenon, causedby the oscillating edge modes that convolute in a destruc-tive manner for certain values of the separation. Theexpression in Eq. (C7) agrees with numerical diagonal-ization, already when R is a small multiple of ξ . TheCorbino geometry has been studied for small disk size,in which case the same zero energy criterion is found inthe limit 2 mv F µ (cid:29) ∆ by imposing Dirichlet boundaryconditions on the wavefunctions directly. Appendix D: Surface-Bulk Scattering with ScreenedDisorder
In Ref. 44 Fermi’s Golden rule is used to find the scat-tering rate from the surface ( S ) to the bulk ( B ) of aTI in the presence of static and dilute point impurities.Here, we take Ref. 44 as a starting point (the reader isreferred to this reference for further details) and apply the formalism to the case of long range scatterers. Asdescribed in the main text we study scattering from sur-face initial wave vector k = ( k (cid:107) ,
0) (in practice we let k (cid:107) = ( k x , k (cid:48) = ( k (cid:48)(cid:107) , k (cid:48) z ).The energy of the incoming surface state determines theFermi energy, ξ S, k (cid:107) = v F k (cid:107) ≡ (cid:15) F . Assuming low energyelastic scattering, and making use of the continuum limit (cid:80) k (cid:48) → V (2 π ) (cid:82) d k (cid:48) , we find the scattering rateΓ impSB ( (cid:15) F ) = 2 π (cid:88) k (cid:48) ,n (cid:48) | g imp k − k (cid:48) | (cid:0) | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | + | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | (cid:1) δ ( ξ B,n (cid:48) ,k (cid:48)(cid:107) − (cid:15) F ) (D1) ≈ V (2 π ) (cid:88) n (cid:48) :min { ξ B,n (cid:48) ,k (cid:48)(cid:107) } <(cid:15) F k (cid:48)(cid:107) |∇ ξ B,n (cid:48) ,k (cid:48)(cid:107) | (cid:90) ∞ d k (cid:48) z (cid:90) π d ϕ | g imp k (cid:107) − k (cid:48)(cid:107) ,k (cid:48) z ,ϕ | (cid:0) | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | + | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | (cid:1) . In going from the first to the second line we integratedover k (cid:48)(cid:107) , which in the last line is then defined implicitlyby ξ B,n (cid:48) ,k (cid:48)(cid:107) = (cid:15) F . Above, ϕ is the polar angle between theincoming and the outgoing momentum projected onto the k x k y -plane, i.e. k (cid:48) x = k (cid:48)(cid:107) cos ϕ and k (cid:48) y = k (cid:48)(cid:107) sin ϕ . The F ’sare defined as convoluted overlaps between the surfaceand the bulk states, F S, k (cid:107) ; B, k (cid:48) , n (cid:48) = (cid:10) Ψ S, k (cid:107) (cid:12)(cid:12) e ik (cid:48) z z (cid:12)(cid:12) Ψ B, k (cid:48) , n (cid:48) (cid:11) , (D2) F S, k (cid:107) ; B, k (cid:48) , n (cid:48) = (cid:10) Ψ S, k (cid:107) (cid:12)(cid:12) e ik (cid:48) z z (cid:12)(cid:12) Ψ B, k (cid:48) , n (cid:48) (cid:11) , (D3)where two bulk bands 1 n (cid:48) and 2 n (cid:48) are degenerate. Thefour-component wavefunctions in these expression arefound by solving the BdG equations for the 3D TI ex-actly at k (cid:107) = 0 and then applying perturbation theoryto leading order in the wave vector (see Appendix C inRef. 44 where the full expressions are listed in Eq. (C1)–(C11)). This procedure also yields the dispersion rela-tions to leading order in the parallel wave vector. In Fig.4 (a) and (b) the surface and some of the bulk statesare visualised in a thin film with model parameters asfor Bi Se , except for the bulk gap which is set to∆ TI = 0 . g imp k − k (cid:48) is here an ensemble averagedCoulomb potential due to screened charge impurities.Assuming randomly distributed impurities with zeromean, the coupling should be proportional to | g imp k − k (cid:48) | ∝ n ( Ze / ( (cid:15) (cid:15) r )) (cid:0) | k − k (cid:48) | + k T F (cid:1) − , where (cid:15) r is therelative permittivity, k T F is the Thomas-Fermi wave vec-tor, and n is the average dopant concentration. Incylindrical coordinates the coupling is expressed as | g imp k (cid:107) − k (cid:48)(cid:107) ,k (cid:48) z ,ϕ | = n V (cid:32) Ze (cid:15) (cid:15) r k (cid:107) (cid:33) (D4) × (cid:34)(cid:16) − k (cid:48)(cid:107) /k (cid:107) (cid:17) + (cid:0) k (cid:48) z /k (cid:107) (cid:1) + r s + 4 k (cid:48)(cid:107) k (cid:107) sin ϕ (cid:35) − , where r s = k T F /k (cid:107) was introduced. In the main textwe assume that the screened charges have Z = 1. Withthis coupling established we define the differential crosssection for the screened Coulomb scattering asd σ ( n (cid:48) )long ( (cid:15) F )d k (cid:48) z ≡ (D5) (cid:90) π d ϕ | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | + | F S, k (cid:107) ; B, k (cid:48) , n (cid:48) | (cid:20)(cid:16) − k (cid:48)(cid:107) k (cid:107) (cid:17) + (cid:16) k (cid:48) z k (cid:107) (cid:17) + r s + 4 k (cid:48)(cid:107) k (cid:107) sin ϕ (cid:21) . This function is shown in Fig. 4 (c) for r s = 0 . (cid:15) r = 100. Note that in the case of pointscatterers, the differential cross section is obtained simplyby replacing the denominator in Eq. (D5) by 1.Using Eq. (D4) and (D5) in (D1) assembles the expres-sion for the scattering rate in the main text, Eq. (23). Wenote that the rates as shown in figure in Fig. 4 (d) areseen to be in good agreement with the Ref. 44.
1. Effect of Spatially Confined TransverseWavefunction
When applying the formula of the electronic scatteringrate in Eq. (23) to the case of the surface Majorana weneglect the transverse profile of the surface wavefunction(inset of Fig. 4 (a)). Here, we argue why this is a validapproximation in our parameter regime.Assume for simplicity that the two surface gaps areof similar magnitude ∆ ≈ M (cid:29) µ , so that thethe transverse wavefunction is confined over a the scale σ y ∼ ξ = v F / ∆ . By the uncertainty principle, this con-finement leads to an uncertainty in k (cid:48) y of σ k (cid:48) y ∼ ξ − . Tosimulate this uncertainty we draw random additions to k (cid:48) y = k (cid:48)(cid:107) sin ϕ from a normal distribution with the stan-dard deviation above and zero mean for each scattering3direction ϕ . The resulting scattering rates display an av-erage absolute deviation (cid:104)| ε |(cid:105) (averaged over the Fermienergy range in Fig. 4 (d)) from the sharp k (cid:48) y curve ofroughly (cid:104)| ε |(cid:105) ≈ ( σ k (cid:48) y / max { k (cid:48) y } ) . Since max { k (cid:48) y } (cid:29) ξ − in the energy range and with the proximity gaps we con-sider, the effect of uncertainty in k (cid:48) y , and hence confine-ment in the y -direction, can be ignored. Appendix E: Scattering with Acoustic Phonons
Consider the toy model coupling electrons to (acoustic)phonons through a deformation potential, H ep = (cid:88) q , q ,s,s (cid:48) M s,s (cid:48) q , q c † q + q ,s c q ,s (cid:48) (cid:16) a q + a †− q (cid:17) . (E1)For small momenta, the electron-phonon coupling goes as | M ( q ) | ∼ q for surface phonons. Acoustic phononshave linear dispersion, ω ( q ) = c R | q | and are ex-pected to dominate the coupling Hamiltonian at low temperatures. The scattering rate follows once againfrom Fermi’s Golden rule, Γ ph i → f = 2 πν f ( E ) | (cid:104) f | H ep | i (cid:105)| ,where ν f ( E ) is the final density of states. Recall also thatstates with a linear dispersion in two dimensions have ν ( E ) ∝ E . In a superconducting system, there will be anon-zero amplitude for the creation of a phonon with thecost of annihilating two quasiparticles. For illustrativepurposes, we consider quasiparticles excitations of the( s -wave) Bardeen-Cooper-Schrieffer (BCS) ground state | Ω (cid:105) , | i ; σ , σ (cid:105) = γ † k ,σ γ † k ,σ | Ω (cid:105) , (E2) | f (cid:105) = a † q | Ω (cid:105) , (E3) | Ω (cid:105) = (cid:89) k ( u k + v k c † k , ↑ c †− k , ↓ ) | (cid:105) . (E4)Above, the quasiparticle creation operators are γ † k , + = u ∗ k c † k , ↑ − v ∗ k c − k , ↓ and γ † k , − = u ∗ k c †− k , ↓ + v ∗ k c k , ↑ . The scat-tering element is found to be (cid:104) f | H ep | i ; σ , σ (cid:105) = (cid:16) M ↑ ,s ( σ ) σ k , − σ k − k δ q ,σ k + k − M ↓ ,s ( σ ) σ k , − σ k + k δ q ,σ k − k (cid:17) v ∗− k u ∗ k − (cid:16) M ↑ ,s ( σ ) σ k , − σ k − k δ q ,σ k + k − M ↓ ,s ( σ ) σ k , − σ k + k δ q ,σ k − k (cid:17) v ∗− k u ∗ k . (E5)Above, we introduced the symbol s (+) = ↑ and s ( − ) = ↓ .The scattering amplitude depends only on the momen-tum transfer for small energies, in which case the couplingabove vanishes identically. This is presumably becausethe average charge of the quasiparticles becomes zero inthe low energy limit.A careful analysis for the electron decay rate aloneyields a Γ ph ∼ T law well below the Bloch-Gr¨uneisentemperature at the Fermi surface. Away from the Fermi surface, by biasing the electrons at a small voltage V , thedecay rate is finite at T = 0 and goes as Γ ph ∼ ( eV ) . In-serting the exact prefactor (for Bi Te ), by following thesteps in Ref. 65, we find a decay rate in the peV rangewhen V (cid:39) µ V. This means that the electron lifetimedue to acoustic phonon scattering τ ph = 1 / Γ ph is in thems range, making this effect negligible at T = 0 at lowvoltage. ∗ [email protected] J. Alicea, Rep. Prog. Phys. , 076501 (2012). N. Read and D. Green, Phys. Rev. B , 10267 (2000). N. B. Kopnin and M. M. Salomaa, Phys. Rev. B , 9667(1991). G. E. Volovik, J. Exp. Theor. Phys. , 609 (1999). V. Gurarie and L. Radzihovsky, Phys. Rev. B , 212509(2007). C. Kallin and J. Berlinsky, Rep. Prog. Phys. , 054502(2016). L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407 (2008). F. Qu, F. Yang, J. Shen, Y. Ding, J. Chen, Z. Ji, G. Liu,J. Fan, X. Jing, C. Yang, and L. Lu, Sci. Rep. , 339(2012). A. D. K. Flinck, C. Kurter, E. D. Huemiller, Y. S. Hor,and D. J. V. Harlingen, Appl. Phys. Lett. , 203101 (2016). S. Sasaki, M. Kriener, K. Segawa, K. Yada, Y. Tanaka,M. Sato, and Y. Ando, Phys. Rev. Lett. , 217001(2011). J.-P. Xu, C. Liu, M.-X. Wang, J. Ge, Z.-L. Liu, X. Yang,Y. Chen, Y. Liu, Z.-A. Xu, C.-L. Gao, D. Qian, F.-C.Zhang, and J.-F. Jia, Phys. Rev. Lett. , 217001 (2014). G. Koren, T. Kirzhner, E. Lahoud, K. B. Chashka, andA. Kanigel, Phys. Rev. B , 224521 (2011). P. Zareapour, A. Hayat, S. Y. F. Zhao, M. Kreshchuk,A. Jain, D. C. Kwok, N. Lee, S. Cheong, Z. Xu, A. Yang,G. D. Gu, S. Jia, R. J. Cava, and K. S. Burch, Nat.Commun. , 1056 EP (2012). V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P.A. M. Bakkers, and L. P. Kouwenhoven, Science ,1003 (2012). S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon,J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani,Science , 602 (2014). A. R. Akhmerov, J. Nilsson, and C. W. J. Beenakker,Phys. Rev. Lett. , 216404 (2009). L. Fu and C. L. Kane, Phys. Rev. Lett. , 216403 (2009). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010). Q. L. He, L. Pan, A. L. Stern, E. C. Burks, X. Che, G. Yin,J. Wang, B. Lian, Q. Zhou, E. S. Choi, K. Murata, X. Kou,Z. Chen, T. Nie, Q. Shao, Y. Fan, S.-C. Zhang, K. Liu,J. Xia, and K. L. Wang, Science , 294 (2017). Y. Ando, J. Phys. Soc. Jpn. , 102001 (2013). F. Ortmann, S. Roche, S. O. Valenzuela, L. W.Molenkamp, and I. Aguilera,
Topological insulators : fun-damentals and perspectives (Wiley-VCH, Weinheim, Ger-many, 2015). S. Wolgast, C¸ . Kurdak, K. Sun, J. W. Allen, D.-J. Kim,and Z. Fisk, Phys. Rev. B , 180405 (2013). X. Zhang, N. P. Butch, P. Syers, S. Ziemak, R. L. Greene,and J. Paglione, Phys. Rev. X , 011011 (2013). J. G. Checkelsky, Y. S. Hor, M.-H. Liu, D.-X. Qu, R. J.Cava, and N. P. Ong, Phys. Rev. Lett. , 246601 (2009). B. Skinner, T. Chen, and B. I. Shklovskii, J. Exp. Theor.Phys. , 579 (2013). K. T. Law, P. A. Lee, and T. K. Ng, Phys. Rev. Lett. ,237001 (2009). J. Li, G. Fleury, and M. B¨uttiker, Phys. Rev. B , 125440(2012). P. Fendley, M. P. Fisher, and C. Nayak, Ann. Phys. ,1547 (2009), july 2009 Special Issue. W. Bishara and C. Nayak, Phys. Rev. B , 155304 (2009). B. Rosenow, B. I. Halperin, S. H. Simon, and A. Stern,Phys. Rev. Lett. , 226803 (2008). B. Rosenow, B. I. Halperin, S. H. Simon, and A. Stern,Phys. Rev. B , 155305 (2009). H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C.Zhang, Nat , 438 (2009). R. S. Akzyanov, A. L. Rakhmanov, A. V. Rozhkov, andF. Nori, Phys. Rev. B , 075432 (2015). M. Cheng, R. M. Lutchyn, V. Galitski, and S. Das Sarma,Phys. Rev. B , 094504 (2010). C.-K. Chiu, D. I. Pikulin, and M. Franz, Phys. Rev. B ,165402 (2015). H. Beidenkopf, P. Roushan, J. Seo, L. Gorman, I. Drozdov,Y. S. Hor, R. J. Cava, and A. Yazdani, Nat. Phys. ,10.1038/nphys2108. C. Caroli, P. D. Gennes, and J. Matricon, Phys. Lett. ,307 (1964). G. M¨oller, N. R. Cooper, and V. Gurarie, Phys. Rev. B , 014513 (2011). T. Mizushima and K. Machida, Phys. Rev. A , 053605(2010). A. L. Rakhmanov, A. V. Rozhkov, and F. Nori, Phys.Rev. B , 075141 (2011). R. S. Akzyanov, A. V. Rozhkov, A. L. Rakhmanov, andF. Nori, Phys. Rev. B , 085409 (2014). C. Nayak, S. H. Simon, A. Stern, M. Freedman, andS. Das Sarma, Rev. Mod. Phys. , 1083 (2008). This can be proven formally; the proof is a straightfor-ward generalization of Eq. (18) and (19) when identifyingthe structure of the sub blocks in T . Blocks characterizingtransport between different leads always give a vanishingcontribution to the average current, related to the obser- vation right below Eq. (20). If two leads are not causallyconnected by the chiral arm, the corresponding block inthe T matrix is zero. Finally, no matter which arms the N leads are coupled to, the drain current is reduced by thefactor (cid:81) Nj =1 r j , where r j is the local reflection amplitudeof lead j . K. Saha and I. Garate, Phys. Rev. B , 245418 (2014). W.-Y. Shan, H.-Z. Lu, and S.-Q. Shen, New J. Phys. ,043048 (2010). D. Culcer, E. H. Hwang, T. D. Stanescu, andS. Das Sarma, Phys. Rev. B , 155457 (2010). D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier,J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong,A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor,R. J. Cava, and M. Z. Hasan, Nature , 0028 (2009). H. Li, T. Zhou, J. He, H.-W. Wang, H. Zhang, H.-C. Liu,Y. Yi, C. Wu, K. T. Law, H. He, and J. Wang, Phys. Rev.B , 075107 (2017). N. Borgwardt, J. Lux, I. Vergara, Z. Wang, A. A. Taskin,K. Segawa, P. H. M. van Loosdrecht, Y. Ando, A. Rosch,and M. Gr¨uninger, Phys. Rev. B , 245149 (2016). Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Ban-sil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan,Nat. Phys. , 389 (2009). J. A. Sobota, S. Yang, J. G. Analytis, Y. L. Chen, I. R.Fisher, P. S. Kirchmann, and Z.-X. Shen, Phys. Rev. Lett. , 117403 (2012). M. Neupane, S.-Y. Xu, Y. Ishida, S. Jia, B. M. Fregoso,C. Liu, I. Belopolski, G. Bian, N. Alidoust, T. Durakiewicz,V. Galitski, S. Shin, R. J. Cava, and M. Z. Hasan, Phys.Rev. Lett. , 116801 (2015). The prefactor of Γ impSB in Eq. (23) has a k − (cid:107) dependency.There is, however, also a less significant dependency on k (cid:107) in d σ ( n (cid:48) )long ( (cid:15) F ) / d k (cid:48) z , which can be inferred from Fig. 4 (b)(recall that (cid:15) F ∝ k (cid:107) ) to be very close to linear. S. K. Kushwaha, I. Pletikosi, T. Liang, A. Gyenis, S. H.Lapidus, Y. Tian, H. Zhao, K. S. Burch, J. Lin, W. Wang,H. Ji, A. V. Fedorov, A. Yazdani, N. P. Ong, T. Valla, andR. J. Cava, Nat. Commun. , 11456 (2016). R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwen-hoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, (2017),arXiv:1707.04899 [cond-mat.supr-con]. D. J. Kim, S. Thomas, T. Grant, J. Botimer, Z. Fisk, andJ. Xia, Sci. Rep. , 3150 EP (2013). S. Lee, X. Zhang, Y. Liang, S. W. Fackler, J. Yong,X. Wang, J. Paglione, R. L. Greene, and I. Takeuchi,Phys. Rev. X , 031031 (2016). H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C.Zhang, Nat. Phys. , 438 (2009). M. B¨uttiker, Phys. Rev. B , 12485 (1992). J. D. Sau, S. Tewari, R. M. Lutchyn, T. D. Stanescu, andS. Das Sarma, Phys. Rev. B , 214509 (2010). M. Cheng, R. M. Lutchyn, V. Galitski, and S. Das Sarma,Phys. Rev. Lett. , 107001 (2009). C. Knapp, T. Karzig, R. Lutchyn, and C. Nayak, (2017),arXiv:1711.03968 [cond-mat.mes-hall]. B. Rosenstein, I. Shapiro, and B. Y. Shapiro, EPL ,57002 (2013). J. M. Ziman,
Electrons and Phonons , 5th ed. (Oxford Uni-versity Press, Walton Street, Oxford, OX2 6DP, 1960). S. Giraud and R. Egger, Phys. Rev. B , 245322 (2011). V. Parente, A. Tagliacozzo, F. von Oppen, and F. Guinea,Phys. Rev. B88