Spin diffusion equation in superconductors in the vicinity of T_{\rm c}

# Spin diffusion equation in superconductors in the vicinity of Tc

Takuya Taira Department of Physics, Okayama University, Okayama 700-8530, Japan    Masanori Ichioka Research Institute for Interdisciplinary Science, Okayama University, Okayama 700-8530, Japan Department of Physics, Okayama University, Okayama 700-8530, Japan    So Takei Department of Physics, Queens College of the City University of New York, Queens, New York 11367, USA    Hiroto Adachi Research Institute for Interdisciplinary Science, Okayama University, Okayama 700-8530, Japan Department of Physics, Okayama University, Okayama 700-8530, Japan
July 15, 2019
###### Abstract

We microscopically derive the spin diffusion equation in an -wave superconductor in the vicinity of the superconducting transition temperature . Applying the general relation between the relaxation function and the response function to the present spin diffusion problem, we examine how the spin relaxation time and the spin diffusion coefficient are renormalized in the superconducting state. The analysis reveals that, below , both the spin relaxation time and the spin diffusion coefficient are increased, resulting in an enhancement of the spin diffusion length. The present result may provide an explanation for the recent observation of an enhanced spin pumping signal below in a Py/Nb/Pt system that is free from the coherence peak effect.

## I Introduction

The interplay of the spin current and superconductivity has been attracting much attention, forming a new research field of superconducting spintronics Linder15 (); Eschrig15 (). Historically, the interplay of spin and superconductivity has been studied for a long time, and it is well known that the static perturbation having the odd time-reversal symmetry, such as spin, affects the thermodynamic properties of an -wave superconductor Anderson59 (); Abrikosov-Gorkov (); Maki-review (); Balatsky-review (). By contrast, the spin current is a flow of spin and thus has the even time-reversal symmetry Rashba06 (); Nagaosa08 (). Therefore, the interplay of spin current and superconductivity Bell08 (); Yang10 (); Hubler12 (); Quay13 (); Wakamura14 (); Ohnishi14 (); Yao18 (); Umeda18 () has an intriguing issue from the viewpoint of extending Anderson’s theorem Anderson59 () to dynamics, i.e., it requires to clarify the effects of dynamic and even time-reversal symmetry perturbation on the superconductivity.

When investigating such spin current physics, the most basic theoretical apparatus is the spin diffusion equation. Since after Bloembergen Bloembergen49 () discussed the diffusion of nuclear spins in a crystalline solid, this equation has been applied to a number of different physical situations Torrey56 (); Walker71 (); Aronov76 (); Silsbee79 (); vanSon87 (); Johnson88 (); Valet93 (); Hershfield97 (); Schmidt00 (); Rashba00 (); Yu02 (); Takahashi03 (). The advent of the spin transfer torque concept Slonczewski96 (); Berger96 () has further reminded us of the importance of the spin diffusion equation. Indeed the real space profile of the spin current is described by this equation, which can now be experimentally measured by the so-called lateral spin valve technique Jedema01 (). Therefore, in order to investigate the interplay of spin current and superconductivity, it is of great importance to know the spin diffusion equation in the superconducting state. Despite its importance, however, only a little is known about the spin diffusion equation in the superconducting state Yafet83 (); Yamashita02 ().

In this work, we microscopically derive the spin diffusion equation in the superconducting state on the basis of the weak-coupling BCS model with -wave pairing and with impurity spin-orbit scattering. Nevertheless, the derivation is not an easy task as one can infer from the fact that the superconductivity involves a nontrivial many-body ground state, and that we need to deal with the nonequilibrium problem in discussing the spin current. To overcome these difficulties, we limit ourselves to a temperature region near , and derive the spin diffusion equation by perturbation with respect to the superconducting gap .

Our microscopic derivation of the spin diffusion equation is enabled by the use of the general relation between the relaxation function and the response function, which in the present context is translated into the relation between the spin diffusion equation and the dynamic spin susceptibility. This relation is well known in the nonequilibrium statistical mechanics Kubo-text (); Chaikin-text (). Employing the formulation of Inoue17 (), where the dynamic spin susceptibility in the superconducting state was calculated, we derive the spin diffusion equation microscopically. Then, we show that the spin relaxation time as well as the spin diffusion coefficient are increased below , leading to an increment of the spin diffusion length in the superconducting state.

