Power-law Temperature Dependence of the Penetration Depth in a Topological Superconductor due to Surface States

TL;DR

This paper investigates how surface states in a 3D topological superconductor cause a distinctive T^3 power-law temperature dependence in the magnetic penetration depth, offering a potential signature of topological superconductivity.

Contribution

It introduces a theoretical model accounting for surface Majorana states' paramagnetic current, predicting a unique T^3 dependence in the penetration depth at low temperatures.

Findings

Surface states induce a T^3 temperature dependence.

Bulk remains fully gapped, distinguishing from nodal superconductors.

Power-law behavior can serve as an indicator of topological superconductivity.

Abstract

We study the temperature dependence of the magnetic penetration depth in a 3D topological superconductor (TSC), incorporating the paramagnetic current due to the surface states. A TSC is predicted to host a gapless 2D surface Majorana fluid. In addition to the bulk-dominated London response, we identify a $T^3$ power-law-in-temperature contribution from the surface, valid in the low-temperature limit. Our system is fully gapped in the bulk, and should be compared to bulk nodal superconductivity, which also exhibits power-law behavior. Power-law temperature dependence of the penetration depth can be one indicator of topological superconductivity.

Equations

$${\nabla }^2\text{A}(z)=-\frac{4\pi }{c}\text{J}(z),$$

$$\text{J}(z)={\text{J}}_{\mathsf{ext}}(z)-\frac{1}{c}{\int }_0^{\infty }d{z}^{\prime }\Pi (z,z^{\prime })\text{A}(z^{\prime }).$$

$$-2B_0=$$

$$+\frac{4\pi }{c^2}\int \frac{d{Q}_z}{2\pi }{\tilde{\Pi }}_{1,R}^{xx}(0,0,0;q_z,-Q_z)\tilde{A}(Q_z),$$

$${\tilde{\Pi }}_{1,R}^{xx}(\Omega =0,q_x=0,q_y=0;q_z,-Q_z)$$

$$=-i\int dtdxdy{\int }_0^{\infty }d{z}_{1}d{z}_2{e}^{-i{q}_z{z}_1+i{Q}_z{z}_2}$$

$$\times ⟨\left[J_1^x(t,x,y,z_1),J_1^x(0,0,0,z_2)\right]⟩\theta (t)$$

$$H=\frac{1}{2}{\int }_{\text{k}}{\chi }^{†}(\text{k})\hat{h}(\text{k})\chi (\text{k}),$$

$$\hat{h}(\text{k})={\tilde{\epsilon }}_k{\hat{\sigma }}^3+\Delta \hat{\text{s}}\cdot \text{k}{\hat{\sigma }}^2,{\tilde{\epsilon }}_k=\frac{k^2}{2m}-\mu ,$$

$${\tilde{\Pi }}_{1,R}^{xx}(\Omega =0,\text{q}=0)=-\frac{\beta }{6}{\left(\frac{e}{m}\right)}^2{\int }_{\text{k}}{k}^2{\mathrm{sech}}^2\left(\frac{\beta E_k}{2}\right),$$

$${\psi }_{\text{k}}^{\mathsf{s}}(z)=\frac{e^{-z/l_{\mathsf{coh}}}}{\sqrt{N_{\text{k}}^{\mathsf{s}}}}\sin \left(z\sqrt{{\text{lambdabar}}_F^{-2}-l_{\mathsf{coh}}^{-2}-k^2}\right)|{\psi }_{\text{k}}^0⟩,$$

$$H_{\mathsf{s}}={\int }_{\text{r}}\left\{\frac{1}{2}{\eta }^{†}(\text{r})\left[-i\Delta (\hat{s}\wedge \nabla )\right]\eta (\text{r})-\frac{1}{c}\text{A}\cdot \text{J}\right\},$$

$$\text{J}(\text{r})=\frac{e}{4m}{\eta }^{†}(\text{r})i\overleftrightarrow{\nabla }\eta (\text{r}).$$

$${\tilde{\Pi }}_{1,R,\mathsf{s},\mathsf{s}}^{xx}(0,0,0;q_z,-Q_z)$$

$$\simeq -\mathcal{C}{\left(\frac{e}{m}\right)}^2\frac{(k_{B}T{)}^3}{{\Delta }^4}\Theta (0,q_z)\Theta (0,-Q_z),$$

$$\Theta (k,q_z)\equiv {\int }_0^{\infty }dz{e}^{-i{q}_{z}z}{\psi }_{\text{k}}^{\mathsf{s}†}(z){\psi }_{\text{k}}^{\mathsf{s}}(z)$$

$$B_{y}(z)=B_{0}\left\{\operatorname{e}^{-z/\lambda_{L}^{{\scriptscriptstyle{(0)}}}}-\varrho(T)\bigg{[}\partial_{z}G(z)\bigg{]}G(0)\right\},$$

