Theoretical investigation of the applicability of the Meservey-Tedrow technique to the surface states of topological insulators

The spin polarization of topological surface states is of high interest for possible applications in spintronics. At present, the only technique capable to measure the surface state spin texture is spin and angle resolved photoemission spectroscopy (SARPES). However, values reported by SARPES differed strongly. An established technique to measure the spin polarization of ferromagnetic materials is the so-called Meservey-Tedrow technique, which is based on spin dependent tunneling from a superconducting electrode to a ferromagnet. Here, we theoretically investigate how the Meservey-Tedrow technique can be adapted to topological insulators. We demonstrate that with a specific device geometry it is possible to determine the in-plane component of the spin polarization of topological surface states. More complex device geometries can access the full momentum dependence of the spin polarization. We also show that it is possible to extract the spin-flip scattering rate of surface electrons with the same devices.

The spin polarization of topological surface states is of high interest for possible applications in spintronics. At present, the only technique capable to measure the surface state spin texture is spin and angle resolved photoemission spectroscopy (SARPES). However, values reported by SARPES differed strongly. An established technique to measure the spin polarization of ferromagnetic materials is the so-called Meservey-Tedrow technique, which is based on spin dependent tunneling from a superconducting electrode to a ferromagnet. Here, we theoretically investigate how the Meservey-Tedrow technique can be adapted to topological insulators. We demonstrate that with a specific device geometry it is possible to determine the in-plane component of the spin polarization of topological surface states. More complex device geometries can access the full momentum dependence of the spin polarization. We also show that it is possible to extract the spin-flip scattering rate of surface electrons with the same devices.

