Disorder and Zeeman coupling induced gap-filling states in the nodeless chiral superconducting Bi/Ni bilayer system

TL;DR

This paper investigates how impurities and Zeeman coupling induce gap-filling states in a fully gapped chiral superconductor in Bi/Ni bilayers, potentially explaining observed zero-bias anomalies.

Contribution

It demonstrates that both impurities and Zeeman fields can generate in-gap states in a helical d-wave superconductor, offering a new understanding of experimental tunneling anomalies.

Findings

Impurities create states within the superconducting gap.

Zeeman coupling can produce Bogoliubov Fermi surfaces.

Results explain zero-bias anomalies in tunneling measurements.

Abstract

Motivated by the recently discovered time-reversal symmetry-breaking superconductivity in epitaxial Bi/Ni bilayer system with transition temperature $T_c\approx 4.2$K and the observation of zero-bias anomaly in tunneling measurements, we show that gap-filling states can appear in the fully gapped $d_{xy}\pm id_{x^2-y^2}$ superconducting states. We consider a model of helical electron states with d-wave pairing. In particular, we show that both magnetic and non-magnetic impurities can create states within the superconducting gap. Alternatively, we also show that the coupling of the electron spins to the in-plane Zeeman field provided by nickel can also create gap-filling states by producing Bogoliubov Fermi surfaces. Our findings may explain the origin of zero-bias anomaly observed in the point-contact tunneling measurements.

Figures

Equations

$$H_0^{\prime }=\sum _{\mathbf{k},\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}\sigma }^{†}{\left[v_F|\mathbf{k}|(\sin {\theta }_{\mathbf{k}}{\sigma }_y+\cos {\theta }_{\mathbf{k}}{\sigma }_x)-\mu \right]}_{\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}{\sigma }^{\prime }},$$

$$d_{\mathbf{k}\lambda }^{†}=\frac{1}{\sqrt{2}}\left(c_{\mathbf{k}\uparrow }^{†}+\lambda e^{i{\theta }_{\mathbf{k}}}{c}_{\mathbf{k}\downarrow }^{†}\right),$$

$$H_{sc}=\sum _{\mathbf{k}}\Delta e^{-i2{\theta }_{\mathbf{k}}}{d}_{\mathbf{k}+}^{†}{\tilde{d}}_{\mathbf{k}+}^{†}+h.c.,$$

$$H_0=\sum _{\mathbf{k}}{\psi }_{\mathbf{k}}^{†}\left[(v_F|\mathbf{k}|-\mu ){\tau }_3+\Delta \cos 2{\theta }_{\mathbf{k}}{\tau }_1+\Delta \sin 2{\theta }_{\mathbf{k}}{\tau }_2\right]{\psi }_{\mathbf{k}},$$

$$H=H_0+H_{\text{imp}}.$$

$$H_{\text{imp}}=\sum _{i=1}^{N_{\text{imp}}}\sum _{\sigma }{c}_{i\sigma }^{†}\left({\mu }_i{\delta }_{\sigma {\sigma }^{\prime }}+J{\mathbf{S}}_i\cdot {\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}\right)c_{i{\sigma }^{\prime }}$$

$$H_{\text{imp}}=\frac{1}{V}\sum _{i,\sigma ,\mathbf{k}{\mathbf{k}}^{\prime }}{e}^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{R}}_i}{c}_{\mathbf{k}\sigma }^{†}\left({\mu }_i{\delta }_{\sigma {\sigma }^{\prime }}+J{\mathbf{S}}_i\cdot {\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}\right)c_{{\mathbf{k}}^{\prime }\sigma },$$

$$H_{\text{imp}}=\frac{1}{2V}\sum _{i,\lambda {\lambda }^{\prime },\mathbf{k}{\mathbf{k}}^{\prime }}{d}_{\mathbf{k}\lambda }^{†}{\bar{\mathcal{H}}}_{\text{imp}}(\mathbf{k},{\mathbf{k}}^{\prime })d_{{\mathbf{k}}^{\prime }{\lambda }^{\prime }},$$

$${\bar{\mathcal{H}}}_{\text{imp}}(\mathbf{k},{\mathbf{k}}^{\prime })=e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{R}}_i}\left({\mu }_i{v}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}+J{\mathbf{S}}_i\cdot {\mathbf{m}}_{\lambda ,\lambda }^{\mathbf{k},{\mathbf{k}}^{\prime }}\right)$$

$$v_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}\equiv 1+\lambda {\lambda }^{\prime }{e}^{-i{\theta }_{\mathbf{k}}+i{\theta }_{{\mathbf{k}}^{\prime }}}$$