$$\varrho (T)\equiv \left[2^4{3}^3\zeta (3)\pi \right]{\left({\lambda }_L^{(0)}\right)}^6{l}_{\mathsf{coh}}^{-1}{\text{lambdabar}}_F^{-4}{t}^3$$

$$\lambda (T)\equiv \frac{1}{B_0}{\int }_0^{\infty }dz{B}_y(z)={\lambda }_L^{(0)}+\varrho (T){\left[G(0)\right]}^2.$$

$$\!\!\!\!\!\begin{aligned} G_{\mathrm{I}}(z)\simeq&\,\frac{\lambdabar_{F}^{2}}{2\left(\lambda_{L}^{{\scriptscriptstyle{(0)}}}\right)^{2}\left[4\left(\lambda_{L}^{{\scriptscriptstyle{(0)}}}\right)^{2}+\lambdabar_{F}^{2}\right]}\bigg{\{}-\left(\lambda_{L}^{{\scriptscriptstyle{(0)}}}\right)^{2}e^{-z/\lambda_{L}^{{\scriptscriptstyle{(0)}}}}\\ &\,+\left[2\left(\lambda_{L}^{{\scriptscriptstyle{(0)}}}\right)^{2}+\lambdabar_{F}^{2}\sin^{2}\left(\frac{z}{\lambdabar_{F}}\right)\right]e^{-2z/l_{\mathsf{coh}}}\bigg{\}}.\end{aligned}\!\!\!\!\!$$

$$G_{\mathrm{II}}(z)\simeq \frac{{\text{lambdabar}}_F^2{l}_{\mathsf{coh}}}{16{\left({\lambda }_L^{(0)}\right)}^4}\left[{\lambda }_L^{(0)}{e}^{-z/{\lambda }_L^{(0)}}-l_{\mathsf{coh}}{e}^{-2z/l_{\mathsf{coh}}}\right].$$

$$\hat{h}(\text{k},-i{\partial }_z){\psi }_{\epsilon }(z)=\epsilon {\psi }_{\epsilon }(z),$$

$${\psi }_{\text{k}}^{\mathsf{s}}(z)=\frac{1}{\sqrt{N_{\text{k}}^{\mathsf{s}}}}{e}^{-z/l_{\mathsf{coh}}}\sin \left(z\sqrt{{\text{lambdabar}}_F^{-2}-l_{\mathsf{coh}}^{-2}-k^2}\right)\left[\begin{array}{c}1\\ -i{e}^{i{\varphi }_{\text{k}}}\\ 1\\ i{e}^{i{\varphi }_{\text{k}}}\end{array}\right],$$

$$N_{\mathrm{k}}^{\mathsf{s}}=l_{\mathsf{coh}}\frac{\left({\text{lambdabar}}_F^{-2}-l_{\mathsf{coh}}^{-2}-k^2\right)}{\left({\text{lambdabar}}_F^{-2}-k^2\right)}.$$

$$-{\hat{M}}_{\mathsf{P}}^{-1}{\hat{h}}^{\mathsf{T}}(-\text{k}){\hat{M}}_{\mathsf{P}}$$

$${\hat{M}}_{\mathsf{T}}^{-1}{\hat{h}}^{\ast }(-\text{k}){\hat{M}}_{\mathsf{T}}$$

$$-{\hat{M}}_{\mathsf{S}}\hat{h}(\text{k}){\hat{M}}_{\mathsf{S}}$$

$${\psi }_{\text{k}}^{\mathsf{s}†}(z_1){\hat{\sigma }}^1{\psi }_{\text{k}}^{\mathsf{s}}(z_2)={\psi }_{\text{k}}^{\mathsf{s}†}(z_1){\hat{s}}^3{\psi }_{\text{k}}^{\mathsf{s}}(z_2)=0.$$

$$S$$

$${\Pi }_1^{xx}(\tau ,\text{r},{\text{r}}^{\prime })=$$

$$=$$

$${\tilde{\Pi }}_1^{xx}(i{\Omega }_n=0,\text{q}=0)=$$

$$=$$

$${\Pi }_1^{xx}(\tau ,\text{r}-{\text{r}}^{\prime },z,z^{\prime })=-\frac{1}{2}{\left(\frac{e}{m}\right)}^2{T}^2\sum _{{\omega }_1,{\omega }_2}{\int }_{{\text{k}}_1,{\text{k}}_2}{e}^{-i({\omega }_1-{\omega }_2)\tau +i({\text{k}}_1-{\text{k}}_2)\cdot (\text{r}-{\text{r}}^{\prime })}(k_1^x)(k_2^x)\mathsf{Tr}\left[\hat{G}(i{\omega }_1,{\text{k}}_1;z,z^{\prime })\hat{G}(i{\omega }_2,{\text{k}}_2;z^{\prime },z)\right],$$

Full text

