Introduction

In spintronics, spin current plays an essential role to transfer the information associated with the spin degrees of freedom. Therefore the generation of spin current is an important issue, and several methods have been proposed and experimentally verified1. Various methods, such as the spin-polarized current injection from the ferromagnet2,3,4, spin battery5,6,7,8,9, and spin Hall effect10,11,12,13,14,15,16, have been employed to generate the spin current. It has been proposed also that the spin current second order in the electric field can be generated in noncentrosymmetric systems with spin-orbit interaction17,18. The spin currents mentioned above are either carried by itinerant electrons through dissipative electric currents, or by localized magnetic moments through exchange interactions. However, there can also be non-dissipative spin current in itinerant-electron systems. In this case, it is an equilibrium spin current without dissipation19. Although such a spin current is detectable according to Sonin20, it does not contribute to the transport property in a set up where the spin current can flow in and out.

On the other hand, superconducting spintronics is an emerging field attracting recent interest21,22,23,24,25,26,27,28,29,30,31,32,33,34. Although the spin degrees of freedom are usually quenched in singlet superconductors, the triplet component can be finite in noncentrosymmetric superconductors, ferromagnet-superconductor hybrids and Josephson junctions, or odd-parity superconductors, where the spins become (partially) active. Therefore, it is an intriguing issue if one can generate a spin current in superconductors with zero or small dissipation. Actually, the spin supercurrent has been discussed in He3 related to the internal degrees of freedom39. Recently, Leurs et al.40 have reexamined the spin supercurrent in spin-orbit coupled systems, and classified it to the coherent and non-coherent parts, only the latter of which contributes to the continuity equation of the spin density and its generation or manipulation is the focus of superconducting spintronics.

The superconducting magnetoelectric effect35,36,37,38 is also a result of the spins of Cooper pairs being active. For example, the spin density (or magnetization) induced by a supercurrent has been discussed by Edelstein35 in the case where the Fermi energy is large compared to the spin splitting by the Rashba interaction. Although this is justified in many situations, there are also systems where the spin-orbit splitting is comparable or even larger than the Fermi energy, e.g., the interface between LaAlO\({}_{3}\) and SrTiO\({}_{3}\)41, where the electron density can be tuned by gating. Therefore, it is desirable to cover the wide range of parameters, e.g., chemical potential, the strength of Rashba interaction, and temperature.

Here, we study spin density and spin current produced by a superconducting current in a two-dimensional superconductor with Rashba spin-orbit interaction for wide parameter regions. We investigate in detail the properties of spin current in this system and find that the spin current generation efficiency is comparable to or even larger than that of normal state when normalized with the charge current density. The spin current we discuss here corresponds to the non-coherent part in the classification by Leurs et al.40 and hence contributes to the spin accumulation, unlike the equilibrium spin current in normal systems19. Furthermore, it is different from those in other superconducting spintronic setups involving ferromagnet-superconductor hybrids where the spin degrees of freedom of Cooper pairs become active due to the presence of ferromagnets. In contrast, what we obtain here is a bulk spin supercurrent induced by charge supercurrent without time-reversal symmetry breaking of the ground state. The spin degrees of freedom become active because of the spin-orbit coupling and the flow of a Cooper pair condensate. Similar to normal spin currents, such a spin supercurrent may be detected by connecting to a material that shows inverse spin Hall effect and converts the injected spin currents into a transverse voltage. The study of the bulk spin supercurrent carried by Cooper pairs is a significant step towards non-dissipative spintronics.

Results

Model

A two-dimensional superconductor (SC) with Rashba spin-orbit interaction (SOI) can be described by the following Hamiltonian,

$${\hat{H}}_{0}=\sum _{{\bf{k}}}\left[\left(\frac{{\hslash }^{2}{k}^{2}}{2m}-\mu\right) {c}_{s,{\bf{k}}}^{\dagger }{c}_{s,{\bf{k}}}+\alpha ({{\boldsymbol{\sigma }}}_{ss^{\prime} }\times \hslash {\bf{k}})\cdot \hat{z}{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}} +\Delta {c}_{\uparrow ,{\bf{k}}}^{\dagger }{c}_{\downarrow ,-{\bf{k}}}^{\dagger }+h.c.\right],$$
(1)

where \(\alpha\) is the Rashba SOI strength, \(\Delta\) is the SC order parameter and \({\boldsymbol{\sigma }}=\hat{{\bf{x}}}{\sigma }^{x}+\hat{{\bf{y}}}{\sigma }^{y}+\hat{{\bf{z}}}{\sigma }^{z}\) is a spatial vector of Pauli matrices. The constants \(m\) and \(\hslash\) are the electron mass and the Planck constant respectively. Summation over implicit indexes is assumed.

Without superconductivity, the normal-state band structure is schematically shown in Fig. 1a. When \(\mu \gg {E}_{{\rm{R}}}\), \({E}_{{\rm{R}}}=m{\alpha }^{2}/2\) being the Rasbha band splitting energy, the Fermi level is far above the band touching point. If \(\mu =0\), the Fermi level cuts through this point and the inner Fermi surface shrinks to a point. The band bottom is reached when \(\mu =-{E}_{{\rm{R}}}\).

Fig. 1
figure 1

Schematic pictures of the bands and the spin current. a The schematic normal state band structure of materials with Rashba splitting energy \({E}_{{\rm{R}}}\). \(E\) is the energy eigenvalue and \(k\) is the magnitude of the wave vector. The band edge is reached if the chemical potential \(\mu =-{E}_{{\rm{R}}}\). b Cartoon picture showing the generation of magnetization \({M}_{y}\) and spin current \({J}_{x}^{y}\) by a charge supercurrent \({J}_{x}\) in a two-dimensional superconductor (SC) with Rashba spin-orbit interaction

Full size image