$$({\mathbf{m}}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}{)}_x\equiv \lambda e^{-i{\theta }_{\mathbf{k}}}+{\lambda }^{\prime }{e}^{i{\theta }_{{\mathbf{k}}^{\prime }}}$$

$$({\mathbf{m}}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}{)}_y\equiv i\lambda e^{-i{\theta }_{\mathbf{k}}}-i{\lambda }^{\prime }{e}^{i{\theta }_{{\mathbf{k}}^{\prime }}}$$

$$({\mathbf{m}}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}{)}_z\equiv 1-\lambda {\lambda }^{\prime }{e}^{-i{\theta }_{\mathbf{k}}+i{\theta }_{{\mathbf{k}}^{\prime }}}$$

$$H_{\text{imp}}=\frac{1}{2V}\sum _{i,\mathbf{k}{\mathbf{k}}^{\prime }}{\psi }_{\mathbf{k}}^{†}{\mathcal{H}}_{\text{imp}}(\mathbf{k},{\mathbf{k}}^{\prime }){\psi }_{{\mathbf{k}}^{\prime }},$$

$${\mathcal{H}}_{\text{imp}}(\mathbf{k},{\mathbf{k}}^{\prime })=e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{R}}_i}\left({\mu }_i{V}^{\mathbf{k},{\mathbf{k}}^{\prime }}+J{\mathbf{S}}_i\cdot {\mathbf{M}}^{\mathbf{k},{\mathbf{k}}^{\prime }}\right).$$

$$V^{\mathbf{k},{\mathbf{k}}^{\prime }}=\left(\begin{array}{cc}{v}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{v}}_{+,+}^{\mathbf{k},{\mathbf{k}}^{\prime }}\end{array}\right),$$

$${\mathbf{M}}^{\mathbf{k},{\mathbf{k}}^{\prime }}=\left(\begin{array}{cc}{\mathbf{m}}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{\mathbf{m}}}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}\end{array}\right),$$

$${\tilde{v}}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}=(-\lambda {\lambda }^{\prime }{e}^{-i{\theta }_{\mathbf{k}}+i{\theta }_{{\mathbf{k}}^{\prime }}})v_{{\lambda }^{\prime },\lambda }^{-{\mathbf{k}}^{\prime },-\mathbf{k}},$$

$${\tilde{\mathbf{m}}}_{\lambda ,{\lambda }^{\prime }}^{\mathbf{k},{\mathbf{k}}^{\prime }}\equiv -\lambda {\lambda }^{\prime }{e}^{-i{\theta }_{\mathbf{k}}+i{\theta }_{{\mathbf{k}}^{\prime }}}{\mathbf{m}}_{{\lambda }^{\prime },\lambda }^{-{\mathbf{k}}^{\prime },-\mathbf{k}}.$$

$$G(\mathbf{k},i{\omega }_n{)}^{-1}=G^{(0)}(\mathbf{k},i{\omega }_n{)}^{-1}-\Sigma (\mathbf{k},i{\omega }_n),$$

$$\Sigma (\mathbf{k},i{\omega }_n)$$

$$+\sum _{{\mathbf{k}}^{\prime }}⟨{\mathcal{H}}_{\text{imp}}(\mathbf{k},{\mathbf{k}}^{\prime })G({\mathbf{k}}^{\prime },i{\omega }_n){\mathcal{H}}_{\text{imp}}({\mathbf{k}}^{\prime },\mathbf{k})⟩,$$

$$\frac{1}{V}{⟨u_i{u}_j{e}^{-i\mathbf{q}\cdot R_i}{e}^{-i{\mathbf{q}}^{\prime }\cdot R_j}⟩}_{\text{imp}}=n_{\text{imp}}\overline{u^2}{\delta }_{i,j}{\delta }_{\mathbf{q},{\mathbf{q}}^{\prime }},$$

$$\frac{1}{V}{⟨S_i^{\alpha }{S}_j^{\beta }{e}^{-i\mathbf{q}\cdot R_i}{e}^{i{\mathbf{q}}^{\prime }\cdot R_j}⟩}_{\text{imp}}=\frac{1}{3}{n}_{\text{imp}}S(S+1){\delta }_{ij}{\delta }_{\alpha \beta }{\delta }_{\mathbf{q}+{\mathbf{q}}^{\prime },0},$$

$${\tilde{\omega }}_n={\omega }_n+\Gamma \frac{{\tilde{\omega }}_n}{2\sqrt{{\tilde{\Delta }}_n^2+{\tilde{\omega }}_n^2}}$$

$${\tilde{\Delta }}_n=\Delta +\alpha \Gamma \frac{{\tilde{\Delta }}_n}{2\sqrt{{\tilde{\Delta }}_n^2+{\tilde{\omega }}_n^2}}$$