In order to compare the theory with experiments, we next apply the present result to a quite recent experiment measuring the spin transport through a superconducting Nb Jeon18 (). In Jeon18 (), the spin pumping in a Py/Nb/Pt system was studied, and it was found that the spin transport across the Nb layer from the Py to Pt layers is enhanced below . Then it is argued that the enhancement of the spin pumping in Nb below is due to a possible formation of triplet Cooper pairs. We analyze the corresponding physical situation using the result derived in the present paper. Then, we provide an alternative explanation for the enhancement in the spin transport below and close to without relying on the spin transport mediated by triplet Cooper pairs.

This paper is organized as follows. In the next section, we briefly review the derivation of the important relation between the spin diffusion equation and the dynamic spin susceptibility that constitutes the basis of the present work. In Sec. III, we introduce our model and the Green’s function, the latter of which is expanded in powers of the superconducting gap . In Sec. IV, using the relation between the spin diffusion equation and the spin susceptibility that is obtained in Sec. II, we derive the spin diffusion equation both in the normal state and in the superconducting state near . In Sec. V, relying on the result derived in the present paper, we calculate the spin pumping signal in a situation of Ref. Jeon18 () in the vicinity of . Finally, in Sec. VI we discuss and summarize our results. We use units throughout this paper.

## Ii Relation between spin diffusion equation and spin susceptibility

In this section, we briefly review the derivation of the key relation between the spin diffusion equation and the dynamic spin susceptibility. In Sec. IV, this relation will be used to microscopically derive the spin diffusion equation.

We begin with the following spin diffusion equation for the spin density :

 (∂∂t−Ds∇2+1τs)⟨s(r,t)⟩=0, (1)

where is the spin diffusion coefficient and is the spin relaxation time. Following Chaikin-text (), we consider a physical situation [Fig. 1(a)] where a spatially non-uniform static external magnetic field is applied to the system at time with which the spin density is in equilibrium. Then, the external magnetic field is turned off at , and the spin density diffuses for to achieve the spatially uniform equilibrium state. The external perturbation describing this situation can be written as

 hq(t)=Θ(−t)e0+thq, (2)

where and are respectively the step function and positive infinitesimal constant.

Consider first the time evolution of for . Using the momentum representation as well as introducing the Laplace-Fourier transformation , we obtain the solution to the spin diffusion equation (1):

 ⟨sq(ω)⟩=Rq(ω)⟨sq(t=0)⟩, (3)

where is the spin density at time , and

 Rq(ω)=1−iω+Dsq2+τ−1s (4)

is the relaxation function of the spin diffusion equation. Applying the result of static linear response to the initial spin density, i.e., , we obtain

 ⟨sq(ω)⟩=Rq(ω)χq(0)hq, (5)

where is the static spin susceptibility.

Consider next the dynamic linear response to the external magnetic field represented by Eq. (2). The basic idea is that the time dependence of the external magnetic field can be decomposed into the static field and the step-function field [Fig. 1(b)]. Since the static magnetic field defines the static spin susceptibility and the step-function field defines the dynamic (retarded) spin susceptibility , we can relate in Eq. (5) with and . Indeed, for the external magnetic field of Eq. (2), we have the relation

 ⟨sq(t)⟩=∫0−∞dt′χRq(t−t′)e0+t′hq. (6)

Then, going into the frequency space and using the spectral representation

 χRq(ω)=∫∞−∞dω′πImχRq(ω′)ω′−ω−i0+, (7)

Eq. (6) can be rewritten as

 ⟨sq(ω)⟩=∫∞−∞dω′iπImχRq(ω′)hq(ω′−ω)(ω′−i0+). (8)

After using of the partial fraction decomposition

 1(ω′−ω)(ω′−i0+)=1ω(1ω′−ω−1ω′−i0+), (9)

we arrive at

 ⟨sq(ω)⟩=1iω(χRq(ω)−χq(0))hq, (10)

where we again used the spectral representation (7).

