Flux flow spin Hall effect in type-II superconductors with spin-splitting field

We predict the very large spin Hall effect in type-II superconductors whose mechanism is drastically different from the previously known ones. We find that in the flux-flow regime the spin is transported by the spin-polarized Abrikosov vortices moving under the action of the Lorenz force in the direction perpendicular to the applied electric current. Due to the large vortex velocities the spin Hall angle can be of the order of unity in realistic systems based on the high-field superconductors, superconductor/ferromagnet hybrid structures or the recently developed superconductor/ferromagnetic insulator proximity structures. We propose the realization of high-frequency pure spin current generator based on the periodic structure of moving vortex lattices. We find the patterns of charge imbalance and spin accumulation generated by moving vortices, which can be used for the electrical detection of individual vortex motion. The new mechanism of inverse flux-flow spin Hall effect is found based on the driving force acting on the vortices in the presence of injected spin current which results in the generation of transverse voltage.

expansion as described in main text. Close to the upper critical field, H c2 , we expand non-stationary GFĝ R = g 1 g 2 −g + 2 −g 1 in orders of |∆| and obtain linear equations where Π ± = ∇ ± 2ieA A A. The diagonal part of the GF can be found from the normalization condition g 2 1 = 1 + g 2 g + 2 . We apply gauge A A A = y y yH c2 x − x x xEt and use Abrikosov vortex lattice solution ∆ = b 0 e −2ieEtx ∑ n C n e inp(y−v L t) L (x − nx 0 ), see notations in main text, to get where q = DeH c2 . This solution can be rewritten for separate spin up or down band, see Eq.
(3) in the main text, and v L = E = 0 should be substituted to obtain stationary limit,ĝ R 0 .

