Topological Defects in Topological Insulators and Bound States at Topological Superconductor Vortices.

The scattering of Dirac electrons by topological defects could be one of the most relevant sources of resistance in graphene and at the boundary surfaces of a three-dimensional topological insulator (3D TI). In the long wavelength, continuous limit of the Dirac equation, the topological defect can be described as a distortion of the metric in curved space, which can be accounted for by a rotation of the Gamma matrices and by a spin connection inherited with the curvature. These features modify the scattering properties of the carriers. We discuss the self-energy of defect formation with this approach and the electron cross-section for intra-valley scattering at an edge dislocation in graphene, including corrections coming from the local stress. The cross-section contribution to the resistivity, ρ, is derived within the Boltzmann theory of transport. On the same lines, we discuss the scattering of a screw dislocation in a two-band 3D TI, like Bi1-xSbx, and we present the analytical simplified form of the wavefunction for gapless helical states bound at the defect. When a 3D TI is sandwiched between two even-parity superconductors, Dirac boundary states acquire superconductive correlations by proximity. In the presence of a magnetic vortex piercing the heterostructure, two Majorana states are localized at the two interfaces and bound to the vortex core. They have a half integer total angular momentum each, to match with the unitary orbital angular momentum of the vortex charge.

Introduction

Boundaries in a topological insulator (TI) host Dirac electrons propagating with a linear dispersion in energy [1]. On the other hand, it was recognized long ago that the low energy electronic properties of a graphene sheet, which is a weak topological material, is a semimetal and can be described close to the neutrality point with the Dirac Hamiltonian [2].

In the recent past, the charge carrier mobility in a single graphene layer has been extensively investigated [3,4]. Graphene resistivity was experimentally found to be inversely proportional to concentration n of the charge carriers, which means that their mobility is almost independent of n [5-7]. This uncommon behavior cannot be attributed to short-range potentials, due to defects on the scale of the lattice parameter, a. Indeed, these defects are believed to provide only a small additional contribution to the resistivity proportional to the impurity concentration [8,9]. Instead, the contribution due to charged impurities, acting as long-range Coulomb scatterers [10], seems to fit with the experimental finding in the case of graphene sheets attached to substrates [8,11-13]. In the case of suspended graphene [14], the matter is still unsettled and has motivated an extended search for other scattering mechanisms. Particularly, corrugations and ripples have been found, and the systematic investigation of their effect on the electrical and optical properties are far from being clarified [3,15-18]. Suspended graphene should not be so influenced by the charge impurities, except for, possibly, trapped clusters in the wake of the corrugations [19,20]. Other sources of scattering, topological point-like lattice defects, can induce only little stretching locally, but contribute significantly to the reduction of the mobility. We find that this is the case of an edge dislocation centered on a pentagon-heptagon defect, which does not generate any curvature in the graphene sheet. Furthermore, wedge disclinations due to isolated pentagons or heptagons in the lattice structure could be considered, although they are believed to cost higher formation energy.

On the other hand, the analysis of the conductance properties of boundary states at the surface of a three-dimensional (3D) TI, like Bi2Se3, Bi2Te3 and Bi1−Sb, is still in its infancy. Experimentally, it is difficult to tune the Fermi energy within the gap of the 3D TI to study the mobility of the Dirac states located at the boundaries. Besides, in most cases, mesoscopic samples, like TI flakes or TI nanowires, are plagued by impurity bands within the insulating gap. Therefore, it is difficult to isolate the contribution of the Dirac electrons from the bulk contribution [21-25]. Measures of the mobility as derived from Shubnikov–de Haas oscillations or Hall resistance give different results. Aharonov–Bohm oscillations of

Dirac electrons at the surface of nanowires have been measured [26-28]. To the best of our knowledge, the influence of topological defects, such as the screw dislocations or the wedge disclinations on the transport properties, has not been investigated yet.

In a previous paper [29], we considered the contributions to the resistivity, due to elastic scattering by a smooth curved bump at the surface of a 3D TI [30]. We considered the unrelaxed situation and attacked the problem directly in the continuum low energy limit, by introducing the defect as a change of the metric experienced by the electrons [18,31]. In the unrelaxed lattice, the local Lorentz invariance is assumed to be still conserved. The defect itself is effectively described as an Aharonov–Bohm flux pinned at its center.

In the case of graphene, we see that the contribution of a collection of isolated edge dislocation to resistivity in the Boltzmann semiclassical limit is found to be ∝ 1/n. Scattering of an edge dislocation also depends on the orientation of the Burgers vector and, therefore, on the angle of the incoming electron wave, but this feature does not introduce any additional k dependence, if some averaging over randomness is implied. Besides, the resistivity of the edge dislocation is found to be finite, close to the neutrality point. We find that any additional contribution to the resistivity due to weak local strain vanishes close to the neutrality point with positive powers of k. Here, k is the Fermi wavevector and ℓ is the size of the perturbation, which is much larger than the lattice spacing, a. Although the graphene monolayer is characterized by two Dirac cones centered at the two valley points, K and , in the Brillouin zone [32], elastic scattering potentials, smooth on the lattice scale, such as ripples or point-like defects, are likely not to change the isospin and the valley of the scattering electrons, nor their energy spectrum drastically, unless the strain produces gaps or zero energy states [20,33]. We point out that intravalley scattering produced by edge dislocations can give a relevant contribution to the resistivity, especially in the absence of sizable scattering by charges or resonant impurities. The question arises: how large is their formation energy in graphene? We will address this point in Section 3.3.