A decade after the widespread infiltration of topology into quantum materials research, the search for electronically correlated topological phases beyond the fractional quantum Hall effect remains an urgent, but still largely unfulfilled quest. Topological superconductivity Qi_Zhang_Review ; Sato_Ando_Review is sought as a platform for Majorana fermion zero modes Jason_Review and topological quantum computation Q_computer1 . Majorana fermions could be detected by various means, including tunneling spectroscopy MF_STM1 ; MF_STM2 ; MF_STM3 ; MF_STM4 ; MF_STM5 , the Josephson effect 4pi_JE1 ; 4pi_JE2 ; MF_STM1 , as well as spin and optical responses spin_resp_TSC ; optical_resp_TSC .

Only a handful of materials have emerged as bulk topological superconductor (TSC) candidates. Topological superconductivity has long been suspected in Sr2RuO4, although a consensus on chiral $p$-wave order is yet to be reached Sato_Ando_Review ; Mackenzie-Maeno03 ; Kallin12 ; Kallin16 . Time-reversal invariant TSCs could serve as solid-state analogs of the topological superfluid phase in liquid helium (3He-$B$) Ludwig08 ; Roy08 ; Volovik09 ; Qi09 ; Sato09 ; Sato10 ; Machida16 . In the absence of a magnetic field, the predicted hallmark of such a bulk TSC is a gapless, two-dimensional (2D) Majorana fermion surface fluid. It has been argued that the odd-parity “$\hat{\bm{\mathrm{s}}}\cdot\bm{\mathrm{k}}$” pairing of 3He-$B$ VolovikBook ; BernevigHughesBook could naturally arise in doped Dirac semimetals or topological insulators CuBiSe_pairing ; CuBiSe_Anisopairing ; Wray10 ; here $\hat{\bm{\mathrm{s}}}$ and $\bm{\mathrm{k}}$ respectively denote the spin-operator and momentum vectors. Alternately, doped Weyl semimetals have been shown to be natural platforms for topological superconductivity Hosur14 .

There is now substantial experimental evidence for nematicity in the superconducting doped topological insulators (Cu,Nb)xBi2Se3 Zheng16 ; Yonezawa17 ; Qiu17 ; Tao18 (see also Wen18 ), possibly indicative of odd-parity pairing. However, gapless Majorana fermions have not been conclusively detected Ando11 ; Chu13 ; Stroscio13 ; Tao18 . Recently, signatures consistent with topological superconductivity were also found in doped $\beta$-PdBi2, but a conclusive detection remains elusive Kolapo18 . In CuxBi2Se3, only a small percentage of the exposed crystal surface was found to exhibit signatures of superconductivity in STM Tao18 , highlighting the possibility that in inhomogeneous TSCs, there is no guarantee that Majorana fermions will appear at the physical surface of the sample.

It is therefore natural to seek global probes of topological superconductivity. Meissner effect penetration depth measurements in NbxBi2Se3 Smylie16 ; Smylie17 and in the half-Heusler compounds YPtBi Paglione18 , YPdBi and TbPdBi Spinu18 exhibit power-law temperature suppression, which is interpreted as evidence for non-$s$-wave, bulk nodal superconductivity Tinkham . In the case of NbxBi2Se3, the results were interpreted as indicative of nodal odd-parity bulk pairing Smylie16 ; Smylie17 , while the YPtBi results were attributed to either an exotic nodal-line “septet” pairing scenario Paglione18 ; Brydon16 , or $d+id$ bulk Weyl superconductivity Brydon16 , smeared by disorder Roy19 .

In this Letter, we show that power-law temperature ($T$) dependence of the penetration depth $\lambda(T)$ can arise at arbitrarily low $T$ in a TSC with a fully gapped bulk, due to the paramagnetic response of the gapless Majorana surface fluid. This requires a nontrivial calculation employing a TSC model with a physical surface-vacuum boundary, and the result involves the convolution of the surface paramagnetic and bulk-dominated diamagnetic responses. It cannot be obtained from the surface Hamiltonian alone. Thus, the observation of non-exponential behavior in $\lambda(T)$ does not necessarily indicate bulk nodal pairing, and can be one diagnostic for screening possible TSCs. We also show that the magnetic field in a fully gapped TSC is sensitive to the spatial profile of the surface states. This can be used to differentiate gapless states from the surface versus those from the bulk. By contrast, STM is sensitive only to the sample surface density of states, while the finite energy resolution of ARPES might prevent the direct detection of Majorana fermions in a low-temperature superconductor. While power-law temperature dependence in the penetration depth can arise from multiple mechanisms Giannetta06 ; Cooper97 , it can serve as one possible indicator for bulk topological superconductivity. In general, one wants as many independent tests as possible in order to identify TSCs.

