Superconducting Diode Effect in Topological Hybrid Structures

Abstract

Currently, the superconducting diode effect (SDE) is being actively discussed, due to its large application potential in superconducting electronics. In particular, superconducting hybrid structures, based on three-dimensional (3D) topological insulators, are among the best candidates, due to their having the strongest spin–orbit coupling (SOC). Most theoretical studies on the SDE focus either on a full numerical calculation, which is often rather complicated, or on the phenomenological approach. In the present paper, we compare the linearized and nonlinear microscopic approaches in the superconductor/ferromagnet/3D topological insulator (S/F/TI) hybrid structure. Employing the quasiclassical Green’s function formalism we solve the problem self-consistently. We show that the results obtained by the linearized approximation are not qualitatively different from the nonlinear solution. The main distinction in the results between the two methods was quantitative, i.e., they yielded different supercurrent amplitudes. However, when calculating the so-called diode quality factor the quantitative difference is eliminated and both approaches result in good agreement.

1. Introduction

The field of superconducting electronics is an important area in the research and development of hybrid quantum devices with lower power consumption. Superconducting hybrid structures, consisting of a superconductor and non-superconducting material (normal metal N, ferromagnet F, etc.), operate by means of the proximity effect. This effect can be described as a leakage of the superconducting correlations into the adjacent non-superconducting layer [1,2,3,4,5,6,7,8,9,10,11]. Superconductor/ferromagnet (S/F) structures have been proposed in many nanoelectronic applications, such as the following: memory devices [12], quantum and classical logic devices [12,13], artificial neural networks [14], detectors and bolometers [15], nanorefrigerators [16,17] and spin-valves [18]. Placing two-dimensional (2D) S/F structures on the surface of a 3D topological insulator, a material with strong spin–orbit coupling, may add new functionality and create superconducting diode, see Figure 1.

The superconducting diode effect (SDE) is an active area of research because of its great application potential in the fields of superconducting electronics and spintronics. Generally, the SDE is observed in two-dimensional superconducting systems with broken inversion and time reversal symmetries [19]. While the former usually implies the presence of the spin–orbit field, the latter can be achieved by the exchange field from the ferromagnet, or by exposing the system to an external magnetic field. These conditions allow for the possibility of superconducting helical state realization. The helical state is characterized by the order parameter modulated in the direction transverse to the field by the phase factor $exp(i{\mathit{q}}_{\mathbf{0}}·\mathit{r})$, where ${\mathit{q}}_{\mathbf{0}}$ is the Cooper pair momentum [20]. In this case, the critical supercurrents are different in the directions parallel to and anti-parallel to the momentum ${\mathit{q}}_{\mathbf{0}}$. Significant advances have been made since the experimental discovery of the diode effect by Ando et al. [21]. There have been numerous reports involving both experimental [21,22,23,24,25,26,27] and theoretical studies [28,29,30,31,32,33,34,35,36] of the SDE. Hybrid SDE devices deserve special attention [33,37,38]. In such structures, the ingredients for the SDE effect are brought together by the proximity effect. For instance, the S/F/TI hybrid structure is a promising platform for realization of the superconducting diode [38]. It should be noticed that placing the S/F structures on the surface of a 3D topological insulator leads to a number of striking phenomena of magnetoelectric nature [39,40,41,42,43,44]. Moreover, new electronic states have been predicted to appear in such structures, including magnetic monopoles [45] and Majorana fermions [46,47,48,49]. It has also been predicted that the presence of helical magnetization in the F layer leads to nonmonotonic dependence of the critical temperature on the F layer width in S/F/TI structures [50].

The majority of the existing theoretical studies on the SDE focus either on the microscopic numerical calculations [30,32,51,52] or on the phenomenological approach [29,33]. In this work, we consider both linear and nonlinear approaches to calculate the SDE in the hybrid S/F/TI structure. We use the microscopic quasiclassical Green’s function formalism in the diffusive regime. We compare the results obtained by linear and nonlinear methods and discuss their ranges of applicability.

2. Materials and Methods

In this section, we present the model under consideration. The system is described by the following Hamiltonian model:

$$H=H_0+H_F+H_S,$$