Similarly, the presence of screw dislocations of Burgers vector can be envisaged in a 3D TI. In the present work, we apply the analytical methods used to discuss an edge dislocation in the graphene sheet to a screw dislocation in a two-band 3D TI. In this case, helical bound states can propagate along its axis. The Hamiltonian can be expanded linearly in close to a time reversal invariant momentum (TRIM), M. Two gapless modes of opposite helicity are present, provided the constraint, , is satisfied [34]. The reference case is the alloy, Bi1−Sb, while TIs like Bi2Se3 and Bi2Te3 cannot fulfill the constraint, because their 3D Dirac point is located at Γ. We exhibit the wavefunctions of the bound states and their properties in Section 5.

Superconductive proximity induced at an interface between a 3D TI and a conventional superconductor has been attracting a lot of interest recently, due to the expectation that Majorana bound states (MBS) could exist under appropriate conditions [35-37]. In the presence of a magnetic field, the screw dislocation could host a vortex line with little energy cost. A vortex line piercing a 3D TI sandwiched between two even-parity superconductors could bind two Majorana quasiparticle excitations at the opposite interfaces. After the pioneering work of Fu and Kane, this feature has been discussed by various authors [38-40] in great detail. Some of us, by solving the Bogolubov–de Gennes equations analytically in the long wavelength limit, exhibited the wavefunctions of the MBS in the two-band 3D TI model and discussed how the MBS depends on the parity of the order parameter of the superconductor inducing the proximity [41]. These results are briefly summarized in Section 6. When proximity involves odd-parity pairing, the modes appearing in the superconducting gap opened in the Dirac dispersion are charged surface Andreev bound states. They originate from interfacial circular states of definite chirality, centered at the vortex singularity and decay with damped oscillations away from the interfaces of the TI film. The case of d-wave-induced pairing in quasi-one-dimensional TI nanostructures was also discussed [42,43].

In Section 2, we sketch our approach by introducing some generalities about the change of coordinates in the presence of defects and curvature, to make the paper self contained [44]. In Section 3, we discuss the edge dislocation in a graphene sheet and compare the contributions to the cross-section coming from the topological defect when the lattice is unrelaxed with the one due to the stress-induced relaxation. The latter turns out to be negligible with respect to the former in the low electron density limit. A detailed account of the derivation of the phase shifts for the scattering of an Aharonov–Bohm (A–B) flux can be found in Appendix A. The contribution to the resistivity is also calculated, in the relaxation time approximation, as well as the self-energy of the defect, which requires Green’s function, which is derived in Appendix B. In Section 4, we introduce a two-band model for the 3D TI and refer to Appendix C for the calculation of the wavefunctions of the Dirac dispersed states localized at the boundary. Section 7 concludes the paper with a summary and few final remarks.

Dirac Electrons on a Free Surface in the Long Wavelength Limit

The long wavelength dynamics of Dirac electrons propagating on a flat two-dimensional boundary at energies close to the neutrality point of the Dirac cone is described by the Dirac equation:

where a = 0,1,2, and the Dirac matrices are γ0 = −iσz, γ1 = σ, γ2 = −σ. In this Section, we consider the generalization of the algebra of the Dirac matrices describing a flat surface, to include topological defects and possible curvatures, in the absence of strain. According to the Equivalence Principle [18], given a frame to be referred to as the “curved frame”; from now on, it is possible to introduce a local flat, x, frame at each point. The components of the Jacobian matrix for the transformation from the coordinates, x, defined on the whole manifold and the local parametrization, are the tetrads, e :

The inverse of the tetrads, e , are defined through the orthogonality relation e e = δ . Our aim is to substitute the Minkowski metric, η, of the flat space with the metric, g, which is singular in the case of a topological defect. We follow the convention to use Latin letters and overlined numbers, a, b, …, 1̄,2̄, to refer to the local frame, as opposed to the Greek letters, which are refer to the curved frame. The tetrads satisfy the completeness relation, e e = δ . They are linked to the metric according to:

The Dirac matrices γ = e γ satisfy the anticommutation relation:

On formulating the Lorentz covariance of the Dirac equation locally, we replace the derivatives with the covariant derivatives: where ∑ are the generators of the spinorial representation of the Lorentz group and are expressed in terms of the commutators as ∑ = i/2[γ, γ]. The connection coefficients, Γ , are given by Γ = e . The massless Dirac equation on curved space then reads:

Making explicit the rotation of the Dirac matrices and the covariant derivatives, the stationary part of Dirac equation on curved space time is:

In full generality, defining the components of the affine connection as Γ = (∂) e, the components of the torsion are:

while the Riemann tensor is defined as:

Equation (8) is clearly zero if the tetrads are regular functions, i.e., if they have continuous second derivatives, due to the Schwartz lemma. If the Riemann tensor also vanishes, the transformation, x(x), is just a change of coordinates of the flat space. Tetrads that are not regular could give curvature, torsion or both. An unrelaxed topological defect does not introduce any curvature, so that only the rotation of the Dirac matrices has to be accounted for.

The Edge Dislocation in Graphene

An edge dislocation in the graphene sheet, centered on the origin, is obtained by cutting the plane; let us say in correspondence of the negative x1 half axis and by adding a half line of carbon atoms [45]. In the case of a hexagonal lattice, the line of atoms could be interpreted as an armchair row. In fact, the whole honeycomb lattice could be recovered by arranging arrays of armchair rows in a square structure. The dislocation produces a pentagon-heptagon pair in the origin of the frame, as shown in Figure 1. In the unrelaxed structure, the other hexagons are not affected. This kind of defect does not mix the valleys. Smooth potentials are likely not to change the isospin and the valley of the scattering electrons, nor their energy spectrum, drastically, unless the strain produces gaps or zero energy states [4,33].

