The quasiclassical theory for interfaces with spin--orbit coupling

TL;DR

This paper develops a new quasiclassical boundary condition framework that explicitly includes spin--orbit coupling, enabling the analysis of superconducting hybrid structures with heavy metal layers and ferromagnets, revealing long-range triplet superconductivity.

Contribution

It introduces a novel set of boundary conditions for quasiclassical theory that incorporate spin--orbit coupling, advancing the understanding of superconducting interfaces in spintronics.

Findings

Derived boundary conditions accounting for spin--orbit coupling.

Demonstrated emergence of long-range triplet superconductivity.

Applied theory to superconductor-ferromagnet systems.

Abstract

In recent years, substantial progress has been made regarding the application of quasiclassical theory on superconducting hybrid structures. This theoretical framework is reliant on a proper set of boundary conditions in order to describe multilayered systems. With the advent of the field of superconducting spintronics, systems which combine heavy metal layers, in which there is large spin--orbit coupling, with ferromagnets have received a great deal of attention, due to their potential for generating long range triplet superconductivity. In contrast to interfaces of strongly spin polarized materials, which are well understood, a quasiclassical theory for interfaces in systems where there is significant spin--orbit coupling does not yet exist. After reviewing the quasiclassical theory for interfaces, we here solve this problem by deriving a new set of boundary conditions which take spin--orbit coupling explicitly into account. We then go on to apply these boundary conditions to a superconductor-ferromagnet (SF) bilayer and an SFS Josephson weak link, demonstrating the emergence of long range triplet superconductivity in these systems.

Equations

$$\hat{n}\cdot {\check{g}}_1\nabla {\check{g}}_1=T\left[{\check{g}}_1,{\check{g}}_2\right]+T_{\alpha }\left[{\check{g}}_1,{\check{\mathbf{\sigma }}}_{||}{\check{g}}_2{\check{\mathbf{\sigma }}}_{||}\right]$$

$$+T_{\alpha }^{\prime }\left[{\check{g}}_1,{\check{\mathbf{\sigma }}}_{||}{\check{g}}_1{\check{\mathbf{\sigma }}}_{||}\right]+i\sqrt{T_{\alpha }^{\prime \prime }T}\left[{\check{g}}_1,\left\{{\check{g}}_2\left[{\check{\mathbf{\sigma }}}_{||},{\check{g}}_2\right],{\check{\mathbf{\sigma }}}_{||}\right\}\right]$$

$$+i\sqrt{T_{\alpha }^{\prime \prime }T}\left[{\check{g}}_1,\left\{{\check{g}}_2,{\check{\mathbf{\sigma }}}_{||}\right\}{\check{g}}_1{\check{\mathbf{\sigma }}}_{||}+{\check{\mathbf{\sigma }}}_{||}{\check{g}}_1\left\{{\check{g}}_2,{\check{\mathbf{\sigma }}}_{||}\right\}\right].$$

$$\check{G}=\left(\begin{array}{cc}{\hat{G}}^R& {\hat{G}}^K\\ 0& {\hat{G}}^A\end{array}\right),{\hat{G}}^R=\left(\begin{array}{cc}{G}^R& F^R\\ {\left(F^R\right)}^{\ast }& {\left(G^R\right)}^{\ast }\end{array}\right),$$

$$[\epsilon {\hat{\rho }}_3-\frac{1}{2m}{\left(-i{\nabla }_0\check{I}-\check{\mathbf{A}}({\mathbf{r}}_0)\right)}^2+\mu \check{I}-\check{\Sigma }({\mathbf{r}}_0)$$

$$-V_{\text{imp}}({\mathbf{r}}_0)\check{I}]\check{G}({\mathbf{r}}_0,{\mathbf{r}}_n;\epsilon )=\delta ({\mathbf{r}}_0-{\mathbf{r}}_n)\check{I},$$

$$\check{G}({\mathbf{k}}_0,{\mathbf{k}}_n)={\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_n)$$

$$+\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}\int \frac{d{\mathbf{k}}_2}{(2\pi {)}^3}{\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_1)V_{\text{imp}}({\mathbf{k}}_1-{\mathbf{k}}_2)\check{G}({\mathbf{k}}_2,{\mathbf{k}}_n),$$

$$\check{G}({\mathbf{k}}_0,{\mathbf{k}}_n)=\int d{\mathbf{r}}_0\int d{\mathbf{r}}_n\check{G}({\mathbf{r}}_0,{\mathbf{r}}_n)e^{-i{\mathbf{k}}_0\cdot {\mathbf{r}}_0+i{\mathbf{k}}_n\cdot {\mathbf{r}}_n}.$$

$$\check{G}({\mathbf{k}}_0,{\mathbf{k}}_n)=\int d\mathbf{R}\check{G}(\mathbf{k},\mathbf{R})e^{-i\Delta \mathbf{k}\cdot \mathbf{R}},$$

$$\check{G}(\mathbf{k},\mathbf{R})=\int d\mathbf{r}\check{G}(\mathbf{r},\mathbf{R})e^{-i\mathbf{k}\cdot \mathbf{R}}.$$

$$\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}$$

$$-\frac{1}{2m}{\left[\check{\mathbf{A}}({\mathbf{k}}_0-{\mathbf{k}}_1)\right]}^2-\check{\Sigma }({\mathbf{k}}_0-{\mathbf{k}}_1)]{\check{G}}_0({\mathbf{k}}_1,{\mathbf{k}}_n)$$

$$=\check{I}\delta ({\mathbf{k}}_0-{\mathbf{k}}_n).$$

$$\check{G}({\mathbf{k}}_0,{\mathbf{k}}_n)={\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_n)$$

$$+\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}\int \frac{d{\mathbf{k}}_2}{(2\pi {)}^3}{\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_1){\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_1,{\mathbf{k}}_2)\check{G}({\mathbf{k}}_2,{\mathbf{k}}_n),$$

$${\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_0,{\mathbf{k}}_n)=V_{\text{imp}}({\mathbf{k}}_0-{\mathbf{k}}_n)+$$

$$\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}\int \frac{d{\mathbf{k}}_2}{(2\pi {)}^3}{V}_{\text{imp}}({\mathbf{k}}_0-{\mathbf{k}}_1){\check{G}}_0({\mathbf{k}}_1,{\mathbf{k}}_2){\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_2,{\mathbf{k}}_n).$$

$${\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_0,{\mathbf{k}}_n)=\int \frac{d\mathbf{q}}{(2\pi {)}^3}|U({\mathbf{k}}_0-\mathbf{q}){|}^2\check{G}(\mathbf{q},\mathbf{q}-\Delta \mathbf{k}).$$

$${\check{\Sigma }}_{\text{imp}}(\mathbf{k},\mathbf{R})=\int \frac{d\mathbf{q}}{(2\pi {)}^3}|U(\mathbf{k}-\mathbf{q}){|}^2\check{G}(\mathbf{q},\mathbf{R}).$$