Now we compare Eqs. (5) and (10), and obtain the key equation:

 Rq(ω)=1iω(χRq(ω)χq(0)−1), (11)

which relates the spin diffusion equation with the dynamic spin susceptibility. Moreover, we use Eq. (4) and arrive at the following identity:

 −iω+Dsq2+1τs=iω(χRq(ω)χq(0)−1), (12)

where the small , small limit should be taken. Eq. (12) means that spin diffusion equation can be derived microscopically by calculating the dynamic spin susceptibility .

## Iii Model & Green’s function near Tc

In this section, we first introduce our model Hamiltonian and define the corresponding Green’s function which is expanded in powers of the superconducting gap . Then, following Inoue17 (), we develop a formalism to calculate the dynamic spin susceptibility up to , by making use of the Green’s function technique. As expected, the present result is in line with that of Ref. Inoue17 () in the limit of small . Indeed, as we will see in Appendix A, our calculation can correctly reproduce the appearance of the coherence peak predicted for the spin pumping into superconductors Inoue17 ().

A key in our approach is that the vertex corrections are fully taken into account in calculating the dynamic spin susceptibility. It is well known that for a microscopic calculation of transport quantities, we need to keep a certain relation between the single-particle and two-particle Green’s functions in order to satisfy the conservation laws Baym-Kadanoff61 (). In the present formulation, this corresponds to the relation between the selfenergy function and the vertex function, which is usually referred to as Ward identity Nozieres-text (). This relation should be kept in the diffusion problem as well Maleev75 (), where it is known that only after taking account of the vertex corrections can we derive the correct diffusion equation that satisfies the charge conservation.

As a starting model, we consider the following BCS Hamiltonian for an impure -wave superconductor Inoue17 ():

 H = ∑pc†pξpcp−∑p(Δ∗cp↓c−p↑+Δc†−p↓c†p↑)+Himp,

where is the electron creation operator for spin projection and , is the kinetic energy measured from the chemical potential, and is the superconducting gap. The last term in Eq. (LABEL:eq:BCS01),

 Himp = ∑p,p′c†p^Up,p′cp′, (14)

describes the scattering by impurities. Here,

 ^Up,p′=u0(p−p′)+iuso(p−p′)^σ⋅(p×p′) (15)

is the impurity potential, where is the Pauli matrices, and are the amplitudes of the momentum and spin-orbit scatterings, respectively. Here and hereafter, we use “hat” notation to represent an operator in the spin space, .

As in Inoue17 () we use the representation for the impurity-averaged Green’s function,

 ˇGp(iεn)=(Gp(iεn),Fp(iεn)i^σyF†p(iεn)i^σy,G†p(iεn)), (16)

where is a fermionic Matsubara frequency with being an integer, and we denote a matrix in the particle-hole (Nambu) space by “check” accent. Then, the Matsubara susceptibility is obtained from

 χq(iων)=−T2∑εn∫pTr[^σz^Υp,p−(iεn,iεn−)], (17)

where is a bosonic Matsubara frequency with being an integer, and we have introduced the shorthand notation . Here, is the component of the matrix function defined by

 ˇΥp,p−(iεn,iεn−)=ˇGp(iεn)ˇΛq(iεn,iεn−)ˇGp−(iεn−), (18)

where and . In the above equation, the vertex function arises from the effects of impurity ladder, and it satisfies the following equation:

 ˇΛq(iεn,iεn−) = ^σz+nimp∫p′ˇVp,p′ˇGp′(iεn) (19) ×ˇΛq(iεn,iεn−)ˇGp′−(iεn−)ˇVp′,p,

where the effect of impurities is described by

 ˇVp,p′=(^Up,p′,00,^Utp′,p), (20)

with being the transpose of a matrix in the spin space.

In order to solve Eq. (19), we introduce the representation Gorkov64 ():

 ˇΛq(iεn,iεn−)=(Λ1^σz,Λ2^σx−Λ†2^σx,Λ†1^σz). (21)

In the normal state where , since and commute as they are diagonal in the Nambu space, the relation and holds. In the superconducting state, if the static limit () is concerned, since and coincide and thus commute, we have the important symmetry

 Λ1=Λ†1,Λ2=Λ†2. (22)