I. INTRODUCTION
Since their discovery about a decade ago 1-5 , topological insulators (TIs) have attracted great interest in the field of spintronics [6][7][8][9][10][11][12][13] . This interest originates from the presence of topologically protected surface states in an otherwise insulating bulk gap. Particularly, these surface states have a Dirac cone like dispersion and, due to strong spin orbit coupling, spin and momentum are locked, i.e. electrons propagating in opposite directions possess opposite spin.
Theoretical calculations for three-dimensional TIs like Bi 2 Se 3 predict that the spin is orthogonal to the momentum and lies mainly in the surface plane, with a small out-of-plane component in some materials due to the hexagonal deformation of the Fermi surface for larger wavevectors 14 . While this is well confirmed by spin-and angle-resolved photoemission spectroscopy (SARPES) experiments 15 , the measurement of the degree of spinpolarization is still problematic. While most theoretical calculations yield values around 50% − 65% [16][17][18] , the values reported by SARPES measurements differ strongly between ∼ 45% and 100% [18][19][20][21] . The problem with this technique is that, depending on the photon energy and the photon polarization, the spin of the photoelectrons is different from that of the original electrons in the topological surface states. 18,22 However, a detailed knowledge about the surface state spin-texture is important, because it can be crucial for the efficiency of spintronic devices based on topological insulators.
An established technique to measure the spin polarization of ferromagnetic materials is a method first used by Meservey and Tedrow 23,24 . In this technique spindependent tunnelling from a superconductor to a ferromagnet is used to determine the degree of spin polarization. Meservey, Tedrow, and Fulde showed that in thin superconducting aluminum (Al) films the quasiparticle states split in a strong parallel magnetic field B and the The energy of electrons with parallel and antiparallel spin orientation with respect to the applied magnetic field is shifted in opposite directions with respect to the chemical potential µ (dashed purple and dotted green). In a tunnel junction the total differential conductance is a superposition of these two spectra depending on the polarization of the other material. For an unpolarized material one finds an equal superposition (blue), while for a polarized material the degree of polarization can be obtained from the asymmetry of the peak heights observed in the differential tunnelling conductance.
BCS energy spectra of spin-up and spin-down electrons are thereby shifted by ±µ B B with respect to the original spectrum (see Fig. 1). 25 In junctions made out of Al, an insulating barrier (I), and a third material, this allows to measure the spin polarization of the third material from the conductance of the junction as the spin of tunnelling electrons is conserved. 23 If the third material is a material in which all spins are oriented either parallel or antiparallel to the magnetic field, e.g. a ferromagnet, the spectrum is a simple polarization-dependent superposition of the two shifted spectra like in Fig. 1 The Meservey-Tedrow technique cannot be applied directly to a topological insulator surface, because the spin direction rotates around the Fermi surface. Therefore, the total polarization seen in the tunneling conductance would average out to zero. In order to still allow a determination of the spin polarization we suggest to take advantage of the spin-momentum locking of the topological surface states. The main idea can be understood from the device geometry shown in Fig. 2. The superconducting Al electrode (light blue) is placed on top of the topological insulator surface (dark blue) and a thin insulating barrier (yellow). Two metallic electrodes (grey) are placed to the left and right of the Al electrode with a distance d between them. A tunnelling current is fed through the Al electrode and extracted on the metallic electrodes, which should both be on the same potential. Due to the helical spin texture and the Dirac cone like dispersion of the topological surface states, the mean spin polarization of electrons propagating in one direction is opposite to that of electrons propagating in the opposite direction. As a result, the current extracted at the right metallic electrode will be dominated by one spin direction and the current extracted at the left metallic electrode by the opposite spin direction. From the imbalance of the differential conductance one can then determine the degree of spin polarization of the topological surface states.
In the following we will provide theoretical calculations of the differential conductance of such a device and demonstrate that an extraction of the spin polarization is feasible. 26 In contrast to most ferromagnetic materials the density of states of the topological surface states varies strongly as a function of energy, which needs to be taken into account in the calculations. Also, the tunnelling conductance will depend on the spin-flip mean free path of the electrons in the topological surface states. In the presence of spin-flip scattering an electron that initially propagates in positive x-direction can be scattered and eventually appear at the opposite electrode in negative x-direction. For this reason the distance d between the electrodes should be chosen smaller than a few spin-flip mean free paths to be able to extract the spin polarization. Experimental values for the spin-flip mean free path range between 200 nm in disordered Bi 2 Se 3 27 and a few microns in HgTe 28 . In section III we will discuss this dependence on the spin-flip mean free path and show that it provides the possibility to extract the spinflip mean free path from such devices.