$$[\epsilon {\hat{\rho }}_3-\frac{1}{2m}{\left(-i\left(i\mathbf{k}+\frac{1}{2}{\nabla }_R\right)\check{I}-\check{\mathbf{A}}(\mathbf{R})\right)}^2+\mu \check{I}-\check{\Sigma }(\mathbf{R})$$

$$-{\check{\Sigma }}_{\text{imp}}(\mathbf{k},\mathbf{R})]\otimes \check{G}(\mathbf{k},\mathbf{R})=\check{I},$$

$$\check{g}({\mathbf{k}}_F,\mathbf{R})=\frac{i}{\pi }\int d{\xi }_k\check{G}(\mathbf{k},\mathbf{R}),$$

$$i{\mathbf{v}}_F\cdot \tilde{\nabla }\check{g}({\mathbf{k}}_F,\mathbf{r})$$

$$+\left[\epsilon {\check{\rho }}_3+\check{\Sigma }({\mathbf{k}}_F,\mathbf{r})-{\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_F,\mathbf{r}),\check{g}({\mathbf{k}}_F,\mathbf{r})\right]=0.$$

$$D\nabla \cdot {\check{g}}_d\nabla {\check{g}}_d+i\left[\epsilon {\rho }_3-\check{\Sigma },{\check{g}}_d\right]=0,$$

$$\check{G}({\mathbf{k}}_0,{\mathbf{k}}_n)={\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_n)$$

$$+\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}\int \frac{d{\mathbf{k}}_2}{(2\pi {)}^3}{\check{G}}_0({\mathbf{k}}_0,{\mathbf{k}}_1)\check{T}({\mathbf{k}}_1,{\mathbf{k}}_2){\check{G}}_0({\mathbf{k}}_2,{\mathbf{k}}_n),$$

$$\check{T}({\mathbf{k}}_0,{\mathbf{k}}_n)=\check{V}({\mathbf{k}}_0-{\mathbf{k}}_n)$$

$$+\int \frac{d{\mathbf{k}}_1}{(2\pi {)}^3}\int \frac{d{\mathbf{k}}_2}{(2\pi {)}^3}\check{V}({\mathbf{k}}_0-{\mathbf{k}}_1){\check{G}}_0({\mathbf{k}}_1,{\mathbf{k}}_2)\check{T}({\mathbf{k}}_2,{\mathbf{k}}_n).$$

$$\check{T}({\mathbf{k}}_0,{\mathbf{k}}_n)={\check{T}}_0({\mathbf{k}}_0,{\mathbf{k}}_n)e^{-i({\mathbf{k}}_0-{\mathbf{k}}_n)\cdot {\mathbf{R}}_0},$$

$$\check{T}(\mathbf{k},\mathbf{R})={\check{T}}_0(\mathbf{k},\mathbf{k})\delta (\mathbf{R}-{\mathbf{R}}_0).$$

$$i{\mathbf{v}}_F\cdot \tilde{\nabla }\check{g}({\mathbf{k}}_F,\mathbf{R})+\left[\epsilon {\check{\rho }}_3+\check{\Sigma }({\mathbf{k}}_F,\mathbf{R})-{\check{\Sigma }}_{\text{imp}}({\mathbf{k}}_F,\mathbf{R}),\check{g}({\mathbf{k}}_F,\mathbf{R})\right]$$

$$=\left[{\check{t}}_0({\mathbf{k}}_F,\mathbf{R}),\check{g}({\mathbf{k}}_F,\mathbf{R})\right]\delta (\mathbf{R}-{\mathbf{R}}_0),$$

$${\check{t}}_0({\mathbf{k}}_F,\mathbf{R})={\check{V}}_0(0)+N_0\int \frac{d{\Omega }_q}{4\pi }{\check{V}}_0({\mathbf{k}}_F-{\mathbf{q}}_F)\check{g}({\mathbf{q}}_F,\mathbf{r}){\check{t}}_0({\mathbf{q}}_F,\mathbf{R}).$$

$$\check{T}({\mathbf{k}}_0,{\mathbf{k}}_n)={\check{T}}_s({\mathbf{k}}_{0,||},{\mathbf{k}}_{n,||})e^{-i(k_{0,\perp }-k_{n,\perp })R_0},$$

$${\check{T}}_s({\mathbf{k}}_{0,||},{\mathbf{k}}_{n,||})$$

Full text

I Introduction

The quasiclassical approximation serene_pr_1983 ; belzig_sm_1999 ; rammer_rmp_1986 is a versatile tool with which complex quantum mechanical problems can be simplified to such an extent that they become numerically solvable. The main assumption of this approximation is that the relevant quantities under study vary on length scales which are much larger than the Fermi wavelength so that the computationally challenging shorter length scale oscillations may be integrated out. This makes the quasiclassical approximation particularly well suited for superconductors, in which the superconducting condensate may remain correlated over mesoscopic distances. Indeed, the most general theoretical framework used to describe superconductors, the Green function technique, requires the solution of the Gor’kov equation AGD , which is too cumbersome in all but a few select problems. Instead, progress can be achieved with its quasiclassical equivalents, the Eilenberger eilenberger_zp_1968 and the Usadel equation usadel_prl_1970 , which govern quassiclassical Green functions describing only the envelopes of the original propagators, and remain the only viable solution method for many problems of practical interest.

When non-superconducting materials are attached to a superconductor to form superconducting hybrid structures, an interesting phenomenon occurs. In such systems, superconducting correlations may leak into the adjoining non-superconducting materials, so that they too attain superconducting properties. This is known as the proximity effect. The study of such systems necessarily involves the proper treatment of interfaces between materials. However, while the bulk properties of superconductors are easily described within quasiclassical theory, interfaces between materials is another matter entirely. In the vicinity of an interface, the governing Hamiltonian changes abruptly, which invalidates the use of the quasiclassical approximation. The consequence of this is that the quasiclassical Green functions feature a discontinuous jump at interfaces, the size of which is impossible to determine within quasiclassical theory—additional information is needed. This jump was first computed for ballistic SN structures using a full microscopic description of the interface zaitsev_jetp_1984 , thereby giving a set of boundary conditions linking the two materials. These boundary conditions were generalized with the use of a projection operator method yip_jltp_1997 ; eschrig_prb_2000 ; shelankov_prb_2000 . Alternative derivations, where the interface is treated perturbatively via a $T$ matrix approach ashauer_jltp_1986 ; nagai_jltp_1988 ; graser_prb_2007 have also been proposed. In the diffusive limit, boundary conditions may be arrived at by connecting the momentum independent diffusive Green functions far away from the interface to ballistic Green functions present in a region immediately surrounding the interface kupriyanov_jetp_1988 ; nazarov_sm_1999 .