Because we are interested in the small region, we keep this symmetry and focus on a small correction to the spin diffusion equation in the normal state.

Using Eqs. (21) and (22), we can transforms Eq. (19) into a set of linear equations for and :

 (Λ1Λ2) = (10)+(A,−BC,D)(Λ1Λ2), (23)

where the coefficients , , , are calculated from Eq. (19) as

 A = g∫p{Gp(iεn)Gp−(iεn−)+Fp(iεn)F†p−(iε−)}, (24) B = g∫p{Gp(iεn)F†p−(iεn−)+Fp(iεn)Gp−(iεn−)}, (25) C = g∫p{Gp(iεn)Fp−(iεn−)−Fp(iεn)G†p−(iε−)}, (26) D = g∫p{Gp(iεn)G†p−(iεn−)−Fp(iεn)Fp−(iε−)}. (27)

Here, , and is the density of states per one spin projection. The scattering rate arising from the vertex corrections is given by

 Γ(−)=12(1τ0−13τso), (28)

where and are respectively the momentum relaxation time and spin-orbit relaxation time, which are assumed to satisfy and given by and . Here, is the number density of impurities, and means the average over the Fermi surface. Substituting the solution of Eq. (23) into Eq. (17), the Matsubara spin susceptibility is calculated to be

 χq(iων) = −πN(0)Γ(−)T∑εnA(1−D)−BC(1−A)(1−D)+BC. (29)

Since we are interested in the region immediately below , we expand the normal Green’s functions in powers of  Gorkov60 ():

 Gp(iεn) = G(0)p(iεn)+δGp(iεn)+O(Δ4), (30) G†p(iεn) = G†(0)p(iεn)+δG†p(iεn)+O(Δ4), (31)

where

 G(0)p(iεn)=1i˜εn−ξp,G†(0)p(iεn)=−1i˜εn+ξp, (32)

are the Green’s functions per one spin projection in the normal state, with

 η=1+Γ(+)|εn|, (33)

and the scattering rate is given by

 Γ(+)=12(1τ0+1τso). (34)

Now we perform a perturbation expansion with respect to . Following the standard procedure (e.g., Werthamer-review ()), we find that and are given by

 δGp(iεn) = (G(0)p(iεn))2{iΓ(+)sgn(εn)2|εn|2−η2G†(0)p(iεn)}Δ2, (35) δG†p(iεn) = (G†(0)p(iεn))2{η2G(0)p(iεn)−iΓ(+)sgn(εn)2|εn|2}Δ2.

In a similar way the anomalous Green’s function, , which is equal to in the Meissner state, is expanded as

 Fp(iεn)=δFp(iεn)+O(Δ3), (37)

where is given by

 δFp(iεn)=G†(0)p(iεn)ΔηG(0)p(iεn). (38)

## Iv Derivation of the Spin diffusion equation

In this section, we derive the spin diffusion equation microscopically. First, we focus on the normal state and demonstrate that our procedure correctly reproduces the well-known expression for the spin diffusion equation in the normal state. Next, we extend the derivation to the superconducting state near , and show that the spin diffusion length is increased below .

### A Normal state

Let us first calculate the Matsubara susceptibility in the normal state, by using Eq. (29). Noting that in the normal state, the Matsubara susceptibility is calculated to be

 χq(iων)=−πN(0)Γ(−)T∑εnA(0)1−A(0), (39)

where

 A(0)=g∫pG(0)p(iεn)G(0)p−q(iεn−iων) (40)

is the normal state value of . The summation over Matsubara frequency in Eq. (39) requires a care, because the analytic behavior of changes discontinuously depending on the sign of . Therefore, it is convenient to divide the susceptibility into two parts AGD ():

 χq(iων)=χ′q(iων)+χ′′q(iων), (41)

where comes from the “non-dissipative” region , i.e.,

 (42)

whereas comes from the “dissipative” region , i.e.,

 χ′′q(iων)=−πN(0)Γ(−)T∑0<εn<ωνA(0)1−A(0). (43)

On the complex plane where the Matsubara frequency is located at (Fig. 2), is given by the sum over regions and , whereas is given by the sum over region .