The superconductivity order parameter considered here is spin-singlet s-wave. In systems with Rashba spin-orbit interaction and superconductivity, s-wave may be mixed with p-wave and a general theory should include all of them. It has been shown that the actual amplitudes of the two different pairing potential depend on the form of the interaction42,43. When only even-parity interaction is considered, the p-wave order parameter vanishes. We assume pure s-wave order parameter in our calculation for simplicity.

Although we do not assume any spin-triplet order parameters in the Hamiltonian, the anomalous Green’s function \({F}(\omega ,{\bf{k}})=({\psi }_{s}+{\bf{d}}\cdot {\boldsymbol{\sigma }})i{\sigma }_{y}\), which carries the full information about the Cooper pairs, has both spin-singlet part \({\psi }_{s}\) and spin-triplet components \({\bf{d}}\) when spin-orbit interaction is present43. As a result, the Cooper pairs become spin-active due to SOI even though the pairing potential is purely s-wave. When time-reversal symmetry is broken, such as by the supercurrent \({J}_{x} \sim {q}_{x}\) in our case, the spin of the Cooper pairs become polarized and the polarization is given by \({\bf{S}} \sim i{\bf{d}}\times {{\bf{d}}}^{* }\) (more discussion in Supplementary Note 1). Since the pairs carry both polarized spin and momentum, they generate spin currents. This is illustrated schematically in Fig. 1b.

To investigate the spin and charge properties, it is convenient to introduce an SU(2) gauge field \({{\bf{W}}}^{\nu }{\sigma }^{\nu }\)40,44,45 along with the electromagnetic field \({\bf{A}}\). Electrons coupled with both fields are described by the Hamiltonian

$$\hat{H}=\sum _{{\bf{k}}}\left\{{\left[\frac{1}{2m}{\left(\hslash {\bf{k}}+e{\bf{A}}-\frac{\hslash }{2}{{\bf{W}}}^{\nu }{\sigma }^{\nu }\right)}^{2}-\mu \right]}_{ss^{\prime} }{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}}+\Delta {c}_{\uparrow ,{\bf{k}}}^{\dagger }{c}_{\downarrow ,-{\bf{k}}}^{\dagger }+h.c.\right\}.$$
(2)

Apparently, the Rashba term in \({\hat{H}}_{0}\) breaks SU(2) symmetry. It is actually equivalent to a constant gauge field (neglecting a numeric constant) \(\bar{{\bf{W}}}=({\bar{W}}_{y}^{x}{\sigma }^{x},{\bar{W}}_{x}^{y}{\sigma }^{y},0)\), with \({\bar{W}}_{y}^{x}=-{\bar{W}}_{x}^{y}=2m\alpha /\hslash\)46. In this formulism, the spin current operator is conveniently defined as \({\hat{{\bf{J}}}}^{\nu }=\partial \hat{H}/\partial {{\bf{W}}}^{\nu }\). Or, starting from the free energy, the spin current can be obtained as \({{\bf{J}}}^{\nu }=-\partial {\mathcal{F}}/\partial {{\bf{W}}}^{\nu }\).

When a uniform supercurrent is flowing (say along \(x\)-direction), the SC order parameter acquires a constant phase gradient, i.e. \(\Delta ={\Delta }_{0}\exp [2i{q}_{x}x]\). Or, by a gauge transformation, it is equivalent to a constant vector potential \({A}_{x}=\hslash {q}_{x}/e\), which is effective only for SCs. Before showing the results, it is helpful to discuss the symmetry properties of the involved physical quantities. The spin current (spin density) is odd (even) under spatial inversion but even (odd) under time-reversal operation, while charge current, or \({q}_{x}\), is odd under both of them. Consequently, the lowest order of the spin current is \({J}^{\nu } \sim {q}_{x}^{2}\) and that of spin density is \({M}_{\nu } \sim {q}_{x}\). Both of them must be odd functions of \(\alpha\).

Spin density and spin current

At zero temperature, the free energy \({\mathcal{F}}\) is just the ground state energy. For arbitrary Rashba SOI strength and chemical potential, the spin polarization is along \(y\)-direction and the density is (derivation in the Methods section)

$${M}_{y}=\frac{m\alpha {q}_{x}}{4\pi }g\left(\frac{\Delta }{{E}_{{\rm{R}}}},\frac{\mu +{E}_{{\rm{R}}}}{{E}_{{\rm{R}}}}\right),$$
(3)
$$g(d,u)=\frac{1}{2}-\frac{1}{2}{\int _{0}^{1}}\frac{{x}^{2}-u}{\sqrt{{d}^{2}+{\left({x}^{2}-u\right)}^{2}}}dx,$$
(4)

where we have defined the Rashba splitting energy \({E}_{{\rm{R}}}=m{\alpha }^{2}/2\). This is a generalized SC Edelstein effect35 in large-Rashba systems. The \(\mu\)-dependence of the spin density \({M}_{y}\) is shown in Fig. 2a with various values of the pairing potential. When \(\mu \gg {E}_{{\rm{R}}}\gg \Delta\), we get \({M}_{y}(\mu \gg {E}_{{\rm{R}}})=m\alpha {q}_{x}/4\pi\), in which the dependence on \(\mu\) is absent. When \(\mu < 0\), the spin density gradually decreases as \(\mu\) goes down.

Fig. 2
figure 2

Spin density and spin current density for a given Cooper pair wave vector \({q}_{x}\). a The spin density \({M}_{y}\) as a function of the chemical potential \(\mu\). The Rashba splitting energy \({E}_{{\rm{R}}}\) is used as the reference energy scale. (The explicit dependence on the spin-orbit interaction strength is shown in Supplementary Note 6.) Note that the band edge in the normal state is reached when \(\mu =-{E}_{{\rm{R}}}\) below which the normal state is insulating. The inset shows the dependence of the same quantity on the order parameter \(\Delta\) for chosen values of \(\mu\). The quantity \(m\) is the electron mass and \(\alpha\) is the Rashba strength. b The spin current as a function of chemical potential with various values of \(\Delta\). Peaks appear near the band edge. The inset shows how the peak position \({\mu }_{\text{p}eak}\) and peak height \({j}_{\,\text{m}\,ax}^{y}/{q}_{x}^{2}\) change as we increase \(\Delta\). In all the results of this paper, we have set the electron charge \(e=1\) and the Planck constant \(\hslash =1\)