For hybrid structures involving strongly polarized magnetic materials, the interfacial boundary conditions generally become spin active. Such boundary conditions have been formulated heuristically using a tunnelling Hamiltonian bergeret_prb_2012 . Another, more fundamental, approach is to connect the two sides of the interface by means of scattering or transfer matrices, in which the spin dependence of the scattering processes are taken into account millis_prb_1988 ; tokuyasu_prb_1988 ; cottet_prb_2009 . Both the projection operator and the $T$ matrix method have also been successfully generalized to handle spin active interfaces in ballistic systems cuevas_prb_2001 ; eschrig_prl_2003 ; kopu_prb_2004 ; fogelstrom_prb_2000 ; zhao_prb_2004 . The latter method was also applied to diffusive systems in Ref. eschrig_njp_2015 , in which a completely general theory for boundary conditions in spin polarized hybrid structures was derived.

When hybrid structures are made from superconducting and ferromagnetic materials, the proximity effect allows for the coexistence of both magnetic and superconducting correlations. This produces a number of interesting effects. One of the more fascinating effects is perhaps the appearance of triplet superconductivity, due to the Zeeman splitting endowing the Cooper pairs with a net momentum buzdin_rmp_2005 . In homogeneous ferromagnets, the triplet Cooper pairs remain in the spinless state, and hence the magnetization has a strong depairing effect. In ferromagnets where the magnetization is inhomogeneous, on the other hand, the triplet Cooper pairs may be converted to an equal spin state bergeret_prl_2001 ; bergeret_rmp_2005 , giving them a net spin. In such a state, the Cooper pairs are insensitive to a parallel magnetization, and may therefore persist for long distances into proximitized ferromagnets, a phenomenon known as long range triplet superconductivity. In addition, since these correlations are spin polarized, they may carry spin currents, a realisation which has strongly contributed to the field of superconducting spintronics linder_np_2015 ; eschrig_rpp_2015 .

Another avenue towards triplet superconductivity is by introducing layers with spin–orbit coupling gorkov_prl_01 ; annunziata_prb_12 ; bergeret_prb_2014 . In such systems, spin polarized supercurrents may be generated even in homogeneous ferromagnets, which is advantageous from an experimental point of view. However, while the boundary conditions for spin polarized systems is well understood, a theory for interfaces in which spin–orbit coupling is prominent does not yet exist. In this paper, we seek to remedy this by deriving a new set of boundary conditions which take spin–orbit coupling into account, using the general framework of Ref. eschrig_njp_2015 . This result consequently allows for a proper treatment of spin-orbit coupled interfaces in quasiclassical theory, which is in principle relevant for any heterostructure, since interfaces break inversion symmetry, but particularly so for heterostructures with heavy metal interlayers.

We here state our main analytical result and comment briefly on the qualitative physical meaning of each term in the boundary condition. Below, we shall derive this result rigorously and provide a clear description of how the various terms arise, and the meaning of each of the symbols.

[TABLE]

The usual Kupriyanov-Lukichev term kupriyanov_jetp_1988 is the one proportional to $T$. The term proportional to $T_{\alpha}$ represents a correction to the usual tunneling boundary condition from the spin-orbit coupling part of the tunneling matrix. The term $T_{\alpha}^{\prime}$ represents spin-dependent phase-shifts occurring due to spin-orbit coupling at the interface and thus exists even in the absence of any tunnelling. The terms $\sqrt{T_{\alpha}^{\prime\prime}T}$ are higher order corrections which exist only in the presence of spin-independent tunneling, spin-orbit coupled tunneling, and spin-orbit coupled reflection.

The quasiclassical theory for boundary conditions is quite intricate, involving many details, so to ensure complete clarity of the ensuing derivation, we include in Section II a review of the treatment of interfaces within quasiclassical theory. In Section III we formulate the boundary conditions in the presence of spin–orbit coupling, and in Section IV we apply the new boundary conditions to example problems, in order to demonstrate their use.

II Review of the quasiclassical theory for interfaces

In this section, a review of quasiclassical boundary conditions will be given, starting with a brief excursion into general quasiclassical theory.

II.1 Quasiclassical equations of motion

All physical observables of interest may be expressed in terms of Green functions, and we use the Keldysh formalism, in which the Green functions take the form of $8\times 8$ matrices in spin$\,\otimes\,$Nambu$\,\otimes\,$Keldysh space, defined as

[TABLE]

where $\hat{G}^{A}=-\hat{\rho}_{3}\left(\hat{G}^{R}\right)^{\dagger}\hat{\rho}_{3}$, with $\hat{\rho}_{3}=\text{diag}\left(+1,+1,-1,-1\right)$, and $\hat{G}^{K}=\hat{G}^{R}\hat{h}-\hat{h}\hat{G}^{A}$ for a given distribution function $\hat{h}$ belzig_sm_1999 . Furthermore, $G^{X}$ and $F^{X}$, with $X\in\left\{R,A,K\right\}$ are $2\times 2$ matrices in spin space. We assume time translation invariance, and Fourier transform in the relative time coordinate, so that we may write $\check{G}=\check{G}(\bm{r}_{0},\bm{r}_{n};\varepsilon)\equiv\check{G}(\bm{r}_{0},\bm{r}_{n})$, where $\varepsilon$ is the quasiparticle energy. The equation of motion may then be written as

[TABLE]

where $\check{I}$ is the $8\times 8$ identity matrix, $\check{\bm{A}}(\bm{r}_{0})$ is the vector potential, $V_{\text{imp}}(\bm{r}_{0})$ is the impurity potential, and $\check{\Sigma}(\bm{r}_{0})$ encompasses any other local self energies, such as the superconducting gap, or an exchange field. The impurity potential in the conventional way—as a perturbation in momentum space. Towards that end, Eq. 3 is reformulated as

[TABLE]

where the Fourier transform of the Green function is defined as

[TABLE]

Note that by changing variables to $\bm{R}=(\bm{r}_{0}+\bm{r}_{n})/2$ and $\bm{r}=\bm{r}_{0}-\bm{r}_{n}$, the Fourier transform may also be written as

[TABLE]

with $\bm{k}=(\bm{k}_{0}+\bm{k}_{n})/2$ and $\Delta\bm{k}=\bm{k}_{0}-\bm{k}_{n}$, and

[TABLE]

Eq. 7 is known as the mixed representation of $\check{G}$, as it involves both the center of mass position $\bm{R}$ and the center of mass momentum $\bm{k}$. The unperturbed Green function $\check{G}_{0}(\bm{k}_{0},\bm{k}_{n})$ satisfies

[TABLE]

By making use of standard techniques for diagram summation AGD , taking into consideration that $\check{G}_{0}(\bm{k}_{0},\bm{k}_{n})$ depends on two momentum indices, Eq. 4 can be written as

[TABLE]

with the impurity self energy $\check{\Sigma}(\bm{k}_{0},\bm{k}_{n})$ defined as

[TABLE]