(1)

where

$$\begin{array}{c}\hfill H_0=\int d^{2}r{\Psi }^{†}\left(\mathit{r}\right)\left[-i\alpha ({\nabla }_{\mathit{r}}\times \widehat{z})\mathbf{\sigma }-\mu +V\left(\mathit{r}\right)\right]\Psi \left(\mathit{r}\right),\end{array}$$

(2)

$$\begin{array}{c}\hfill H_F=-\int d^{2}r{\Psi }^{†}\left(\mathit{r}\right)\left[\mathit{h}\mathbf{\sigma }\right]\Psi \left(\mathit{r}\right),\end{array}$$

(3)

$$\begin{array}{c}\hfill H_S=\Delta \left(\mathit{r}\right){\Psi }_{\uparrow }^{†}\left(\mathit{r}\right){\Psi }_{\downarrow }^{†}\left(\mathit{r}\right)+{\Delta }^{*}\left(\mathit{r}\right){\Psi }_{\downarrow }\left(\mathit{r}\right){\Psi }_{\uparrow }\left(\mathit{r}\right).\end{array}$$

(4)

Here, ${\Psi }^{†}\left(\mathit{r}\right)=({\Psi }_{\uparrow }^{†}\left(\mathit{r}\right),{\Psi }_{\downarrow }^{†}\left(\mathit{r}\right))$ is the creation operator of an electron at the 3D TI surface, $\widehat{z}$ is the unit vector normal to the surface of TI, $\alpha$ is the Fermi velocity of electrons at the 3D TI surface and $\mu$ is the chemical potential. $\mathbf{\sigma }=({\sigma }_x,{\sigma }_y,{\sigma }_z)$ is a vector of Pauli matrices in spin space and $\mathit{h}=(h_x,h_y,0)$ is an in-plane exchange field, which is assumed to be nonzero only at $x<0$. The superconducting pairing potential $\Delta$ is nonzero only at $x>0$. Therefore, effectively, the TI surface states are divided into two parts: one at $x<0$ possesses $h\ne 0$, which can be called “ferromagnetic”, while the other corresponds to $x>0$ with $\Delta \ne 0$, which can be called “superconducting”. Below we use subscripts f and s to denote quantities related to the appropriate parts of the TI surface. The potential term $V\left(\mathit{r}\right)$ includes the nonmagnetic impurity scattering potential $V_{imp}={\displaystyle \sum _{{\mathit{r}}_i}}{V}_i\delta (\mathit{r}-{\mathit{r}}_i)$, which is of a Gaussian form $\langle V\left(\mathit{r}\right)V\left({\mathit{r}}^{\prime }\right)\rangle =(1/\pi \nu \tau )\delta (\mathit{r}-{\mathit{r}}^{\prime })$ with $\nu =\mu /\left(2\pi {\alpha }^2\right)$.

The superconductivity and in-plane exchange field is assumed to be proximity-induced, due to adjacent superconducting and ferromagnetic layers. Thus, we can imagine the system as a planar hybrid structure, consisting of a superconductor, S, and a ferromagnetic layer, F, on top of a three-dimensional topological insulator, TI, as shown schematically in Figure 1. The role of the TI surface is to provide strong spin–orbit coupling which produces a full spin–momentum locking effect. In this case, only one helical band crossing the Fermi energy is present. We employ the quasiclassical Green’s function formalism in the diffusive regime. In principle, Green’s function matrices have two degrees of freedom that are particle–hole and spin. In our model, the spin structure is characterized by a projection onto the conduction band:

$${\check{g}}_{s,f}({\mathit{n}}_F,\mathit{r},\epsilon )={\widehat{g}}_{s,f}(\mathit{r},\epsilon )\frac{(1+{\mathit{n}}_{\perp }\mathbf{\sigma })}{2},$$

(5)

where, ${\widehat{g}}_{s\left(f\right)}$ is the spinless Green’s function matrix in the particle–hole space in the superconducting (ferromagnetic) part of the 3D TI layer, ${\mathit{n}}_F={\mathit{p}}_F/p_F=(n_{F,x},n_{F,y},0)$ is a unit vector directed along the quasiparticle trajectory and ${\mathit{n}}_{\perp }=(n_{F,y},-n_{F,x},0)$ is a unit vector perpendicular to the quasiparticle trajectory, and directed along the quasiparticle spin, which is locked to the quasiparticle momentum.