II. MODEL
To provide a realistic calculation of the tunnelling conductance of an Al/I/TI junction, we use a tight binding model for the Bi 2 Se 3 class of materials as given in Ref. 9. This tight binding model has been derived from bandstructure calculations up to third order in momentum k by Liu et al 14 using k · p theory.
In k-space the Hamiltonian can be written in terms of the Dirac Γ matrices Γ 1,2,3,4,5 = (τ 1 ⊗ σ 1 , τ 1 ⊗ σ 2 , τ 2 ⊗ I 2×2 , τ 3 ⊗ I 2×2 , τ 1 ⊗ σ 3 ), where the Pauli matrices τ i and σ i act in orbital and spin space, respectively. It takes into account four dominant bands at the Fermi surface and reads (1) On a bilayer hexagonal lattice the tight-binding parameters can be defined as follows: 9,29 The corresponding model parameters for Bi 2 Se 3 , which we consider here as an example, are derived from Liu et al 14 using the atomic distances a = 4.14Å and c = 28.64 15Å 596eV, R 1 = 0.713eV, and R 2 = −1.597eV. The parameter C 0 adds only a tiny energy shift and will therefore be neglected in the following.
As regards the aluminum electrode, only its density of states near the Fermi level is of importance for the tunnelling current. Thus, for simplicity we describe the superconducting aluminum by Hamiltonian Eq. (1) with only one orbital and parameters C Al The parameter C 0 was chosen such that the center of the Al band fits the Fermi level of the TI. In this way the normal state density of states of the Al band is nearly constant within the bulk gap of the TI. The BCS density of states in the superconducting state is described by the Dynes formula, which accounts for a finite lifetime broadening Γ in the aluminum: 31 Here, E is energy, µ B the electron magnetic moment, B the applied magnetic field, ∆ = 0.35meV the superconducting gap of Al, and Γ = 0.03meV. In order to calculate the transition rate of the junction, we Fouriertransform the Hamiltonian Eq. (1) into real space onto its lattice in z-direction, i.e. in the direction perpendicular to the junction. Periodic boundary conditions are used in the in-plane directions which allows to keep the in-plane momenta k x and k y as good quantum numbers. Using Fermi's golden rule for the transition rate from an initial state |m into a final state |n , the total tunnelling current from the Al film into the TI, which takes into account tunnelling processes at finite temperature in both directions, is given by for a bias voltage U between the Al electrode and the TI.
Here, m and n number the unperturbed eigenstates of the Al film and the TI, respectively. The Fermi function describes the occupation of these states at finite temperature. The insulating barrier is modeled by a tunneling Hamiltonian of the form where d † kx,ky,α.σ creates an electron in orbital α with spin σ in the top layer of the topological insulator and c kx,ky,σ destroys an electron in the bottom layer of the aluminum.
The differential conductance (DC) which is usually measured in experiments is obtained from the derivative of I with respect to U III. CALCULATIONS In this section we derive a method to calculate the polarization of the topological surface states from a given tunnel spectrum. As pointed out above we cannot simply take the differential conductance of an Al/I/TI junction but have to exploit the locking between spin and propagation direction. Concerning our calculations the propagation direction of an electron can be obtained via its group velocity v ∝ ∂E ∂k . We first consider the device geometry shown in Fig. 2. Other geometries will be discussed in Appendix A.
In order to obtain useful approximation formulas for the differential conductance we start from an analytical approximation of the topological surface states, which is valid in the vicinity of the Dirac node. The results of this approximate calculation will be compared with full numerical results below. When we expand Hamiltonian Eq. (1) up to second order in k, we can derive an analytical expression for the four components of the surface state wave function: 9 Here, ± is for the upper and lower Dirac cone, ϕ is the in-plane polar angle of the momentum, −1 ≤ p ≤ 1 is the degree of spin-polarization of the surface states, and the orientation of the spin is given by the phase e −i(ϕ− π 2 ) , i.e. it is always rotated by π 2 with respect to ϕ. The corresponding eigenenergies only depend on the magnitude k = k 2 x + k 2 y of the in-plane momentum 9 (9) These surface states describe an isotropic Dirac cone, i.e. they neglect the hexagonal deformation, which is of third order in k. The position of the Dirac node is at the energy E 0 = − C1M0 M1 = 0.234eV. In our model for the Al electrode we simply have two degenerate eigenstates for each pair of k x and k y . The spatial dependence in z-direction is a superposition of an incoming and a reflected wave and thus given by sin zk z . An appropriate linear combination of these eigenstates then leads to states with a specific spin polarization. For a spin polarization within the x-y surface plane, this linear combination is given by where ϕ Al is the in-plane polar angle of the polarization with respect to the k x -axis and for given momenta k x and k y and energy E. Using Eq. (8) and (10) with z = 1 for the bottom layer of the aluminum, the transfer matrix element of the junction can be calculated where C B is the hopping matrix element of the barrier. This can then be inserted into Eq. (7) multiplied with the shifted BCS density of states. To get the DC with respect to the electrode in positive x-direction we follow Ref. 9 and introduce a function f (ϕ), which gives the probability that an electron initially propagating under an angle of ϕ ends up at this electrode. f (ϕ) depends on the geometry of the device and can also be used to include effects like spin scattering in the TI. For the simple case shown in Fig. 2, where all electrons with a positive group velocity component in x direction end up at that electrode it is simply 9 More complex cases are discussed in Appendix A. In the following we assume that f (ϕ) is an even function, i.e. f (ϕ) = f (−ϕ), which is always satisfied when mirror symmetry with respect to the x-z-plane holds.
Let us choose the applied magnetic field to point in the direction of ϕ Al . Then, the differential conductance Eq. 7 consists of two contributions coming from the electrons with spin either parallel or antiparallel to the magnetic field. For the contribution with parallel spin we find const. Tˆk Here, the DC becomes a product of a ϕ-integral, which depends on the geometry of the device and the relative polarization of the Al film with respect to the TI, and the term G − (T, U ), which contains the densities of states of the two materials. Analogously, we find for electrons with spin oriented antiparallel The total DC of the junction is then Fig. 3 as a function of bias voltage U for T = 0.4K, B = 3T, ϕ Al = π 2 , p = 1, and a chemical potential of µ = 0.2eV. The chemical potential is chosen somewhat below the Dirac node, i.e. at an energy where the hexagonal deformation of the Dirac-cone is small. For comparison, we also show the result of a numerical calculation of the differential conductance Eq. (7). This calculation was based on the full Hamiltonian Eq. (1) on a hexagonal lattice with 50 layers along z ([001]). The momenta k x and k y were uniformly distributed over the first Brillouin zone with a discretization of 2 √ 3 2π N and N = 48000, corresponding to a sample width of about 20µm. In spite of the simplifications made in the analytical approximations, there is only a small deviation of the two DC curves. Please note that one can discern four peaks in the DC curves, even though the surface states are fully polarized here (p = 1). The physical reason for this is the spin texture of the topological surface states: even though the current from the Al electrode to the M electrode in positive x-direction is dominated by spin up electrons, there is a small but finite probability that a spin down electron can tunnel into a topological surface state with positive group velocity v x in x-direction.