It is customary to employ the self-consistent Born approximation, in which Eq. 9 is truncated at second order in $V_{\text{imp}}$, and the replacement $\check{G}_{0}\to\check{G}$ is made on the right hand side. The latter is equivalent to including a much greater number of diagrams in the self energy. The impurity potential is assumed to consist of a large number of randomly distributed identical impurities, so that $V_{\text{imp}}(\bm{r})=\sum_{j}U(\bm{r}-\bm{r}_{j})$. After impurity averaging, ignoring terms which only depend on $U(0)$ since these simply renormalize the chemical potential, the self energy becomes

[TABLE]

In the mixed representation, arrived at by Fourier transforming in the relative momentum $\Delta\bm{k}$, this expression becomes particularly simple,

[TABLE]

By comparing with the full Fourier transform, given in Eq. 5, we may identify $\bm{R}$ as the center of mass position. Eq. 3 takes the following form in the mixed representation,

[TABLE]

where the operator $\otimes$ indicates the Moyal product, which is defined as $A(\bm{k},\bm{R})\otimes B(\bm{k},\bm{R})=e^{i(\nabla^{(A)}_{R}\cdot\nabla^{(B)}_{k}-\nabla^{(A)}_{k}\cdot\nabla^{(B)}_{R})/2}A(\bm{k},\bm{R})B(\bm{k},\bm{R})$. If the spatial variation of $A(\bm{k},\bm{R})$ and $B(\bm{k},\bm{R})$ is slow, one may approximate $A(\bm{k},\bm{R})\otimes B(\bm{k},\bm{R})\simeq A(\bm{k},\bm{R})B(\bm{k},\bm{R})$.

Eq. 12 may be further simplified by introducing the quasiclassical approximation, wherein the rapid oscillations of the Green function are integrated out rammer_rmp_1986 ,

[TABLE]

where $\xi_{k}=\frac{1}{2m}(\bm{k}^{2}-k_{F}^{2})$. In Eq. 13 there is an implicit assumption that the Green function $\check{G}(\bm{k},\bm{R})$ is strongly peaked at the Fermi level $k_{F}$, so that only the angular dependence of the momentum $\bm{k}$ appears in the quasiclassical Green function $\check{g}(\bm{k}_{F},\bm{R})$. This is satisfied as long as the spatial variation of the self-energies appearing in $\check{G}$ is sufficiently slow. The quasiclassical approximation may not be applied to Eq. 12 directly, as it contains both constant terms and terms proportional to $\xi_{k}$. These terms can be removed by employing the so-called “left–right” trick eilenberger_zp_1968 ; larkin_jetp_1969 ; serene_pr_1983 , where one instead considers the difference between Eq. 12 and its adjoint, thereby cancelling out the problematic terms. Doing so leads to the Eilenberger equation eilenberger_zp_1968 ,

[TABLE]

Eq. 14 is accompanied by a normalization condition on the quasiclassical Green function, $\check{g}^{2}=\check{I}$.

In the limit of large concentrations of impurities, the effect of frequent scatterings may be included by averaging over momentum direction. This defines a diffusive Green function $\check{g}_{d}=\langle\check{g}\rangle$, and its governing equation of motion, the Usadel equation belzig_sm_1999 ; rammer_rmp_1986 ; usadel_prl_1970

[TABLE]

where $D$ is the diffusion constant.

II.2 Distinguished impurities

We next consider a case where there is an additional impurity $\check{V}(\bm{r}_{0})$ present, which may in some way be distinguished from the averaged impurities described by $\check{\Sigma}_{\text{imp}}$. We further assume that this impurity is strongly localized at some position. This means that impurity averaging is not possible. Even so, the quasiclassical formulation of the equation of motion for such a system may be arrived at by perturbation theory thuneberg_prb_1984 . Indeed, if any interference between the averaged and the localized impurity is neglected, the integral equation for the Green function once again takes the form of Eq. 4, where $\check{V}$ replaces $V_{\text{imp}}$ as the perturbing potential, and $\check{G}_{0}$ is the Green function for a system where $\check{\Sigma}_{\text{imp}}$ is included, but where $\check{V}=0$. By repeated iteration of this equation, it is seen that it may be written in the form

[TABLE]

with the $T$ matrix defined as

[TABLE]

If the distinguished impurity is localized at a position $\bm{R}_{0}$, the Fourier transformed of the impurity potential is given as $\check{V}(\bm{q})=\check{V}_{0}(\bm{q})e^{-i\bm{q}\cdot\bm{R}_{0}}$, where $\check{V}_{0}(\bm{q})$ is a slowly varying function of $\bm{q}$. By inserting this into Eq. 17, it is seen that the $T$ matrix can be written as

[TABLE]

where $\check{T}_{0}(\bm{k}_{0},\bm{k}_{n})$ has the exact same form as Eq. 17, with the replacements $\check{V}\to\check{V}_{0}$ and $\check{T}\to\check{T}_{0}$. It is thus a slowly varying function of $\bm{k}_{0}$ and $\bm{k}_{n}$, so that we may approximate $\check{T}_{0}(\bm{k}_{0},\bm{k}_{n})\simeq\check{T}_{0}(\bm{k},\bm{k})$, where $\bm{k}$ is the center of mass momentum. In the mixed representation, the $T$ matrix then becomes

[TABLE]

The equation of motion in the mixed representation becomes identical to Eq. 12, but with the addition of a term $\check{T}(\bm{k},\bm{R})\otimes\check{G}_{0}(\bm{k},\bm{R})$, which is once again approximated by a product. Following the same steps used in deriving Eq. 14 then gives

[TABLE]

where $\check{t}_{0}(\bm{k}_{F},\bm{R})$ is the quasiclassical version of the $T$ matrix, given as

[TABLE]

II.3 Interface

An interface, which is a plane in three dimensions, may be treated as an impurity that is localized along a specific direction, having translation invariance along the two orthogonal directions. This implies that a similar perturbation expansion as was discussed in the previous section may be applied also in this case. However, the interface may not be constructed from an ensemble of point impurities satisfying Eq. 20 thuneberg_prb_1984 . This is because i) the pointlike nature of the impurity was explicitly made use of in Eq. 18, and ii) interference between different impurities is neglected. Instead, we follow Ref. buchholtz_zp_1979 and consider a model surface of the form $\check{V}(\bm{r}_{0})=\check{V}_{0}\delta[\hat{n}\cdot(\bm{r}_{0}-\bm{R}_{n})]$, where $\hat{n}$ is the normal vector of the surface, and $\bm{R}_{n}$ is a point on the surface. To simplify the notation, we define $\hat{n}\cdot\bm{r}_{0}=r_{\perp}$, and $\hat{n}\cdot\bm{R}_{n}=R_{0}$. Furthermore, we have that $\bm{r}_{0}=r_{\perp}\hat{n}+\bm{r}_{||}$. Insertion into Eq. 17 allows us to define

[TABLE]

with $k_{j,\perp}$ and $\bm{k}_{j,||}$, respectively, the orthogonal and parallel components of momentum $j$, with respect to the surface. Moreover,