We consider “minimal” TSCs in class DIII with winding number $\nu=1$, possessing a single surface Majorana cone. We show that the leading correction to the London response due to the presence of Majorana surface states scales as $T^{3}$. The same temperature dependence is predicted to arise in the suppression of the mass flow of superfluid 3He-$B$ through a channel, due to surface currents Wu13 . We calculate explicitly the magnetic field profile inside the slab, which incorporates new features introduced by the surface states. For type I TSCs, the field penetrates much deeper than the London depth into the bulk, with the scale set by the coherence length. For type II TSCs the field is modulated in a shallow region near the surface, and then decays at deeper depths according to the London length, but with an enhanced field amplitude.

Model.—We consider a superconducting slab filling up the $z>0$ half-space, with an external magnetic field $\textbf{B}_{\mathsf{ext}}=B_{0}\,\hat{y}$ as shown in Fig. 1. We assume the field is weak enough such that $B_{0}<B_{c_{1}}$ where $B_{c_{1}}$ is the lower critical field of a type II superconductor. Under this assumption, the effect of vortex formation can be neglected and a linear response treatment is valid. The total current $\bm{\mathrm{J}}$ and the vector potential $\bm{\mathrm{A}}$ satisfy the static Maxwell’s equation

[TABLE]

where $\textbf{A}(z)=A(z)\,\hat{x}$, $\bm{\mathrm{B}}(z)=\partial_{z}A(z)\,\hat{y}$, $\textbf{J}(z)=J(z)\,\hat{x}$, and

[TABLE]

Here, $\textbf{J}_{\mathsf{ext}}(z)$ is a fictitious current generating the external magnetic field, and $\Pi(z,z^{\prime})$ is the current-current correlation function capturing the linear response of the TSC. The above equations can be expressed as

[TABLE]

where $\lambda_{L}^{{\scriptscriptstyle{(0)}}}$ is the London penetration depth and

[TABLE]

is the retarded paramagnetic current-current correlation function. Here, $J_{1}^{x}(t,x,y,z)$ is the paramagnetic current flowing along the $x$ direction and $\theta(t)$ is the Heaviside step function. The first term on the right-hand-side of Eq. (3) represents the diamagnetic London response, while the second term is the paramagnetic response from both the bulk and surface states. The above framework is general and the magnetic field in the slab is determined once the current-current correlation function is specified. In what follows we consider a clean system. Weak nonmagnetic disorder is not expected to modify the low-temperature response of the surface response (it is strongly irrelevant Ghorashi19 ) or of the fully gapped bulk.

Although the low-temperature, $T^{3}$-dependence of the penetration depth derived below depends only on the low-energy dispersion of the 2D Majorana surface fluid, here we consider a microscopic model for both the bulk and surface modes of the TSC in order to completely specify the problem. Solid state models analogous to 3He provide a fertile playground to study topologically nontrivial superconductivity Sato_Ando_Review ; VolovikBook ; Volovik09 . “Solid-state 3He-$A$” would correspond to a Weyl superconductor, which has nodal Weyl points in the bulk that connect to a surface Majorana arc. In this work, we consider “solid-state 3He-$B$” VolovikBook ; BernevigHughesBook , with isotropic $p$-wave pairing of spin-1/2 electrons, represented by the following Bogoliubov-de Gennes Hamiltonian

[TABLE]

where

[TABLE]

and where $\int_{\textbf{k}}\equiv\int\frac{d^{3}\textbf{k}}{(2\pi)^{3}},$ $\hat{\textbf{s}}$ and $\bm{\hat{\sigma}}$ respectively denote Pauli matrices acting in the spin and particle-hole spaces, and $\chi(\textbf{k})\equiv\left[c(\textbf{k}),\hat{s}^{2}\left[c^{\dagger}(-\textbf{k})\right]^{\mathsf{T}}\right]$ is the four-component Balian-Werthammer spinor BW63 . The latter satisfies the reality (“Majorana”) condition $\chi^{\dagger}(\textbf{k})=i\chi^{\mathsf{T}}(-\textbf{k})\,\hat{M}_{\mathsf{P}},$ where $\hat{M}_{\mathsf{P}}=\hat{s}^{2}\hat{\sigma}^{2}$ defines particle-hole symmetry for $\hat{h}(\bm{\mathrm{k}})$. In Eq. (6), $\mu$ is the chemical potential and $\Delta$ is the superconducting order parameter amplitude. For $\mu>0$ the above Hamiltonian has winding number $\nu=1$. The retarded paramagnetic current-current correlation function due to the bulk is

[TABLE]

where $E_{k}=\sqrt{\tilde{\varepsilon}_{k}^{2}+\Delta^{2}k^{2}}$ is the eigenenergy of $\hat{h}(\textbf{k})$ and $\beta=T^{-1}$ is the inverse temperature. In the low-temperature limit, the bulk paramagnetic response is exponentially suppressed and the diamagnetic London response dominates. The London depth is given by $\lambda_{L}^{{\scriptscriptstyle{(0)}}}=\sqrt{mc^{2}/(4\pi e^{2}n)},$ where $n$ is the charge number density.