This is shown in
In the following we demonstrate that the polarization p of the topological surface states can be reliably inferred from such differential conductance curves. In Ref. 24 Tedrow and Meservey derived a formula to extract the polarization of ferromagnets from the four peak heights of the DC curve. Here, we adapt that formula to the present case and show that the polarization p can be obtained from it within a good approximation. If the density of states of the TI and Al film were constant as a function of energy except for the BCS density of states, we could express G − (T, U ) and G + (T, U ) in terms of the unsplit DC G (T, U ). 24 Defining F ± =´π −π dϕf (ϕ) (1 ± p sin ϕ Al cos ϕ) the differential conductance g at some arbitrary bias voltage x = U − µ e could be written is the splitting of spin up and spin down densities of states. One now evaluates g(x) at four bias voltages ±x and ±(x − 2b), where x ∼ ∆+µ B B e is chosen at the peak position of the outermost peak. The values of the four conductances g 1 to g 4 (from left to right) are then given as Assuming that G (T, x) is a symmetric function of x, this set of equations can be solved and the polarization p can be obtained from the four conductances g 1 to g 4 leading to the formula Here, the factor accounts for the geometry of the device. For f (ϕ) given by Eq. (13) we have γ = π 2 . Since we required that we can write the splitted DC curves Eq. (16) and (17) in terms of the unsplit curve G (T, U ) in the derivation of Eq. (25), Eq. (25) becomes inaccurate if the densities of states of spin-up and spin-down electrons in the superconductor are not the same function of energy. This is the case if spin-orbit scattering in the superconductor is present. However, concerning thin aluminum films, spin-orbit scattering is small and the deviation from Eq. (25) should therefore be negligible. 24 In the case of significant spin-orbit scattering in the superconductor, one can still fit the DC curves to obtain the polarization, if the energy dependence of the separate spin densities of states are known.
In the above calculation we assumed the normal state densities of states to be constant as a function of energy. Since this is not the case for TIs we have to remove this energy dependence from the measured DC in order to calculate the polarization p. This can be done by fitting the DC curve with a low order polynomial h (U ), after cutting out the part with a high influence of the BCS density of states, as shown in Fig. 3 (dashed line). Afterwards one then multiplies G (T, U, ϕ Al ) with 1 h(U ) and analyzes this corrected DC curve. The corrected DC curves are shown in Fig. 4. While there was a small difference between our analytical and numerical calculations in Fig. 3, the curves now agree very well and thereby support the general validity of this procedure and of Eq. (25). Because of the broadening of the BCS density of states, the exact peak maxima in general do not fulfill the symmetry requirement around µ. To reduce the error on the polarization value, the positions for the g i have to be chosen such that the slope of the DC at these positions is small. Here, we choose them such that the largest outer and largest inner peak, which are either g 1 and g 3 or g 4 and g 2 , are met exactly. The other two positions are then automatically given by the symmetry requirement around µ. Applying this to the analytical DC curve in Fig. 4 (black crosses) yields p ≈ 0.9996, consistent with the model parameter of p = 1. The deviation from the absolute value of the actual polarization of 100% is only 0.04 percentage points. It is however the inverse of the actual polarization for k x < 0 at the lower Dirac cone. The numerical DC curve yields p ≈ 1.0038, with only a slightly larger deviation. Applying the same formula, with the same γ factor, to the DC of the opposite metallic electrode, we find p ≈ −1.0038, i.e. the same value with opposite sign.
Equation (25) is based on the topological surface states given by Eq. (8) where the spin is always perpendicular to the in-plane momentum. As this may not always be the case it is useful to rewrite Eq. (25) such that it does not depend on the direction of the magnetic field but instead on the angular difference ∆ϕ Al of magnetic field and surface state polarization. When we account for the counterclockwise rotation of the spin in the lower Dirac cone as well, the new model independent formula reads Provided that the γ factor is calculated for a function f (ϕ) that gives the probability with respect to the electrode in positive direction, it always yields the correct polarization value for surface electrons propagating to the considered electrode. In practice the magnetic field should be oriented parallel to the polarization of the TI (∆ϕ Al = 0), i.e. such that the DC becomes maximal. By this, on the one hand we get the orientation of the surface state spin and on the other hand maxi- When we apply the above scheme to analytical DC curves for various polarization values, the difference ∆p = p ex − p between extracted polarization and actual polarization is always of the same order of magnitude, as shown in Fig. 5. There also seems to be no correlation between ∆p and p.
In Fig. 6 we show the polarization extracted from DCs of the full Hamiltonian as a function of chemical potential µ. We see no systematic deviation from p = −1 . This is not intuitively clear considering the energy dependent hexagonal deformation of the Fermi surface along with an out of plane tilt of the spin polarization. To analyze the influence of a hexagonal deformation we therefore assume a correction of the surface state approximation, where the spin is tilted out of the surface plane with alternating sign with a threefold period: Here, 0 ≤ q ≤ 1 is the absolute value of the out of plane polarization. The spin expectation values of this improved approximation are given by n x = ψ ± |Σ x | ψ ± = ± 1 − q 2 cos 2 3ϕ p sin ϕ (29) n y = ψ ± |Σ y | ψ ± = ∓ 1 − q 2 cos 2 3ϕ p cos ϕ (30) n z = ψ ± |Σ z | ψ ± = ±q cos 3ϕ (31) With this new approximation, the ratio of the DC peaks remains basically the same, with a new factor γ =´π −π dϕf (ϕ) π −π dϕf (ϕ) 1 − q 2 cos 2 3ϕ cos ϕ which depends on q. So it becomes possible to calculate the deviation from Eq. (27) for specified values of q and estimate possible errors. For the case shown in Fig. 6 a maximal out-of-plane polarization of q ≈ 0.17 at the lower edge of the bulk gap yields a relative deviation γ γ ≈ 1.008 for DC measurements along the x-direction. This is only of the order of the measurement accuracy. Along the y-direction (replace cos 3ϕ with cos 3ϕ − π 2 ), this ratio is even slightly smaller γ γ ≈ 1.007, because the out-of-plane polarization is zero along the y-axis and states close to the axis contribute strongest to the DC. In the extreme case of f (ϕ) = δ (ϕ), which would be valid for a two dimensional device, these deviations become somewhat larger, but are still small. At the lower end of the bulk gap, a deviation of γ γ ≈ 1.015 could be expected. Note, however that in this case, the out-ofplane polarization strongly depends on the measurement direction. Along some crystal axes it reaches the maximum value, while for others it completely vanishes. So, by varying the measurement direction, one can get rid of the out-of-plane spin component in order to access the in-plane component. It is however unlikely that the outof-plane component can be determined from how γ varies as a function of ϕ. A device that is capable of measuring along certain crystal axes is presented in appendix A 3.
From Eq. (27) one sees that the polarization p depends linearly on the geometrical factor γ. Instead of measuring p for a given value of γ, one can alternatively also measure γ for a given p. This grants access to another physical variable of the topological insulator: the spinflip mean free path, which is crucial to applications in spintronics. The spin-flip mean free path ξ is the average path after which an electron has lost information on its original spin and hence also on its propagation direction. If there is spin-scattering in the TI, the measured apparent spin-polarization of the surface states will depend on the length of the path between the Al and M electrode and the spin-flip mean free path ξ thereby enters into the distance dependent γ factor. For the device in Fig. 2 we derived the probability distribution f (ϕ) accounting for a finite spin-flip mean free path ξ in Appendix A 1. From this expression γ can be calculated numerically. For a known polarization p one can then simply calculate ξ by fitting γ to Eq. (27). If p is unknown, it is still possible to calculate ξ from how γ varies with distance d. However, this is more inaccurate as it requires multiple devices with different distances d. To reduce errors, the different devices should at least be prepared on the same sample, since ξ and p may depend on the quality of the TI material. When we define the experimentally accessible quantity p is of the form and identical for devices with different distances d i . Then, we can solve for the ratio of two γ i : The spin-flip mean free path ξ can be fitted to this ratio as shown in Fig. 7. In order to get an accurate result, d or, in the case of an unknown polarization, the distances d 1 , d 2 and |d 2 − d 1 | should be in the same order of magnitude as ξ.