[TABLE]

where

[TABLE]

In the mixed represantion, Eq. 22 simply becomes

[TABLE]

with

[TABLE]

and

[TABLE]

To find the quasiclassical version of Eq. 27, we insert the inverse of Eq. 13, namely $\check{G}(\bm{k},\bm{R})=-i\pi\check{g}(\bm{k}_{F},\bm{R})\,\delta\left(\frac{1}{2m}(\bm{k}^{2}-k_{F}^{2})\right)$. Performing the integral over $q_{\perp}$ gives

[TABLE]

with $v_{n}=\frac{\bm{k}_{F}\cdot\hat{n}}{m}$,

[TABLE]

and $\bm{k}_{\pm}=\pm\sqrt{k_{F}^{2}-\bm{k}_{||}^{2}}\hat{n}+\bm{k}_{||}$. This means that the quasiclassical $T$ matrix for an interface,

[TABLE]

only depends on the average of Green functions whose normal component of the momentum direction point, respectively toward and away from the interface. The equation of motion, takes the same form as Eq. 20, with the replacements $\check{t}_{0}\to\check{t}_{s}$ and $\delta(\bm{R}-\bm{R}_{0})\to\delta(R_{\perp}-R_{0})$.

II.4 Formulation of boundary conditions

We next want to consider an interface between two different materials. This is done by expanding Hilbert space into two domains, which represent the two sides of the interface. The Green function in this space can be written as

[TABLE]

where the subscripts “1” and “2” indicate the two materials. The matrices $\check{g}_{12}$ and $\check{g}_{21}$ contain creation and annihilation operators on both sides of the interface. These quantities are drone amplitudes that do not have physical meaning, and since they will be eliminated from the theory, we do not specify them further. The interface itself is described as an infinite surface located at some position, and which mediates tunnelling between its two sides. Such a potential may be described as

[TABLE]

The surface is treated as a perturbation, and the $T$ matrix is given by Eq. 30.

Without the presence of the interface potential, there is no coupling between the two sides. The unperturbed Green function, $\breve{g}_{0}$, therefore takes the form

[TABLE]

Note that Eq. 33 satisfies a generalized version of Eq. 14, given as

[TABLE]

where

[TABLE]

Similarly, $\breve{g}$ satisfies

[TABLE]

To find a relationship between $\breve{g}$ and $\breve{g}_{0}$, Eq. 36 may be integrated along a small interval surrounding $R_{0}$. For a trajectory (as determined by $\bm{k}_{F}$) crossing the interface, this leads to

[TABLE]

where we henceforth define $\hat{n}$ to be an outwards pointing surface normal. While Eq. 31 is defined everywhere in space, its diagonal elements, $\check{g}_{11}$ and $\check{g}_{22}$, only make sense physically in, respectively, in material 1 and 2, i.e., on opposite sides of the interface. Without loss of generality, we choose material 1 to be the active material, that is, the material for which we formulate the boundary conditions. This means that $\hat{n}$ points from material 1 to material 2. If $\breve{g}(R_{0}^{-})$ is located in material 1, we thus need to eliminate $\breve{g}(R_{0}^{+})$, which is located in material 2. This can be done by making use of a generalized normalization condition, given by

[TABLE]

These conditions are clearly satisfied in the special case of $\breve{g}=\breve{g}_{0}$. That they are valid also in the more general case, can be seen by considering the original derivation of Shelankov shelankov_jltp_1985 , introducing the interface and taking care to only evaluate points away from the interface, so that $\breve{g}$ and $\breve{g}_{0}$ both satisfy Eq. 34—albeit with different boundary conditions buchholtz_jltp_2002 ; rainer_book_1988 . In material 1 the Green function which describes a particle on a trajectory towards the interface satisfies $\text{sgn}(\bm{k}_{F}\cdot\hat{n})=+1$. We label these Green functions as $\check{g}_{11}^{i}$ and $\check{g}_{0,1}^{i}$, indicating that they are incoming with respect to the interface. Similarly, Green functions where $\text{sgn}(\bm{k}_{F}\cdot\hat{n})=-1$ are labelled as outgoing; $\check{g}_{11}^{o}$ and $\check{g}_{0,1}^{o}$. By evaluating Eqs. 38 and 39 immediately adjacent to, and on opposite sides of the interface, and inserting Eq. 37, the following boundary conditions may be derived at the interface,

[TABLE]

Due to the form of the interface potential $\breve{V}_{0}$, $\breve{t}_{s}$ is in general dense in 1–2-space. However, a closed solution for the $\check{t}_{s,11}$ element may be found by iterating Eq. 30 once cuevas_prb_2001 ,

[TABLE]

Note that $\check{\bar{g}}_{0,j}=\frac{1}{2}\left(\check{g}_{0,j}^{i}+\check{g}_{0,j}^{o}\right)$, for a given side of the interface $j$. Eq. 42 may easily be solved for $\check{t}_{s,11}$, giving

[TABLE]

To summarize the progress so far, we have found boundary conditions for the quantities $\check{g}_{11}^{i}$ and $\check{g}_{11}^{o}$, given in Eqs. 40 and 41, expressed entirely in terms of $\check{g}_{0,j}^{i/o}$. While these unperturbed Green functions exist everywhere in space, they are only physically valid solutions on their respective sides of the interface. Furthermore, as discussed in Ref. eschrig_njp_2015 , they can be easily modified to describe a system with an impenetrable interface by having them satisfy the condition

[TABLE]

where $\check{S}$ is a scattering matrix. For spin independent scattering, $\check{S}=\check{I}$. The point is that rather than considering the interface as a perturbation to omnipresent Green functions, we may exploit the fact that we only evaluate $\check{g}_{0,j}$ in one region of space to redefine them to represent a system with an impenetrable interface, simply by imposing Eq. 44. This provides significant benefits. Ultimately, the goal is to derive boundary conditions for diffusive systems, which are governed by diffusive Green functions $\check{g}_{d,j}$. Even for such a system, there will always be a ballistic zone immediately surrounding an interface, since a particle travelling away from the interface will traverse a distance on the order of the mean free path before encountering its first impurity. Diffusive systems with an interface are therefore governed by ballistic Green functions $\check{g}_{0,j}^{i/o}$, with $\check{\Sigma}_{\text{imp}}=0$, close to the interface, and by $\check{g}_{d,j}$ far away from the interface. In between is an asymptotic matching region known as the isotropization zone, within which the momentum dependence of the ballistic Green functions are averaged out by repeated impurity scatterings. Matching the Green functions of the two regimes is only possible if the following conditions are satisfied millis_prb_1988 ; nazarov_sm_1999 ; cottet_prb_2009 ,

[TABLE]

Using Eqs. 45, 46, 47 and 48 and Eq. 44, the ballistic Green functions at the interface may be expressed in terms of their diffusive counterparts as

[TABLE]