To consider the response from the surface, we replace $k_{z}\rightarrow-i\partial_{z}$ in $\hat{h}(\textbf{k})$ [Eq. (6)] and solve for the Majorana surface states $\psi_{\textbf{k}}^{\mathsf{s}}(z)$ with eigenenergies $\pm\Delta|\textbf{k}|$. Here and in what follows, $\textbf{k}=(k_{x},k_{y})$ specifies the momentum transverse to the interface. With hard wall boundary conditions at $z=0$, we obtain surface wave functions of the form SM

[TABLE]

where $N_{k}^{\mathsf{s}}$ is a normalization constant and $\ket{\psi_{\textbf{k}}^{0}}$ is a spin-momentum-locked spinor in (spin)$\otimes$(particle-hole) space. The two length scales in Eq. (8) are the (reduced) Fermi wavelength $\lambdabar_{F}\equiv 1/\sqrt{2m\mu}$ and the coherence length $l_{\mathsf{coh}}\equiv 1/m\Delta$.

We can see how the magnetic field couples to the surface fluid by incorporating a vector potential $\bm{\mathrm{A}}$ into Eq. (5), and then projecting onto the low-energy surface states. The result is

[TABLE]

where $\eta=\eta_{s}$ is the two-component surface Majorana fermion operator ($s\in\{\uparrow,\downarrow\}$), $\hat{s}\wedge\bm{\nabla}\equiv\hat{s}^{1}\partial_{y}-\hat{s}^{2}\partial_{x},$ and the surface paramagnetic current operator

[TABLE]

Here $\textbf{r}=(x,y)$ and $\bm{\overleftrightarrow{\nabla}}\equiv\bm{\overrightarrow{\nabla}}-\bm{\overleftarrow{\nabla}}$. Eq. (9) assumes that the field $\bm{\mathrm{B}}(z)$ and the vector potential $\bm{\mathrm{A}}(z)$ (in London gauge) both reside in the $(x,y)$ plane. Zeeman coupling to a nonzero component $B_{z}$ would induce a Majorana mass, gapping out the surface fluid Qi_Zhang_Review , but this is prevented by bulk Meissner screening. On the other hand, a very strong in-plane field $A_{x,y}$ could “overtilt” the surface Majorana cone, creating a surface Fermi pocket. The latter should be included in the diamagnetic current Guinea19 , but we exclude this situation here by restricting to linear response.

In the low-temperature limit, the paramagnetic current-current correlation function due to the surface state fluid evaluates to SM

[TABLE]

$\mathcal{C}\equiv\left[3^{2}\zeta(3)/(2^{3}\pi)\right]$, $\zeta(z)$ is the Riemann Zeta function, and where

[TABLE]

is the Fourier transformed probability density of the surface states along the $z$-direction. Unlike the paramagnetic response from the bulk, the one from the surface [Eq. (11)] has a non-trivial $T^{3}$ power-law dependence at low temperature. Two factors of temperature arise from the form of the paramagnetic current operator (a derivative) in Eq. (10), while the third stems from the surface density of states of the Majorana fluid. One should also consider the surface-bulk cross terms when evaluating the paramagnetic current-current correlator appearing in Eq. (3). However, these cross terms exhibit higher-power temperature-dependence at low $T$, and are thus subleading SM . We neglect these surface-bulk contributions in the following.

Results.—Taking only the diamagnetic and surface paramagnetic responses into account, which is valid at low temperature as discussed above, we can formally invert the integral equation Eq. (3) and solve for the vector potential (and hence the magnetic field) profile inside the slab. To leading order in temperature, the final result is SM

[TABLE]

where

[TABLE]

is a temperature-dependent length, with $t\equiv k_{B}T/\Delta k_{F}$ being the dimensionless temperature. Here, $\Delta k_{F}$ is the energy gap of the $p$-wave TSC. In Eq. (13), the function $G(z)$ is a temperature-independent, real-valued function encoding the convolution of the bulk and surface responses. It is a dimensionless function only of $z$ and of the three lengths $\{\lambdabar_{F},\lambda_{L}^{{\scriptscriptstyle{(0)}}},l_{\mathsf{coh}}\}$. $G(z)$ emerges when we invert the integral equation and Fourier transform the quantities back to real space SM .

The first term in Eq. (13) describes the Meissner screening due to the diamagnetic London response, while the second term is the correction due to the Majorana surface fluid. The correction term depends on the length $\varrho(T)$ that encodes the $T^{3}$ dependence, while its spatial dependence is captured by $G(z)$. The function $G(z)$ decays exponentially for large $z$; its spatial extent is governed by the maximum of $\{\lambda_{L}^{{\scriptscriptstyle{(0)}}},l_{\mathsf{coh}}/2\}$, assuming that $\lambdabar_{F}$ is the shortest scale. This leads to different qualitative type I and II behaviors. Nevertheless, to characterize the overall spatial extent of the magnetic field, we can define the effective penetration depth of the system via Tinkham