IV. SUMMARY AND CONCLUSIONS
We have studied a generalization of the Meservey-Tedrow method to topological insulators. Starting from an analytical approximation of TI surface states, we showed how quantum tunneling from a superconductor into these surface states can be used to measure their spin polarization. In contrast to the application to ferromagnets, one has to measure the tunneling current with respect to different spatial directions and take into account that the density of states in TI surface states is strongly energy dependent. Considering these aspects we derived formulas that allow easy calculation of the in-plane spin polarization from measured tunneling spectra, where the geometry of the device enters as a single factor. As there is no chemical potential dependence of the in-plane polarization, the tunneling spectra seem to be insensitive to the out-of-plane polarization of the surface states, at least for the device geometries considered here. When spin-flip scattering is included in the calculation of the geometrical factor, it can be measured as well and hence, if measured prior to the spin polarization, can increase the accuracy of the calculated spin polarization. where Θ is the Heaviside step function, l and h are the length and height of the two small electrodes and d the distance between them. For given values of d, l and h, γ can be calculated numerically. If x h is large for all x, Eq. (A4) may be approximated by f (ϕ) = δ (ϕ), i.e. γ = 1. The DC then effectively becomes that of a two dimensional device, where all electrons reaching the metallic electrode have the same momentum direction and spin. In principle, using multiple devices with different orientations of the electrodes with respect to the lattice of the TI the k-dependent spin structure of the surface states can be studied. It is also possible to split the U-shaped electrode into multiple parts. In this way measurements in multiple directions can be done with a single device.