The quasiclassical boundary conditions are completed by computing the matrix current directed at the interface. In the ballistic zone it is defined as

[TABLE]

where the integration measure is an angular average over a hemisphere. In the diffusive zone the matrix current is given as nazarov_sm_1999

[TABLE]

where $\sigma_{j}$ is the normal state conductivity of material $j$. The matrix current is conserved across the isotropization zone cottet_prb_2009 , and hence Eq. 52 may be equated with Eq. 51, which in turn is determined from Eqs. 40 and 41 and Eqs. 49 and 50, thus giving the complete boundary conditions solely in terms of the diffusive Green functions.

For spin independent scattering we have $\check{S}=\check{I}$, which gives $\check{g}_{0,j}^{i}=\check{g}_{0,j}^{o}=\check{g}_{d,j}$. This means that the difference between Eqs. 40 and 41 reduces to $\check{g}^{i}_{11}-\check{g}^{o}_{11}=\frac{1}{i|v_{n}|}\left[\check{t}_{11},\check{g}_{d,1}\right]$. By using Eqs. 43, 51 and 52 one may, with some algebra kopu_prb_2004 , produce Nazarovs boundary conditions nazarov_sm_1999 ,

[TABLE]

where the coupling constant $\tau(\theta)$ is given as

[TABLE]

with $\upsilon=V_{0}/v_{F}$, and $\theta$ is the angle of incidence with respect to the interface. $\tau(\theta)$ is shown in Fig. 1a), where it is seen that it attains its maximum value for normal incidence ($\theta=0$), and goes to zero for trajectories parallel to the interface ($\theta=\frac{\pi}{2}$), which is intuitively reasonable.

When $\tau(\theta)$ is small, e.g., in the tunnelling limit, we may neglect its contribution to the denominator of Eq. 53. This gives the Kupriyanov–Lukichev boundary conditions kupriyanov_jetp_1988 ,

[TABLE]

where $T$ is the angular average of the coupling constant, and is given as

[TABLE]

$T$ is shown for increasing $V_{0}$ in Fig. 1b). For $V_{0}=0$, there is no coupling between the two sides. The system reduces to its unperturbed state in which the interface is impenetrable, and hence $T=0$. For large values of $V_{0}$, $T$ also goes to zero. The reason for this is that when the barrier potential increases, incoming particles are more likely to be reflected than transmitted, which means that the two sides eventually become decoupled also in this limit.

III Interfaces with spin–orbit coupling

We will now consider an interface to a material within which spin–orbit coupling plays a prominent role, for instance a heavy metal. This means that the transmission probability will depend on both the spin of the incoming particle and its angle of incidence. We model this with a Rashba-like tunnelling coupling,

[TABLE]

where $\hat{n}_{\alpha}$ is a unit vector indicating the direction in which the symmetry is broken, and $w$ is a spin independent contribution to the tunnelling between the two sides separated by the interface, while $w_{\alpha}$ is the strength of the spin–orbit coupling contribution to the tunnelling. As the barrier region increases in width, both $w$ and $w_{\alpha}$ decrease toward zero.

III.1 Scattering matrix

As a step towards formulating the boundary conditions, we need to find the scattering matrix for a system where the interface is impenetrable. However, even though the probability for transmission through the interface is zero, an incoming particle may still penetrate the barrier for some distance before being reflected. While inside the scattering region, the particle will experience the spin–orbit coupling and accumulate a spin dependent phase, the size of which depends on the incoming angle of incidence. To find the scattering matrix for this process, we assume that sufficiently close to the interface, the interface potential is large enough to dominate all other self energies millis_prb_1988 . This means that we may follow the procedure of Ref. tokuyasu_prb_1988 and describe the interface as a step function potential in a free electron gas. For simplicity, we use a cylindrical coordinate system in which the $z$ axis is aligned with the surface normal. The Hamiltonian for such a system is given as

[TABLE]

where $\alpha$ expresses the strength of the spin–orbit interaction at the interface, $\mu=\varepsilon_{F}$ is the chemical potential, and $\theta(z)$ is the Heaviside step function. We consider particles with energies close to the Fermi energy, meaning excitation energies $\varepsilon\simeq 0$, which is the relevant energy regime for quasiclassical theory. The wave function therefore satisfies $H\psi\simeq 0$. The parameter $\varepsilon_{g}$ expresses an energy gap, and is included to ensure that the wave function is evanescent in the barrier region. Since the interface is assumed to be perfectly smooth, the momentum parallel to the interface, $\bm{k}_{||}$, is conserved during the scattering process. Hence, we may use the ansatz $\psi(\bm{r})=e^{i\bm{k}_{||}\cdot\bm{r}_{||}}\phi(z)$. Eq. 58 then takes the form

[TABLE]

For $z<0$, the solution of Eq. 59 is given as

[TABLE]

where $k_{\perp}=\sqrt{k_{F}^{2}-\bm{k}_{||}^{2}}$. For $z>0$, we get

[TABLE]

where $\varphi$ is the azimuthal incidence angle. The momentum $q_{\pm}$ is given as

[TABLE]

The scattering matrix is found by relating the coefficients $A_{1,2}$ to $B_{1,2}$ via the matrix equation $B=SA$. This is done by enforcing continuity of $\phi$ and $\phi^{\prime}$ at $z=0$, and leads to

[TABLE]

Since the scattering matrix is unitary, $SS^{\dagger}=I$, it may be parametrized as

[TABLE]

where $\gamma$ is a spin mixing angle and $\hat{e}$ is a unit vector. By inspection, we see that we may define ${\hat{e}\cdot\bm{\sigma}=-\sigma_{x}\sin\theta\sin\varphi+\sigma_{y}\sin\theta\cos\varphi}$. To find an expression for $\gamma$, we define the constants $a=\sqrt{\frac{\varepsilon_{g}}{\varepsilon_{F}}}$ and $b=\frac{2m\alpha}{k_{F}}$. The former expresses the strength of the barrier potential, and the latter the strength of the spin–orbit coupling. Furthermore, we have $|\bm{k}_{||}|=k_{F}\sin\theta$, where $\theta\in[0,\pi/2]$ is the polar angle of incidence. We assume that the barrier potential is strong, so that an incoming particle only penetrates a short distance into the scattering region before being reflected, and hence $a\gg b$. From Eq. 64 we may then identify the spin mixing angle as

[TABLE]

with $K=\frac{b}{a^{3}}$ when we assume that $a\gg 1$, corresponding to a large band gap in the insulator. The parameter $\beta$ describes an overall phase, which is inconsequential for the boundary conditions, and hence we set $\beta=0$. In coordinate free form, the scattering matrix therefore takes the form,

[TABLE]

As before, $\hat{n}_{\alpha}$ is a unit vector that is either parallel or antiparallel to the interface normal, depending on the direction of symmetry breaking. In Nambu-space, the scattering matrix becomes

[TABLE]