$$T_c(n_i)=T_c(0)-\frac{\pi }{8}(1-\alpha )\Gamma .$$

$$N(\omega )=-\frac{1}{\pi }N(0)\mathrm{Im}\left[\underset{i{\omega }_n\to \omega +i{0}^{+}}{\lim }\frac{-i{\tilde{\omega }}_n}{\sqrt{{\tilde{\omega }}_n^2+{\Delta }^2}}\right].$$

$$H_0=\sum _{\mathbf{k},\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}\sigma }^{†}{H}_{\sigma {\sigma }^{\prime }}(\mathbf{k})c_{\mathbf{k}{\sigma }^{\prime }}$$

$$H_{\sigma {\sigma }^{\prime }}(\mathbf{k})={\left[\frac{|\mathbf{k}{|}^2}{2m}{\delta }_{\sigma {\sigma }^{\prime }}+v_R|\mathbf{k}|(\sin {\theta }_{\mathbf{k}}{\sigma }_y+\cos {\theta }_{\mathbf{k}}{\sigma }_x)-\mu \right]}_{\sigma {\sigma }^{\prime }},$$

$$H_{sc}=\sum _{\mathbf{k},\lambda }{\Delta }_{\lambda }{e}^{-i2{\theta }_{\mathbf{k}}}{d}_{\mathbf{k}\lambda }^{†}{\tilde{d}}_{\mathbf{k}\lambda }^{†}+h.c.,$$

$$H=\sum _{\mathbf{k}\lambda }{\psi }_{\mathbf{k}}^{†}\left(\frac{{\sigma }_0+\lambda {\sigma }_3}{2}\right)\otimes \left[\left(\frac{|\mathbf{k}{|}^2}{2m}+\lambda v_R|\mathbf{k}|-\mu \right){\tau }_3+{\Delta }_{\lambda }(\cos 2{\theta }_{\mathbf{k}}{\tau }_1+\sin 2{\theta }_{\mathbf{k}}{\tau }_2)\right]{\psi }_{\mathbf{k}}.$$

$$V^{\mathbf{k},{\mathbf{k}}^{\prime }}=\left(\begin{array}{cccc}{v}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& v_{+-}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{v}}_{+,+}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\tilde{v}}_{+,-}^{\mathbf{k},{\mathbf{k}}^{\prime }}\\ v_{-+}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& v_{--}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{v}}_{-+}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\tilde{v}}_{--}^{\mathbf{k},{\mathbf{k}}^{\prime }}\end{array}\right),$$

$${\mathbf{M}}^{\mathbf{k},{\mathbf{k}}^{\prime }}=\left(\begin{array}{cccc}{\mathbf{m}}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\mathbf{m}}_{+-}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{\mathbf{m}}}_{++}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\tilde{\mathbf{m}}}_{+-}^{\mathbf{k},{\mathbf{k}}^{\prime }}\\ {\mathbf{m}}_{-+}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\mathbf{m}}_{--}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0\\ 0& {\tilde{\mathbf{m}}}_{-+}^{\mathbf{k},{\mathbf{k}}^{\prime }}& 0& {\tilde{\mathbf{m}}}_{--}^{\mathbf{k},{\mathbf{k}}^{\prime }}\end{array}\right).$$

Full text

I Introduction

Interplay between topology and electronic band structures in insulators, superconductors, and metals has received enormous attention and interest in recent years and has become a central issue in the condensed matter physics Hasan and Kane (2010); Hasan and Moore (2011); Qi and Zhang (2011). The nontrivial topology encoded in the bulk band structures, usually characterized by integer numbers, yield nontrivial consequences for the electronic states living on the boundary of the system: the topologically protected gapless surface states. Helical edge states on the one-dimensional boundary of HgTe Bernevig et al. (2006); König et al. (2007), Dirac electrons on the surface of three-dimensional topological materials Bi2Se3, Bi2Te3Zhang et al. (2009); Chen et al. (2009), Majorana fermions at the open ends of a one-dimensional topological superconductorKitaev (2001), and Femi arcs in topological Weyl semimetalsWan et al. (2011) such as TaAsLv et al. (2015a, b); Xu et al. (2015) are a few known examples.

