Theory of Josephson Current on a Magnetically Doped Topological Insulator

TL;DR

This paper develops a theoretical framework to analyze how the Josephson current and local density of states in magnetically doped topological insulators change across topological and metal-insulator phase transitions, highlighting the emergence of chiral Majorana modes.

Contribution

It provides an analytical calculation of the Josephson current and density of states considering topological phase transitions in magnetically doped topological insulators.

Findings

Josephson current magnitude increases when surface states become metallic.

Temperature dependence of Josephson current changes drastically at the topological transition.

Chiral Majorana modes influence the superconducting properties at the phase transition.

Abstract

Proximity induced superconducting states in the surface of magnetically doped topological insulators can host chiral Majorana modes. We consider a Josephson junction in that system with changing the chemical potential, which drives a topological phase transition in the induced superconducting states as well as a metal-insulator transition in the surface states. The local density of states and the Josephson current are analytically calculated by McMillan's Green's function method in terms of the Andreev reflection coefficient. We show that although the magnitude of the Josephson current is greatly enhanced when the surface state changes from insulating to metallic, its temperature dependence drastically changes at the topological phase transition point, reflecting the appearance of the chiral Majorana modes

Equations

$$H(\mathbf{x},{\partial }_{\mathbf{x}})=\left(\begin{array}{cc}h({\partial }_{\mathbf{x}})+M& i{\sigma }_y\bar{\Delta }(x)\\ -i{\sigma }_y\bar{\Delta }(x)& -h^{\ast }({\partial }_{\mathbf{x}})-M^{\ast }\end{array}\right)$$

$$\left[E-H(\mathbf{x},{\partial }_{\mathbf{x}})\right]G(\mathbf{x},{\mathbf{x}}^{\prime })=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })$$