[TABLE]

The second term in Eq. (15) is the change of the penetration depth due to the surface states. This term is always positive, meaning that the magnetic field can penetrate deeper into the slab for any $T>0$, due to the surface Majorana fluid. It is instructive to roughly estimate the order of magnitude for such correction. For CuxBi2Se3, which is in the extreme type II regime, we substitute typical experimental data $l_{\mathsf{coh}}\sim 10$ nm, $\lambdabar_{F}\sim 1$ nm, and $\lambda_{L}^{{\scriptscriptstyle{(0)}}}\sim 1$ $\mu$m Ando12 ; Tao18 to obtain $\delta\lambda(T)\sim 0.1\lambda_{L}^{{\scriptscriptstyle{(0)}}}t^{3}$.

The two physical quantities $B_{y}(z)$ and $\lambda(T)$ we focus on inherit the $T^{3}$ dependence from the surface current-current correlation function [Eq. (11)]. Similar power-law-dependence is observed in bulk nodal superconductors. In contrast to those systems, the model we considered is fully gapped in the bulk, and thus the Majorana surface states are responsible for the gapless excitations.

Although the exact expression for $G(z)$ is complicated SM , it takes relatively simple forms in the strong type-I and type-II limits. For a type-I TSC ($l_{\mathsf{coh}}\gg\lambda_{L}^{{\scriptscriptstyle{(0)}}}$),

[TABLE]

Since the coherence length sets the depth of the surface fluid [Eq. (8)], the latter allows a much deeper penetration of the field in the type-I limit than the bulk London depth. The slower decay is modulated by a Friedel oscillation. Representative field profiles are shown in Fig. 2(a).

For a type-II TSC ($l_{\mathsf{coh}}\ll\lambda_{L}^{{\scriptscriptstyle{(0)}}}$),

[TABLE]

In this case the Friedel oscillating terms are subleading, and can be neglected. Note that the first term in Eq. (17) dominates the second. In this case the spatial field penetration is governed by the London depth, but the correction in Eqs. (13) and (17) effectively enhances the field amplitude. Representative field profiles are indicated in Fig. 2(b).

Conclusion.—Our calculations suggest an alternative way to search for Majorana surface states in TSCs by measuring the change in the penetration depth $\delta\lambda(T)$. ARPES is often employed as the key tool to detect smoking gun signatures of topology in quantum materials. Unfortunately, TSC candidates typically have a small gap, making surface states difficult for ARPES to resolve Wray10 . A necessary, but not sufficient condition for the existence of a gapless Majorana surface fluid can be the signature power-law dependence $\delta\lambda(T)\sim T^{3}$. This can be probed by means of tunnel-diode oscillator techniques Smylie16 ; Smylie17 ; Paglione18 ; Spinu18 . While $\delta\lambda(T)$ is exponentially suppressed in conventional, topologically trivial superconductors, a power-law $\delta\lambda(T)\sim T^{\alpha}$ is also expected in superconductors with nodal bulk excitations, where $\alpha$ depends on whether the nodes are points or lines Tinkham ; Sigrist_Review . For example, a power-law with $\alpha=1$ has been observed in high-$T_{c}$ $d$-wave superconductors Giannetta06 .

One way to distinguish the origin of a penetration-depth power law (Majorana surface states versus bulk nodes) is via a specific-heat measurement versus temperature. The specific heat due to the 2D surface of a fully gapped 3D TSC is negligible. In addition to $\delta\lambda(T)\sim T^{\alpha}$, a TSC with a fully gapped bulk (nodal superconductor) should therefore demonstrate exponential suppression (power-law dependence) in specific heat Sigrist_Review .

The superconducting doped topological insulators (Cu,Nb)xBi2Se3 Ando12 ; Smylie16 ; Smylie17 and the half-Heusler compound YPtBi Paglione18 ; Bay12 are all strongly type II. Power-law dependence of $\lambda(T)$ was observed in NbxBi2Se3 Smylie16 ; Smylie17 , YPtBi Paglione18 , and in YPdBi and TbPdBi Spinu18 . It would be interesting to assess whether any of these could be attributed to the presence of Majorana surface states. In particular, it is worthwhile to note that while specific-heat measurements suggest a fully gapped bulk in CuxBi2Se3 and NixBi2Se3Yonezawa17 ; Qiu17 , penetration depth data for NixBi2Se3 shows a power-law dependence Smylie16 ; Smylie17 . In the case of half-Heusler compounds such as YPtBi, it has been suggested that optical-phonon-mediated pairing could favor a fully gapped TSC state with winding number $\nu=3$ Savary17 , and this should induce novel, cubic-dispersing Majorana surface fermions Fang15 ; Roy19 . In this case, we would expect a very slow $\delta\lambda(T)\sim T^{1/3}$ dependence for a clean, cubically-dispersing Majorana surface fluid. However, the results due to surface states with $\nu\geq 3$ might be strongly modified by quenched disorder Roy19 ; Foster14 .