Full size image

In Fig. 2a, one readily notes that the spin density remains finite when \(\Delta\) approaches zero. On the other hand, if \(\Delta =0\), i.e. in a normal state, spin density due to the vector potential must vanish since one can trivially gauge out \({q}_{x}\). Thus, there is a discontinuity in \({M}_{y}/{q}_{x}\) at \(\Delta \to 0\). However, as we will see later, \({M}_{y}\) actually vanishes when \(\Delta\) decreases to zero since the supercurrent \({J}_{x}\) (and thus \({q}_{x}\)) also approaches zero, as discussed in Supplementary Note 2. Also, this discontinuity is absent at finite temperature, as shown in Supplementary Note 3. The \(\Delta\)-dependence of the function \(g\) is shown in the inset of Fig. 2a, where a finite value is obtained at \(\Delta \to 0\) for \(\mu /{E}_{{\rm{R}}}> -\!\!1\). In the limit \(\Delta /{E}_{{\rm{R}}}\to 0\), Eq. (3) simply becomes

$${M}_{y}=\frac{m\alpha {q}_{x}}{4\pi }\left\{\begin{array}{cc}1,&{\rm{when}}\ \mu \ > \ 0,\\ \sqrt{\mu /{E}_{{\rm{R}}}+1},&{\rm{when}}\ 0\ > \ \mu \ > \ -\!{E}_{{\rm{R}}},\\ 0,&{\rm{when}}\ \mu \ < \ -\!{E}_{{\rm{R}}}.\end{array}\right.$$
(5)

Due to the form of the Rashba SOI, and as indicated by the direction of the spin density obtained above, the only non-vanishing spin current induced by the supercurrent in \(x\)-direction is \({{\bf{J}}}^{y}\). In the limit where the chemical potential is high and the order parameter is small, i.e. \(\mu \gg {E}_{{\rm{R}}}\gg \Delta\), the spin current density is (derivation in the Methods section)

$${{\bf{j}}}_{+}^{y}=\alpha \hslash {q}_{x}^{2}\hat{{\bf{x}}}/16\pi ,$$
(6)

and all other spin current components vanish. This expression is independent of \(\mu\) and \(\Delta\), similar to the spin density discussed previously.

When \(-{E}_{{\rm{R}}} \, < \, \mu \, < \, 0\) and \(\Delta\) is small, the spin current density can be written as

$${{\bf{j}}}_{-}^{y}=\frac{\alpha \hslash {q}_{x}^{2}}{32\pi }f\left(\frac{\Delta }{{E}_{{\rm{R}}}},\frac{\mu +{E}_{{\rm{R}}}}{{E}_{{\rm{R}}}}\right)\hat{{\bf{x}}},$$
(7)
$$f(d,u)=-5\sqrt{u}+{\int _{0}^{1}}\frac{{d}^{2}(4{x}^{2}+3)}{{[{d}^{2}+{({x}^{2}-u)}^{2}]}^{3/2}}dx.$$
(8)

Particularly, when \(d=\Delta /{E}_{{\rm{R}}}\to 0\), the function \(f(d,u)\) simplifies to \(f(d\to 0,u)=3/\sqrt{u}-\sqrt{u}\). Combined with Eq. (7), this leads to the same spin current density as Eq. (6) when \(u=1\) (corresponding to the band crossing point \(\mu =0\)). Also, the spin current increases because of the increase of the density of states when the chemical potential goes down until it reaches the band edge at \(\mu \to -{E}_{{\rm{R}}}\) where it diverges, as shown by the dashed curve in Fig. 2b. The two Fermi surfaces have the same (opposite) spin textures and the contributions from the two sub-bands add up (tend to cancel) when \(\mu\) is negative (positive). While the spin current is similar to spin density for large \(\mu\) in the sense that they both keep constant, it shows different behavior when \(\mu \, < \, 0\). Note that the \(\mu\)-dependence of spin current resembles the density of states of the normal Rashba band, which also diverges at the band edge.

For general parameters, the spin current of the SC is calculated numerically and shown in Fig. 2b. Because of finite \(\Delta\), the aforementioned divergence of the spin current at the band edge is smoothed out while the spin current at \(\mu \, > \, 0\) is hardly changed even the pairing potential reaches \(0.5{E}_{{\rm{R}}}\). The change of the peak position (\({\mu }_{\text{p}eak}\)) and peak height (\({j}_{\,\text{m}\,ax}^{y}/{q}_{x}^{2}\)) as functions of \(\Delta\) are shown in the inset of Fig. 2b. It turns out that the peak position shifts almost linearly from the band edge when the order parameter increases while the height of the peak decreases in a nonlinearly.

Bose-Einstein condensation regime

When the chemical potential is below the band edge, i.e. \(\mu \, < -\!{E}_{{\rm{R}}}\), the normal electronic state has no Fermi surface and SC cannot happen in the weak coupling limit of the Bardeen-Cooper-Schrieffer (BCS) theory. Assuming that the SC exists, it must be in BEC regime where electrons are tightly bond and the critical temperature is determined by the Bose-Einstein condensation (BEC) temperature47,48. As shown in Fig. 2, spin density and spin current quickly drops to zero when \(\mu \, < -\!{E}_{{\rm{R}}}\) if \(\Delta\) is small. However, when the pairing potential is large and a BEC superconductor is achieved, they become finite. Especially, when \(-u\gg d\) and \(-u\gg 1\), direct expansions of  Eq. (3) and Eq. (7) respect to \(d\) and \({u}^{-1}\) give \(M_y \approx \frac{m \alpha q_x}{32\pi}{d}^{2}/{u}^{2}\) and \({{\mathbf{j}}}_{\text{BEC}}^{y}=\alpha\hslash{q}_{x}^{2}{d}^{2}{(12\pi)}^{-1}|u{|}^{-3}\hat{{\mathbf{x}}}\). Both the spin density and spin current show power-law decay as \(\mu\) decreases. Note that the discontinuity at \(\Delta \to 0\) does not appear when \(\mu \, < -\!{E}_{{\rm{R}}}\).

Efficiency

To relate our calculation to experiments, it is useful to convert the wave vector \({q}_{x}\) to the supercurrent density \({j}_{x}\). To calculate \({j}_{x}\), we note that (charge) current operator in a Rashba system is \({\hat{J}}_{x}=e\int {d}^{2}{\bf{k}}{[\frac{\hslash }{m}({k}_{x}+{q}_{x})-\alpha {\sigma }^{y}]}_{ss^{\prime} }{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}}.\) At \(T=0\), the usual paramagnetic term vanishes while the diamagnetic term remains. The last term proportional to \(\alpha\) contributes a supercurrent of \(e\alpha {M}_{y}/\hslash\). So the zero temperature supercurrent density can be written as