$$\displaystyle G(x,x^{\prime},k_{y})=\left\{\begin{array}[]{lll}\hat{\psi}_{-}(x)\hat{\tilde{\psi}}^{\mathsf{T}}_{+}(x^{\prime})&\textrm{for}&x<x^{\prime}\\ \hat{\psi}_{+}(x)\hat{\tilde{\psi}}^{\mathsf{T}}_{-}(x^{\prime})&\textrm{for}&x>x^{\prime}\end{array}\right.,$$

$${\hat{A}}_1=\left(\begin{array}{c}i\\ {\zeta }_1{e}^{i{\theta }_1}\\ -{\zeta }_1{e}^{i{\theta }_1}{\gamma }_{-}\\ i{\gamma }_{+}\end{array}\right),{\hat{A}}_2=\left(\begin{array}{c}i{\zeta }_2{e}^{i{\theta }_2}{\gamma }_{-}\\ -{\gamma }_{+}\\ 1\\ i{\zeta }_2{e}^{i{\theta }_2}\end{array}\right),{\hat{A}}_3=\left(\begin{array}{c}i{e}^{i{\theta }_1}\\ -{\zeta }_1\\ {\zeta }_1{\gamma }_{-}\\ i{e}^{i{\theta }_1}{\gamma }_{+}\end{array}\right),{\hat{A}}_4=\left(\begin{array}{c}i{\zeta }_2{\gamma }_{-}\\ e^{i{\theta }_2}{\gamma }_{+}\\ -e^{i{\theta }_2}\\ i{\zeta }_2\end{array}\right),$$

$${\zeta }_j$$

$$e^{\pm i{\theta }_j}$$

$${\gamma }_{\pm }$$

$$Z_{j\pm }$$

$${\Omega }_z$$

$${\Delta }_z$$

$$k_j$$

$$G^{\mathrm{R}}(x,x^{\prime },k_y)=\{\begin{array}{l}{g}_1\left[{\hat{A}}_1{\hat{B}}_3^{\mathsf{T}}{e}^{i{k}_1(x-x^{\prime })}+a_1{\hat{A}}_4{\hat{B}}_3^{\mathsf{T}}{e}^{i{k}_{2}x-i{k}_1{x}^{\prime }}+b_1{\hat{A}}_3{\hat{B}}_3^{\mathsf{T}}{e}^{-i{k}_1(x+x^{\prime })}\right]\\ +g_4\left[{\hat{A}}_2{\hat{B}}_4^{\mathsf{T}}{e}^{-i{k}_2(x-x^{\prime })}+a_2{\hat{A}}_3{\hat{B}}_4^{\mathsf{T}}{e}^{-i{k}_{1}x+i{k}_2{x}^{\prime }}+b_2{\hat{A}}_4{\hat{B}}_4^{\mathsf{T}}{e}^{i{k}_2(x+x^{\prime })}\right](x>x^{\prime })\\ g_1\left[{\hat{A}}_3{\hat{B}}_1^{\mathsf{T}}{e}^{-i{k}_1(x-x^{\prime })}+{\tilde{a}}_1{\hat{A}}_3{\hat{B}}_4^{\mathsf{T}}{e}^{-i{k}_{1}x+i{k}_2{x}^{\prime }}+{\tilde{b}}_1{\hat{A}}_3{\hat{B}}_3^{\mathsf{T}}{e}^{-i{k}_1(x+x^{\prime })}\right]\\ +g_4\left[{\hat{A}}_4{\hat{B}}_2^{\mathsf{T}}{e}^{i{k}_2(x-x^{\prime })}+{\tilde{a}}_2{\hat{A}}_4{\hat{B}}_3^{\mathsf{T}}{e}^{i{k}_{2}x-i{k}_1{x}^{\prime }}+{\tilde{b}}_2{\hat{A}}_4{\hat{B}}_4^{\mathsf{T}}{e}^{i{k}_2(x+x^{\prime })}\right](x<x^{\prime })\end{array}.$$

$$G^{\mathrm{R}}(x,x+0,k_y)-G^{\mathrm{R}}(x,x-0,k_y)=i{v}_{\mathrm{F}}^{-1}\left(\begin{array}{cc}{\sigma }_y& 0\\ 0& -{\sigma }_y\end{array}\right),$$

$$g_1$$

$$g_4$$

$$\rho (x,E)$$

$$=\int \frac{d{k}_y}{2\pi }\rho (x,k_y,E)$$

$$=-\frac{1}{\pi }\int \frac{d{k}_y}{2\pi }\mathrm{Im}\left[G_{11}^{\mathrm{R}}(x,x,k_y,E)+G_{22}^{\mathrm{R}}(x,x,k_y,E)\right].$$

$$⟨j⟩$$

$$⟨S⟩$$

$$I$$

$$=\frac{\Delta (e^{i{\theta }_1}{\zeta }_2+e^{i{\theta }_2}{\zeta }_1)}{k_1-k_2}\left(\frac{e^{-i{\theta }_1}}{{\zeta }_1\cos {\theta }_1}{a}_1-\frac{e^{-i{\theta }_2}}{{\zeta }_2\cos {\theta }_2}{a}_2\right),$$

$$I=e{k}_{\mathrm{B}}T\sum _{k_y,i{\omega }_n}\frac{{\Delta }_z}{\sqrt{{\omega }_n^2+{\Delta }_z^2}}(a_1-a_2).$$

$${\hat{\psi }}_1(x)$$

$${\hat{\psi }}_2(x)$$

$${\hat{\psi }}_3(x)$$

$${\hat{\psi }}_4(x)$$

$$\tilde{H}({\mathbf{x}}^{\prime },{\partial }_{{\mathbf{x}}^{\prime }})=\left(\begin{array}{cc}h({\partial }_{x^{\prime }},-{\partial }_{y^{\prime }})+M^{\ast }& i{\sigma }_y\Delta (x^{\prime })\\ -i{\sigma }_y\Delta (x^{\prime })& -h^{\ast }({\partial }_{x^{\prime }},-{\partial }_{y^{\prime }})-M\end{array}\right)$$

$${\hat{B}}_1=\left(\begin{array}{c}i{e}^{-i{\theta }_1}\\ {\zeta }_1\\ -{\zeta }_1{\gamma }_{-}\\ i{e}^{-i{\theta }_1}{\gamma }_{+}\end{array}\right);{\hat{B}}_2=\left(\begin{array}{c}i{\zeta }_2{\gamma }_{-}\\ -e^{-i{\theta }_2}{\gamma }_{+}\\ e^{-i{\theta }_2}\\ i{\zeta }_2\end{array}\right);{\hat{B}}_3=\left(\begin{array}{c}i\\ -{\zeta }_1{e}^{-i{\theta }_1}\\ {\zeta }_1{e}^{-i{\theta }_1}{\gamma }_{-}\\ i{\gamma }_{+}\end{array}\right);{\hat{B}}_4=\left(\begin{array}{c}i{\zeta }_2{e}^{-i{\theta }_2}{\gamma }_{-}\\ {\gamma }_{+}\\ -1\\ i{\zeta }_2{e}^{-i{\theta }_2}\end{array}\right).$$

$${\hat{\tilde{\psi }}}_1(x^{\prime })$$

$${\hat{\tilde{\psi }}}_2(x^{\prime })$$

$${\hat{\tilde{\psi }}}_3(x^{\prime })$$

$${\hat{\tilde{\psi }}}_4(x^{\prime })$$

$$G^{\mathrm{R}}(x,x^{\prime })=\{\begin{array}{ll}{\alpha }_1{\hat{\psi }}_1(x){\hat{\tilde{\psi }}}_3^{\mathsf{T}}(x^{\prime })+{\alpha }_2{\hat{\psi }}_1(x){\hat{\tilde{\psi }}}_4^{\mathsf{T}}(x^{\prime })+{\alpha }_3{\hat{\psi }}_2(x){\hat{\tilde{\psi }}}_3^{\mathsf{T}}(x^{\prime })+{\alpha }_4{\hat{\psi }}_2(x){\hat{\tilde{\psi }}}_4^{\mathsf{T}}(x^{\prime })& (x>x^{\prime })\\ {\beta }_1{\hat{\psi }}_3(x){\hat{\tilde{\psi }}}_1^{\mathsf{T}}(x^{\prime })+{\beta }_2{\hat{\psi }}_4(x){\hat{\tilde{\psi }}}_1^{\mathsf{T}}(x^{\prime })+{\beta }_3{\hat{\psi }}_3(x){\hat{\tilde{\psi }}}_2^{\mathsf{T}}(x^{\prime })+{\beta }_4{\hat{\psi }}_4(x){\hat{\tilde{\psi }}}_2^{\mathsf{T}}(x^{\prime })& (x<x^{\prime })\end{array}.$$

$${\alpha }_1{\tilde{c}}_3+{\alpha }_2{\tilde{d}}_4={\beta }_1{c}_3+{\beta }_2{d}_4$$

$${\alpha }_3{\tilde{d}}_3+{\alpha }_4{\tilde{c}}_4={\beta }_3{d}_3+{\beta }_4{c}_4$$

$${\alpha }_3{\tilde{c}}_3+{\alpha }_4{\tilde{d}}_4={\alpha }_3{\tilde{c}}_3+{\alpha }_4{\tilde{d}}_4$$

$$({\alpha }_1{\tilde{c}}_3+{\alpha }_2{\tilde{d}}_4)a_1$$

$$({\alpha }_3{\tilde{d}}_3+{\alpha }_4{\tilde{c}}_4)a_2$$

$$({\alpha }_1{\tilde{c}}_3+{\alpha }_2{\tilde{d}}_4)b_1$$

$$({\alpha }_3{\tilde{d}}_3+{\alpha }_4{\tilde{c}}_4)b_2$$

Full text

I Introduction

It is known that the surface state of a topological insulator (TI) is expressed by the massless Dirac electron with spin-momentum locking Hasan and Kane (2010). This state is protected by time reversal symmetry and can be gapped by magnetization. Also, if we put a superconductor on the surface of a TI, the energy gap opens due to the induced pair potential in this surface state Fu and Kane (2008). It has been predicted that chiral Majorana modes are generated at the boundary between the ferromagnetic insulator (FI) and the superconductor (SC) in FI/SC junctions on the surface of TIs Fu and Kane (2009). Then, the tunneling effect and the Josephson effect via chiral Majorana modes in FI/SC junctions have been theoretically studied Akhmerov et al. (2009); Linder et al. (2009); Tanaka et al. (2009); Linder et al. (2010); Lu et al. (2015). There are also several theoretical works about charge transport in superconducting junctions on a TI Tkachov and Hankiewicz (2013); Potter and Fu (2013); Vali and Khouzestani (2015); Vasenko et al. (2017a). Besides these researches, odd-frequency pairings have been predicted to emerge in such systems Black-Schaffer and Balatsky (2013); Burset et al. (2015); Cayao and Black-Schaffer (2017); Vasenko et al. (2017b); Bobkova et al. (2016).

Up to now, Josephson junctions on the surface of TIs Veldhorst et al. (2012); Williams et al. (2012); Sacépé et al. (2011); Cho et al. (2013); Galletti et al. (2014); Sochnikov et al. (2015); Molenaar et al. (2014); Charpentier et al. (2017) have been reported, where Bi2Se3, doped Bi2Se3, Bi2Te3, and strained HgTe were used as the TIs. The observed nonsinusoidal current-phase relation of the Josephson current is consistent with the existence of helical edge modes which is a Kramers pair of Majorana modes Sochnikov et al. (2015). However, it is still difficult to simultaneously stack a ferromagnet and a superconductor on a TI, and FI/SC junctions on the surface of TIs, which supports chiral edge modes, have not been experimentally realized yet.

A possible alternative to the FI/SC junction is a superconducting junction on a magnetically doped TI. There are several theoretical works on systems with coexisting magnetization and superconducting pairing in the same spatial region Burset et al. (2015); Goudarzi et al. (2017); Khezerlou et al. (2017, 2018); Wang (2018). Chiral Majorana modes are predicted to appear in such systems. In this paper, we consider the simplest configuration created just by fabricating superconducting islands on a magnetically doped TI as shown in Fig. 1, and discuss how the existence of the chiral Majorana modes affects the Josephson current. Recently, there is a report on the experimental realization of Josephson junctions on magnetically doped TIs Yano et al. (2019, 2018). In magnetically doped TIs, the energy dispersion of the surface Dirac electron has a gap opening when the direction of the magnetization is perpendicular to the surface and its magnitude is larger than the chemical potential of surface Dirac electron $\mu$. In other words, by tuning the chemical potential, metallic (insulating) surface state is realized for $\mu>|m_{z}|$ ($|m_{z}|>\mu$) for $\mu>0$ where $|m_{z}|$ is the magnitude of the out of plane magnetization. At the same time, if we put a superconductor with pair potential $\Delta$ on the magnetically doped TI, topological superconducting state is generated for $|m_{z}|<\sqrt{\mu^{2}+\Delta^{2}}$ Alicea (2012). Our focus is on the case realized for $\mu<|m_{z}|<\sqrt{\mu^{2}+\Delta^{2}}$ where the chiral Majorana mode is localized at the edge of superconducting region on the TI. It is a challenging issue to clarify how the Josephson current flows via chiral Majorana modes in this regime since it is in the insulating phase in the normal state.

Returning to the history of the investigation of the Josephson current Golubov et al. (2004), the Josephson current shows a clear dependence on the presence or absence of zero energy surface Andreev bound states (ZESABSs) Kashiwaya and Tanaka (2000); Löfwander et al. (2001). In the junctions with low transmittance at the interface, the temperature dependence of the Josephson current obeys the Ambegaokar-Baratoff relation in the absence of ZESABSs Ambegaokar and Baratoff (1963). The maximum Josephson current $I_{c}$ saturates with the decrease of temperature $T$, and it is inversely proportional to $R_{N}$, where $R_{N}$ is the resistance in the normal state Ambegaokar and Baratoff (1963). On the other hand, in the presence of ZESABSs like the case of $d$-wave superconductor junctions, $I_{c}$ is proportional to $1/T$ in an intermediate temperature regime and proportional to $1/\sqrt{R_{N}}$ at low temperature Tanaka and Kashiwaya (1996, 1997); Barash et al. (1996). In superconducting junctions on magnetically doped TIs, one can create or annihilate Majorana modes by tuning the chemical potential on the surface of the TI. Then, we expect drastic change of the magnitude and the temperature dependence of the Josephson current at the topological phase transition point since Majorana modes are nothing but ZESABSs appearing at the edge of the superconducting regions.

There has been recent development of a theoretical methodology to calculate the Josephson current on a TI based on Green’s function of the two-dimensional (2D) Dirac electron system Lu and Tanaka (2018). This formulation is an extension of McMillan’s formulation originally developed for the free electron model with parabolic dispersion McMillan (1968). In this formalism, the Josephson current is expressed by Andreev reflection coefficient like Furusaki-Tsukada formula Furusaki and Tsukada (1991), which has been developed for $s$-wave superconductor junctions. By an analytic continuation of the Andreev reflection coefficient obtained in the regime of metallic normal state ($|m_{z}|<\mu$), it is possible to calculate the Josephson current even when the normal state is insulating ($|m_{z}|>\mu$). Thus, it is timely to calculate the Josephson current using this formulation in superconducting junctions on magnetically doped TIs.

In this paper, we consider a Josephson junction on the surface of a magnetically doped TI. We analytically calculate Green’s function on the surface of the TI available for various values of the chemical potential, which drives a topological phase transition in the induced superconducting states and a metal-insulator transition in the surface states. Depending on the relative magnitude of $m_{z}$ to $\mu$, there are three phases: (I) trivial superconductivity and insulating in the normal state for $|m_{z}|>\sqrt{\mu^{2}+\Delta^{2}}$, (II) topological superconductivity and insulating in the normal state for $\sqrt{\mu^{2}+\Delta^{2}}>|m_{z}|>\mu$, and (III) topological superconductivity and metallic in the normal state for $\mu>|m_{z}|$ (see Fig. 2). As expected, the magnitude of the Josephson current in Phase (III) is much larger than those in Phases (I) and (II) since the normal state is metallic (insulating) for the former (latter). Besides, we find that the temperature dependence of the Josephson current drastically changes at the phase boundary between Phases (I) and (II), namely, at the topological phase transition point. We focus on the enhancement of the Josephson current at low temperature and evaluate it by the ratio $I_{c}(T=0)/I_{c}(T=0.2T_{c})$, where $I_{\mathrm{c}}$ is the maximum Josephson current as a function of temperature $T$ and $T_{\mathrm{c}}$ is the superconducting transition temperature. The ratio is almost unity in Phase (I), shows rapid increase in Phase (II), and remains large in Phase (III). We attribute the enhancement of the ratio to the chiral Majorana modes appearing in Phases (II) and (III).

The organization of this paper is as follows. In Sec. II, we introduce the model and formulation. In Sec. III, we show the local density of states. In Sec. IV, we show the current-phase relation and temperature dependence of the Josephson current. In Sec. V, we summarize our results.

II Model and formulation

We consider a Josephson junction of two spin-singlet $s$-wave superconductors deposited on the surface of a magnetically doped TI shown in Fig. 1. The superconductors provide pair potential in the surface state of the TI via proximity effect while the magnetic moment produces the Zeeman coupling through exchange interaction. Therefore, the Hamiltonian of our interest can read

[TABLE]

where $h(\partial_{\bm{x}})=v_{\mathrm{F}}(-i\partial_{y}\sigma_{x}+i\partial_{x}\sigma_{y})-\mu$ in the diagonal part is the surface Dirac Hamiltonian with the Fermi velocity $v_{\mathrm{F}}$, the vector differential operator $\partial_{\bm{x}}=(\partial_{x},\partial_{y})=(\partial/\partial x,\partial/\partial y)$, and the chemical potential $\mu$. $M=\bm{m}\cdot\bm{\sigma}$ is the Zeeman coupling with $\bm{\sigma}$ being the Pauli matrix vector acting on the spin space. The magnetic moment are finite and fixed in all the region, and we assume it points perpendicular to the surface, i.e., $\bm{m}=(0,0,m_{z})$. Without loss of generality, we assume $m_{z}>0$. The off-diagonal part in Eq. (1) is the proximity induced pair potential of the spin-singlet $s$-wave form. To consider the junction system depicted in Fig. 1 we employ $\bar{\Delta}(x)=\Delta(\Theta(-x)+e^{-i\varphi}\Theta(x-L))$ with the interval of the junction $L$ and the phase difference between two superconductors $\varphi$. Although, in general, the induced pair potential depends on the detail of the hybridization and states which the superconductor couples to, we approximate the pair potential has a step-function profile $\Theta(x)$ with a constant magnitude $\Delta$. We choose $\Psi=(c_{\bm{k},\uparrow},c_{\bm{k},\downarrow},c_{-\bm{k},\uparrow}^{\dagger},c_{-\bm{k},\downarrow})^{\mathsf{T}}$ as the basis.

Figure 2 illustrates the phase diagram of the system. The Zeeman term produces a magnetic gap in the surface state. Phases in the normal state or in the region uncovered with the superconductors are shown in Fig. 2 (a). It becomes metallic when the Fermi energy lies above the magnetic gap and becomes insulating otherwise. Figure 2 (b) represents a trivial-topological phase transition in the region covered with the superconductor. According to Ref. Fu and Kane, 2008, topologically nontrivial superconducting states, which are equivalent to a spinless $p+ip$ superconductor, are induced on the surface state coupled to a superconductor. When $0<m_{z}<\sqrt{\Delta^{2}+\mu^{2}}$, the induced pair potential mediates time reversal broken $p$-wave superconductivity Alicea (2012). The bulk gap closes along the solid line in the figure, and in the trivial regime, the gap is dominated by the time reversal breaking Zeeman coupling. We superimpose these two diagrams and obtain Fig. 2 (c). It consists three phases. Along the arrow in the figure: (I) Josephson junction of topologically trivial superconducting states with an insulating interval; (II) Josephson junction of topological superconductors with an insulating interval; and (III) Josephson junction of topological superconductors with a metallic interval. In general, topological states host gapless edge modes at their boundary. In the present case, Phases (II) and (III) support Majorana edge modes. Since time reversal symmetry is explicitly broken by the Zeeman coupling, the edge modes are chiral, i.e., a chiral Majorana mode appears at each edge of the two superconductors facing each other across the junction. However, since the interval of the junction is metallic in Phase (III), the Majorana modes cannot be localized at the edge of the superconductors, and, instead, spread into the interval area. There are also number of propagating modes in the interval which can mediate the current across the junction. In this case, the contribution from Majorana modes to the current becomes less emphasized. The signature of the presence of Majorana modes is strongly manifested in Phase (II). These are explicitly shown in the following sections.

We will follow McMillan’s formula McMillan (1968); Kashiwaya and Tanaka (2000) for constructing Green’s function to calculate the density of states and the Josephson current. The calculation is performed along Ref. Lu and Tanaka, 2018, where they discuss a similar problem but without magnetization in the superconducting regions. We here briefly summarize the strategy. In general, Green’s function satisfying

[TABLE]

is constructed by the outer product of a solution of the homogeneous differential equation $\left[E-H(\bm{x},\partial_{\bm{x}})\right]\hat{\psi}(\bm{x})=0$ and that for the conjugate equation $\left[E-\tilde{H}(\bm{x}^{\prime},\partial_{\bm{x}^{\prime}})\right]\hat{\tilde{\psi}}(\bm{x}^{\prime})=0$, where $\tilde{H}(\bm{x},\partial_{\bm{x}})\equiv H^{\mathsf{T}}(\bm{x},-\partial_{\bm{x}})$ is the adjoint operator of $H(\bm{x},\partial_{\bm{x}})$ with $\mathsf{T}$ denoting the matrix transpose. In the present case with the translational symmetry along the $y$ direction, the Fourier transform of $G(\bm{x},\bm{x}^{\prime})$ with respect to $y$ is given by

[TABLE]

where $\hat{\psi}_{\pm}(x)\ (\hat{\tilde{\psi}}_{\pm}(x^{\prime}))$ is the eigenstate of $H|_{\partial_{y}=ik_{y}}$ ($\tilde{H}|_{\partial_{y}=-ik_{y}}$) that satisfies the boundary condition at $x=\pm\infty$. For the retarded (advanced) Green’s function $G^{\mathrm{R(A)}}(x,x^{\prime},k_{y})$, we impose the out-going (in-coming) boundary condition, i.e., we adopt right-going (left-going) electron solution and left-going (right-going) hole solution at $x=+\infty$ and left-going (right-going) electron solution and right-going (left-going) hole solution at $x=-\infty$. When there are several independent solutions for each of $\hat{\psi}_{\pm}(x)$ and $\hat{\tilde{\psi}}_{\pm}(x)$, the Green’s function is written as a linear combination of the outer products of all possible pairs of eigenstates. The coefficients of the outer products will be determined so as to satisfy the boundary condition at $x=x^{\prime}$. This boundary condition ensures that the obtained Green’s function satisfies Eq. (2). For the case of a Josephson junction, by expressing $\hat{\psi}_{\pm}(x)$ and $\hat{\tilde{\psi}}_{\pm}(x^{\prime})$ using the Andreev-reflection amplitude, we analytically obtain the Green’s function and the resulting local density of states and the Josephson current in terms of these parameters.

Now we move to the concrete calculation for our system given by Eq. (1). To obtain $\hat{\psi}_{\pm}(x)$, we first solve the eigen-solution of $H|_{\partial_{y}=ik_{y}}$. There are four propagating modes with energy $E$; right-going (in the $x$ direction) electron mode $\hat{A}_{1}e^{ik_{1}x+ik_{y}y}$, and right-going hole mode $\hat{A}_{2}e^{-ik_{2}x+ik_{y}y}$, left-going electron mode $\hat{A}_{3}e^{-ik_{1}x+ik_{y}y}$, and left-going hole mode $\hat{A}_{4}e^{ik_{2}x+ik_{y}y}$. The explicit form of the eigenstates of Eq. (1) in the $x<0$ region are given as

[TABLE]

where $\zeta_{j}$, $\theta_{j}$ ($j=1,2$), and $\gamma_{\pm}$ are defined by

[TABLE]

with

[TABLE]

To obtain normalizable wavefunctions, we need to take $\mathrm{Im}k_{1}\geq 0$ and $\mathrm{Im}k_{2}\leq 0$. $\Delta_{z}$ is introduced for the short notation, and it also gives an effective gap in the quasiparticle spectrum when $m_{z}<\mu^{2}/\sqrt{\Delta^{2}+\mu^{2}}$. The lowest unoccupied band in the spectrum has double-minima at $k_{\mathrm{min}}=\pm\sqrt{\mu^{2}-m_{z}^{2}(1+\Delta^{2}/\mu^{2})}$ with the gap size $\Delta_{z}$. By increasing $m_{z}$, two minima approach each other, and at $m_{z}=\mu^{2}/\sqrt{\Delta^{2}+\mu^{2}}$ they join up at $k=0$. Then the energy gap becomes $(m_{z}^{2}+\Delta^{2}+\mu^{2}-2m_{z}\sqrt{\Delta^{2}+\mu^{2}})^{1/2}$, which finally closes at $m_{z}=\sqrt{\Delta^{2}+\mu^{2}}$ (this is the topological phase transition point in Fig. 2 (b)). At $m_{z}>\mu$, $\Delta_{z}$ becomes pure imaginary, i.e., $\Delta_{z}^{2}<0$ in the definition of $\Omega_{z}$.

By solving the scattering problem with the Josephson junction, we obtain two independent solutions that satisfy the out-going boundary condition at $x=+\infty$. Their analytic forms for $x<0$ is given by $\hat{\psi}_{1}(x)=\hat{A}_{1}e^{ik_{1}x}+a_{1}\hat{A}_{4}e^{ik_{2}x}+b_{1}\hat{A}_{3}e^{-ik_{1}x}$ and $\hat{\psi}_{2}(x)=\hat{A}_{2}e^{-ik_{2}x}+a_{2}\hat{A}_{3}e^{-ik_{1}x}+b_{2}\hat{A}_{4}e^{ik_{2}x}$, where $a_{1(2)}$ and $b_{1(2)}$ are the Andreev-reflection and normal reflection-amplitudes for electron (hole), respectively. In a similar manner, we calculate two independent solutions for each of $\hat{\psi}_{-}(x)$, $\hat{\tilde{\psi}}_{+}(x^{\prime})$, and $\hat{\tilde{\psi}}_{-}(x^{\prime})$, obtaining the retarded Green’s function for $x,x^{\prime}<0$ as

[TABLE]

Here, $\hat{B}_{i}\,(i=1,2,3,4)$, $\tilde{a}_{1,2}$, and $\tilde{b}_{1,2}$ are the counterparts of $\hat{A}_{i}$, $a_{1,2}$, and $b_{1,2}$, respectively, in the conjugate process described by $\tilde{H}|_{\partial_{y}=-ik_{y}}$. More details are given in Appendix A. In the derivation of Eq. (14), we have used the boundary condition at $x=x^{\prime}$:

[TABLE]

obtained by integrating Eq. (2) over $x^{\prime}$ in the section $[x-0,x+0]$, and $g_{1}$ and $g_{4}$ are given by

[TABLE]

III Local density of states

The local density of states (LDOS) is calculated by

[TABLE]

LDOS at the edge of the left superconductor resolved with $k_{y}$, $\rho(x=0,k_{y},E)$, gives the dispersion of modes. The normal and Andreev coefficients appearing in Green’s function in Eq. (14) are numerically obtained. They are plotted in Fig. 3 for several sets of $m_{z}$ and the length of the interval $L$. The phase difference between two superconductors and the chemical potential are set to $\varphi=\pi$ and $\mu/\Delta=5.0$, respectively. In Fig. 3, the panels belong to Phases (I), (II), and (III) from left to right. There are no states observed at zero energy in Fig. 3 (a) and (d) whereas there exist modes with linear dispersion at low energy regime in Figs. 3(b), (c), and (e). The presence of in-gap states is the consequence of topological property. From Fig. 3 (b), it seems two modes propagate to the $+y$ and $-y$ direction. This is because we have two edges of superconductors at $x=0$ and $x=L$, and at each edge one mode is localized. The tail of the wavefunction of the mode at $x=L$ is still large enough at $x=0$ when the length of the interval is comparable with or less than the coherence length $\xi$. This point is explicitly confirmed by looking at the LDOS in an SN junction shown in Figs. 3 (d-f), which is equivalent to taking $L\to\infty$ in our setting. Now, only the mode propagating to the $+y$ direction remains and the other part disappears (Fig. 3 (e)) which means we have one chiral Majorana mode localized at $x=0$ Fu and Kane (2008); Tanaka et al. (2009). The LDOS of the SN junction also reveals a clear distinction between Phases (II) and (III). When the junction is metallic, there are number of modes at the Fermi energy (see Fig. 3 (f)). The chiral Majorana modes potentially localized at the edges are no longer distinguishable. It is hybridized with continuous modes in the junction area and spreads into there.

IV Josephson current

We move on to the Josephson current. The analytic continuation $E+i\delta\to i\omega_{n}$ transforms the retarded Green’s function to the Matsubara Green’s function, where $\omega_{n}$ is given by $\omega_{n}=2\pi k_{B}T(n+1/2)$ with an integer $n$. The observables are expressed in the sum over the Matsubara frequency. The Josephson current is derived to be $I=\langle j\rangle+\int^{x}_{0}dx\langle S\rangle$ with the current operator $j=iev_{\mathrm{F}}(c_{\downarrow}^{\dagger}c_{\uparrow}-c_{\uparrow}^{\dagger}c_{\downarrow})$ and the source term operator $S=ie\Delta(c_{\downarrow}c_{\uparrow}\allowbreak-c_{\uparrow}^{\dagger}c_{\downarrow}^{\dagger}\allowbreak-c_{\uparrow}c_{\downarrow}\allowbreak+c_{\downarrow}^{\dagger}c_{\uparrow}^{\dagger})$, where $c_{s}^{\dagger}$ ($c_{s}$) is the Fourier transform of the creation (annihilation) operator $c_{\bm{k},s}$ of an electron with spin $s=\uparrow,\downarrow$ with respect to $x$. They can be calculated with respect to the Green’s function in the superconductor in $x<0$ as

[TABLE]

After calculation represented in Appendix B, the Josephson current is shown to be

[TABLE]

where we use $k_{1}^{2}-k_{2}^{2}=4i\mu\sqrt{\omega_{n}^{2}+\Delta_{z}^{2}}$ to get the second equation. In the latter expression, it is obvious that there is no singularity even in the case with $\Delta_{z}$ being a complex number since $\mathrm{Im}k_{1}>0$ and $\mathrm{Im}k_{2}<0$ as we mentioned.

When the chemical potential is the dominating energy scale, $\mu\gg|E|,\Delta,m_{z}$, so that the quasiclassical approximation is valid, Eq. (24) can be approximated in a simple form of

[TABLE]

This is an extension of Furusaki-Tsukada formula applied for the Josephson current on magnetically doped topological insulators.

The current-phase relation calculated from Eq. (24) is shown in Fig. 4. The temperature dependence of the pair potential is assumed to be $\Delta(T)=\Delta_{0}\tanh[1.74\sqrt{T_{\mathrm{c}}/T-1}]$. The current is normalized by $\Delta_{0}/eR_{\mathrm{N}0}$ where $R_{\mathrm{N}0}$ is the resistivity of the junction without magnetization in the normal state (We cannot use the resistivity of the junction at $T>T_{c}$ for the normalization, because the normal state is insulating at $m_{z}>\mu$.) The current-phase relation in Phase (I) has an ordinary sinusoidal form and little temperature dependence below $T/T_{\mathrm{c}}<0.1$ (Fig. 4 (a)). This is consistent with Ambegaokar-Baratoff relation, $\sin\varphi\tanh[\Delta(T)/2k_{\mathrm{B}}T]$, obtained in a junction with low transmittance Ambegaokar and Baratoff (1963). In Phase (II) (Fig. 4 (b)), the behavior of the current-phase relation is nearly in the form of a sinusoidal function. When the temperature goes down, however, the maximum Josephson current grows and its position is slightly away from $\varphi=\pi/2$. These properties are different from the naively expected Ambegaokar-Baratoff relation. Based on the result of LDOS (Fig. 3 (b)), these behavior can be attributed to the chiral Majorana modes arising in the superconducting gap. On the other hand, in Phase (III) (Fig. 4 (c)), the maximum current flows when the phase difference is close to $\pi$. The behavior is similar to superconductor-short normal metal-superconductor junctions with fully transparent interface, known as Kulik-Omelyanchuk relation Kulik and Omelyanchuk (1977), or unconventional superconductor-short junctions with low transmission in the presence of ZESABSs Tanaka and Kashiwaya (1996, 1997). Figure 5 shows the temperature dependence of the maximum value of the Josephson current $I_{\mathrm{c}}$. The maximum current saturates near zero temperature in all cases. It should be noted again that the surface state without superconductor attached is insulating in Phases (I) and (II). It is then natural that the saturated value in Phases (III) is much larger than those in Phases (I) and (II). If we closely look at the temperature dependence, however, we notice that the saturation occurs relatively higher temperature when the superconducting state is trivial, whereas the maximum current keeps getting larger till close to zero temperature when the superconducting state becomes topological.

To extract these increase at low temperature, we take the ratio between the current at zero temperature and at a reference temperature, which we choose as $T=0.2T_{\mathrm{c}}$. Figure 6 shows the ratio as a function of the chemical potential. By increasing the chemical potential, the system goes through two transition points around $\mu/\Delta_{0}=1.75$ and $2.0$. When the surface state is insulating, i.e., with small $\mu$, the current shows saturation below the reference temperature. The ratio significantly increases when the superconducting regions become topological and chiral Majorana modes appear. It keeps increasing until the system is in the metallic phase and then decreases converging to a constant. To emphasize the contribution of Majorana modes, we also calculate the ratio in the one dimensional space, which is illustrated in Fig. 7. This case corresponds to taking only the contribution from $k_{y}=0$ in 2D case. Now the peak locates in Phase (II), which indicates that the increase of the Josephson current is mostly attributed to the presence of Majorana modes not to continuous modes in metallic surface states. The shift of the peak in 2D space can be understood as follows. As shown in Fig. 3, the energy of chiral Majorana modes depends on the wave number parallel to the junction, and it is on $E=0$ only when $k_{y}=0$. Therefore, at low temperature comparable to the finite energy of chiral Majorana modes, the contribution of Majorana modes to the Josephson current is overwhelmed by the contribution from the continuous modes.

V Summary

In this paper, we have studied charge transport of the Josephson junction on the surface of magnetically doped TIs. We have analytically obtained the Green’s function on the surface of TIs available for the various chemical potential, which drives a topological phase transition in the induced superconducting states and a metal-insulator transition in the surface states of 2D Dirac electron. There are three parameter regions: (I) trivial superconductivity with insulating normal state for $|m_{z}|>\sqrt{\mu^{2}+\Delta^{2}}$, (II) topological superconductivity with insulating normal state for $\sqrt{\mu^{2}+\Delta^{2}}>|m_{z}|>\mu$, and (III) topological superconductivity with metallic normal state for $\mu>|m_{z}|$. It is remarkable that non-zero Josephson current is obtained even for (I) and (II) though it is rather small. By focusing on the ratio $I_{c}(T=0)/I_{c}(T=0.2T_{c})$, we have distinguished three regions clearly. This ratio is almost unity in Phase (I), suddenly increases in Phase (II), and remains enhanced in Phase (III). Figure 6 tells that the existence of chiral Majorana states can be detected by an experimental observation of a prominent increase in the Josephson current at low temperature in the insulating phase (Phase (II)). The current available experimental technique would allow us to tune the chemical potential widely enough to go across these phase boundaries. Thus, the observation of the temperature dependence of the Josephson current provides detection of Majorana modes.

In this paper, we have calculated Josephson current in the ballistic regime. In the actual surface state of TIs, the impurity scattering effect also influences on the charge transport. It is interesting to clarify the possible odd-frequency spin-triplet $s$-wave pairing which is robust against the impurity scattering Tanaka and Golubov (2007); Tanaka et al. (2012). In addition to this point, one of the remarkable finding of this paper is the presence of Josephson current in the regime where electronic property of surface state of TIs in the normal state is insulating. There are several studies on charge transport in superconducting quantum anomalous Hall system where the normal state is insulating. One can see our model as the proximity induced superconducting state in an anomalous quantum Hall state. The magnetic gap in the surface Dirac state of a TI provides a half quantized Chern number $\sigma=1/2$, which have been intensively investigated both in theories Qi et al. (2010); Ii et al. (2011, 2012); Yamakage and Sato (2014); He et al. (2014) and experiments He et al. (2017). Besides the setup shown in Fig. 1, our system may be realized at the interface between a 2D TI and an FI with a superconducting junction deposited on it, or same configuration at hinge states in higher order TIs Benalcazar et al. (2017).

Acknowledgements.

The authors would like to thank R. Yano, S. Kashiwaya, B. Lu, and P. Burset for fruitful discussion. This work was supported by JST CREST Grant No. JPMJCR16F2, and JSPS KAKENHI Grant Nos. JP15H05851, JP15H05853, JP15K17726, JP18H01176, and JP19H01824.

Appendix A McMillan’s Green’s function method

Here we give all the detail of the calculation for constructing McMillan’s Green’s function.

A.1 Scattering solutions

Solving the scattering problem with the Josephson junction results in the following four independent wavefunctions in the $x<0$ region Burset et al. (2015):

[TABLE]

The former (latter) two satisfy the out-going boundary condition at $x=+\infty$ ($x=-\infty$). $c_{3}$ ($d_{4}$) and $d_{3}$ ($c_{4}$) are the transmission coefficients for electron and hole, respectively, when electron (hole) is injected from the right. On the other hand, the eigenstates of the conjugate processes obtained by solving

[TABLE]

are in the form of $\hat{B}_{1}e^{ik_{1}x^{\prime}-ik_{y}y^{\prime}},\hat{B}_{2}e^{-ik_{2}x^{\prime}-ik_{y}y^{\prime}},\hat{B}_{3}e^{-ik_{1}x^{\prime}-ik_{y}y^{\prime}},\hat{B}_{4}e^{ik_{2}x^{\prime}-ik_{y}y^{\prime}}$ with

[TABLE]

The scattering solutions satisfying the out-going boundary condition at $x=+\infty$ ($x=-\infty$), denoted by $\hat{\tilde{\psi}}_{1}(x^{\prime})$ and $\hat{\tilde{\psi}}_{2}(x^{\prime})$ ($\hat{\tilde{\psi}}_{3}(x^{\prime})$ and $\hat{\tilde{\psi}}_{4}(x^{\prime})$), are given in the $x^{\prime}<0$ region as

[TABLE]

where $\tilde{c}_{3}$ ($\tilde{d}_{4}$) and $\tilde{d}_{3}$ ($\tilde{c}_{4}$) are the transmission coefficients for electron and hole, respectively, when electron (hole) is injected from the right.

A.2 Determining the coefficients

Green’s function is constructed by making the linear combination of the outer products of the above obtained solutions, resulting in McMillan (1968); Lu and Tanaka (2018)

[TABLE]

The boundary condition of Green’s function at $x=x^{\prime}$ (Eq. (15)) gives following set of equations

[TABLE]

The first two equations give the definitions for $g_{1}$ and $g_{4}$ in the main text, respectively. Combining these equations, the Green’s function is derived in the form of Eq. (14).

Appendix B Josephson current

The calculation of Josephson current is done with evaluating the following quantities by using Green’s function defined in the left superconducting region ($x<0$). The current satisfies the continuity equation

[TABLE]

where $\rho$, $j$, and $S$ are the charge density, the current, and the source term operators, respectively. They are defined as

[TABLE]

As we defined in the main text, $c_{s}^{\dagger}$ ($c_{s}$) is the Fourier transform of the creation (annihilation) operator $c_{\bm{k},s}$ of an electron with spin $s=\uparrow,\downarrow$ with respect to $x$. The expectation values of them at $x$ are evaluated in terms of Green’s function with the summation over Matsubara frequency and the wave number along the $y$ direction:

[TABLE]

[TABLE]

and

[TABLE]

obtaining

[TABLE]

The same calculation with the advanced Green’s function should be performed for the summation $\omega_{n}<0$ and it turns out to be the same as the above equation. Therefore, the Josephson current is obtained as Eq. (24) in the main text. In the quasiclassical limit, the variables appearing in the eigenstates are approximated as

[TABLE]

resulting in Eq. (25) in the main text.