In this Letter, calculations were performed for a clean system. On one hand, low-energy Majorana surface states with winding number $\nu=1$ are not affected by non-magnetic impurities Foster14 ; Ghorashi19 ; weak nonmagnetic disorder should therefore not alter the cubic power law $\delta\lambda(T)\sim T^{3}$ predicted here as $T\rightarrow 0$. The power can possibly be modified by strong (resonant) impurity scattering Giannetta06 . On the other hand, magnetic impurities can strongly perturb the surface states, and even gap them out. Moreover, magnetic impurities may independently induce power-law dependence in $\delta\lambda(T)$, by altering the magnetic permeability Giannetta06 ; Cooper97 .

We thank Andriy Nevidomskyy for useful discussions. T.C.W. and M.S.F. acknowledge support by the Welch Foundation Grant No. C-1809, by NSF CAREER Grant No. DMR-1552327, and by the U.S. Army Research Office Grant No. W911NF-17-1-0259. M.S.F. thanks the Aspen Center for Physics, which is supported by the NSF Grant No. PHY-1607611, for its hospitality while part of this work was performed. P.H. and H.K.P. were supported by the Department of Physics, College of Natural Sciences and Mathematics at the University of Houston. H.K.P acknowledges support from IRCC, IIT Bombay (RD/0518-IRCCSH0-029).

I I. Majorana surface states

We solve for the surface states by converting $k_{z}\rightarrow-i\partial_{z}$ in Eq. (6),

[TABLE]

where $\bm{\mathrm{k}}=(k_{x},k_{y})$ specifies the momentum transverse to the vacuum-TSC interface at $z=0$ (Fig. 1). For this model with a hard wall boundary condition at $z=0$, the surface Majorana fluid has an exactly linear dispersion relation $\varepsilon=\Delta k$, corresponding to the surface wave function [Eq. (8)]

[TABLE]

where $e^{i\phi_{\textbf{k}}}\equiv(k_{x}+ik_{y})/k$ and the normalization constant

[TABLE]

The four-component spinor in Eq. (S2) is expressed in the [spin ($s$)]$\otimes$[particle-hole ($\sigma$)] basis such that $\hat{s}^{3}\rightarrow\operatorname{diag}(1,-1,1,-1)$ and $\hat{\sigma}^{3}\rightarrow\operatorname{diag}(1,1,-1,-1)$.

The Bogoliubov-de Gennes Hamiltonian $\hat{h}(\bm{\mathrm{k}})$ in Eq. (6) has the following particle-hole (P), time-reversal (T), and chiral (S) symmetries:

[TABLE]

where $\hat{A}^{\mathsf{T}}$ is the matrix transpose of $\hat{A}$. The chiral symmetry is a product of time-reversal and particle-hole; since the latter is an automatic consequence of fermion antisymmetry, chiral is equivalent to time-reversal.

The negative-energy surface eigenstate with momentum $\bm{\mathrm{k}}$ is the chiral transform of Eq. (S2), $\hat{\sigma}^{1}\,\psi_{\textbf{k}}^{\mathsf{s}}(z)$. Positive- and negative-energy surface states are bi-locally orthogonal (due to the spin-momentum–locked spinors),

[TABLE]

II II. Paramagnetic Current-Current Correlation Function

II.1 A. Correlation function from the bulk alone

The imaginary time action corresponding to Eq. (5) is

[TABLE]

where $\omega_{n}$ denotes a fermionic Matsubara frequency. The imaginary time paramagnetic current-current correlation function is

[TABLE]

where we have made use of the translational invariance in the bulk. Using the Green’s function $\hat{G}(i\omega,\textbf{k})=\left[-i\omega+\hat{h}(\textbf{k})\right]^{-1},$ the Fourier transformed correlation function is

[TABLE]

where $E_{k}=\sqrt{\tilde{\varepsilon}_{k}^{2}+\Delta^{2}k^{2}}$. The retarded version is given by Eq. (7) in the main text.

II.2 B. Correlation functions for the semi-infinite slab

To consider the effect of the surface, we must retain the $(z,z^{\prime})$-dependence of the Green’s function. Eqs. (S2) and (S3) are replaced by

[TABLE]

where $\textbf{r}=(x,y)$ and $\textbf{k}=(k_{x},k_{y})$ are 2D vectors parallel to the interface. Then

[TABLE]

We assume a generic eigenstate decomposition for $\hat{G}$,

[TABLE]

where the sum runs over all positive- and negative-energy bulk and surface states of $\hat{h}$ in Eq. (S1), so that