$$\begin{array}{ll}{j}_{x}&=e\hslash {q}_{x}{n}_{e}(u)/m-2e\alpha {M}_{y}/\hslash \\ &=[\hslash {n}_{e}(u)/m-g(d,u){E}_{{\rm{R}}}/(\pi \hslash)]e{q}_{x}.\end{array}$$
(9)

The quantity \({n}_{e}\) is the electron density. Eq. (9) may be regarded as the generalized London equation for Rashba systems.

Since \({j}_{x} \sim {q}_{x}\), we define two coefficients

$$\gamma ={M}_{y}/{j}_{x},\eta ={j}^{y}/{j}_{x}^{2},$$
(10)

which denote the efficiency of spin density and spin-current generation, respectively, for given charge current density. They are shown in Fig. 3. For large \(\mu\), both the spin density and spin current decrease in power-law, \(\gamma \sim {(\mu +{E}_{{\rm{R}}})}^{-1}\), \(\eta \sim {(\mu +{E}_{{\rm{R}}})}^{-2}\). When \(\mu\) goes below zero, both \(\gamma\) and \(\eta\) increase. Interestingly, if \(\Delta\) is small compared to Rashba splitting, \(\gamma\) stays constant for \(\mu \, < \, 0\) until it reaches the band edge, below which it decreases again. For spin current, the ratio \(\eta\) keep increasing before \(\mu\) reaches the band edge. Below that, \(\eta\) decreases very slowly. That means the efficiency of spin current generation in the BEC regime is very high.

Fig. 3
figure 3

Spin density \({M}_{y}\) and spin current density \({j}^{y}\) for a given charge current density \({j}_{x}\). a The logarithm of spin density generation efficiency, defined as, \(\gamma ={M}_{y}/{j}_{x}\) as a function of the chemical potential \(\mu\) for various values of the superconductivity order parameter \(\Delta\). The black dashed curve is for the normal state. b The logarithm of spin current generation efficiency, defined as \(\eta ={j}^{y}/{j}_{x}^{2}\), as a function of the chemical potential \(\mu\) for various values of the \(\Delta\). The dashed colorful curves are for the quantity \({\eta }_{0}={j}_{0}^{y}/{j}_{x}^{2}\) where \({j}_{0}^{y}=({j}_{x}/e)({M}_{y}/{n}_{e})\) denotes a direction product of the particle current \(({j}_{x}/e)\) and the spin polarization per particle \(({M}_{y}/{n}_{e})\). \({n}_{e}\) denotes the electron density and \(e\) is the electron charge. The same color corresponds to the same value of \(\Delta\). The insets are the log–log plot of \(\eta\) emphasizing the regime with negative \(\mu\)

Full size image

Relation between spin density and spin current

As shown above, the functional behaviors of the spin density and spin current are similar in some parameter regions but very different in others. To further clarify the mechanism of the spin-current generation, we investigate its relation to the spin polarization.

When spin polarization is induced by supercurrent, we may intuitively expect a spin current

$${j}_{0}^{y}=({M}_{y}/{n}_{e})({j}_{x}/e)={\eta }_{0}{j}_{x}^{2},$$
(11)

which is just the particle current times the spin polarization of each particle. Using Eqs. (9)–(11) and Eq. (3), we get \({\eta }_{0}=\gamma e/{n}_{e}\) and \(\gamma =\frac{m\alpha \hslash }{4e{E}_{{\rm{R}}}}\frac{g(d,u)}{\pi {\hslash }^{2}{n}_{e}(u)/(m{E}_{{\rm{R}}})-g(d,u)}.\) In Fig. 3b, by comparing \(\eta\) (solid curves) with \({\eta }_{0}\) (dashed curves in corresponding colors), we find that the differences between them are not small in general. Thus, Eq. (11) does not fully describe the origin of the spin current because the spin and momentum degrees of freedom cannot be separated in systems with SOI.

Comparison to normal states

In normal dissipative Rashba systems, spin density49 and spin current17,18 can also be generated when an external electric field \({\bf{E}}\) is applied. In such a case, the coefficients \(\gamma\) and \(\eta\) are shown by the black dashed curves in Fig. 3. It turns out that the spin polarization efficiency for SC state and that for normal state are almost the same if \(\Delta\) is small. For spin current generation, \(\eta\) for the SC state is smaller than the normal state in most of the parameter regime. Only when \(\mu\) is very close to the band edge and \(\Delta\) is small, the SC state shows a higher efficiency. The comparison of the normal state and the SC state becomes clearer in the log-log plots as shown in the insets of Fig. 3. When \(\mu \, < -\!\!{E}_{{\rm{R}}}\), it is not possible to pass a current in the normal state at \(T=0\). However, there could be BEC SCs in this limit, which can carry supercurrent and gives nonzero spin densities and spin currents.