with $\hat{\bm{\sigma}}=\text{diag}\left(\bm{\sigma}\;,\;\bm{\sigma}^{*}\right)$. Note that the spin mixing angle is antisymmetric in $\bm{k}$, e.g., $\gamma(\bm{k})=-\gamma(-\bm{k})$. Finally, the scattering matrix is diagonal in Keldysh space, e.g., $\check{S}=\text{diag}\left(\hat{S}\;,\;\hat{S}\right)$.

III.2 Boundary conditions

Finding the correct boundary conditions has now become a matter of identifying the terms in Eq. 51. To achieve this, we include only the lowest order tunnelling contributions,

[TABLE]

where $\check{\zeta}_{k}=(\hat{n}_{\alpha}\times\hat{k})\cdot\check{\sigma}\check{\rho}_{3}$. Furthermore, the spin mixing angle $\gamma$ is assumed to be small, and we therefore keep terms in the scattering matrix only up to second order. This gives

[TABLE]

In addition, we approximate

[TABLE]

and similarly for $\left[\check{S}\check{g}_{d,j}\check{S}^{\dagger}+\check{g}_{d,j}\right]^{-1}$, leading to

[TABLE]

By performing the first order expansion in this way, we ensure that in the limit of spin independent scattering, $\check{S}=\check{I}$, we get $\check{g}_{0,j}^{i}=\check{g}_{0,j}^{o}=\check{g}_{d,j}$. By inserting Eq. 70 we get

[TABLE]

Finally, we compute the full Green functions from Eqs. 40 and 41, and find the boundary conditions from Eqs. 51 and 52. Note that due to the angular averaging, all odd terms in $\check{\zeta}_{k}=(\hat{n}_{\alpha}\times\hat{k})\cdot\check{\sigma}\check{\rho}_{3}$ cancel, and thus we remove them immediately. Furthermore, the spin–orbit coupling is assumed to stem from the interface, which means that $\hat{n}_{\alpha}$ is parallel to $\hat{n}$. Since only even orders of the former appears, we may set $\hat{n}_{\alpha}=\hat{n}$. The matrix current is then given as

[TABLE]

with $d\Omega=\sin\theta\;d\theta\;d\phi$, $v_{n}=v_{F}\cos\theta$, and

[TABLE]

where we have neglected terms of order $w^{2}\gamma^{2}$. From the $\phi$ integration we find that, for an arbitrary matrix $\check{M}$,

[TABLE]

where $\check{\bm{\sigma}}_{||}=[\check{\bm{\sigma}}-\hat{n}\left(\hat{n}\cdot\check{\bm{\sigma}}\right)]\check{\rho}_{3}$, i.e., only spin directions parallel to the interface contribute to the boundary conditions.

The $\theta$ integration of the spin independent term in Eq. 75 diverges. However, when $\alpha=0$, we know that including all orders of the $T$ matrix, given in Eq. 43, yields a finite expression—namely, Eq. 53. This means that the divergence appears when the $T$ matrix is truncated to give Eq. 69. The interpretation of this is that microscopic analytical expressions for the coupling constants due to $w$, $w_{\alpha}$, and $\gamma$ cannot be found within the present theory, and they instead become input parameters. After the $\theta$ integration we therefore get

[TABLE]

where $T_{\alpha}^{\prime\prime}=\frac{1}{2}T_{\alpha}T_{\alpha}^{\prime}$. The parameter $T$ may be identified by comparing with Eq. 55, and is hence given by Eq. 56. The parameter $T_{\alpha}$ arises from the spin-orbit coupling part ($w_{\alpha})$ of the tunnelling potential in Eq. (57) whereas $T_{\alpha}^{\prime}$ arises from the interfacial spin-orbit coupling ($\alpha$) giving rise to a spin-mixing angle in Eq. (66). Eq. 77 is the main result of this paper.

A special case is worth commenting on. In the absence of any tunnelling, as is the case for a superconductor interfaced by a spin-orbit coupled insulator, only $T_{\alpha}^{\prime}$ is non-zero, whereas all other terms vanish, giving the boundary condition

[TABLE]

This boundary condition could thus be used to look for possible bound states induced at the interface of a superconductor and a spin–orbit coupled insulator. The fact that the $T_{\alpha}^{\prime}$ term exists despite the absence of a material to tunnel into on the other side of the interface is clear from the fact that this term only depends on $\check{g}_{1}$. In this sense, it may be thought of as the spin-orbit coupled equivalent of the spin-dependent phase-shift term $G_{\phi}$ previously discussed in the context of ferromagnetic insulators cottet_prb_2009 ; eschrig_njp_2015 .

IV Applications

In the following we will apply the boundary conditions derived in Eq. 77 to a set of example problems. In particular, we will consider superconducting hybrid structures in which a non-superconducting material is proximitized to a superconductor. The boundary conditions are assumed to represent a thin intermediary layer of a material with strong spin–orbit coupling. We neglect the inverse proximity effect, in which case the superconductors are approximately described by the Bardeen-Cooper-Schrieffer (BCS) bulk Green function,

[TABLE]

where $\theta=\arctan\frac{\Delta}{\varepsilon}$ for a given quasiparticle energy $\varepsilon$ and gap size $\Delta$, and $\phi$ is the superconducting phase. In other words, we solve the Usadel equation only in the non-superconducting material. We note that $\left[\check{g}_{\text{BCS}},\check{\bm{\sigma}}_{||}\right]=0$, and hence the boundary conditions simplify to

[TABLE]

For simplicity we will in the following set $T_{\alpha}=T_{\alpha}^{\prime}$.

IV.1 SF bilayer

As a first example, we consider a bilayer consisting of a superconductor and a ferromagnet, as shown in Fig. 2. The exchange field in the ferromagnet is directed at an angle $\theta$ from the plane of the interface, with a strength of $|\bm{h}|=2\Delta$. In this case the Usadel equation takes the following form,

[TABLE]

When both spin–orbit coupling and a magnetization is present in the system, this may lead to long range spin triplet superconducting correlations. Indeed, such correlations were found in a previous work which considered a similar system jacobsen_prb_2015 . There, a ferromagnet with Rashba spin–orbit coupling was attached to a superconductor via conventional spin independent boundary conditions. Long range spin triplet correlations were then observed as a zero energy peak in the density of states, the size of which depended upon the angle of the exchange field, $\theta$. Here, we seek to explore whether similar results emerge when the sole contribution to the spin–orbit coupling stems from the boundary conditions.

To quantify the presence of long range spin triplet correlations, we compute the density of states, which is given as

[TABLE]

where $N_{0}$ is the density of states at the Fermi level, and $g_{\sigma\sigma}(\bm{R},\varepsilon)$ are spin components of the normal Green function. The presence of long range spin triplets can be inferred from density of states at $\varepsilon=0$, at which point Eq. 82 may be expressed in terms of the contributions from the anomalous Green function $f=\left(f_{s}I+\bm{f}_{t}\cdot\bm{\sigma}\right)i\sigma_{y}$. With this particular parametrization, where the scalars $f_{s}$ and $\bm{f}_{t}$ give the singlet and triplet parts of $f$, respectively, Eq. 82 takes the form