Beside all, especial attentions have been paid to realization of topological superconductors promising a platform for topological quantum computationsNayak et al. (2008), a central paradigm in building a quantum computer. At the heart of material realization lies the nontrivial pairing wave functions of Cooper pairs of electrons in the vicinity of the Fermi surface. While the phonon-based mechanism for superconductivity leads to the conventional s-wave superconductor, the intrinsic spin or charge density fluctuations due to electron correlations may result in more complicated structure for pairing wave functionsLee et al. (2006); Scalapino (2012) for instance in cuprates, Sr2RuO4 Kallin and Berlinsky (2016); Luke et al. (1998),UPt3Schemm et al. (2014). Alternatively the nontrivial pairing structures can be induced by an ingenious combination of more conventional materials. A famously celebrated structure has been introduced by Fu and Kane in Ref.[Fu and Kane, 2008] where a conventional s-wave superconductor proximitized to the surface of topological insulators (sSc-TI interface) is shown to support Majorana zero-energy states in vortex cores. Yet, the surface of topological insulators can be replaced with more conventional two-dimensional electron gas with strong Rashba spin-orbit coupling, where the spin degeneracy is effectively lifted, and addition of a magnetic or Zeeman field turn the conventional induced s-wave superconductivity to a topological superconductorSau et al. (2010a, b).

The recently discovered superconductivity with rather high critical temperature $T_{c}\approx 4.2$K in epitaxial Bismuth-Nickel (Bi/Ni) bilayer system Xin-Xin et al. (2015) provided another example of a superconducting state with nontrivial pairing. The bilayer system is schematically shown in Fig. 1 where the surface termination Bi(110) is exposed to vacuum and Ni is a ferromagnet with in-plane magnetization. The optical measurements of Kerr signal from the Bi(110) side clearly shows that the superconducting states breaks the time-reversal symmetry (TRS) Gong et al. (2017). Not only that, cooling below $T_{c}$ in a weak magnetic field applied to the sample in both directions and subsequently measuring resistivity on warming up and in zero field shows that the pairing states must break the TRS spontaneously. The theoretical model proposed in Ref.[Gong et al., 2017] uses the dominantly spin-orbit coupled electronic states localized on the surface of Bi(110) and shows that the pairing symmetry should be $d_{xy}\pm id_{x^{2}-y^{2}}$, a chiral topological superconductor characterized by Chern numbers $c=\pm 2$. The complex structure of the superconducting order parameter indicates that the superconductivity in Bi/Ni system is fully gapped, an observation which is in agreement with the recently measured optical conductivity using time-domain THz spectroscopy Chauhan et al. (2019). The latter measurements further show that the superconductivity develops through the bulk of Bi with possible $p\pm ip$ pairing Chao (2018).

Beside the transport and optical measurements discussed above, a zero-bias anomaly has also been observed in the point-contact Andreev reflection Xin-Xin et al. (2015) which may seemingly contradict with having a fully gapped superconductor. This is exactly the motivation of our work in this paper. We ask the following specific question: how do the states filling the gap arise in a topologically nodless superconductor yielding a zero-bias anomaly? We give an affirmative answer to this question by introducing two possible scenarios. The gap-filling states appear as a result of (i) magnetic and non-magnetic impurities distributed randomly throughout, and/or (ii) emergence of a Bogolon Fermi surface due to the coupling of the electron spins to the in-plane magnetic moments as shown by arrows in Fig. 1. In the second scenario we treat the disorder effects using the Abrikosov-Grokov formalismAbrikosov and Gor’Kov (1960) and ignore the inhomogeneity which might result in enhancement of superconductivityAndersen et al. (2006); Gastiasoro et al. (2013); Gastiasoro and Andersen (2018); Gastiasoro et al. (2016); Rømer et al. (2012). In particular, we show that there is a critical disorder strength beyond which the gap closes. The first scenario was first introduced for the sSC-TI interface and the gapless states emerge in the presence of a Zeeman fieldYuan and Fu (2018). We extend this finding to the chiral superconductors which may explain the observed zero-bias anomaly in Bi/Ni system.

The paper is organized as follows. We first discuss the disorder-induced states in Sec II and then the effects of in-plane Zeeman coupling is discussed in Sec III, and Sec. LABEL:onclusion concludes.

II disorder-induced gap-filling states in chiral superconductors

The electron states at the Bi(110) surface is strongly spin-momentum locked due breaking of inversion symmetry near the surfaceHofmann (2006) with a large hole pocket located around the center of surface Brillouin zone. Though our main motivation is to understand the origin of zero-bias anomaly in Bi/Ni system, in this section we consider two cases relevant to surface states of topological insulators and an 2D electron gas in the presence of Rashba spin-orbit coupling. We also assume that the magnetic fluctuations of Ni provide the pairing glue as detailed in Ref. [Gong et al., 2017]. We study the effect of disorder on the superconducting states.

II.1 Surface of Topological Insulator

The effective theory for the surface states of a topological insulators is proportional to the 2D Hamiltonian $\bm{\sigma}\times\mathbf{k}\cdot\hat{z}$, where $\mathbf{k}=(k_{x},k_{y})$ is the wave vector and $\bm{\sigma}=(\sigma_{x},\sigma_{y})$ are Pauli matricesHasan and Kane (2010). After a rotation it can be brought to the canonical Dirac form $\bm{\sigma}\cdot\mathbf{k}$, or equivalently,