Although a dissipative normal system seems actually better at generating spin current than its SC counterpart for most of the parameter space, we should note that, in the normal state, the spin current, as well as the charge current, is generated by an electric field and suffers from dissipation and heating, making it impossible to maintain long-range coherence. On the other hand, the spin current in the SC state has very low dissipation (if any) because it is purely due to the flow of a supercurrent. In other words, it is the property of the Cooper pair condensate and long-range coherence is guaranteed, making it a spin supercurrent27,50,51,52.

It should be emphasized that the spin current in a coherent system does not necessarily correspond to the coherent spin current as defined by Leurs et al.40 as mentioned previously. In the formalism of non-Abelian gauge fields, the spin current operator can be decomposed into two parts, \({\hat{{\mathbf{J}}}}^{a}={\hat{{\mathbf{J}}}}_{\,\text{NC}}^{a}+{\hat{{\mathbf{J}}}}_{\,\text{C}\,}^{a}\). The first term \({\hat{{\mathbf{J}}}}_{\,\text{NC}}^{a}=\hat{{\mathbf{J}}}{\hat{S}}^{a}\) is a direct product of the charge current operator and the spin operator and is called non-coherent. The second term \({\hat{{\mathbf{J}}}}_{\,\text{C}\,}^{a}\) involves spin entanglement, according to Leurs at al., and is called coherent. In our calculation, the spin current operator is actually \({\hat{{\mathbf{J}}}}_{\,\text{NC}}^{a}\). Such spin supercurrent carried by Cooper pairs does induce spin transport and may be used to create new spintronic devices with long-range coherence.

Discussion

We have shown that a supercurrent in an superconductor with Rashba SOI can induce spin density and spin current. Let’s take the interface between a LaAlO\({}_{3}\) and a SrTiO\({}_{3}\) as an example and estimate the realistic magnitude of these effects. The critical current density (2D) is \(\sim \!\!1{0}^{-3}{\mathrm{A/cm}}\), effective mass \({m}^{* }=1.5{m}_{e}\) and the super-fluid density \({n}_{s} \sim 1{0}^{13} \, {\mathrm{cm}}^{-2}\)53. Assuming \({E}_{{\rm{R}}} \sim 10 \, {\mathrm{meV}}\) (\(\alpha \sim 1{0}^{5} \, {\mathrm{m}}/ \!{\mathrm{s}}\)) and replacing the electron density \({n}_{e}\) by the super-fluid density \({n}_{s}\), the second term in Eq. (9) has the same order of magnitude as the first term when the Fermi level is high, i.e. or \(u\gg 1\). For a super-current density of \({j}_{x} \sim 1{0}^{-5} \, {\mathrm{A/cm}}\), which is about \(1/100\) of the critical current density, the corresponding \({q}_{x} \sim 1{0}^{3} \, {\mathrm{m}}^{-1}\). Then, according to Eqs. (6) and (3), the spin current density \({j}^{y} \sim 1{0}^{8}\hslash \cdot {\mathrm{s}}^{-1}{\mathrm{cm}}^{-1}\) and spin density \({M}_{y} \sim 1{0}^{7}\hslash \cdot {\mathrm{cm}}^{-2}\). Converted to charge current by replacing \(\hslash /2\) by \(e\), such a spin current corresponds to an charge current density of \(\sim \!\!0.01 \, {\mathrm{nA/cm}}\). This is rather small. However, If the Fermi level is decreased, say by gating41,54,55,56, and it is still superconducting, the spin current can be enhanced by several orders of magnitude as shown in Fig. 3b. In that case, it should be detectable by inverse spin Hall effect57, by connecting the superconductor with a light-emitting diode58, or by other methods. Our results may be generalized to other spin-orbit coupled systems where general properties of the Edelstein effect have been recently studied38.

Another important issue is the effect of finite temperature which affects the results through both the Fermi distribution of the quasi-particles and change of \(\Delta\). In the following, we assume that the function \(\Delta (T)\) is determined by the BCS theory. Fig. 4a shows the spin density as a function of temperature for different values of \(\mu\). Generally, the spin density starts to grow linearly when \(T \, < \, {T}_{c}\), and it saturates when \(T\to 0\). It is exactly the same feature as the superfluid density \({\rho }_{s}\) in standard BCS theory. This is no surprise since Eq. (9) indicates that the spin density is part of the supercurrent. When the Fermi level moves downward, the temperature dependence remains the same. In Fig. 4b, we show the temperature dependence of the spin current. Further discussion about temperature effect can be found in Supplementary Note 3. The temperature dependence in the case of normal state is shown in Supplementary Note 4 for comparison. The spin current generation efficiency as a function of chemical potential is obtained for finite temperature in Supplementary Note 5 which shows that our conclusions are still valid.

Fig. 4
figure 4

Temperature effect. a The spin density (with a given Cooper pair wave vector \({q}_{x}\)) as a function of temperature for various values of the chemical potential \(\mu\). The Rashba splitting energy \({E}_{{\rm{R}}}\) is used as the energy scale. The superconductivity critical temperature is set to be \({T}_{\text{c}}=0.1{E}_{{\rm{R}}}/{k}_{{\rm{B}}}\) with \({k}_{{\rm{B}}}\) being the Boltzmann constant. b The spin current density (with a given \({q}_{x}\)) as a function of temperature for various values of the chemical potential \(\mu\)

Full size image