In our theoretical analysis, we consider the diffusive limit, in which the superconducting coherence length is given by the expression ${\xi }_s=\sqrt{D_s/2\pi T_{cs}}$, where $D_s$ is the diffusion coefficient and $T_{cs}$ is the critical temperature of the bulk superconductor (we assume $\hslash =k_B=1$) and the elastic scattering length $\ell \ll {\xi }_s$. We also neglect the nonequilibrium effects in the structure [53,54,55].

In the following we outline the nonlinear and linear equations to calculate the SDE effect in the system under consideration.

2.1. Nonlinear Usadel Equations

The quasiclassical Usadel equation for spinless Green’s functions is [56,57]:

$$D\widehat{\nabla }\left(\widehat{g}\widehat{\nabla }\widehat{g}\right)=\left[{\omega }_n{\tau }_z+i\widehat{\Delta },\widehat{g}\right].$$

(6)

Here, D is the diffusion constant, ${\tau }_z$ is the Pauli matrix in the particle–hole space, $\widehat{\nabla }X=\nabla X+i\left(h_x{\widehat{e}}_y-h_y{\widehat{e}}_x\right)\left[{\tau }_z,\widehat{g}\right]/\alpha$. The gap matrix $\widehat{\Delta }$ is defined as $\widehat{\Delta }=\widehat{U}i{\tau }_x\Delta \left(x\right){\widehat{U}}^{†}$, where $\Delta \left(x\right)$ is a real function and transformation matrix $\widehat{U}=exp\left(iqy{\tau }_z/2\right)$. The finite center of mass momentum q takes into account the helical state. The Green’s function matrix is also transformed into $\widehat{g}=\widehat{U}{\widehat{g}}_q{\widehat{U}}^{†}$. Inclusion of the magnetization component $h_y$ produces no quantitative effect on either the supercurrent in y direction of the bilayer or on the critical temperature in the S part. It only enters the solution $f_f$ as a phase factor $\exp \left(2i{h}_{y}x/\alpha \right)$ [50,56]. Thus, we did not take it into consideration in our model, and define $h_x=h$.

In the hybrid structure under consideration, the helical state appears in the system in the following way. In our system, the Zeeman field and superconducting region are spatially separated, so that the helical state is realized via the proximity effect through the S/F interface. The helical state in the system is also characterized by the order parameter with a spatially inhomogeneous phase. However, the supercurrent density is not uniform in the structure. The total current across the hybrid structure is equal to zero. Thus, we call this state the hybrid helical state. More detailed analysis of the hybrid helical state can be found in [38].

To facilitate the solution procedures of the nonlinear Usadel equations we employ $\theta$ parametrization of the Green’s functions [58]:

$${\widehat{g}}_q=\left(\begin{array}{cc}cos\theta & sin\theta \\ sin\theta & -cos\theta \end{array}\right).$$

(7)

Substituting the above matrix into the Usadel Equation (6), we obtain the following in the S part of the TI surface $x>0$:

$${\xi }_s^2\pi T_{cs}\left[{\partial }_x^2{\theta }_s-\frac{q^2}{2}sin2{\theta }_s\right]={\omega }_{n}sin{\theta }_s-\Delta \left(x\right)cos{\theta }_s,$$

and in the F part $x<0$:

$${\xi }_f^2\pi T_{cs}\left[{\partial }_x^2{\theta }_f-\frac{q_m^2}{2}sin2{\theta }_f\right]={\omega }_{n}sin{\theta }_f,$$

(8)

where ${\xi }_f=\sqrt{D_f/2\pi T_{cs}}$, and $D_f$ is the diffusion coefficient of the ferromagnetic layer. $q_m=q+2h/\alpha$ and $X_{s\left(f\right)}$ means the value of X in the S(F) part of the TI surface, respectively. The self-consistency equation for the pair potential reads:

$$\Delta \left(x\right)ln\frac{T_{cs}}{T}=\pi T\sum _{{\omega }_n}\left(\frac{\Delta \left(x\right)}{|{\omega }_n|}-2sin{\theta }_s\right).$$

(9)

We supplement the above equations with the following boundary conditions at the S/F interface ($x=0)$ [59]:

$$\begin{array}{c}\hfill {\gamma }_B\frac{\partial {\theta }_f}{\partial x}{|}_{x=0}=sin\left({\theta }_s-{\theta }_f\right),\end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{{\gamma }_B}{\gamma }\frac{\partial {\theta }_s}{\partial x}{|}_{x=0}=sin\left({\theta }_s-{\theta }_f\right),\end{array}$$

(11)

where $\gamma ={\xi }_s{\sigma }_f/{\xi }_f{\sigma }_s$, ${\gamma }_B=R{\sigma }_f/x{i}_f$, and ${\sigma }_{s\left(f\right)}$ is the conductivity of the S (F)layer. The parameter $\gamma$ determines the strength of suppression of superconductivity in the S lead near the interface compared to the bulk. No suppression occurs for $\gamma =0$, while strong suppression takes place for $\gamma \gg 1$. The parameter ${\gamma }_B$ is a dimensionless parameter, describing the transparency of the S/F interface [59,60,61]. For example, the imperfect transparency can take place due to the Fermi velocity mismatch or presence of the potential barrier at the S/F interface. To complete the boundary problem, we also set boundary conditions at the free edges:

$$\frac{\partial {\theta }_f}{\partial x}{|}_{x=-d_f}=0,\frac{\partial {\theta }_s}{\partial x}{|}_{x=d_s}=0.$$

(12)

In order to calculate the superconducting current we utilize the expression for supercurrent density:

$${\mathbf{J}}_{s\left(f\right)}=\frac{-i\pi {\sigma }_{s\left(f\right)}}{4e}T\sum _{{\omega }_n}Tr\left[{\tau }_z{\widehat{g}}_{s\left(f\right)}\widehat{\nabla }{\widehat{g}}_{s\left(f\right)}\right].$$

(13)

Performing the unitary transformation U, the current density transforms as follows:

$$\begin{array}{c}\hfill j_y^s\left(x\right)=-\frac{\pi {\sigma }_{s}q}{2e}T\sum _{{\omega }_n}{sin}^2{\theta }_s,\end{array}$$

(14)

$$\begin{array}{c}\hfill j_y^f\left(x\right)=-\frac{\pi {\sigma }_n}{2e}\left[q+\frac{2h}{\alpha }\right]T\sum _{{\omega }_n}{sin}^2{\theta }_f.\end{array}$$

(15)

The total supercurrent flowing via the system along the y-direction can be calculated by integrating the current density of the total width of the S/F bilayer $d_f+d_s$:

$$I={\int }_{-d_f}^0{j}_y^f\left(x\right)dx+{\int }_0^{d_s}{j}_y^s\left(x\right)dx.$$

(16)

2.2. Linear Usadel Equations

At the limit when $T\approx T_c$, the Usadel equations (6) can be linearized, since the normal Green’s function is close to unity, i.e., ${\widehat{g}}_q\approx {\tau }_z+\theta \left(x\right){\tau }_x$.

In the superconducting S layer ($0<x<d_s$), the linearized Usadel equation for the spinless amplitude ${\theta }_s$ reads [56,57,58,62]:

$${\xi }_s^2\pi T_{cs}\left({\partial }_x^2-q^2\right){\theta }_s-{\omega }_n{\theta }_s+\Delta =0.$$

(17)

In the ferromagnetic region of the TI, the linearized Usadel equation takes the form:

$${\partial }_x^2{\theta }_f=\left[\frac{{\omega }_n}{{\xi }_f^2\pi T_{cs}}+q_m^2\right]{\theta }_f.$$

(18)

The solution of Equation (18) is found in the form:

$${\theta }_f=C\left({\omega }_n\right)cosh{k}_q\left(x+d_f\right),$$

(19)

where

$$\begin{array}{c}\hfill k_q=\sqrt{\frac{|{\omega }_n|}{{\xi }_f^2\pi T_{cs}}+q_m^2}.\end{array}$$