[TABLE]

where $\theta_{\mathbf{k}}$ is the polar angle of $\mathbf{k}$ plus $\pi/2$. This Hamiltonian can be diagonalized by the following transformation

[TABLE]

where $\lambda=\pm$ is helicity and label the energy bands. Motivated by earlier work on the TR breaking superconductivity in Bi/Ni system and the theoretical proposal of $d\pm id$ in this system Gong et al. (2017), let us formulate the effect of impurity scattering for such a superconducting state. Using the time reversal (TR) operator $\Theta=i\sigma_{y}\mathcal{K}$, with $\mathcal{K}$ as complex conjugation operator, the TR transformation of creation operator $d^{{\dagger}}_{\mathbf{k}\lambda}$ becomes $\tilde{d}^{{\dagger}}_{\mathbf{k}\lambda}=\Theta d^{{\dagger}}_{\mathbf{k}\lambda}\Theta^{-1}=\lambda e^{-i\theta_{\mathbf{k}}}d^{{\dagger}}_{-\mathbf{k}\lambda}$ which can be used to construct the pairing interaction from time-reversed partners as

[TABLE]

We have only included the pairing between the electron states with positive helicity which is justified when the chemical potential is much larger than the pairing energy scale, i.e. the paring Hamiltonian is projected on the Fermi surface. In terms of the Nambu spinor $\psi_{\mathbf{k}}=(d_{\mathbf{k}+},\tilde{d}^{{\dagger}}_{\mathbf{k}+})^{T}$ we can write the superconducting Hamiltonian as

[TABLE]

where Pauli matrices $\tau$’s act within the particle-hole space. On top of this clean pairing Hamiltonian we add a disorder term $H_{\text{\text{imp}}}$ such that the total Hamiltonian becomes,

[TABLE]

The $H_{\text{imp}}$ is given by

[TABLE]

which describes scattering from magnetic and non-magnetic impurities. The random scalar potential is given by the random variables $\mu_{i}$, while the randomness in the magnetic impurities is determined by the random orientation of the local spin $\mathbf{S}_{i}$. The exchange coupling $J$ is essentially determined by the hybridization of conduction and impurity electrons and the strength of the on-site Hubbard term which is assumed to be non-random. The average value of $\mu_{i}$ and $\mathbf{S}_{i}$ is zero, but their standard deviations are non-zero and together with concentration of impurities determine the strength of interactions.

By Fourier transform to momentum space the $H_{\text{imp}}$ becomes

[TABLE]

and in band basis takes the following form

[TABLE]

where

[TABLE]

where $v_{\lambda,\lambda^{\prime}}^{\mathbf{k},\mathbf{k^{\prime}}}$ and $\mathbf{m}^{\mathbf{k},\mathbf{k^{\prime}}}$ are

[TABLE]

The compact representation of the above Hamiltonian in the the Nambu space is,

[TABLE]

where

[TABLE]

Here by assumption of large positive chemical potential we focus on $\lambda=+1$ matrix elements, i.e.

[TABLE]

where we have introduced,

[TABLE]

In order to study the effect of disorder, we use the language of Green’s function. The disordered and clean Green function are defined as $G(\mathbf{k},i\omega_{n})=(i\omega_{n}-H)^{-1}$ and $G^{(0)}(i\omega_{n})=(i\omega_{n}-H_{0})^{-1}$, respectively, where $i\omega_{n}$’s are Fermionic Matsubara frequencies. According to Dyson equation, these two are related by

[TABLE]

where in the self-consistent Born approximation up to second order, the self energy is given byYuan and Fu (2018); Bruus and Flensberg (2004)

[TABLE]

where $\left<\cdots\right>$ means the ensemble average is taken over disorder configurations. Following the approach presented in Refs. [Ozfidan et al., 2016-Ito et al., 2012] for treatment of disordered system, we assume that the excitations can be described in terms of an effective clean single-particle Hamiltonian $\tilde{H}_{0}$, which is similar to the original $H_{0}$ with $\Delta$ now depends on $\omega_{n}$, i.e. we make the following replacement: $\Delta\to\Delta_{n}$ in $H_{0}$. Similarly for the Green function the spectral re-arrangements due to disorder are taken into account by $i\omega_{n}\to i\tilde{\omega}_{n}$ where $i\tilde{\omega}_{n}$ is a function of $i\omega_{n}$ as described below. The effective Hamiltonian $\tilde{H}_{0}$ satisfies $G(i\omega_{n})=(i\tilde{\omega}_{n}-\tilde{H}_{0})^{-1}$. The Eq. (13) and  (14) can now be cast into an algebraic self-consistency between $i\tilde{\omega}_{n}$ and $\Delta_{n}$; See Eqs. (17)-(18).