One may also note that this system is a paradigm for the study of topological superconductors. In fact, an out-of-plane Zeeman field can drive it into a topological phase with Majorana edge modes if \(\mu \sim 0\). Then the following questions become interesting: (1) What is the spin current contributed by the edge modes, if any? (2) How are topological phase transitions affected by the supercurrent? In fact, the second question has recently been discussed but in a quasi-one-dimensional system where it is found that a supercurrent shifts the topological phase transitions59. Similar things may happen for two-dimensional systems as well. In this way, one may control topological phases by a supercurrent, achieving new methods for coherent quantum-device manipulation.

To conclude, we have studied the spin and spin current generation by the supercurrent in two-dimensional noncentrosymmetric superconductor with Rashba spin-orbit interaction. The spin degrees of freedom are partially activated due to the noncentrosymmetry even when on-site pairing between up and down spins is considered. When the chemical potential is below the band crossing point and approaching to the Band edge, i.e., in dilute electron density limit, the large enhancement of the spin supercurrent generation occurs although the spin density is small. The efficiency of the spin supercurrent does not decrease so much even when the chemical potential is below the band edge, i.e., in the BEC limit of the superconductivity. Furthermore, the carrier density can be controlled there by gating, and the chemical potential dependence can be studied experimentally. These studies will pave a route toward the dissipationless superconducting spintronics, where the transfer of spin information without the energy loss is possible through the long-range quantum coherence of the system.

Methods

Non-Abelian gauge fields and Rashba spin-orbit interaction

The Hamiltonian with SU(2) gauge symmetry of free electrons in magnetic field is

$${\hat{H}}_{N}=\sum _{{\bf{k}}}{c}_{{s}_{1},{\bf{k}}}^{\dagger }{H}_{{s}_{1},{s}_{2}}({\bf{k}}){c}_{{s}_{2},{\bf{k}}},$$
(12)
$$H({\bf{k}})=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}}-\bar{g}{{\bf{W}}}^{\nu }{\sigma }^{\nu })}^{2}+{\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu ,$$
(13)

where \({\bf{A}}\) is the U(1) gauge potential of electromagnetic field and \({{\bf{W}}}^{\nu }\) is the SU(2) gauge potential. The constant \({\mu }_{{\rm{B}}}\) denotes the Bohr magneton and \({\sigma }^{\nu =1,2,3}\) are Pauli matrices. The coefficient \(\bar{g}\) is a coupling constant to be assigned later. Rashba SOI corresponds to the existence of the following SU(2) gauge field46,

$$\bar{g}{\bar{W}}_{i}^{\nu }=\frac{e\hslash }{m{c}^{2}}{\epsilon }_{ij\nu }{E}_{j}.$$
(14)

Note that we have introduced the speed of light \(c\) and the Levi-Civita symbol \({\epsilon }_{ij\nu }\). The Hamiltonian with this gauge field can be rewritten as

$$H({\bf{k}})=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}})}^{2}+\frac{1}{2m}({\bar{W}}_{j}^{\nu }{\sigma }^{\nu }{\bar{W}}_{j}^{\gamma }{\sigma }^{\gamma }){\!}-\frac{1}{m}[(\hslash {k}_{i}+e{A}_{i}){\bar{W}}_{i}^{\nu }{\sigma }^{\nu }]+{\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }{\!}-\mu ,$$
(15)
$$=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}})}^{2}+\frac{e\hslash {E}_{j}}{{m}^{2}{c}^{2}}(\hslash {k}_{i}+e{A}_{i}){\epsilon }_{ji\nu }{\sigma }^{\nu }+{\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu +\frac{1}{2m}{\bar{W}}_{j}^{\nu }{\bar{W}}_{j}^{\nu },$$
(16)
$$=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}})}^{2}+\frac{e\hslash }{{m}^{2}{c}^{2}}{\bf{E}}\cdot (\hslash {\bf{k}}+e{\bf{A}})\times {\boldsymbol{\sigma }}+{\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu +\frac{1}{m}{\left(\frac{e\hslash }{m{c}^{2}}\right)}^{2}| {\bf{E}}{| }^{2}.$$
(17)

The second term describes the Rashba SOI. So, the Rashba SOI will be included by assuming a constant SU(2) gauge field \({\tilde{{\bf{W}}}}^{\nu }\) and the SU(2) gauge symmetric Hamiltonian with certain Rashba strength becomes

$${H}_{{\rm{R}}}({\bf{k}})=\frac{1}{2m}{[\hslash {\bf{k}}+e{\bf{A}}-\bar{g}(\delta {{\bf{W}}}^{\nu }+{\bar{{\bf{W}}}}^{\nu }){\sigma }^{\nu }]}^{2}$$
(18)
$$+ \, {\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu ,$$
(19)
$$=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}}-\bar{g}{\delta {\bf{W}}}^{\nu }{\sigma }^{\nu })}^{2}+\alpha (\hslash {\bf{k}}+e{\bf{A}})\times {\boldsymbol{\sigma }}$$
(20)
$$+ \, {\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu +\frac{\bar{g}^{2}}{m}\delta {W}_{j}^{\nu }{\bar{W}}_{j}^{\nu }+\frac{1}{m}{\left(\frac{e\hslash }{m{c}^{2}}\right)}^{2}| {\bf{E}}{| }^{2}.$$
(21)

The last two terms are just constant energy shift due to the electric field of the Rashba SOI, which we ignore hereafter. Then we have

$${H}_{{\rm{R}}}({\bf{k}})=\frac{1}{2m}{(\hslash {\bf{k}}+e{\bf{A}}-\bar{g}{\delta {\bf{W}}}^{\nu }{\sigma }^{\nu })}^{2}$$
(22)
$$+ \, \alpha (\hslash {\bf{k}}+e{\bf{A}})\times {\boldsymbol{\sigma }}+{\mu }_{{\rm{B}}}{B}^{\nu }{\sigma }^{\nu }-\mu .$$
(23)

In order for the non-Abelian gauge field \({{\bf{W}}}^{\nu }\) to couple with the spin, we set

$$\bar{g}=\hslash /2.$$
(24)

Derivation of spin density

We calculate the spin density induced by the supercurrent, or by the vector potential \({\bf{A}}=(\hslash {q}_{x}/e)\hat{{\bf{x}}}\) in Eq. (2), using perturbation method. The perturbation Hamiltonian \(\hat{H}^{\prime}\) including both the test fields (\({{\bf{W}}}^{\nu }\) and \({\bf{B}}\)) and the external field (\({\bf{A}}\)) is

$$\hat{H}^{\prime} =\delta {\hat{H}}_{A1}+\delta {\hat{H}}_{A2}+\delta {\hat{H}}_{W}+\delta {\hat{H}}_{AW}+\delta {\hat{H}}_{B},$$
(25)

with

$$\delta {\hat{H}}_{A1}=\sum _{{\bf{k}}}\left[-\frac{\hslash {\bf{k}}\cdot e{\bf{A}}}{m}{\delta }_{ss^{\prime} }+\alpha e{({A}_{x}{\sigma }^{y})}_{ss^{\prime} }\right]{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}},$$
(26)
$$\delta {\hat{H}}_{A2}=\sum _{{\bf{k}}}\left[-\frac{{e}^{2}| {\bf{A}}{| }^{2}}{2m}{\delta }_{ss^{\prime} }\right]{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}},$$
(27)
$$\delta {\hat{H}}_{W}=\sum _{{\bf{k}}}\left[\frac{-{\hslash }^{2}}{m}\right.{\bf{k}}\cdot \delta {{\bf{W}}}^{\nu }{\sigma }_{ss^{\prime} }^{\nu }$$
(28)
$$+ \, \alpha \hslash (\delta {W}_{x}^{y}-\delta {W}_{y}^{x}){\delta }_{ss^{\prime} }\bigg]{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}},$$
(29)
$$\delta {\hat{H}}_{AW}=\sum _{{\bf{k}}}\left[\frac{e\hslash }{m}{\bf{A}}\cdot \delta {{\bf{W}}}^{\nu }{\sigma }_{ss^{\prime} }^{\nu }\right]{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}},$$
(30)
$$\delta {\hat{H}}_{B}=\sum _{{\bf{k}}}{\mu }_{{\rm{B}}}{\bf{B}}\cdot {{\boldsymbol{\sigma }}}_{ss^{\prime} }{c}_{s,{\bf{k}}}^{\dagger }{c}_{s^{\prime} ,{\bf{k}}},$$
(31)