(20)

Here, $C\left({\omega }_n\right)$ is found from the boundary conditions. Using boundary conditions (10), the problem is written in a closed form with respect to the Green’s function $f_s$. At $x=0$ the boundary conditions can be written as:

$${\xi }_s\frac{\partial {\theta }_s\left(0\right)}{\partial x}=W^q\left({\omega }_n\right){\theta }_s\left(0\right),$$

(21)

where,

$$W^q\left({\omega }_n\right)=\frac{\gamma }{{\gamma }_B+A_{qT}\left({\omega }_n\right)},A_{qT}\left({\omega }_n\right)=\frac{1}{k_q{\xi }_f}coth{k}_q{d}_f.$$

(22)

In general, the boundary condition (21) can be complex. However, in the considered system, $A_{qT}$ is real. Hence, the condition (21) coincides with its real-valued form.

Then, we write the self-consistency equation for $\Delta$ considering only positive Matsubara frequencies:

$$\Delta ln\frac{T_{cs}}{T}=\pi T\sum _{{\omega }_n>0}\left(\frac{2\Delta }{{\omega }_n}-2{\theta }_s\right),$$

(23)

The Usadel equation in the superconducting part is given by:

$${\xi }_s^2\left(\frac{{\partial }^2{\theta }_s}{\partial x^2}-{\kappa }_{qs}^2{\theta }_s\right)+\frac{\Delta }{\pi T_{cs}}=0.$$

(24)

Within the linearized Usadel equations, the supercurrent is also calculated in the self-consistent manner, using Equations (21)–(24).

Single-Mode Approximation

In the framework of the single-mode approximation, the solution in S is introduced in the form [63,64]:

$${\theta }_s(x,{\omega }_n)=f\left({\omega }_n\right)cos\left(\Omega \frac{x-d_s}{{\xi }_s}\right),$$

(25)

$$\Delta \left(x\right)=\delta cos\left(\Omega \frac{x-d_s}{{\xi }_s}\right).$$

(26)

The solution presented above automatically satisfies boundary condition (12) at $x=d_s$. Substituting expressions (25) and (26) into the Usadel equation for ${\theta }_s$ (24) yields:

$$\begin{array}{c}\hfill f\left({\omega }_n\right)=\frac{\delta }{{\omega }_n+{\Omega }^2\pi T_{cs}+q^2{\xi }_s^2\pi T_{cs}}.\end{array}$$

(27)

Employing the single-mode approximation and using the solution in the TI layer, as well as the boundary conditions, we find the total supercurrent flowing through the system in the y direction. We obtain the following expression:

$$\begin{array}{c}\hfill I=-\frac{\pi }{4e}T\sum _{{\omega }_n}{f}^2\left({\omega }_n\right)\left[{\sigma }_{s}q{C}_s+\frac{{\sigma }_n}{{\beta }^2}{cos}^2\left(\Omega \frac{d_s}{{\xi }_s}\right)q_m{C}_f\right],\\ \hfill C_s=\left(d_s+\frac{{\xi }_s}{2\Omega }sin\left(2\Omega \frac{d_s}{\xi }\right)\right),\\ \hfill C_f=\left(d_f+\frac{1}{2k_q}sinh\left(2k_q{d}_f\right)\right),\end{array}$$

(28)

where the coefficient $\beta$ is defined as:

$$\begin{array}{c}\hfill \beta ={\gamma }_B{k}_q{\xi }_{f}sinh{k}_q{d}_f+cosh{k}_q{d}_f,\end{array}$$

(29)

and $\Omega$ is calculated from the boundary condition for the single-mode approximation (21). In the following section, we present the results of the supercurrent calculations for both linear (self-consistent and single-mode) and nonlinear approaches.

3. Results

In this section, we present the results of the calculations based on the model presented above. For simplicity we set ${\xi }_s={\xi }_f=\xi$.