We next evaluate and in the normal state. Regarding , it is well known that this quantity is dominated by the static susceptibility

 χ′q(iων)≈N(0), (44)

where the result is insensitive to the effects of impurities Fulde68 (). Regarding , since in the “dissipative” region is given by with being the diffusion coefficient, we obtain

 χ′′q(iων) = −N(0)ων2Γ(+)A(0)1−A(0), (45)

where in moving from Eq. (43) to Eq. (45), the summation over is replaced by multiplication by a factor , i.e., the number of Matsubara frequencies in the region , as the summand is independent of  Fukuyama85 ().

As the final step, we perform the analytic continuation . Then, is calculated to be

 A(0)=Γ(−)Γ(+)(1−−iω+Dq22Γ(+)), (46)

from which we obtain the retarded spin susceptibility,

 χRq(ω)=N(0)Dq2+(4/3τso)−iω+Dq2+(4/3τso), (47)

due to Eq. (45). Substituting the above equation into Eq. (12), we obtain the spin diffusion coefficient

 D(n)s=v2Fτtr3, (48)

and the spin relaxation time

 τ(n)s=34τso, (49)

where the superscript means the value in the normal state.

### B Superconducting state near Tc

We now calculate the spin diffusion equation in the vicinity of the superconducting transition, and show that the spin diffusion length is increased below . Since we are interested in an analytic result valid near , we expand the dynamic spin susceptibility up to .

First, we consider the “non-dissipative” component, , which is dominated by the uniform susceptibility derived in Abrikosov62 (). The result can be expanded up to (see Appendix B),

 χ′q(iων) = N(0)⎧⎨⎩1−2πT∑εn>0Δ2ε3n+23τso⎫⎬⎭+O(Δ4) = N(0){1−(7ζ(3)4−5ζ(4)4ε0τso)Δ2ε20}+O(Δ4),

where is the lowest fermionic Matsubara frequency,

 ε0=πT. (51)

Next, we consider the “dissipative” region . To calculate up to , we expand the coefficients in powers of :

 A = A(0)+δA+O(Δ4), (52)

where is defined by Eq. (40). In the above equation, is given by

 δA = g∫p{δG(p)G(0)(p−)+G(0)(p)δG(p−) (53) +δF(p)δF†(p−)+δF(p)δF†(p−)},

where we have introduced the four-vector notation , namely, , , etc. Recalling that the normal state values of and vanish, we expand and as follows:

 B = δB+O(Δ3), (54) C = δC+O(Δ3), (55)

where and are both of the first order with respect to , which are explicitly given by

 δB = g∫p{G(0)(p)δF(p−)+δF(p)G(0)(p−)}, (56) δC = g∫p{G(0)(p)δF(p−)−δF(p)G†(0)(p−)}, (57)

where is defined below Eq. (27).

Finally, the coefficient in the “dissipative” region is expanded as

 D = δD+O(Δ4), (58)

where, as we will see below, only the fact that is important and the detailed expression of is unnecessary.

With these apparatus, let us calculate . Substituting Eqs. (52), (54), (55), and (58) into Eq.(29), and expand the numerator and the denominator up to , we obtain

 χ′′q(iων)=−πN(0)Γ(−)T∑0<εn<ωνS(εn) (59)

where the summand is defined by

 S(εn)=A(0)+[δA−δBδC]1−A(0)−[δA−δBδC]. (60)

To proceed further, we first recall that the summand is independent of in the normal state, and that we are considering a small correction to the normal state, of the order of . Under this circumstances, the summand depends very weakly on . Therefore we regard the summand as nearly independent of in the “dissipative” region , and approximate

 T∑0<εn<ωνS(εn)≈Tων2πTS(ε0), (61)

where is defined by Eq. (51). Note that this approximation is justified in the case with sizable impurity spin-orbit scattering where the unperturbed term of the denominator in Eq. (60), , acquires a nonvanishing value even in the limit . Then, we obtain the expression

 χ′′q(iων)=−N(0)ων2Γ(−)A(0)+[δA−δBδC]1−A(0)−[δA−δBδC], (62)