Since we have assumed that the mean value of $\mu_{i}$ and $\mathbf{S}_{i}$ is zero, the first term in Eq. (14) gives zero. However their standard deviations are non-zero and are given as follows: Bruus and Flensberg (2004)

[TABLE]

Obviously there is no correlation between electrostatic and magnetic scattering forces, $\langle\mu_{i}\mathbf{S}_{j}\rangle=0$. Therefore the non-zero contribution of the second term in Eq. (14) arises from the above fluctuations. By assumption of constant density of states (DOS) $N(0)$ at the Fermi surface the integral in the second term of Eq. (14) can be performed giving rise to the following equations for $i\tilde{\omega}_{n}$ and $\Delta_{n}$:

[TABLE]

where $\Gamma=\left(\bar{u^{2}}+S(S+1)J^{2}\right)\pi n_{\text{imp}}N(0)$ defines the disorder strength. Note that the assumption of constant $N(0)$ around the Fermi level is justified for parabolic bands, as well as highly doped Dirac cone. For chiral $d\pm id$ superconductor, $\alpha=0$. This is because the $e^{2i\theta_{\mathbf{k}}}$ factors in the pairing enforce the angular integration to become zero. Therefore in $d\pm id$ case, the pairing potential $\tilde{\Delta}_{n}$ does not change in the presence of disorder, in sharp contrast to s-wave pairing. For the case of s-wave pairing, and within a model that scatters single particles only (not the Cooper pairs), $\alpha=+1$Ito et al. (2012) ($-1$Ozfidan et al. (2016)) for non-magnetic (magnetic) disorders. This leads to an important difference between the chiral $d\pm id$ superconductivity and s-wave pairing. In the s-wave case, the fact that $\tilde{\Delta}_{n}$ depends on $n$, implies that for magnetic impurities ($\alpha=-1$) the $T_{c}$ will be suppressed much faster than the non-magnetic case ($\alpha=+1$). In the chiral $d\pm id$ case, the $\alpha$ is zero anyhow, and to that extent, magnetic and non-magnetic impurities have comparable effect on the $T_{c}$. As a function of impurity concentration, the transition temperature reads asOzfidan et al. (2016); Ito et al. (2012)

[TABLE]

Therefore the chiral superconductivity is fragile against impurity scattering, irrespective of the magnetic nature of scattering center.

The DOS can be found by momentum integration after analytical continuation $i\omega_{n}\to\omega+i0^{+}$ of the Green function yielding

[TABLE]

One way to perform the analytical continuation is to use the Padé approximation to represent $-i\tilde{\omega}_{n}/\sqrt{\tilde{\omega}_{n}^{2}+\Delta^{2}}$ by a function $F(i\omega_{n})=Q(i\omega_{n})/P(i\omega_{n})$, where $Q$ and $P$ are polynomials of $i\omega_{n}$, from which its analytical continuation $F(\omega+i0^{+})$ can be evaluated. Alternative way is to perform the analytical continuation before solving the self-consistency equations (17) and (18) Fathi and Jafari (2010). We have checked that these two approaches give identical results.

In Fig. 2, we have plotted the DOS for $\Delta=0.2$ as a function of energy for different disorder strength $\Gamma$ indicated in the legend. As can be seen, beyond a certain critical disorder strength the superconducting gap is completely filled and we will have a gapless superconductor. The value of DOS at $\omega=0$ determines whether the superconductor is gapless or gapped. Therefore in left panel of Fig. 3 we probe the density of states for $\omega\to 0$ in the plane of pairing and disorder strength. This gives us a border that separate gapless and gapped superconducting phases. To emphasize the difference between the behavior of chiral and non-chiral superconductors against disorder, in right panel of Fig. 3 for a fixed disorder strength of $\Gamma=0.5$, we have generated a map of the DOS at $\omega=0$ in the plane of $\alpha$ and $\Delta$. The chiral superconductor corresponds to $\alpha=0$. The s-wave superconductors with magnetic (non-magnetic) scatterers correspond to $\alpha=-1$ ($\alpha=+1$). This comparison indicates that the chiral superconductor is capable of producing gapless superconductivity much easier than s-wave superconductor. This can be a possible explanation for the observed in-gap features in the Bi/Ni system Xin-Xin et al. (2015).

II.2 Electron Gas with Rashba Coupling