In Figure 2, we compare $I\left(q\right)$ dependencies calculated by linear and nonlinear approaches. Both of these curves were calculated in a numerical self-consistent approach. At equilibrium, when no external supercurrent was applied, the superconducting system chose the state with non-zero $q_0$. This state corresponded to zero total supercurrent (in the y direction), i.e., the condition $I\left(q_0\right)=0$ was satisfied, which can be seen in Figure 2. In the system, q cannot be varied directly, but one can apply an external supercurrent through the structure in the range between the critical supercurrents $I_c^{+}$ and $I_c^{-}$, where $I_c^{+(-)}$ is the critical supercurrent for the positive (negative) direction. Applying a certain value of the supercurrent leads to the superconducting state with a corresponding value of $q\ne q_0$. When current $\left|I\right|$ was applied in the range between $I_c^{+}$ and $|I_c^{-}|$ the system was still in a superconducting state for $+I$ (since there was a superconducting state with a corresponding q), and the system was in the normal state for $-I$ (there was no superconducting state with q that corresponded to the applied current $-I$ ).

We noticed that the linearized solution resulted in higher values of the critical currents. As expected, the linearized approach did not capture nonlinearities in the current behavior as a function of q.

When calculating the critical temperature in the hybrid structure, it is common to use the single-mode approach within the linearized Usadel equations. We tested the application possibility of the single-mode method in calculating the supercurrent. The main disadvantages of Equations (25) and (26) are that these expressions disregard the dependencies of the amplitude $\delta$ on the parameter q. Moreover, the amplitude of the pair potential cannot be obtained within the solution provided by the single-mode, i.e., $\delta$ remains as a fitting parameter. In Figure 3, we compare the full nonlinear approach and the single-mode approximation. We observed that the supercurrent derived by the single-mode was in fairly good agreement with the nonlinear method in the vicinity of the equilibrium value of $q=q_0$. However, for larger values of $q-q_0$, it was clear that the single-mode approach tended to fail, resulting in much larger values of the critical current.

In the calculations we set the temperature $T=0.25T_{cs}$, which may seem to be out of the applicability range of the linearized Usadel equations. However, the true value of the critical temperature in the S part can be substantially lower, due to strong suppression from the adjacent F part of the structure. In our system, the critical temperature did not exceed $0.3T_{c}s$, as seen in Figure 2 and Figure 3.

It is more instructive to discuss the diode quality factor, which is defined in the following way:

$$\eta =\frac{\Delta I_c}{I_c^{+}+\left|I_c^{-}\right|}=\frac{I_c^{+}-\left|I_c^{-}\right|}{I_c^{+}+\left|I_c^{-}\right|}.$$

(30)

However, before doing so, we briefly discuss the temperature dependencies of the supercurrent nonreciprocity for the linear and nonlinear approaches. In Figure 4, the temperature dependence of $\Delta I=I_c^{+}-\left|I_c^{-}\right|$ is demonstrated. Both methods showed approximately the same temperature when the current nonreciprocity vanished. This corresponded to the transition from the superconducting state to the normal state. It can be seen from the figure that there was a substantial quantitative difference between the two methods.

In Figure 5, we demonstrate the SDE quality factor as a function of the exchange field $\xi h/\alpha$ for the nonlinear and linear approaches. We emphasize that the quality factors calculated for both cases were quite similar, despite the fact that $I\left(q\right)$ may be substantially different both quantitatively and qualitatively (Figure 2). This fact can be connected with the definition of the quality factor $\eta$. In brief, since $\eta$ is defined as a ratio between a sum and a difference of the critical currents, it loses information about current values and q dependencies.

4. Discussion

In this work, we calculated the superconducting diode effect using three different methods, which included linear and nonlinear equations.The results obtained above suggest several conclusions. The simplified single-mode approximation for the linearized Usadel equation is only applicable for a qualitative critical temperature calculation. Although the single-mode approach may be used at the vicinity of $(q-q_0)$, it does not capture the possible q dependency of the pair potential and fails at larger $|q-q_0|$. When operating close to the critical temperature, the full solution of the linearized Usadel equation provides adequate results. In particular, $\eta$, calculated via the linearized approach, is in good agreement with the nonlinear case (Figure 5). Nevertheless, in order to get a valid description of the helical state and the SDE in a wide range of parameters one should use fully nonlinear equations.