[TABLE]

where

[TABLE]

For a system that is isotropic (rotationally invariant) in the $(x,y)$ plane parallel to the interface, the double-Fourier transform of the retarded version is

[TABLE]

II.2.1 1. Surface-surface response

The surface eigenstates $\psi_{\bm{\mathrm{k}}}^{\mathsf{s}}(z)$ and $\hat{\sigma}^{1}\,\psi_{\bm{\mathrm{k}}}^{\mathsf{s}}(z)$ [Eqs. (8) and (S2)] respectively have eigenenergies $\pm\Delta|\bm{\mathrm{k}}|$. The surface-surface contribution to Eq. (II.2) is

[TABLE]

where we have used Eq. (S7), and where [Eq. (12)]

[TABLE]

The ultraviolet momentum cutoff $\lambdabar_{F}^{-1}=k_{F}$ in Eq. (S10) is where the surface Majorana band merges with the bulk quasiparticle continuum. For low temperatures $\beta\rightarrow\infty$, we can extend the upper limit of the $k$ integration to infinity, and drop the dependence of $\Theta(k,q)$ on $k$. To leading order in temperature, one obtains Eq. (11) in the main text.

II.2.2 2. Surface-bulk cross terms

The surface-bulk cross term contributions to Eq. (II.2) take the form

[TABLE]

where $\alpha_{\mathsf{s},\mathsf{b}}(\bm{\mathrm{k}};q_{z},-Q_{z})$ is a temperature-independent coefficient encoding the $(q_{z},-Q_{z})$-transformed overlaps between the surface bound Majorana and bulk standing wave quasiparticle states. The summation runs over all surface and bulk eigenstates with eigenenergies $E_{\mathsf{s}}(\bm{\mathrm{k}})$ and $E_{\mathsf{b}}(\bm{\mathrm{k}})$, respectively. At low temperatures, the expression is dominated by small energies, but the mismatch between the gapless surface [$\lim_{k\rightarrow 0}|E_{\mathsf{s}}(\bm{\mathrm{k}})|=\Delta k\rightarrow 0$] and the gapped bulk [$\min|E_{\mathsf{b}}(\bm{\mathrm{k}})|=\Delta k_{F}$] means that the leading temperature dependence of this term is

[TABLE]

which is higher order than the contribution of the surface-surface term.

III III. Low-temperature field penetration

III.1 A. Solution to the integral equation

The kernel in Eq. (11) can be identified as the matrix elements of an outer product

[TABLE]

Here we assume the norm and resolution of the identity,

[TABLE]

The Meissner response in Eq. (3) can then be written as

[TABLE]

Formally, we can invert the operator to obtain

[TABLE]

where

[TABLE]

is a temperature-dependent constant that goes to $1+\bm{\mathit{O}}\left(T^{3}\right)$ as $k_{B}T\rightarrow 0$. Using $l_{\mathsf{coh}}=(m\Delta)^{-1}$, $\lambdabar_{F}=1/\sqrt{2m\mu}$, $\lambda_{L}^{{\scriptscriptstyle{(0)}}}=\sqrt{mc^{2}/(4\pi e^{2}n)}$, $n=k_{F}^{3}/3\pi^{2}$, and defining physical BCS gap as $\Delta_{BCS}=k_{F}\Delta$, we define

[TABLE]

which is the length scale introduced in Eq. (14), after replacing $\Xi(T)\rightarrow 1$ (valid in the low-temperature limit). Then Eq. (S4) becomes

[TABLE]

where the dimensionless function $G(z)$ is given by

[TABLE]

In the weak-pairing BCS limit ($\lambdabar_{F}<l_{\mathsf{coh}}$), this evaluates to

[TABLE]

From Eq. (S7), the magnetic field inside the slab is given by Eq. (13) in the main text. The results in Eqs. (16) and (17) obtain from the type I ($\{\lambdabar_{F},\lambda_{L}^{{\scriptscriptstyle{(0)}}}\}\ll l_{\mathsf{coh}}$) and type II ($\lambdabar_{F}\ll l_{\mathsf{coh}}\ll\lambda_{L}^{{\scriptscriptstyle{(0)}}}$) limits of Eq. (III.1).

III.2 B. Penetration Depth

In the expression for the effective penetration depth given by Eq. (15), the parameter $G(0)$ evaluates to [Eq. (III.1)]

[TABLE]

In the type-I limit ($\{\lambdabar_{F},\lambda_{L}^{{\scriptscriptstyle{(0)}}}\}\ll l_{\mathsf{coh}}$), this simplifies to

[TABLE]

In the opposite type-II limit ($\{\lambdabar_{F},l_{\mathsf{coh}}\}\ll\lambda_{L}^{{\scriptscriptstyle{(0)}}}$), Eq. (S10) instead becomes

[TABLE]