Figure 1.

Edge dislocation appearing as a pentagon-heptagon pair in the perfect lattice of a graphene monolayer.

Let be the Burgers vector pointing down along the x2 –axis, orthogonal to the negative x1-axis. The Burgers vector defines a singular coordinate transformation:

where the branch cut of the inverse tangent is on the negative x1-axis. The tetrads are easily derived from the infinitesimal transformation dx = e , thus obtaining:

Equation (11), as such, describes an unrelaxed configuration. The equilibrium configuration could be restored by adding an effective gauge potential to the Dirac equation, which may further change the spin connection, without influencing, however, the holonomy on the wave function, since the elastic deformation does not add any curvature. This is a way of restating the Saint Venant conditions for the two-dimensional case.

The Dirac matrices γ = e γ can be easily found by inversion of the tetrads:

In the case of an edge dislocation, the Riemann tensor generated by these tetrads, R , vanishes [31], so that the connection on an edge dislocation can be put to zero. On the contrary, the torsion is δ – like, as it encodes the mismatch in the parallel transport occurring at the topological defect [46]:

Since the connection vanishes, the spin connection is also zero. Therefore, only the rotation of Dirac matrices appears in the Dirac equation.

Let be the component of the wavevector referring to the valley point, , and θ the angle of the local momentum, . If is a spinor satisfying the equation :

the solution of the rotated Dirac equation can be written down in full generality as (μ = 0, 1, 2) [47]:

Here, the integral is to be performed along the geodesic path, C, connecting a reference point, , with the actual point, . The linear transformation e = k defines k(r).

This is indeed a solution, because, by substituting the spinor of Equation (15) in the Dirac equation, we get:

where the equality to zero stems from the definition of Φ0 given in Equation (14).

Equation (15) is a general way to take account of the rotation of the Dirac matrices, and it provides the full solution when no spin connection is present. In the case of the edge dislocation, the differential form appearing in the phase of Equation (15), when calculated using the tetrads of Equation (11), is:

The quantities, , represent the components of the momentum in the local inertial frame. The curl of this differential form is zero. Indeed, it is the differential of the accumulated phase given by:

It follows that the integral in Equation (15) is independent of the path, and it is readily calculated. The full spinor solution for electrons having a vector close to the valley point, , includes a fast oscillating plane wave prefactor, , giving:

The edge dislocation produces a vortex-like singularity in the graphene sheet, unless , that is, unless the propagation is along the branch cut. The solution corresponds to a plane wave scattering of a flux line of flux , piercing the sheet at the origin, and acquiring an Aharonov–Bohm phase in circulation. The single valuedness of the wavefunction in circulating around the dislocation and the C3 rotational symmetry of the graphene lattice requires that .

Cross-Section of the Unrelaxed Topological Defect

The phase shifts of a particle incoming with momentum scattered with angular momentum quantum number m by the edge dislocation, , can be easily derived [48,49]. They are (see Appendix A for the details):

Defining the outgoing wave (ϕ ≡ θ − θ):

the scattering amplitude, , is:

The total cross-section, with given by Equation (22), is:

Close to the neutrality point (limit ):

The resistivity, ρ, can be estimated with the Boltzmann relaxation-time approximation as:

where is the density of states at the Fermi level for double spin, but one valley The usual definition of the total relaxation rate is:

The relaxation rate, related to the imaginary part of the self-energy, is expressed in terms of the t – matrix element < k′|t|k > for an outgoing circular wave, of incoming wavevector at the Fermi energy, scattered elastically into a plane wave of wavevector by the extra potential arising in Equation (7), of scattering amplitude . It is:

(both waves are normalized to the square root of the area, πR2). It depends on the energy, ϵ, and on the scattering angle, between the incoming and outgoing wave. The factor provides the cancellation of the backward scattering. The integration over ϕ gives a Bessel function, J(kr), and the integral over r can be approximated as:

Finally:

Now, the sum can be performed. Defining α as the non-integer part of the flux, , with α < 1, we obtain :

The prefactor of Equation (26), , makes the first term in Equation (30) disappear, which is likely to be spurious anyway [48]. It also compensates for the divergence at Θ = 0 of the second term. The final result is:

Equation (31) depends on the incoming direction due to the orientation of the Burgers vector contained in α.

When multiplying this result by the number of dislocations, z, and after averaging over their random distribution, the resistivity of Equation (25) due to the elastic scattering on the edge dislocation turns out to be:

The resistivity is proportional to the density of the defects and inversely proportional to the density of carriers . Lattice relaxation, around the branch point, provided it is not too strong, would contribute to the resistivity with a term that is a higher positive power of k [20] and would not change this result qualitatively. This is proven in the Section 3.2.

Cross-Section Due to the Stress at the Origin of the Edge Dislocation

The singularity point, which is the origin of the branch cut, due to the edge dislocation, is also the center of a strain texture induced by the defect. In this section, we show that the contribution of the strain to the total cross-section of Equation (23), due to the stress at the dislocation, is higher order with respect to the one that contributes to the relaxation time, τ (k), calculated in the previous section. The stress provides an extra term coupling the sublattices, A and B, within the same valley (here a (r), b (r) are fields on the two sublattices) which takes the form [32,50]:

where βκ is the stiffness of the lattice. Here, u (r) is derived from the displacement field around the dislocation center, , which is taken to decay radially as 1/r by approximation:

Explicitly, the potential arising from the strain due to the edge dislocation, away from the core of the defect (r > a) is, in cylindrical coordinates,

and is limited to an area < πℓ2, where ℓ is the mean free path between dislocations. The incoming Dirac electron wave function:

is scattered by the potential. The Green’s function is required, which solves the Dirac equation inclusive of the A-B flux, f, at the origin:

Its spectral representation is given in Appendix B. In the Born approximation, we get:

We keep just the contribution coming from the pole Equation (A20), and we take the direction of the incoming orthogonal to as the reference direction.

Using the decomposition of a plane wave in angular momentum eigenfunctions, Equation (36) becomes:

The integral giving the scattered part of the wave function is (B is a constant, including βκ and a normalization factor):

where I are shorthand for the two contributions to the integral given above. Let us analyze the first one:

These integrals are special cases of the Weber–Schafheitlin integral:

with (α, p) = (n + 3, 1) and (n + 2, 0). Thus, the first integral is equal to −1/2k, while the other is 1/2k. Putting things together:

Similarly, the second contribution yields:

where the same integral as Equation (42) appears, with (α, p) = (n − 1, −3) and (n, −2), so that the final result is:

The dominant contributions to the t – matrix scattering amplitude in the long wavelength limit come from the terms in which the order of the Bessel functions is the lowest possible, i.e.:

Putting aside the consideration that these matrix elements imply incoming waves of relatively high order (i.e., |n| = 2, 4), the largest ones among them lead to integrals whose limiting form is of the kind:

We conclude that their contribution to the cross-section goes as:

which is higher order, when compared with , which was derived in Section 3.1.

Self-Energy of the Dislocation

We now evaluate the formation the self-energy of the dislocation. In the long wavelength limit, we cannot account for the cost of the pentagon-heptagon defect formation, but we can include the strain cost, which is long range, because it decays as 1/r. We will make use of the Hellmann–Feynman theorem [51]:

by adding a perturbation potential, , of the kind of Equation (35) to the unperturbed Hamiltonian :

Gλ is the Green function for our system, with an A-B flux, λf, at the origin. λ is a parameter introduced for convenience, which we integrate out from zero to one:

We have subtracted a reference term in which the dislocation is absent. This is independent of λ and is immaterial, except for the fact that it provides the vanishing of the perturbation, when λ → 0. The prescription, e, is fixed by the requirement that one has to impose t′ → t + 0+ to get the correct ordering in the Green’s functions, which is .

Using the Dyson equation , we obtain:

The strain potential of the edge state is taken from Equation (35). It is defined for r > a, where a is the radius of the core of the dislocation.

According to Equation (A13), . On the other hand, the spectral representation of Gλ, Equation (A25), which includes the half pole only , is:

The correct time ordering requires that, before taking the limit, the variables are exchanged: r′ → r and θ ↔−θ. After multiplying the matrices and taking the trace, we get:

Now the integral over ω. The largest contribution comes from ω → 0 and/or r → 0. Hence, we substitute the factor, ω, with sin(ωa)/a, to regularize the integral. Using dimensionless variables, t = 2ωr, b = a/(2r) → 0, we perform the integration on the real axis:

Plugging all together, the awkward prefactor, e2, disappears, and the final result is (` is the size of the defect):

This is the expected self-energy for a long-range strain potential.

Effective Theory for Boundary States in 3D Topological Insulators

Bismuth-based materials are mostly studied since the prediction that Bi(1)Sb is a strong 3D TI [52], which was confirmed soon after with angle-resolved photoemission spectroscopy (ARPES) [53]. The stoichiometric crystal, Bi2Se3, is the prototype of a class of 3D TIs. This material was predicted to have boundary states with energy dispersion forming a Dirac cone centered at the Γ point located within the insulating gap [54]. The prediction has been experimentally confirmed [55]. The atomic structure of the material consists of the stacking of quintuple layers. While the coupling is strong inside each of the quintuple layers, the coupling between quintuple layers is much weaker, predominantly of the van der Waals type. The non-trivial topology in the band structure of the material stems from the inversion, at the Γ point, of the bands of the p orbital character, arising from Se and Bi atoms, due to the strong spin orbit (SO) coupling. A minimal tight-binding model represents the wavefunction for the valence and conduction bands of opposite parity, on the basis of four p− orbitals, spanning the spin and the orbital degree of freedom, Φ = (|p1 ↑〉, |p2 ↑〉, |p1 ↓〉, |p2 ↓〉) (the upper label, , means “transpose”) as:

where are lattice vectors, in terms of the slowly-varying envelope functions Ψα( ). The symmetries of the crystal can be represented as:

the time reversal symmetry

the inversion symmetry

the three-fold rotational symmetry around the z-axis

Here, is the complex conjugation, σ are the Pauli matrices in the spin space and τ are the Pauli matrices acting in the orbital space.

For states with in the vicinity of the Γ point, the Hamiltonian can be written in the or effective mass approximation as [56]: where v is the Fermi velocity close to the chemical potential, μ. h0 − μ does not contribute with any interesting feature to the model and will be dropped in the following. In the bulk, −i → and with M > 0 and C1, C2 < 0, to implement the inversion of the bulk bands [56]. The relative sign between M and C12 qualifies the insulator as being topologically non-trivial or trivial. ħv M is the bulk gap of the two-band model. The bulk Hamiltonian can be rewritten in compact form as: where , , . Note that these matrices satisfy the Clifford algebra {γ, γ} = η, i.e., they are a set of Dirac matrices. This means that they can be rotated as every well-educated set of Dirac matrices.

Other important symmetries that could be present are the particle-hole symmetry, Ξ, and the chiral symmetry, Γ. is an antiunitary symmetry, such that Ξ H(k)Ξ = −H(−k). The chiral symmetry Γ = i · ΘΞ is a unitary transformation that changes the sign of the Hamiltonian without reversing the sign of . In the presence of a non-zero chemical potential, the Hamiltonian of Equation (58) has neither the particle-hole symmetry, Ξ, nor the chiral symmetry, Γ. As Θ2 = −1, Ξ2 = 0, Γ2 = 0, it belongs to the class, AII, of the Altland–Zirnbauer classification [57].

In Appendix C, we derive the boundary states for a flat surface orthogonal to the z^−axis, with the simplifying assumption that C12 = 0 and that M changes sign at the boundary. This provides localized states at the boundary, which decay exponentially with the same decay length 1/M on both sides of the boundary. The simplifying assumption allows for an analytical treatment of the matching condition [29]. The reduced Hamiltonian for the boundary states conserves their helicity, .

Possibility of Gapless Bound States at a Screw Dislocation in a 3D TI

In a 3D material, a screw dislocation can occur. The Volterra process for the creation of such a defect requires cutting the material along a half plane. Next, the two free surfaces are twisted with a relative displacement along the direction of the plane and glued back, so that the right-hand side is displaced upward and the left-hand side is displaced downward, as shown in Figure 2. It has been shown that, in a 3D TI with inversion symmetry, a screw dislocation of Burgers vector, , harbors two gapless Dirac modes, which propagate helically along the screw axis, but are localized in the radial direction, provided the constraint:

Figure 2.

Volterra process for the screw dislocation (taken from [45]). In the picture, the Burgers vector is one lattice spacing long and points downward.

is fulfilled [34]. Here, M = (ν1G1 + ν2G2 + ν3G3))/2 is a time reversal invariant momentum (TRIM) expressed on the basis of reciprocal lattice vectors (G1, G2, G3). (ν1, ν2, ν3) are the “weak topological invariants”, which contribute in classifying the 3D TI [52]. One starts by deriving the Z2 variables δ(Γ) = ±1 from the parities of the occupied bands. Here, Γ are the time reversal invariant k−vector points in the Brillouin zone. By choosing appropriate gauges for the Bloch functions, it is possible to obtain that a single one, δ(M), takes the opposite sign with respect to the others. In the case of Bi2Se3 discussed above, there is no such non-vanishing M point, and M is set equal to the Γ ≡ (0, 0, 0) point. Hence, the constraint of Equation (60) cannot be fulfilled. As Bi is very close to a band inversion transition at the L point, a good test material that can satisfy the constraint of Equation (60) is the strong TI, Bi1Sb (with x ∼ 0.03), of topological invariants ν0 = 1 and (ν1, ν2, ν3) = (1, 1, 1) [58]. If ~b is one lattice spacing along one of the three directions of G, let us say x^3, the condition is fulfilled.

In the case of the screw dislocation, the Burgers vector is oriented along the defect line, at difference with the edge dislocation, which has the Burgers vector orthogonal to the defect axis. The change of coordinates, describing a screw dislocation with , is:

The coordinates with overlined indexes refer to the local inertial set of coordinates, i.e., the be interpreted as the coordinates before the Volterra process takes place. A former plane existing before the creation of the defect, changes when the defect is created in such a way that the x3 coordinate is displaced by b, when circulating of an angle θ = 2π around the vertical axis. The corresponding tetrads are:

while the inverse tetrads are:

The only non-vanishing component of the torsion is , while the curvature vanishes, as in the case of the edge dislocation. The rotation of the Dirac matrices gives (with ħυ = 1):

Here, is the linearized Hamiltonian in the vicinity of the M point. In the case of Bi1Sb, the Time Reversal Invariant Momentum (TRIM) is one of the three equivalent LTRIMs (which is (1, 1, 1) in the appropriate G vector basis), where the band inversion occurs. The in the approximation, to first order in , has been worked out in [58] and is unitarily equivalent to the Hamiltonian of Equations (58) and (59) [59]. Strictly speaking, when choosing an approximated wave function of the type , where:

is the slowly varying part in cylindrical coordinates, the coordinates are in the flat reference frame, as well as the polar axis (i.e., overlined variables ). The Schrödinger equation for Φ(r, θ) is the Dirac equation with an Aharonov–Bohm flux given by and eigenvalue :

In case , two helical states with gapless linear energy dispersion are bound to the screw dislocation. Indeed, a zero energy solution of Equation (66) can be found, exponentially decaying away from the dislocation core, with the antiperiodic boundary condition:

It follows that the full wavefunction , when expressing in terms of x3, acquires an extra phase , which compensates the antiperiodic boundary condition given above, due to the fact that x3 has to be displaced by b.

We explicitly derive Ψ. According to Equation (66) , the explicit equations to be solved are:

The slowly varying function, Ψ, turns out not to be an eigenstate of the integer angular momentum, mħ. A solution of the system given above is:

The functions, K(κr), are the modified Bessel functions decaying at infinity with a length-scale κ1 to be determined. The only normalizable solution is with n = 0, so that we have:

with the eigenvalue: By posing κ2 = M2, the energy dispersion is gapless and linear: The solution satisfies antiperiodic boundary conditions in the θ variable, as required. Note that the same result is obtained for ℳ (k) = (M + C1k2), provided we choose κ = |M + C1k2|.

The state, which is time reversed with respect to Equation (70), is:

These states have opposite helicity. At zero energy, i.e., with k = 0, Θ commutes with the Hamiltonian, and one can construct eigenstates of zero energy from Equations (70) and (71), which are mutually time reversed. This can be easily done by summing and subtracting Equations (70) and (71) together, in the k → 0 limit:

where:

In the class, AII, of the Altland–Zirnbauer [57] classification, the linear defect in 3D has a corresponding topological invariant, Z2 [60]. The protected gapless modes, which are delocalized along the −axis, can be combined in chiral form, because the left-chiral projector, Π = (1 − Γ)/2, and the right-chiral projector, Π = (1 + Γ)/2, commute with Θ.

Superconductive Proximity in a 3D TI and Bound States at an Axial Vortex

A superconductor in close contact with a normal metal induces Cooper pairing in it by proximity. The bulk states located at the Fermi energy in the metal penetrate the superconductor, even when their energy is below the energy of the superconducting gap, thanks to the Andreev reflection mechanism. The matching at the boundary builds up a superposition of particles and holes in the metal, which are quasiparticles nicknamed “bogoliubons”. A pairing order parameter is induced in the normal metal within a distance from the boundary, which depends on whether transport in the metal is diffusive or “clean” [61]. Proximity does not necessarily imply the opening of a gap in the metal, in particular when the transparency of the boundary is high.

An undoped 3D TI, being a semiconductor, has a Fermi energy located inside the gap separating the bulk bands. Therefore, in principle, bulk quasiparticle states cannot be involved in the proximity. However, the interface of a TI hosts boundary states, whose energy dispersion is the Dirac cone occupying energies within the band gap. The boundary acts as a semimetal of reduced dimension and proximity can take place. However, the properties of the bogoliubons differ from those of the topologically trivial metal, because orbital and spin degrees of freedom are strongly coupled in the Dirac boundary states.

In the presence of a magnetic field, a vortex can be trapped, piercing the heterostructure in which a 3D TI slab is sandwiched between two conventional superconductors. We assume that the S/T I/S heterostructure, forming a slab laying in the x − y plane, with flat boundary surfaces at z = 0, L, is immersed in a magnetic field parallel to the z−direction. The interest for this configuration has been recently growing since the work by Fu and Kane [35], which predicts the possibility of binding a Majorana zero energy quasiparticle state at the vortex in this geometry. This could be the elementary brick for building a completely new architecture for quantum information [62-64]. It is interesting that, while a vortex can be viewed as an axial defect, it is not constrained by features of the lattice symmetry, due to its size. On the other hand, it adds in all cases an orbital angular momentum of π to the bogolubon. This implies that zero energy quasiparticles bound to the vortex core automatically satisfy the required periodic boundary conditions in the azimuthal angle, and the constraint of Equation (60) does not apply. In this section, we report on results obtained about zero energy bound quasiparticle states that can form in the core of the vortex at the boundary with the superconductors, both in the case of even and odd parity proximity [41]. It is unclear whether the pairing induced in the TI is to be expected to have even or odd parity. It has been proposed that, when doped with few percent of Cu, the Bi2Se3 becomes a topological superconductor, undergoing the superconducting phase transition with an odd-parity order parameter [65,66]. For this reason, we consider both possibilities, and we have found very different results in the two cases.

Usually, proximity is described in the Nambu basis, of the kind: ; and theHamiltonian takes the compact Bogoliubov de Gennes (BdG) mean field form:

Superconductive proximity induced by an even parity singlet pairing requires that . The odd-parity triplet pairing satisfies In the case of the two-band 3D TI, an 8 × 8 basis is required for the Hamiltonian of Equation (58), with two types of orbitals, labeled generically with g/u, even/odd by inversion, respectively. The Pauli matrices, σ and s, span the orbital and spin space, respectively, while the Pauli matrices, τ, address the particle-hole sectors. In the even parity, s-wave, singlet case, we define ∆ = 〈ψ〉 = 〈ψ〉 = −〈ψ〉 = −〈ψ〉, assumed to be independent of the orbital type. This gives rise to a pairing mean field Hamiltonian:

In the odd parity case, the pairing is chosen with zero spin projection along the spin quantization axis, which is pinned to the z−axis, normal to the surface of the slab (polar ordering). This choice is motivated by the expectation that the superconductive gap can vanish at the interface (z = 0, L) where the Dirac states are located, if transparency is high. The polar pairing is described by:

where ∆ = 〈ψ〉 = −〈ψ〉 = −〈ψ〉 = 〈ψ〉 is the odd-parity order parameter. The change in signs in the expectation values for ∆ arises from the fact that the operators act on a triplet pair with zero spin projection along z. The choice of zero spin projection allows for a closer comparison with the s-wave case.

We will choose different representation bases for the two symmetries. They are connected to the Nambu basis by unitary transformations, but differ from it. This is convenient, as it can be shown that the induced even and odd-parity proximities give rise to the same matrix form of the model Hamiltonian, when the two different bases are adopted, each for the two different cases. We consider a vortex line of charge q = ±1 piercing the TI, with its axis orthogonal to the boundary surface, and we use cylindrical coordinates oriented along the z−axis.The Hamiltonian is (ħ = 1):

where, outside the vortex core, the Hamiltonian blocks, H, read:

with:

(M > 0 and C1, C2 < 0). We have added the vector potential associated to the vortex, which, far away from the vortex core, takes the form of a pure singular gauge:

(ϕ = hc/2e is the flux unit). The phase factor, e, breaks the TR invariance, which holds when ∆ is real (∆ is purely imaginary).

Equation (77) is appropriate for the even-parity superconducting correlations, when the basis is:

while a unitary transformation, which changes the basis to:

transforms the model to describe the odd-parity pairing, provided ∆ is replaced with ∆.

We now search for zero energy excitations corresponding to quasiparticles bound to the vortex. When proximity induces s-wave , singlet superconductive correlations, Majorana quasiparticles are bound to the vortex. The wavefunction decays exponentially, as exp(−|M||z|), at the boundary, so that two bound states are localized at the interface with each of the two topologically trivial superconductors. The zero energy eigenstates are found by matching solutions inside and outside the vortex core. Its boundary is defined as a circle of radius ξ ∼ ħv /∆. An analytic derivation of the quantum field in the two-band model can be given far away from the vortex core in the limit of μ = 0 (mid-gap Majorana bound states (MBS)) [41]. The two zero energy real fermion fields, localized far apart at the two boundary surfaces of the slab, z ∼ 0+, L, in the inside of the TI, outside the vortex core (r > ξ), take the form:

with λ ∼ |M|. Here, are the modified Bessel functions, so that the excitations are localized in the surface plane. The decay length scale is w1 ∼ ξ. Inside the vortex core, the normalizable solution requires an r−dependent order parameter, ∆(r), together with the corresponding vector potential, A(r), and should be matched with the one given previously at r = ξ.

Let us compare the two Majorana excitations of Equations (83) and (84). They mix both u and g orbitals and involve opposite spin orientations. To understand the result, we note that, in the absence of superconductive proximity, the p and h surface states at the two opposite boundaries, z = 0, L, have exchanged chiralities. In the spinless case, this can be seen also in a simple, reduced model in which chirality is the σ operator:

where ϕ = θ is the phase of the order parameter and ψ is a real spinless fermion. The requirement that the wavefunction decays both radially and inside the TI fixes the eigenvector of σ, with eigenvalue ±1, and the chirality of the state follows. At the two surfaces, the z convergency requires opposite eigenvalues.

Now, let us assume the vortex charge to be q = 1. The Majorana of Equation (83), γ(z ∼ 0), does not have any orbital angular momentum in the z−direction. It includes spin ↑, so that the total angular momentum m = 1/2. The Majorana of Equation (84), γ(z ∼ L), has one unit of orbital angular momentum in the z−direction (because its p/h components have opposite chirality with respect to those of γ(z ∼ 0) and its orbital momentum has to fit with the vorticity q = 1). It includes spin ↓, so that the total angular momentum is again m = 1/2. It looks as if the vorticity of the vortex has been split between the two MBSs having m = 1/2 each. This reminds one of the half quantum vortex of Sr2RuO4 [68-70], as just one spin component acquires an orbital angular momentum. This is remarkable as, notwithstanding the fact that the vorticity adds just orbital angular momentum, MBSs are formed in a spin-dependent combination. On the other hand, in the presence of strong SO coupling, just the total angular momentum projection along matters. This result is consistent with the so-called “thick flux limit” of [71], when the magnetic field penetrates extensively into the bulk below the surface, so that the bulk insulating gap is entirely restored. In that paper, however, the correspondence between the opposite chiralities at the two surfaces is not mentioned.

States localized at a vortex core occur also when Cooper pairing induced by proximity in the Hamiltonian of Equation (77) is odd-parity . The change of basis to the one of Equation (82) changes the nature of the bound states trapped at the vortex singularity. They are not of the kind of the ones given in Equations (83) and (84), because they are not MBSs. Instead, they are zero energy surface Andreev bound states (SABS), decaying in an oscillatory fashion within the TI slab (z > 0) [72]. The state involving just u orbitals takes the form:

Outside the vortex core , both κ and w are complex, so that the function decaying in z has also an exponentially decaying factor . The function of r is a Hankel function, , of complex argument, with an oscillating factor, as well as a decaying exponential factor , where we have chosen the gauge in which ∆ is purely imaginary).

Inside the vortex core , the amplitude is a combination of Hankel functions that converges to zero at the origin (i.e., the point where the order parameter vanishes). In a “hard core” approximation, which does not account for the r-dependence of ∆ and inside the core, the value of is fixed by matching the two solutions of Equation (86) at the core boundary, . Inside the TI slab, the bound state decays exponentially along the vortex axis with a wavevector, . These behaviors qualify the result as a SABS, which, by inspection, is not an MBS, nor an eigenstate of Time Reversal (TR) and decays along the vortex line, with oscillations depending on |∆. It is easy to check that the combinations given here in the asymptotic region out of the vortex core are L-chiral, i.e., they form the combinations: and .

The partner state to the one given in Equation (86) involves the ψ field operators, in i the R—chiral mate combination, i.e.: , . These excitations involve both spin components.

Away from the mid-gap (μ ≠ 0), fermionic excitations can be found, which correspond to circular waves propagating at the interface inward or outward from the vortex singularity and merging into the film by traveling across the slab, along the vortex line, with a radial localization length . These states involve both u and g orbitals and have wavevector [41].

There is no possibility for an MBS to exist, when proximity-induced pairing is odd-parity. The reason can be found in the effective parity of the pairing, which is developed in the TI by proximity. It was shown by Fu and Kane [35] that an even-parity superconductor inducing proximity in a TI acts on Dirac fermions as an effective odd-parity pairing, giving rise to the MBS. Vice versa, an odd-parity proximization can be expected to develop into an effective even-parity pairing, which does not give rise to MBSs.

Summary and Final Remarks

The bulk of 3D TI exhibits an odd number of time reversal invariant wavevectors in the Brillouin zone, in the vicinity of which the Hamiltonian can be expanded in the approximation up to linear terms only. Projection of these 3D Dirac points onto the appropriate surface gives rise to a 2D Dirac cone dispersion for boundary states, which conserve helicity thanks to spin-orbit coupling and are robust with respect to static non-magnetic disorder [73-76]. These properties make boundary states of 3D TI very attractive for applications in spintronics [77]. Conductance of graphene, which is a “weak” 2D TI, due to the presence of the two valleys centered at the two time reversal invariant points, K and K, has been widely analyzed in these years, including weak localization corrections [78]. Topological defects are also being studied in graphene [32,79-82].

In this work, we have shown that an approach similar to the one often used in graphene [83] can be applied to discuss conduction at the boundaries of topological insulators. In the long wavelength limit for Dirac electrons, it is possible to consider a topological defect as a metric distortion in curved space, so that a rotation of the Gamma matrices and a spin connection inherited by the curvature modifies the scattering properties of the carriers [45,84]. Far away from the defect, these modifications can be accounted for by gauge fields. Short-range potentials due to local stress in the lattice can be included in the Born approximation to the Dyson scattering equation. The t-matrix for the scattering of 2D Dirac electrons and edge dislocation has been derived in Section 3.1.

To make a connection with reality, we have interpreted the results as the modeling of one single edge dislocation in graphene for electrons belonging to one valley only, with no inter-valley scattering. In a picture in which defects are considered as non-interacting and very dilute, a Boltzmann approach to conductivity is acceptable, when averaging over their orientation. The contribution to the resistivity is found to be, of course, proportional to the density of defects, but inversely proportional to the density of carriers, in agreement with the measured resistivity. We have also derived the self-energy of the defect connected with the stress involved in the rearrangement of the lattice in Section 3.3. This information could be relevant to estimate the temperature at which the Boltzmann approach breaks down and a defect-mediated phase transition to a disordered phase takes place. As we find a log-dependence of the self-energy on the size of the dislocation, one could surmise that the transition is of Kosterlitz–Thouless type, as in 2D crystal melting. In this case, a temperature scale is determined the stiffness of the lattice, βκ. Assuming , where u is the strain tensor and v is the Fermi velocity [32], the energy of the formation of such a defect would turn out to be of the order of tens of electron volts, leading to various thousand of Kelvin [85]. However, the gauge theory cannot exclude cubic terms, which could drive the phase transition to first order [86,87].

On the same lines, the corresponding topological defects in a 3D TI are the screw dislocations. The odd Dirac point, which is responsible for the material being topologically non-trivial, is the best candidate for hosting the branch point of a screw dislocation. Emphasis is put in our work on the fact that the defect could harbor gapless helical electron states propagating along the dislocation axis. Again, we approach the description of boundary states at a flat surface of a 3D TI and at the defect in the long wavelength limit. We have shown that an analytic form of the bound state wavefunctions can be given easily, because the dislocation acts as an effective flux in the 3D Dirac Hamiltonian. The gapless states exist, provided the constraint on the Burgers vector of Equation (60) is satisfied [34]. TIs, like Bi2Se3 and Bi2Te3, which have the 3D Dirac point at Γ, cannot fulfill the constraint, while the alloy, Bi1Sb, could. The constraint on the orientation of the screw dislocation would influence a statistical approach to the proliferation of defects and thoroughly change the conduction properties, as well as the crystal melting. To our knowledge, this topic has not been addressed yet in the literature.

When a 3D TI is sandwiched between two even-parity superconductors, a vortex piercing the structure can host a zero energy bound state, which is a real fermion field. This is a Majorana bound state. There is great excitement at present for the possibility of revealing Majorana bound states in Josephson junctions between TI, in proximity with superconductors [88] and at vortices [66,89]. In this case, the constraint of Equation (60) does not apply, and the reference to the TRIM, which originates the surface Dirac cone, is unnecessary.

The conserved quantity is total angular momentum along the vortex axis due to SO. The two Majorana at opposite surfaces have a total angular momentum of 1/2 along the vortex axis, and the sum matches the unitary vortex orbital momentum (the “vortex charge”). However, it is remarkable that, because of their opposite chirality, one of them has no orbital angular momentum and spin polarization up, while the other one has spin angular momentum down, to be subtracted from one positive unity of orbital angular momentum. Hence, opposite chiralities imply that the orbital angular momentum of the vortex fixes the spin content of the Majorana fields. It is also interesting that, in the case the induced superconductivity by proximity is of the odd-parity type, the zero energy states are Andreev bound states localized at the free surface, with damped oscillations away from the surface.

An alternative route to realize in a controlled way structures with emerging MBS excitations could be to resort to pertinently designed Josephson junction rings [88] and networks. It has been shown that, by properly designing the network, it is possible to realize topological configurations [90-92], which are expected to enucleate defects hosting MBSs [93,94], similar to the ones predicted at the interfaces between a TI and a superconductor.