where the quantities , , and are now evaluated at . Keep this in mind and using approximation consistent with Eq. (22), after a bit tedious algebra, these quantities are calculated to be

 δA = −Γ(−)Γ(+)Δ2ε20γ+(2−3γ)Dq22Γ(+)(1−γ)2, (63) δB = δC=−iΓ(−)Γ(+)Δε01−Dq22Γ(+)1−γ, (64)

where . Note that should be in the region to obtain physical results (see Appendix A).

We are now in a position to calculate the retarded susceptibility. In the following, we normalize the susceptibility by the density of states in the normal state as

 χRq(ω)=N(0)˜χRq(ω). (65)

Starting from Eqs. (LABEL:eq:Rechi01) and (62), then performing analytic continuation , we obtain

 ˜χRq(ω) = 1+δ˜χ0+iω2Γ(−)A(0)+δE1−A(0)−δE, (66)

where

 δ˜χ0=−(7ζ(3)4−5ζ(4)4ε0τso)Δ2ε20. (67)

In Eq. (66), is given by Eq. (46), and is defined by

 δE=δA−δBδC. (68)

From Eqs. (63) and (64), is calculated to be

 δE = δE0+δE1q2, (69) δE0 = Δ2ε2011−γ, (70) δE1 = −Δ2ε20D2Γ(+)4−3γ(1−γ)2, (71)

where we have used in this small correction term. Using the above equation as well as Eq. (46), we can represent the retarded susceptibility in the following form:

 ˜χRq(ω) = (δE0−δ˜χ0)iω+(1+δ˜χ0)(Dsq2+1/τs)−iω+Dsq2+1/τs, (72)

where we have introduced and . Using Eqs. (70) and (71), and can be expressed as

 Ds ≈ D(n)s(1+Δ2ε204−3γ(1−γ)2), (73) τs ≈ τ(n)s⎛⎝1+Δ2ε202Γ(+)τ(n)s1−γ⎞⎠. (74)

Now we substitute Eq. (72) into Eq. (12) to derive the spin diffusion equation. Then, we obtain the following equation:

 1+δ˜χ01+δE0(−iω+Dsq2+1τs)=0, (75)

from which we find the following spin diffusion equation:

 (∂∂t−Ds∇2+1τs)⟨s(r,t)⟩=0, (76)

where and are defined by Eqs. (73) and (74).

The result thus obtained means that, in the superconducting state near , [Eq. (73)] is identified as the spin diffusion coefficient, and [Eq. (74)] is identified as the spin relaxation time. Note that both quantities are increased upon the superconducting transition. Finally, the spin diffusion length is calculated to be

 λs = √Dsτs (77) ≈ λ(n)s⎛⎝1+Δ2ε20⎡⎣Γ(+)τ(n)s1−γ+4−3γ(1−γ)2⎤⎦⎞⎠,

where . Therefore, the spin diffusion length is also increased below .

## V Application to experiments

In this section, we perform a model calculation to compare the present result with experiments. We apply the theory to a quite recent experiment observing an enhanced spin pumping into the superconducting Nb in a Py/Nb/Pt structure Jeon18 (). In Jeon18 (), the observed enhancement of the spin pumping signal is attributed to a putative spin supercurrent carried by triplet Cooper pairs. Below, relying on the spin diffusion carried by quasiparticles, we analyze the corresponding physical situation and put forward an alternative interpretation of the experiment near without invoking the triplet spin supercurrent.

As a simple model calculation, we consider an F/S/X system as shown in Fig. 3, where F and S are a ferromagnet and an -wave superconductor, respectively, and X is either a vacuum or a perfect spin sink. Then, we assume that the distribution of the spin accumulation in is described by the spin diffusion equation [Eq. (76)]. Then, the spin accumulation is written as

 δs(z)=Asinh(zλs)+Bcosh(zλs), (78)

where , and are unknown coefficients to be determined by the boundary conditions. The spin current inside , , is given by

 Js(z)=−Dsλs{Acosh(zλs)+Bsinh(zλs)}. (79)

Now, following ZhangZhang12 (), we consider the boundary conditions. At the F/S interface, we impose

 Jtots=Jpumps−