[TABLE]

For an exchange field direction indicated by the unit vector $\hat{h}$, the triplet correlation is decomposed into a parallel, $\bm{f}_{||}=\left(\bm{f}_{t}\cdot\hat{h}\right)\hat{h}$, and an orthogonal, $\bm{f}_{\perp}=\bm{f}_{t}-\bm{f}_{||}$, component. The motivation for this decomposition is that the spin expectation value of the triplet Cooper pairs is given as $\left\langle\bm{S}\right\rangle\propto i\bm{f}_{t}(\varepsilon)\times\bm{f}_{t}^{*}(-\varepsilon)$ leggett_rmp_1975 ; linder_np_2015 . The triplet Cooper pairs whose spins are aligned with the exchange field, will not experience a pair breaking effect, in contrast to orthogonal spin alignments, and hence $\bm{f}_{||}$ and $\bm{f}_{\perp}$ expresses the short and long range triplet correlations, respectively. From Eq. 83 it is clear that while the presence of singlet superconducting correlations causes a suppression of the density of states, triplet correlations lead to an increase, and thus a potential for the formation of a zero energy peak.

Fig. 3a) shows the change in the density of states at zero energy as the exchange field is rotated away from the interface, $\delta\nu(\theta)\equiv\nu(\varepsilon=0,\bm{R}_{\text{max}},\theta)-\nu(\varepsilon=0,\bm{R}_{\text{max}},0)$, where $\bm{R}_{\text{max}}$ is the location at which the maximum density of states is found. A modulation of the zero energy peak is found, similar to the results of Ref. jacobsen_prb_2015 . At $\theta=0$ and $\theta=\pi$, i.e., for an exchange field parallel to the interface, the $T_{\alpha}^{\prime}$-dependent terms of Eq. 80 do not contribute, and the boundary conditions reduces to a conventional, spin independent tunnelling barrier. As $\theta$ is increased from zero, so too is the zero energy peak of the density of states, indicating the generation of triplet Cooper pairs. However, as $\theta$ approaches $\pi/2$, a dip is found instead. These results are further elucidated in Fig. 3b), which shows the angular dependence of the singlet and triplet correlations. The largest modulation is clearly seen in the long range triplets, $\bm{f}_{\perp}$, which is nonzero only when the exchange field has both an in-plane and an out-of-plane component with respect to the interface. In other words, it vanishes for $\theta\in\left\{0,\frac{\pi}{2},\pi\right\}$, in agreement with Ref. jacobsen_prb_2015 . Interestingly, it is not purely sinusoidal, but has maxima that are slightly tilted towards $\theta=\pi/2$. A small angular dependence in the singlet, $f_{s}$, and the short range triplet, $\bm{f}_{||}$, is also observed. At $\theta=\pi/2$, $f_{s}$ has a slight increase, whereas $|\bm{f}_{||}|$ decreases. This explains the reduction in the zero energy peak of the density of states at $\theta=\frac{\pi}{2}$.

IV.2 SFS Josephson weak link

When two superconductors are separated by a non-superconducting material, they may form a Josephson weak link. When a phase difference $\Delta\phi$ is induced between the superconductors of such systems, for instance by applying a current bias, dissipationless charge currents will flow between them, mediated by the Cooper pairs present in the nonsuperconducting material due to the proximity effect. It is well known that when the intermediary layer consists of a ferromagnet with an inhomogeneous magnetization, long range triplet Cooper pairs may be generated, which are spin polarized, and hence may carry a dissipationless spin current. It was recently predicted that a spin current may also emerge in homogeneous ferromagnets if thin normal metal layers with strong spin–orbit coupling are added between the ferromagnet and the superconductors jacobsen_scirep_2016 . To achieve this, the spin–orbit coupling was introduced in thin separate layers, coupled to the surrounding layers by tunnelling barriers, within which the Usadel equation was solved. While experimental verification of these results has proven elusive satchell_prb_2018 ; satchell_arxiv_2019 , the theoretical predictions provide an excellent benchmark for the new boundary conditions, as similar results should be obtained when the spin–orbit coupling is introduced as an interface effect. To verify this, we consider the system illustrated in Fig. 4. A homogeneous exchange field is defined in the ferromagnet, with a strength of $|\bm{h}|$, pointing in a direction $\theta$ relative to the transversal direction of the weak link. For this system, we compute the spin supercurrent, which in equilibrium is found from the Green function as

[TABLE]

where it has been assumed that junction is aligned along the $x$ axis, and $I_{s,0}=\frac{N_{0}DW}{8}$, for a given junction width $W$, and $\beta=1/k_{\text{B}}T$, with $T$ the temperature.

The results are given in Fig. 5, which shows the spin current component aligned parallel to the exchange field, in other words the spin current induced by the long range triplets. In Fig. 5a) its dependence on the canting angle $\theta$ is shown. It is noticed that the spin current goes to zero for $\theta=0$ and for $\theta=\frac{\pi}{2}$. This means that an exchange field with both an in-plane and an out-of-plane component is required in order to observe an effect, similarly to the SF bilayer. In Fig. 5b) we show the dependence of the spin current on the phase difference $\Delta\phi$ between the superconductors. It is seen that the current phase relation is approximately sinusoidal, similar to the conventional Josephson effect, indicating that the charge currents have become spin polarized. Finally, we show in Fig. 5c) the maximum spin current as a function of the interface spin–orbit coupling $T_{\alpha}$. For low values of the spin–orbit coupling, the spin current has an approximately parabolic form, but reaches a plateau as $T_{\alpha}$ approaches 0.5, the maximum value investigated in this study. A possible interpretation of this is that we are nearing the edge of the domain of validity for the small angle approximation used in the derivation of the boundary conditions, which requires $T_{\alpha}$ to be small.

V Conclusion

We have derived a new set of boundary conditions for systems in which there is large spin–orbit coupling. This allows the study of, for instance, superconducting hybrid structures with thin heavy metal layers. We demonstrate the use of these boundary conditions by considering an SF bilayer and an SFS Josephson weak link. In both cases we find that whenever the exchange field of the ferromagnet has both an in-plane and an out-of-plane component, long range triplet superconductivity is induced. The findings reported herein are consistent with results found in previous works, where the spin–orbit coupling is approximated by other means.

Acknowledgements.

The authors are grateful to M. Eschrig for valuable input. The authors also wish to thank V. Risinggrd, J. A. Ouassou, and J. R. Eskilt for useful discussions. J. L. acknowledges funding from the Research Council of Norway Center of Excellence Grant Number 262633, Center for Quantum Spintronics, and from the Research Council of Norway Grant No, 240806. J. L. and M. A. also acknowledge funding from the NV-faculty at the Norwegian University of Science and Technology.