which is composed of terms linear and quadratic in \({\bf{A}}\) respectively, a term proportional to non-Abelian gauge field \({{\bf{W}}}^{\nu }\), a term bilinear in \({\bf{A}}\) and \({{\bf{W}}}^{\nu }\), and the term due to the field \({\bf{B}}\). As mentioned previously, the spin density \({M}_{y}\) is linear in \({q}_{x}\) (or \({A}_{x}\)) to the lowest order. Consequently, we should calculate the contribution to the free energy by \(\hat{H}^{\prime}\) up to the second order.

At zero temperature, the free energy is simply the ground state energy.

$${F}_{0}=\sum _{{\bf{k}},\pm }{E}_{h\pm }(k)=\frac{{L}^{2}}{4{\pi }^{2}}\int {E}_{h\pm }(k)kdkd\theta .$$
(32)

The subscript \(h\) labels the hole bands, i.e. \({E}_{h\pm }(k) \, < \, 0\). Here we ignored the energy due to Cooper pairs, assuming that the order parameter is not affected by the perturbation term. The unperturbed energy spectrum of the system is

$${E}_{e\pm }^{0}({\bf{k}})=\sqrt{{\xi }_{\pm }{({\bf{k}})}^{2}+{\Delta }^{2}},$$
(33)
$${E}_{h\pm }^{0}({\bf{k}})=-\sqrt{{\xi }_{\pm }{({\bf{k}})}^{2}+{\Delta }^{2}}.$$
(34)

\({\xi }_{\pm }=\frac{{\hslash }^{2}{k}^{2}}{2m}-\mu \pm \alpha \hslash k\) are the energy spectrum of the normal Rashba bands. The lowest order term of the spin density is from the bilinear term \(\sim {q}_{x}{B}_{i}\). Then, the perturbation correction to the free energy is

$$\delta {F}_{M}=\frac{{L}^{2}\hslash {B}_{y}{q}_{x}}{16\pi } \,\times {\int _{0}^{\infty }}dk\left(\frac{{\xi }_{-}(k)}{\sqrt{{\xi }_{-}{(k)}^{2}+{\Delta }^{2}}}-\frac{{\xi }_{+}(k)}{\sqrt{{\xi }_{+}{(k)}^{2}+{\Delta }^{2}}}\right).$$
(35)

The spin density is

$${M}_{y}=-\frac{1}{{L}^{2}}\frac{\partial \delta {F}_{M}}{\partial {B}_{y}}=\frac{\hslash {q}_{x}}{16\pi }{\int _{0}^{\infty }}dk\left(\frac{{\xi }_{+}(k)}{\sqrt{{\xi }_{+}{(k)}^{2}+{\Delta }^{2}}}-\frac{{\xi }_{-}(k)}{\sqrt{{\xi }_{-}{(k)}^{2}+{\Delta }^{2}}}\right).$$
(36)

After substituting the expressions of \({\xi }_{\pm }\) into the above integral and defining \({k}_{\pm }=k\pm m\alpha /\hslash\), it becomes