In this section we use the methods used in preceding subsection to study the affect of disorder on chiral superconductors in system with multiple Fermi surfaces. For simplicity we consider a two-dimensional electron gas in the presence of a strong Rashba spin-orbit coupling described by the following Hamiltonian

[TABLE]

with

[TABLE]

where the spin-orbit coupling is expressed in the form of a velocity scale $v_{R}$ to make it look similar to the helical surface states considered in the preceding section. This Hamiltonian reduces to the helical Hamiltonian (1) upon replacement $m\to\infty$ and $v_{R}\to v_{F}$. The pairing in the two cases is however different. In the helical metallic case for any chemical potential, there is only one Fermi contour. But in the present Rashba spin-orbit coupled case, the spin-orbit coupling splits them into two Fermi contours. For large enough chemical potential, the two spin-orbit split Fermi contours will have opposite helicity. Therefore the pairing will include both helicities $\lambda=\pm 1$

[TABLE]

which together with Eq. (21) yields the following superconducting Hamiltonian

[TABLE]

where $\psi_{\mathbf{k}}=(d_{\mathbf{k}+},\tilde{d}^{{\dagger}}_{\mathbf{k}+},d_{\mathbf{k}-},\tilde{d}^{{\dagger}}_{\mathbf{k}-})^{T}$ is the Nambu spinor.

The calculation proceeds along the same steps as preceding subsection, with a difference that instead of one helicity ($\lambda=+$ for the surface of TI), now will have both helicities $\lambda=\pm$ and hence for the impurity Hamiltonians in Eqs. (9)-(10) the relevant matrices will be $4\times 4$ which are given by

[TABLE]

and

[TABLE]

Again due to a vanishing angular integration coming from the $e^{2i\theta_{\mathbf{k}}}$, for chiral $d\pm id$ superconductor we obtain $\tilde{\Delta}_{\lambda,n}=\Delta_{\lambda}$ which again corresponds to $\alpha=0$ situation in Eq. (18), albeit with new helicity index acquired due to two Fermi contours. The self consistency equation corresponding to Eq. (17) will now become

[TABLE]

where

[TABLE]

is the DOS for Fermi contour with helicity $\lambda$,

[TABLE]

In this case from Eq. (27) it follows that the total DOS is given by

[TABLE]

which again can be calculated as before, either by Padé analytic continuation, or direct solution of the self-consistency equations slightly above the real frequency axisFathi and Jafari (2010). Density of states for different $\Delta_{+}/\Delta_{-}$ ratio with $\Delta_{+}=0.2$ and $\xi=1.2$ as a function of energy plotted in Fig. 4 for different disorder strength (a) $\Gamma=0$, (b) $\Gamma=0.05$, (c) $\Gamma=0.1$ and (d) $\Gamma=0.15$. In panel (a) where the disorder is zero, both gaps are clearly visible as two coherence peaks. Adding a little bit of disorder, in panel (b), the visible peaks are smeared out and both gaps start to get filled. Indeed within Eq. (27), if the DOS at zero for the two bands would be the same ($\xi=1$) then both gaps would evolve similarly. In the present case corresponding to $\xi=1.2$, both gaps are filled almost in the same manner. In this way, the smaller of the two gaps will be filled in first. This is what happens in panel (c) where the smaller gap is filled in, while the larger gap still survives. Finally in panel (d), both gaps are filled. Transport measurements can clearly indicate the value of $\Gamma$ at which the smaller gap is filled. This is a disorder strength at which zero-bias conductance takes off from zero. The second critical disorder at which the larger gap is filled, shows up as a kink in the trend of zero-bias peak as a function of disorder. The presence of such a kink is an essential difference between the chiral superconductivity in helical states and Rashba spin-orbit coupled systems. If the system is in a regime where the smaller gap is already filled, specific heat measurement will detect the in-gap features and the larger gap.

III Gap-Filing states caused by in-plain Zeeman coupling

In this section we investigate the effect of in-plane magnetic field as another possible origin of mid-gap states. This mechanism was suggested recently Yuan and Fu (2018) in the context of sSc-TI. In this section we investigate this mechanism for the chiral $d\pm id$ superconductivity and also for the $d$-wave pairings. In our set-up shown in Fig. 1 the in-plane Zeeman coupling can be provided by the proximity to Ni. We start by adding an in-plane Zeeman field $V\sigma_{y}$ to Hamiltonian (21)

[TABLE]

where the matrix $\Delta(\mathbf{k})$ is parameterized as follows,

[TABLE]

In following for numerical calculations we set $\left(m,v_{R},\mu\right)=\left(0.3,1,0.5\right)$. To address the surface of topological insulator we let $1/2m\to 0$ and $v_{R}=v_{F}$.