Derivation of Eqs. (4), (5) and (6)
We assume that vortices move with the constant velocity v v v L and consider corrections in the linear-response regime which are realized provided the vortex velocity v v v L is small enough. For this purpose we take into account first-order terms in the gradient expansion of time convolutions as well as the non-equilibrium corrections to the spectral functionsĝ R/A ne and the distribution functionf ne =f − f 0τ0 . As a result, we get the two parts of spin current j j j s = j j j s1 + j j j s2 given by The first term in the r.h.s. of (S6) incorporates corrections to the spectral GF g R/A ne as well as the electric field term which appears from the expansion of the covariant differential operator∂ rX =∇X + eE E E{τ 3 , ∂ εX }/2 and∇X = ∇X − ieA A A[τ 3 ,X]. Here we use the gauge with zero electric potential such that electric field is given by E E E = −∂ t A A A. The second term in the r.h.s. of (S6) comes from the linear-order expansion of time convolution. It contains the equilibrium spectral GF in the moving frameĝ The second part of the spin current j j j s2 is determined by the non-equilibrium distribution function. In the general case the differential operator in the r.h.s. of (S7) contains correction from time convolution expansion so that∂ rf = ∇f + eE E Eτ 3 ∂ ε f 0 .
In order to calculate contribution to spin current due to spectral GF, j j j s1 , one should carefully take into account vortex-motion driven distortions of spectral GF, i.e. corrections to spectral GF,ĝ R ne ∼ v L , linear in vortex velocity. By using ansatzĝ R ≈ĝ R 0 +ĝ R ne and normalization condition we get This relation is used for calculating the first contribution to Eq. (S6), namely where short notation∇ab = (∇a)b + a(∇b) is used. Therefore, the spectral current density defined in (5) reads asĴ J J can be written as where ε ± = ε ± h + iq, I I I ∆ = ∆(Π∂ t ∆) * − ∂ t ∆(Π∆) * + cc and Π = ∇ − 2ieH c2 xy y y. We transform integral of (S10) over energy into the sum over Matsubara frequencies ω n = πT (2n where This results in the Eq. (4). To find the spatial average of (S11) we first notice that due to Abrikosov vortex-lattice solution the components of vector I I I ∆ reads as where x n = x −nx 0 . Stationary vortex-lattice solution is periodic function in y-direction with period L y = 2π/p and in x-direction with period L x = νx 0 , if C n = C n+ν . Here ν = 1 produces square and ν = 2 triangular vortex lattice. Spatial average over vortex lattice reduces to Gaussian integrals so that Combining this result with (S11) we obtain At the same time the average spin density deviation from the normal state induced by the superconducting correlations reads as Low-temperature limits of Eqs. (S13) and (S14) result in the relation (5) and . (S15) The absolute magnitude of the spin current depends on the order parameter amplitude b 0 which is determined by the magnetic field. In the limit of large Ginzburg-Landau parameter we get the usual expression for the average gap function 1, 2 with the only substitution Ψ (k) → ReΨ (k) to get into account the presence of spin splitting. In the limit of low temperatures T T c the amplitude is given by the analytical expression where β L is the Abrikosov parameter. Substituting this expression into (S15) and then the obtained S into (5), we get the spin Hall angle expression (6) by taking into account that near the upper critical field v L = − j/(σ n H c2 ).

Spin conservation
The Keldysh part of Eq. (S1) reads as and its components in the mixed representation are given up to the terms linear in vortex velocity by and E E E = −∂ t A A A is electric field in temporal gauge. We substitute this into Eq. (S17), multiply it by −N 0τ0σ3 /16, take trace and integrate over energy to obtain Since at high energies GF approach the normal-metal ones,ĝ R/A 0 = ±τ 3 , the integrals in the r.h.s. of Eq. (S19) vanish so that in the absence of the scattering processes spin is conserved.

Derivation and solution of Eq. (7)
Non-diagonal part of Usadel Eq. (S1) provides us with kinetic equation for the distribution functionf which parametrizes the Keldysh GFĝ K =ĝ R •f −f •ĝ A . By subtracting the spectral components of Eq.(S1) from Keldysh part one obtains the equation Here we omitted collision integral. Next we introduce four components of distribution function,f = f L + f T 3σ3 + ( f T + f L3σ3 )τ 3 , related to energy, charge and spin imbalances. We need only kinetic equation for f T 3 which can be obtained by using first-order gradient expansion of Eq. (S20) multiplied byσ 3 and traced where energy dependent diffusion coefficients and the spectral charge currents are Close to H c2 , kinetic equation (S21) reduces to (7). By using vortex-lattice solution, the terms in r.h.s of (7) reads as where I ± x = ∆ * ∂ x ∆ ± cc and I ± t = ∆ * ∂ t ∆ ± cc. By integrating (7) over energy and introducing spin accumulation µ s = ∞ −∞ f T 3 dε/2 we obtain close to H c2 ImΨ (1) .
To obtain r.h.s of first line in (S24) we integrated by parts and then transformed energy integral into Matsubara sum. By using vortex-lattice solution, the source term is determined by and I ± t = −iv L I ∓ x . Next we consider overlap between nearest Gaussians only, n − m = ±1. As a result, Eq. (S24) has the form and φ n = arg(C * n C n+1 ). The solution of differential Eq. (S26) has the form µ s = v L 2L 2

Derivation of Eq. (11)
Non-equilibrium corrections to spectral GF are determined by (S8) which in the leading order in |∆| can be written aŝ By using first-order gradient expansion of Keldysh GF and conditions (S28) we obtain in the leading order By performing integration over energy by parts and transforming integral to the summation over Matsubara frequencies, spin accumulation near H c2 becomes After summation this results in Eq. (11). By considering only overlap between nearest-neighbour Gaussians we obtain at low temperatures T → 0 asymptotic expressioñ (S31)

Derivation and solution of Eq. (8)
To calculate electrostatic potential defined by first line in (S29), we first solve kinetic equation for charge imbalance f T . The equation can be obtained by using first-order gradient expansion of Eq. (S20) multiplied byτ 3 and traced Close to H c2 , where order parameter is small, kinetic Eq reduces to Eq. (8) which can be written with the help of vortex-lattice solutions as We integrate r.h.s. of (S33) over energy by parts and transform integrals to Matsubara sum to get for µ = ∞ −∞ f T dε/2 an equation Particular solution of this inhomogeneous equation can be written in the approximation of nearest-neighbour Gaussians as G(x, y) = cos(py + χ n ) Note that last term in G produces linear divergence modulated by weak oscillation due to summation over n in µ. However, linear function is homogeneous solution of Eq. (S34), so that linear divergence in (S35) can be removed by substracting linear homogeneous solution from G.

4/5
Derivation of Eq. (12) By using (S29), full electrostatic potential near H c2 reads as After transforming energy integral to Matsubara sum we obtain Eq. (12). By using (S35) and asymptotic expansions of digamma function, we obtain at low temperatures