$$\begin{array}{ll}{M}_{y}=&\frac{\hslash {q}_{x}}{16\pi }\left[{\int _{\frac{\alpha m}{\hslash }}^{\infty +\frac{\alpha m}{\hslash }}}\frac{\frac{{k}_{+}^{2}{\hslash }^{2}}{2m}-\left(\mu +{E}_{{\rm{R}}}\right)}{\sqrt{{\Delta }^{2}+{\left(\frac{{k}_{+}^{2}{\hslash }^{2}}{2m}-\left(\mu +{E}_{{\rm{R}}}\right)\right)}^{2}}}\ d{k}_{+}\right.\\ &-\left. {\int _{-\frac{\alpha m}{\hslash }}^{\infty -\frac{\alpha m}{\hslash }}}\frac{\frac{{k}_{-}^{2}{\hslash }^{2}}{2m}-\left(\mu +{E}_{{\rm{R}}}\right)}{\sqrt{{\Delta }^{2}+{\left(\frac{{k}_{-}^{2}{\hslash }^{2}}{2m}-\left(\mu +{E}_{{\rm{R}}}\right)\right)}^{2}}}\ d{k}_{-}\right]\end{array}$$
(37)
$$\begin{array}{ll}=&\frac{\hslash {q}_{x}}{16\pi }\frac{\alpha m}{\hslash }\left[{\int _{1}^{\infty +1}}\frac{{x}^{2}-u}{\sqrt{{d}^{2}+{\left({x}^{2}-u\right)}^{2}}}dx\right.\\ &-\left. {\int _{-1}^{\infty -1}}\frac{{x}^{2}-u}{\sqrt{{d}^{2}+{\left({x}^{2}-u\right)}^{2}}}dx\right]\end{array}$$
(38)
$$=\frac{\hslash {q}_{x}}{16\pi }\frac{\alpha m}{\hslash }\left(2-{\int _{-1}^{1}}\frac{{x}^{2}-u}{\sqrt{{d}^{2}+{\left({x}^{2}-u\right)}^{2}}}dx\right)$$
(39)
$$=\frac{m\alpha {q}_{x}}{4\pi }g(d,u),$$
(40)

with \(u=\mu /{E}_{{\rm{R}}}+1\) and \(d=\Delta /{E}_{{\rm{R}}}\).

Derivation of spin current

Similar to the spin density derivation, we use perturbation theory to obtain the correction to the free energy \(\delta {F}_{J}\). The difference here is that we need to go to the third order perturbation due to the fact that the spin current is quadratic in \({q}_{x}\). After lengthy but straightforward calculation, the free energy correction is

$$\delta {F}_{J}={\int _{0}^{\infty }}dk[{F}_{-}(\alpha)-{F}_{-}(-\alpha)],$$
(41)

The general form of the function \({F}_{-}\) is complicated. However, in the limit \(\Delta \to 0\), it becomes

$$\begin{array}{l}{F}_{-}(\alpha) \!\!\!\!\approx\; \displaystyle\frac{{L}^{2}{q}_{x}^{2}{W}_{x}^{y}}{16{\pi }^{2}}\Bigg\{\frac{5\pi {\hslash }^{2}}{8m}{\rm{sgn}}(\xi -k\alpha \hslash)\hfill\\ -\displaystyle\frac{\pi k{\Delta }^{2}{\hslash }^{3}(4k\hslash -\alpha m)}{8{m}^{2}{\left[{\Delta }^{2}+{(\xi -\alpha k\hslash)}^{2}\right]}^{3/2}}\Bigg\}.\end{array}$$
(42)

By change of integral variable from \(k\) to \({k}_{-}\) to \(x\) similar to the previous calculation, the integral can be written as

$$\begin{array}{ll}\delta {F}_{J}=&-\frac{{L}^{2}k{q}_{x}^{2}{W}_{x}^{y}}{16{\pi }^{2}}\Bigg\{-\frac{5\pi {\hslash }^{2}}{4m}({k}_{2}-{k}_{1})\\ &+\pi \alpha \hslash {\int _{0}^{1}}\frac{2{d}^{2}\left({x}^{2}+3/4\right)}{{\left[{d}^{2}+{\left({x}^{2}-u\right)}^{2}\right]}^{3/2}}dx\\ &+\frac{7\pi \alpha \hslash }{4}\left(1+\frac{\mu }{\sqrt{{\Delta }^{2}+{\mu }^{2}}}\right)\Bigg\}.\end{array}$$
(43)

The spin current density is

$$\begin{array}{ll}{j}_{x}^{y}=-\frac{\partial \delta {F}_{J}}{\partial {W}_{x}^{y}}=&\frac{{q}_{x}^{2}}{16{\pi }^{2}}\Bigg\{-\frac{5\pi {\hslash }^{2}}{4m}({k}_{2}-{k}_{1})\\ &+\pi \alpha \hslash {\int _{0}^{1}}\frac{2{d}^{2}\left({x}^{2}+3/4\right)}{{\left[{d}^{2}+{\left({x}^{2}-u\right)}^{2}\right]}^{3/2}}dx\\ &+\frac{7\pi \alpha \hslash }{4}\left(1+\frac{\mu }{\sqrt{{\Delta }^{2}+{\mu }^{2}}}\right)\Bigg\}.\end{array}$$
(44)

We have defined two wave vectors \({k}_{2}\) and \({k}_{1}\) which denotes the Fermi wave number of the outer and inner Fermi surfaces respectively.

$${k}_{2}-{k}_{1}=\left\{\begin{array}{cc}2\alpha m/\hslash ,&{\rm{if}}\,\mu \, > \, 0,\\ \frac{2\sqrt{m\left(2\mu +{\alpha }^{2}m\right)}}{\hslash },&{\rm{if}}-\!\!{E}_{{\rm{R}}} < \, \mu \, < \, 0.\end{array}\right.$$
(45)

When \(\mu\) is large, the integral of the second term is negligible and the first and third terms lead to Eq. (6) of the main text. When \(-{E}_{{\rm{R}}} \, < \,\mu \, < \, 0\), the third term becomes negligible and Eqs. (7)–(8) are obtained.