For the Rashba model with two Fermi contours $\lambda=\pm$, and consequently pairings $\Delta_{\pm}(\mathbf{k})$, the gap functions in Eq. (30) can be written asSamokhin and Mineev (2008); Ghadimi et al. (2018) $\Delta_{1}(\mathbf{k})=\Delta_{+}(\mathbf{k})-\Delta_{-}(\mathbf{k})$ and $\Delta_{0}(\mathbf{k})=\Delta_{+}(\mathbf{k})+\Delta_{-}(\mathbf{k})$. We take a general case for band pairings as $\Delta_{\pm}(\mathbf{k})=(\Delta_{0}\pm\Delta_{1})\left(\zeta_{x^{2}-y^{2}}\cos 2\theta_{\mathbf{k}}+i\zeta_{xy}\sin 2\theta_{\mathbf{k}}\right)$, where $(\zeta_{x^{2}-y^{2}},\zeta_{xy})$ are ${\cal O}(1)$ numerical parameters that encode the relative contribution of $d_{x^{2}-y^{2}}$ and $d_{xy}$ pairings. With this identification for pairing functions, we rewrite the pairing matrix (30) as

[TABLE]

where here $(k_{x},k_{y})=|\mathbf{k}|(\cos\theta_{\mathbf{k}},\sin\theta_{\mathbf{k}})$. In this way the explicit $d\pm id$ structure has been factored out. Now the amplitudes $\Delta_{0}$ and $\Delta_{1}$ do not depend on $\mathbf{k}$. The above pairing matrix has the structure that the superconductivity is nodal unless both $\zeta_{x^{2}-y^{2}}$ and $\zeta_{xy}$ become nonzero, which then breaks the TRS.

In Fig. 5 we present DOS for various cases. The three rows, from top to bottom, correspond to $(\zeta_{x^{2}-y^{2}},\zeta_{xy})=(1,0)$, $(\zeta_{x^{2}-y^{2}},\zeta_{xy})=(0,1)$ and $(\zeta_{x^{2}-y^{2}},\zeta_{xy})=(1,1)$, respectively. The panels in the first column depict the results for the surface of topological insulators, while those in next two columns are for the electron gas with Rashba coupling.

Let us first consider the surface of topological insulator. As pointed out, in this case we need to set $\Delta_{0}=\Delta_{1}=0.5$. The corresponding DOS is shown in panels (a), (b), and (c) in Fig. 5. The first two panels, (a) and (b), correspond to nodal $d_{x^{2}-y^{2}}$ and $d_{xy}$ pairing symmetries. As seen, in the absence of the Zeeman field we have a V-shaped pseudo-gap arising from the node in the pairing function as shown by solid red line. A smallest amount of in-plane Zeeman $V$ gives rise to Bogoliubov Fermi contours, and consequently a finite DOS is created at $E=0$. The most interesting situation occurs when the superconducting state is fully gapped (panel (c)) due to TRS breaking of the $d+id$ pairing wave function. In the absence of the Zeeman coupling there is clear gap in the DOS. However, a rather strong Zeeman coupling can fill up the gap by creating Bogoliubov Fermi contours.

For the electron gas with Rashba coupling we can also tune $\Delta_{0}$ and $\Delta_{1}$ independently. We consider two extreme cases where $(\Delta_{0},\Delta_{1})=(1,0)$ correspond to purely singlet pairing as shown in panels (d), (e), and (f), and $(\Delta_{0},\Delta_{1})=(0,1)$ represent the triplet case in panels (g), (h), and (i). Again it is clearly seen that for nodal cases a tiny Zeeman field creates finite DOS at zero energy, and the fully gapped case needs stronger Zeeman coupling for the gap to be filled with Bogoliubov quasiparticles.

IV Conclusion

This work is mainly motivated by the observation of zero-bias anomaly in the point-contact tunneling measurements in the nodless time-reversal symmetry-breaking superconducting epitaxial Bi/Ni bilayer system. In the previous work by one of the authors and coworkers in Ref.[Gong et al., 2017], it was proposed that the pairing symmetry should be the surface chiral $d\pm id$, in agreement with Kerr measurements and thickness-dependent of transition temperature. In this work we considered the effects of both magnetic and non-magnetic impurities on the superconducting gap. We showed in all cases mid-gap states appear in the gap and start to fill it up. We also showed that the in-plane Zeeman coupling, which is produced by Ni lying underneath Bi, can also create Bogoliubov Fermi contours of quasiparticles with non-vanishing spectral weight within the gap. We believe our finding may explain the anomalous zero-bias signal observed.

V Acknowledgements

M. K. acknowledges the support from the Sharif University of Technology under Grant No. G690208. S. A. J. acknowledges grant number G960214 of the research deputy, Sharif University of Technology.