Extrinsic Spin-Charge Coupling in Diffusive Superconducting Systems

We present a theoretical study of diffusive superconducting systems with extrinsic spin-orbit coupling and arbitrarily strong impurity potential. We derive from a microscopic Hamiltonian a diffusion equation for the quasi-classical Green function, and demonstrate that all mechanisms related to the spin-orbit coupling are expressed in terms of three kinetic coefficients: the spin Hall angle, the spin current swapping coefficient, and the spin relaxation rate due to Elliott-Yafet mechanism. The derived diffusion equation contains a hitherto unknown term describing a spin-orbit torque that appears exclusively in the superconducting state. As an example, we provide a qualitative description of a magnetic vortex in a superconductor with triplet correlations, and show that the novel term describes a spin torque proportional to the vector product between the spectral angular momentum of the condensate and the triplet vector. Our equation opens up the possibility to explore spintronic effects in superconductors with no counterparts in the normal metallic state.


I. INTRODUCTION
Exotic phenomena can occur when two or more materials with different properties are merged into a single hybrid-structure. Of particular interest are hybridstructures consisting of conventional superconductor (S) and ferromagnet (F) in which the interplay between these two quantum states leads to striking effects due to the appearance of odd-triplet superconducting correlations. 1,2 .
The quasiclassical kinetic equation provides a unified description of realistic hybrid-structures, including disorder, interfaces and out-of-equilibrium situations. This approach becomes particular relevant when the hybridstructure is in the diffusive limit. The resulting diffusion equation describing superconducting structures is known as the Usadel equation 3 . Its extraordinary descriptive power has been demonstrated by the agreements between theory and experiments on S-N (N is a normal metal) 4-6 and S-F structures [7][8][9] . As such, the Usadel equation has been extended to discuss spin-relaxation induced by spin-orbit coupling (SOC) and/or magnetic exchange field. 10 . In fact, since the discovery of the Knight shift (paramagnetic response of superconductor), spin-relaxation induced by SOC has been extensively investigated in various superconducting systems [11][12][13][14][15][16] . Besides spin-relaxation 17 , SOC can also lead to fascinating spin-charge coupling phenomena that has spurred intensive research activities in the field of spintronics. Notable spin-charge coupling phenomena like the anomalous Hall effect 18 , the spin Hall effect 19 (SHE) and the spin-galvanic effect [20][21][22][23] have been discussed extensively in metallic systems.
Since the discovery of anomalous Hall effect, the microscopic origin of spin-charge coupling has been receiving everlasting attentions, see Ref. 18 for a historical overview. It is customary to separate the microscopic origin of SOC into the intrinsic and extrinsic type. While intrinsic SOC originates from inversion breaking potential that respects translational symmetry of the underlying material, the extrinsic SOC is generated from disorder potential that breaks translational symmetry such as random impurities.
In a superconducting state, the effect of spin-charge coupling stemming from intrinsic SOC has been investigated extensively. For example, Ref. 24 and 25 discussed the Edelstein effect in non-centrosymmetric superconductors while Ref. 15,26, and 27 discussed the consequences of uniform SOC in superconducting junctions. Two of the present authors also study the effect of intrinsic SOC in S/F hybrid structures [28][29][30] and found that spin-charge coupling is related to the coupling between the singlet and triplet components of the superconducting condensate.
Contrary to intrinsic SOC, spin-charge conversion induced by extrinsic SOC received considerably lesser attention despite its importance in understanding superconducting junctions where disorder is ubiquitous. Recently, some quasiclassical kinetic theories have been put forward to discuss the effect of spin-charge conversion generated by extrinsic SOC 31,32 . In Ref. 31, a quasiclassical kinetic equation is developed using the Born approximation to treat the disorder induced self-energy. Since the Born approximation does not capture skewscattering induced by SOC impurities 18 , the derived equation in Ref. 31 is parameterized by spin Hall an-gle θ that only contains contributions from the so-called side-jump mechanism. In Ref. 32, the skew-scattering is taken into account by evaluating the self-energy at the third Born approximation but the side-jump mechanism is not fully accounted for and the resulting Usadel equation missed a spin torque generated by SOC (i.e. Eq. (3)) which, as we show below, is present in the superconducting state.
In this work, we present a comprehensive and rigorous study of spin-charge conversion in superconductors with extrinsic SOC in the diffusive limit beyond the Born approximation. Our main result generalized the Usadel equation to account for the spin Hall and the spin current swapping effects. Strikingly, we uncover a non-linear and non-local torque induced by SOC that vanishes in the normal metallic state. We demonstrate that all the terms related to spin-charge coupling in the Usadel equation can be parametrized by the spin-relaxation time and two kinetic coefficients -the normal state spin Hall angle θ and the spin swapping 33 coefficient κ. We assume the extrinsic spin-orbit coupling to be weak but treat the impurity scalar scattering exactly without resorting to any finite order Born approximation. In other words, the disorder potential can be arbitrarily large with the (s-wave) scattering length taking any value from minus to plus infinity. On the one hand, this provides a unified picture of extrinsic spin-charge coupling where the comparison between side-jump and skew-scattering is unambiguous. On the other hand, our approach captures resonant scattering induced by SOC disorder [34][35][36] , which may lead to an enhancment of spin-charge conversion, as in the case of graphene decorated with adatoms 37,38 , or materials with Kondo impurities 39,40 .
The article is organized as follow: In Sec. II, we present and discuss our main result -a generalization of the Usadel equation which accounts for extrinsic spin-charge coupling. In Sec. III, we introduce our microscopic model, basis for the Green function, and the general formulation of kinetic theory. In Sec. IV, we describe the evaluation of the collision integral in the presence of SOC and generalize the Usadel equation. Emphasis is placed on careful analysis of spatially non-local selfenergies and the power-counting scheme to correctly capture the side-jump effect and the spin-orbit torque. In Sec. V, we illustrate the effect of spin-orbit torque in the presence of a vortex. We close the article with a outlook in Sec. VI. Technical details concerning the evaluation of self-energy and collision integral are given in the Appendix. Throughout this work, we set = 1 and adopt the rule of summing over repeated indices.

II. MAIN RESULTS
We aim to derive a set of quasiclassical kinetic equations for diffusive superconducting systems with extrinsic spin-orbit coupling. These equations are valid in the regime ξ l p −1 F λ where ξ, l, p −1 F , and λ cor-respond to the superconducting coherence length, meanfree path, Fermi wavelength and the effective Compton wavelength of a material, respectively. The quasiclassical approximation focus on the spatial variation of observables and spectral functions over distances much larger than the Fermi wavelength p −1 F . The diffusive limit further sets ξ l. The SOC, being a relativistic phenomena, is microscopically associated with a material dependent effective Compton wavelength λ. The effective Compton wavelength can be much larger than its vacuum value 41 , but it is still significantly smaller than the Fermi momentum in most materials 42 , i.e. λ 2 p 2 F 1.Note the extrinsic SOC that leads to Mott scattering starts at power λ 2 p 2 F , this is different from an uniform SOC in electron gas where the expansion can occur at αp F where α is the strength of uniform SOC. This justifies a perturbation theory in λ 2 p 2 F while allowing the impurity potential strength to be arbitrarily large. In other words, we sum the entire Born series generated by the impurity potential (see Fig. 2c). This introduces the scattering length a which, within our approach, can take any value −∞ < a < ∞.
Our main result is the generalized Usadel equation describing the effect of extrinsic SOC: Hereǧ =ǧ(r, t 1 , t 2 ) quasiclassical Green function (GF) averaged over momentum at the Fermi surface. It is an 8 × 8 matrix in the Nambu-spin-Keldysh space. The Pauli matrices σ a and τ i span, respectively, the spin and Nambu spaces,∆ is the anomalous superconducting self energy and τ so is the usual spin relaxation time induced by SOC disorder. Nontrivial effects related to the SOC enter the Usadel equation via the generalized 8 × 8 matrix currentJ k and the matrix torqueŤ . The matrix currentJ k flowing in the spatial direction k is defined aš where kja is the total antisymmetric tensor. The first term in the right hand side of Eq. (2) is the standard diffusive current characterized by the diffusion constant D = v F l/3. The second and third terms describe the spin Hall effect (the anticommutator term couples the charge and spin degrees of freedom) and the spin current swapping (the commutator couples different components of the spin flow), respectively. Note that our theory captures the slow variation of all kinetic coefficients, D, θ, κ and τ so on scales larger than p −1 F , as sketched in Fig.  1. The conservation of the generalized currentJ k at interfaces between different materials define a boundary condition for Eq. (1).
Interestingly, we uncover a non-local spin-orbit torquě T in Eq. (1) when we account for the anomalous velocity impurities r Kinetic coe cients FIG. 1. The formulated kinetic theory (Eq. 1) describes the diffusive transport of spin, charge and spectral weight in a superconducting hybrid structures. Importantly, our theory captures spatially varying kinetic coefficients (e.g. spin-Hall angle θ and swap-current coefficient κ) which are important for spin-charge conversion.
induced by SOC disorder 43 . It is given by the following: WhileJ k , after taking corresponding traces, describes charge and spin currents in the normal metallic state,Ť is only finite in the superconducting state where the anomalous components of the GF are non-vanishing. Note that the trace ofŤ over the spin Pauli matrices σ are always zero. It gives a finite contribution if one first multiplies Eq. (1) by τ 3 σ a and then takes the trace, i..e.Ť describes a novel type of spin torque in superconductors. Both the torqueŤ and the generalized currentJ k are parameterized by the spin Hall angle θ and the swapping coefficient κ derived in section IV. They are given by following expressions: where τ is the elastic scattering time. The effective spincharge coupling rates ω 1 and ω 2 can be expressed in terms of components of the single impurity scattering matrix at the Fermi energy: where N F is the density of states at the Fermi energy, and n im is the impurity concentration. To the lowest order in SOC, the coefficient B is real, so that ω 2 is related to ω 1 via the optical theorem ω 2 = p F a ω 1 , where a is the scattering length.
Importantly, the kinetic coefficients θ and κ are exactly those characterizing the coupled spin-charge transport in the normal state. The first and the second terms in Eq. (4) are the renowned skew scattering and the side-jump contributions to the spin Hall angle respectively. In Eq. (5), the first term was identified by Lifshitz and Dyakonov 33 as the swap current coefficient, whereas the second term arises when we consistently include the anomalous velocity induced by SOC. The latter modifies the first term just as side-jump modifies the skew scattering in θ.
As it will become clear later, we shall name the first (second) term in Eq. (5) as the local (nonlocal) swapcurrent coefficient. Similar to the side-jump contribution to the spin Hall conductivity, we found that the "swap-current conductivity" would also have a component that scales independently from the impurity concentration. Note that in the limit of strong scattering potential a → ∞, the skew-scattering dominates over the side-jump mechanism in θ while the nonlocal swap current dominates the local mechanism in κ.

III. MODEL HAMILTONIAN, BASIS AND KINETIC FORMULATION
In this section, we discuss the model Hamiltonian, the basis we use to define the matrix Green functions, and the basic kinetic theory of Green function. The starting mean field Hamiltonian for a superconducting system is The mean-field superconducting order parameter ∆ is local in space, σ 2 is the second Pauli matrix and K αβ (r, −i∂ r ) is the single-particle Hamiltonian given by We consider here a disorder potential V αβ (r, i∂ r ) which contains a spin-independent part, proportional to δ αβ , and a spin-orbit coupling part porportional to σ αβ : Here λ is the material dependent Compton wavelength and we assume λ 2 p 2 F 1. V (r) is a short-range potential induced by N i randomly distributed impurities.
The field operators entering Eq. (6) can be conveniently organized as a spinor: and correspondingly It is customary, see for example Ref. 10 or the chapter by K. Maki in Ref. 44, to define the matrix Green functions as the time ordered correlatorǦ = −i T ΨΨ † of the product between the column bispinor Ψ = (ψ ↑ , ψ ↓ , ψ † ↑ , ψ † ↓ ) and the row bispinor Ψ † = (ψ † ↑ , ψ † ↓ , ψ ↑ , ψ ↓ ). The Gorkov equations for this GF can then be obtained straightforwardly by using the Heisenberg equation of motion of the field operators 45 .
Although the GF defined in the above basis are widely used in the literature, we opt here for another basis which allows for a more intuitive interpretation of different terms in the equation of motion for the GF and subsequently in the kinetic equation. Instead of using the spinors defined in Eq. (10) and (11) we introduce the time-reversal conjugated spinor and construct the GF asǦ = −i TΨΨ † , whereΨ = (Ψ, Ψ c ) andΨ † = (Ψ † , Ψ †c ). A more intuitive form of the equations of motion for such defined GF is related to the fact that it explicitly reflects the superconducting paring between the time-reversal conjugated states. In our basis, the quasiclassical matrix GF that enters Eq.
(1) has the form:ǧ where the hat. quantities are matrices in spin space: One of the advantages of using this basis is that the Pauli matrices only appear multiplying spin-related quantities, in particular, the triplet components of the condensate amplitude f a where a = x, y, z. Furthermore, since the kinetic energy and the impurity potential is time-reversal invariant [i.e. K(r, −i∂ r ) = σ y K(r, i∂ r )σ y ], it is simply proportional to identity in Nambu space τ 0 , when written in our basis.
In contrast, if one uses the GF defined by the basis in Eq. (10)- (11), all components of the anomalous Green's function acquire an additional iσ y factor. Moreover, all spin-dependent fields has to be written using a Nambu diagonal matrix proportional to diag[σ, σ * ] 10 . In Table I, we compare different physical quantities expressed in the basis used in Ref. 10 and our basis. For readers who wish to recover the GF defined in Ref. 10, they can do so by applying the following transformation to our matrix GF defined in Eq.(13):Ǔ †ǧǓ with U = 1 2 (1+iσ y )(1−iτ 3 σ y ) = e i π 4 σ y (1−τ3) . Correspondingly one can use this transformation to transform our Usadel equation, Eq. (1), to the basis used in Refs. 10, 46, and 47.
Having established the basis in which the Green functions are written, we now derive the kinetic equation governing the charge-spin coupling in superconducting systems. The derivation of the quasiclassical kinetic equation from microscopic Hamiltonian can be found in many textbooks 48,49 and reviews 1,50,51 . Here, we provide a brief summary of it and postpone the calculation of the collision integral within the quasiclassical approach to the next section. Given a microscopic Hamiltonian [Eq.
where for abreviation the numbers j = 1, 2, 3 denote the set of spatial and time coordinates r j , t j . HereΣ is the self-energy. In the basis we have chosen to represent the Green's functions the matrix describing the superconducting order parameter readš In order to derive the kinetic equation from Eq. (16), one introduce the Wigner coordinates r = (r 1 + r 2 )/2, s = r 1 − r 2 . Unlike the derivation of the Boltzmann equation 52 , there is no obvious advantage of introducing the Wigner coordinates for the time component. Therefore, we Fourier transform Eq. (16) only with respect to the relative space coordinate s and arrive at the following equation: Note that the (two) time arguments are skipped for brevity and we neglect the small spatial dependence of the gap function. More importantly, we retain the Poission bracket in the right hand side that is usually discarded 53 . It turns out that this term is essential to describe effects associated to charge-spin coupling 31 . Fundamentally, this is because the self-energy describing spin-charge coupling is not only made up by GFs that are local in space. As we shall see in the next section, the presence of SOC generates self-energy terms that depends on nonlocal GFs. Hence, the commutator and the Poisson bracket on the right hand side of Eq. 18 can generate terms that are of the same order in the nonlocality of the GF. Eq. (18) has the form of a kinetic equation; it describes a balance between the driving force (left hand side) and collision integral (right hand side). Note that we have not given any prescription to perform the (quantum) average of the fermion field operators in the Green functions. At zero-temperature, the average is taken over the ground state of a filled Fermi sea. At finite-temperature and thermal equilibrium, the average can be taken over the grand canonical ensemble using the Matsubara formalism. For a system that is out-of-equilibrium, we shall place the time-coordinates onto the Keldysh time contour and promote all 4×4 matrices in Nambu-spin space onto an 8×8 matrix in the Nambu-spin-Keldysh space. Since in all these situations Eq. (18) remains formally unchanged, we use the check symbol,., to denote either the 4×4 matrices in the equilibrium case or 8×8 matrices in the Keldysh formalism.

IV. DERIVATION OF THE USADEL EQUATION
The anitcommutator in Eq. (18) describes important spin-charge coupling also poses a hurdle to continue the derivation of the kinetic equation following standard approach 1,48-51 . In this section, we discuss in detail how we deviate from the standard approach and derive the generalized Usadel equation in the presence of disorder SOC from Eq. (18).
Let us begin by reminding the readers that Eq. (18) still contains superfluous information that is unessential for the description of electronic transport near the Fermi energy. In a superconductor, where the density of states changes dramatically around the Fermi energy, it is customary to simplify Eq. (18) within the quasiclassical approximation. In this approximation, the Fermi-energy is the largest energy scale in the problem and, as mentioned in Sec. II, spatial variations of all observables and spectral functions take place over distances much larger than the inverse of the Fermi momentum. Moreover, the GFs are peaked at the Fermi level and therefore it is convenient to integrate them over the quasiparticle energy (ξ p ) to obtain the so-called quasiclassical Eilenberger Green function:ǧ where n is a unit vector pointing in the direction of the momentum at the Fermi surface. The standard way of deriving the quasiclassical kinetic equation is to integrate Eq. (18) over the quasiparticle energy and to obtain an equation forǧ, the Eilenberger equation 10,54 . This equation is complemented by a normalization conditionǧ 2 = 1. In the present case however, the situation is more complicated and one cannot follow this path straightforwardly. This is because the Poisson bracket, i.e. the anti-commutators on the right hand side of Eq. (18), contains momentum derivatives. They prohibit a straightforward integration over the quasiparticle energy and do not ensure the normalization condition for the GF at this stage.
In order to overcome these difficulties we follow the procedure put forward in Ref. 31 and assume that the system is in the diffusive regime. In this limit the system is almost isotropic in space. We then expandǧ in spherical harmonics and keep only the zeroth and first moments:ǧ (n, r) ≈ǧ(r) + n kǧk (r). (20) Our goal is to obtain a close equation for the zerothmoment Green functionǧ(r), i.e. the Usadel Green function. For this sake, we resort to the following counting scheme of small parameters in the diffusive limit. Let 0 , 1 and 2 describe the characteristic magnitudes of the zeroth moment, first moment and leading non-locality of the quasiclassical GF:ǧ In the diffusive limit, the zeroth moment of the GF is the dominant component and we have the following hierarchy of scales 0 1 2 . The first inequality 0 1 arises from the fact that we are considering the diffusive limit, 1 l/ξ, while the second inequality, 1 2 , arises from the quasiclassical limit l, ξ p −1 F . As mentioned in Sec. II, our theory describes the macroscopic in-homogeneity of the disorder potential and we shall assume that the kinetic coefficients changes on the scale of 2 .
Spin-charge coupling occurs at linear order in λ 2 p 2 F and has contributions from both skew-scattering and sidejump mechanism. The skew-scattering mechanism occurs at 1 ; it does not require spatial nonlocality and can be captured in standard T-matrix calculation with equilibrium/uniform Green function. Unlike skew-scattering mechanism, the side-jump mechanism requires the Green function to be non-uniform in space so it is of the order 1 2 . Since in our power counting scheme 1 2 , in order to catch consistently the side-jump contribution, we also have to retain terms of order 2 1 . These terms are typically discarded in the standard derivation of the Usadel equation without SOC.
At order λ 4 p 4 F , the most dominant contribution to the self-energy is the Elliott-Yafet spin relaxation which occurs at order 0 . Hence, at this order, we shall only retain the 0 term and neglect all other corrections arising from 1 and 2 .
Having established our approximation scheme, we can proceed to compute the equation of motion for the zeroth and first moment from Eq. (18). The resulting equation of motion forǧ(r) andǧ k (r) are given by: In order to lighten the notations, from now and what follows, we shall display the space arguments ofǧ andǧ k only when it is important for discussion. In Eq. (25) we assume that the elastic scattering rate is much larger than the superconducting gap and the typical rate of change ofǧ k i.e. ∂ t1 , ∂ t2 , ∆ τ −1 . The right hand side of Eq. (24) and (25) correspond to the collision integral of the zeroth and first moment GF respectively. They are the essence of extrinsic spin-orbit coupling: The GFs on the right hand side of Eq. (26) and (27) are the Eilenberger Green function,ǧ(n, r). Their arguments are shown explicit to distinguish them from the zeroth moment Green functionǧ, Eq. (20). The angular brackets in Eqs. (26)- (27) stand for integration over the solid angle at the Fermi surface. In deriving Eq. (26), we assume that the Green function is zero at large momentum and used integration by parts to shift the momentum derivative fromǦ toΣ, c.f. the second last term of Eq. (18). Within the quasiclassical approach, the momentum derivative on the self-energy is also evaluated on the Fermi surface, Usually, the self-energy can be expressed as local GF (Σ(p, r) ∝Ǧ(p, r) ), then the second term in Eq. (26) resembles the familiar renormalization of the Fermi velocity due to self-energy v F → v F + ∂ piΣ (p, r). However, the self-energy itself can also be a function of nonlocal GF (Σ(p, r) ∝ ∂ iǦ (p, r) ) in the presence of SOC disorder, so one has to account for the first term in Eq. (26) during the identification of the Fermi velocity renormalization, as described in Sec. IV C. In Eq. (27), we neglect the Poisson bracket (i.e. linear order in gradient) term. This is because, as we will show below, when the linear in SOC self-energy is substituted into the Poisson bracket, it generates terms of the order of 3 1 and 2 1 2 , which are neglected in our approximation scheme. This is not the case for the anticommutator Eq. (26) which has to be kept.
From Eq. (24) and (25), we derive the the Usadel equation as follows. First we evaluate the collision integrals by expanding the self-energyΣ in terms of the small parameter λ 2 p 2 F up to second order (see next subsections): whereΣ (n) ∝ (λp F ) 2n . The zeroth order self-energy describes the usual Drude relaxation. The first and second order describes spin-charge coupling and the Elliott-Yafet spin relaxation process respectively. Then, we substitute all the self-energies in Eq. (29) into Eq. (27) to express the first momentǧ k in terms of the zeroth momentǧ, and obtain the so-called constitutive relation, Eq. (68). From this equation, we can infer the normalization conditioň g 2 = 1. Next, we substitute Eq. (29) into Eq. (26) and identify the anomalous currentJ an k and the spin-orbit torqueŤ . Lastly, we substitute the constitutive relation into Eq. (24) and arrive at the generalized Usadel equation.
We shall now follow the procedure described above and evaluate various self-energies in Eq. (29).
In this section, we evaluate the self-energy induced by spin-independent scattering potential. Throughout this article, the concentration of impurities, n im , is assumed to be small. In this dilute impurity limit, the self-energy depends only on the scattering properties of a single impurity, i.e. non-crossing approximation, see Fig. 2. Let us first neglect SOC and introduce the T-matrix describing the total scattering amplitude for an electron scattered by the scalar part of the impurity potential: where is the Fourier component of the single impurity potential in Eq. (9), and the superscript (0) reflects zeroth order in SOC. HereǦ p (r) is the GF in the Wigner representation that enters the full kinetic equation before quasiclassical approximation, Eq. (18). In Eq. (30), we assume a short-range impurity potential with a dominating s-wave scattering. Under this assumption, one can use the standard renormalization procedure to eliminate the high energy contribution to the momentum integral in Eq. (30) by introducing the physical scattering amplitude t 0 = 2πa/m at zero energy, where a is the zero-energy scattering length. As a result, the T-matrix is expressed in terms of the quasiclassical GFǧ which describes quasiparticle dynamics in the vicinity of the Fermi surface, and the physical scattering amplitude t 0 . The corresponding equation forŤ (0) , The scattering vertex of a) the spin-independent potential and b) the spin-orbit coupling potential acting on the annihilation field operatorΨ in the quasiclassic formulation. Subplot c) depicts the self-energyΣ (0) in Eq. (32). The star symbol represents the impurity density nim while the arrow line represents the Green functions. This "blob" diagram is later used to compose the self-energy of spin-charge coupling, Σ (1) .
which now describes the scattering of quasiparticle at the Fermi surface readš Here N F = mp F /2π 2 is the density of states at the Fermi energy. Note, since the dominant scattering wave induced by the impurity potential is assumed to be s-wave, only the isotropic part of the quasiclassical GF (the first term in Eq. 20) enters theŤ (0) . To leading order in impurity density n im , the self-energy can be expressed in terms ofŤ (0) as follow (cf. Fig. 2c): The quasiclassical GFǧ is a 8 × 8 matrix in spin-Nambu-Keldysh space, which, in the absence of SOC, obeys the normalization conditionǧ 2 (r) = 1. Using this condition, one can solve Eq. (31) explicitly: where we used the identity πN F t 0 = p F a, and introduced the scattering amplitude t F of an electron at the Fermi energy in a normal equilibrium system: It is worth mentioning that since the above expressions are local in space (i.e. they do not involve spatial derivatives), they are also valid for systems with spatially dependent concentration and/or type of impurities as schematically shown in Fig. 1. In other words, the impurity concentration n im , the scattering amplitudes t 0 and t F , and the scattering length a may depend on the (slowly varying) spatial coordinate r.
Unlike t 0 (r), t F (r) is a complex number and its complex phase satisfies the optical theorem: 3. The spin-charge coupling self-energy at linear order in spin-orbit coupling strength λ 2 p 2 F . There are two distinct class of Feynman diagrams: a) SOC vertex at the corner and b) SOC vertex in the middle. This has important consequences as discussed in the last paragraph of Sec. IV A 2. The shaded "blob" vertex stands the full Born series induced by spin-independent part of the impurity potential, while the cross represents the SOC vertex, see Fig. 2.
After substituting Eq. (33) into Eq. (32) and performing the Fourier transform with respect to the difference of the coordinates we arrive at the (zeorth order in SOC) self-energy in the Wigner representation: The elastic (Drude) relaxation time is expressed in terms of t F (r) which can model arbitrarily strong impurity potential: In the next two sections we analyze the SOC scattering at the impurities that leads to the spin-charge coupling and spin-relaxation. We should point out that a small contribution arising from SOC (i.e. λ 4 p 4 F or smaller) are neglected from the Drude relaxation time.

2.Σ (1) : Spin-charge coupling
Let us now include extrinsic SOC perturbatively. To leading order in SOC, we parametrize the renormalized spin-dependent part of the scattering vertex as which has the same form as the SOC term in Eq. (8), but with the bare impurity potential replaced with the zero-energy scattering amplitude t 0 . In particular, this means that the total scattering amplitude for electrons on the Fermi surface in the normal phase is approximated ast pp = t F + iλ 2 σ · (p × p )t 0 , where t F is given by Eq. (34). In other words, in the general form of the Mott scattering T-matrixt pp = A + iσ · (p × p )B/p 2 F , the scalar part A = t F is the full complex scattering amplitude of a scalar potential, while the coefficient B = λ 2 p 2 F t 0 is purely real, which corresponds to the leading perturbative correction due to SOC.
The linear in SOC self-energy,Σ (1) , is constructed from the vertex t so of Eq. (38) and the 8 × 8 zeroth order Tmatrix defined in Eq. (33). This results in two different class of diagrams shown in Fig. 3. The self-energy in Fig. (3a), which we defined asΣ 1a can be analytically read off as follow: For short-range randomly distributed impurities, the correlation functions needed for performing the impurity average in these diagrams are This closely resembles the correlators used in the calculations within the standard Born approximation t 0 (r 2 )t 0 (r 1 ) imp = n im t 2 0 δ(r 1 − r 2 ). However, there are two important differences. Firstly, unlike the scalar t 0 , the full scattering amplitudeŤ (0) (r) is a 8 × 8 matrix that depends on the GFǧ. Secondly, in our general setting we allow for a spatially varying disorder so that botȟ T (0) (r) and t 0 (r) may depend on spatial argument even after impurity average. Next, we substitute Eq. (33) and (40) into Eq. (39) and found that it can be conveniently decompose as follow:Σ HereΣ o (Σ e ) is a series that contains odd (even) power of the isotropic GF,ǧ introduced in Eq. (20). NoteΣ o (Σ e ) corresponds to substituting the first (second) term in Eq. (33) into Eq. (39). To proceed further, we perform the Fourier transform with respect to the relative coordinate to go to the Wigner representation. As detailed in Appendix A, the result takes the following form: Here ω 1 and ω 2 are spatially dependent scattering rates induced by SOC: where B(r) = λ 2 p 2 F t 0 (r) and t F is defined in Eq. (34). Notice that ω 2 is related to ω 1 by the optical theorem ω 2 (r) = ω 1 (r)p F a(r). This form of parameterization is commonly used in discussing the extrinsic spin-charge coupling in normal metallic state, see Ref. 33 and Ref. 38 for the discussion in 3D and 2D respectively. As mentioned earlier, the self-energy in Eq. (42) and (43) have terms that are proportional to different powers of small parameters 1 and 2 . For example, in Eq. (42), the first, second, third and forth terms are of the order of 1 , 2 , 1 2 and 1 2 respectively.
Note that the last terms in Eq. (42) and Eq. (43) become crucial near material boundaries or interface of two materials where the gradient of the scattering rates, ∂ j ω 1 and ∂ j ω 2 are significant. In fact, the current induced spin-orbit torque that drives spin diffusion in metals is precisely induced at the material boundary 55,56 . Within our quasiclassical approach, the spatial variation of the scattering rates has to be over distances much larger than the Fermi wave length.
An important observation is that the spin-dependent corrections to the self-energy, Eqs. (42) and (43), contain terms that are proportional to spatial derivative of the GF. These nonlocal terms give rise to (p F l) −1 correction to the spin-Hall angle and spin swapping coefficient. In the momentum space, these terms stem from the thin shell of the relative momentum near the Fermi surface. Hence, they are not captured in the equation of motion for the diagonal (in momentum space) component of the density matrix 38,57 . Let us now consider the self-energy coming from the diagram presented in Fig. 3b. The corresponding analytic expression takes the form In the second line, we substituted the SOC vertex t so from Eq. (38) and performed the partial integration over r 3 to move the derivative ∂ 3j toǦ(r 1 , r 3 ). To proceed further, we let K = [∂ 3jǦ (r 1 , r 3 )]σ a [∂ 3kǦ (r 3 , r 2 )] as a matrix valued function and note that it is independent of the impurity position. Then, we perform the impurity average in Eq. (46) as follow: In the last line, using the Delta functions, we seť T (0) (r 1 ) =Ť (0) (r),Ť (0) (r 2 ) =Ť (0) (r), and t 0 (r 3 ) = t 0 (r 2 ) = t 0 (r) where r = (r 1 + r 2 )/2. Next, we substitute Eq. (47) into Eq. (46) and perform the spatial integral (of r 3 ) to arrive at a self-energy that is purely local in space: whereΣ with the function K(r) defined as follow: Note thatΣ 1b (r) in Eq. (48) is, in fact, the self-energy in the Wigner representation required for the quasiclassical collision integral. Now we substitute the T-matrix from Eq. (33) and use the following identities to evaluate the function K(r) entering Eq. (49): These identities are derived by using the Wigner representation ofǦ(r, r 3 ) together with the definitions of the quasiclassical GF and its first two moments,ǧ andǧ k . Finally, after straightforward algebra, Eq. (49) reduces to the following form wherẽ In the above expressions we neglect the term proportional to the product ∂ jǧ ∂ kǧ as it is of the order of 2 2 . The latter is beyond our accuracy corresponding to retaining only terms linear in 2 and at most proportional to 2 1 or 1 2 , which exactly corresponds to the terms kept in Eq. (54).
In contrast toΣ 1a in Eqs. (41)-(43), the self-energyΣ 1b in Eqs. (53) does not depend on the external momentum. This has two important consequences for the derivation of the Usadel equation. Firstly, the anticommutator term in the zeroth momentǏ 0 of the collision integral in Eq. (26) vanishes forΣ 1b . Secondly, within our accuracyΣ 1b does not contribute to the first momentǏ k of the collision integral defined by Eq. (27). AsΣ 1b by itself contains terms proportional to 2 1 and 1 2 it brings toǏ k the corrections of the order of 3 1 and 2 1 2 which are irrelevant in our diffusive limit. Exactly the same arguments apply to the part ofΣ 1a that does not depend on the external momentum n. Therefore the self-energyΣ 1b and the nindependent part ofΣ 1a contribute only to the first term inǏ 0 of Eq. (26), while the n-dependent part ofΣ 1a gives nonvanishing contributions to bothǏ 0 andǏ k . We will make use of these properties later in Secs. IV B and IV C.
Before we close this subsection, let us emphasize the need to retain the Poisson bracket in Eq. (26). As shown in Eqs. (42), (43) and (53),Σ 1a andΣ 1b are made up by terms of the order 1 , 2 and 1 2 . When we substitute all of them into Eq. (26) to evaluate the collision integral, terms of order 1 in the self-energy can enter the Poisson bracket and generate terms of order 1 2 since the Poisson bracket involves spatial derivative. This is precisely the result one would get from substituting order 1 2 terms of the self-energy into the commutator in Eq. (26). Hence, in order to derive the Usadel equation in the presence of SOC correctly, it is necessary to retain the Poisson bracket in Eq. (26).

3.Σ (2) : Spin-relaxation
The last self-energy we consider is of the second order in the SOC strength. We only include the leading order term in the counting scheme which describes spin relaxation: Here the Elliott-Yafet spin relaxation time is given by It is worth mentioning that most of the works concerning the effect of extrinsic SOC in superconducting state 12,53,58 focus only on the spin relaxation described by ourΣ (2) while the self-energy related to the spincharge couplingΣ (1) was disregarded.

B. Constitutive relation and the normalization condition
After computing all the self-energies, the next step in the derivation of the Usadel equation is to eliminate the first momentǧ k in Eq. (25) by expressing it in terms of the zeroth momentǧ. In the following, we call the map g →ǧ k the constitutive relation as it is similar to the relation between the current and the density in the usual diffusion theory. By substituting Eq. (20) into Eq. (27) we rewrite the collision integralǏ k entering Eq. (25) more explicitly as followš As we argued in the discussion after Eqs. (53) and (54), the self-energyΣ (1) gives negligible (of the order of 3 1 and 2 1 2 ) contributions to the first term in Eq. (57). Therefore only the Drude self-energyΣ (0) contributes to Σ in Eq. (57). For the same reason, n kΣ in the second term of Eq. (57) is fully determined by the n-dependent part ofΣ 1a ; those n independent part ofΣ 1a contributes to the collision integral at order 3 1 and 2 1 2 . Thus, Eq. (25) leads to the constitutive relation: By substituting Eqs. (36) and (41) into this equation, and rearranging the terms we bring it to the following compact form, whereǍ k is given by the following expression, The first term in the above expression comes from the first (Drude) term in the right hand side of Eq. (58) while the rest corresponds to the second term in Eq. (58) and arises from the parts proportional to the external momentum n k in Eqs. (42) and (43). Importantly, the structure of Eq. (59) suggests that ∂ kǧ 2 = 0 and this allows us to impose the standard normalization condition on the Usadel Green function: To proceed further, we expandǧ k to linear order in SOC:ǧ k =ǧ Because of the normalization conditionǧ 2 = 1, this equation leads to the well-known solutioň where l = v F τ is the mean free path. At the linear in SOC order we substitute the zeroth order solutionǧ (0) k into the terms proportional to ω 1 and ω 2 in Eqs. (59)(60). This generates the following equation for the linear in where the parameters θ and κ are defined as follows The solution that satisfies Eq. (64) is given by, We combine the results of Eqs. (63) and (67) and find the final expression that relatesǧ k toǧ, This equation is the constitutive relation that we need for the derivation of a closed equation for the isotropic GFǧ. Typically the relations of this sort establish a connection between the diffusion "current" and the density. However, due to the presence of SOC, the first moment of the GF,ǧ k , is not the conserved current entering the continuity equation . This is because SOC depends on the particle's momentum and it produces an additional contribution to the current, the so called anomalous current. At the level of the kinetic theory, the anomalous current comes from both the commutator and anticommutator in the collision integral, Eq. (26). It is worth mentioning that θ and κ defined by Eqs. (65) and (66) are not the total spin Hall angle and the spin swapping coefficient. In Eq. (65), the first term is the skew-scattering contribution while the second term is only half of the side-jump contribution to the spin Hall angle. The other half comes from the anomalous current, exactly as it happens at the level of the Born approximation 59 . Similarly, for the spin swapping coefficient κ the second term in Eq. (66) will be doubled due to the anomalous contribution. We will return to this discussion in the next subsection after completing the derivation of the Usadel equation.

C. Usadel equation
Let us substitute the constitutive relation, Eq. (68) into Eq. (24) and obtain the generalized Usadel equation. In the absence of SOC, we recover the usual Usadel equation by setting the right hand side in Eq. (24) to zero, and by using the zeroth order constitutive relation, Eq. (63) on the left hand side. To second order in SOC, the Usadel equation reads: describes the standard (Elliott-Yafet) spin relaxation 12,53,58 and is obtained by substitutingΣ (2) , Eq. (55), into the first term of Eq. (26): The evaluation ofǏ (1) 0 is more cumbersome and requires some care. One has to substitute the first order self-energiesΣ 1a , Eqs. (41)- (43), andΣ 1b , Eqs. (53)-(54), into Eq. (26). Note that both the commutator and the anticommutator 60 terms in Eq. (26) contribute toǏ (1) 0 . After some lengthy algebra detailed in Appendix B, the collision integralǏ (1) 0 can be represented compactly as sum of two distinct contributions, The the 8×8 matrixŤ and the matrix-valued vectorJ an k are defined as follows, Due to the way the above quantities enter the diffusion equation, we identifyJ an k as the anomalous current anď T as the spin-orbit-torque.
It is important to emphasize that all the kinetic coefficients, D, ω 1 , ω 2 , τ, l depend on the spatial coordinate r. Therefore, one cannot redefine the spin-torque and anomalous current by simply absorbing part ofŤ into ∂ kJ an k or vice versa. In other words, the definition of T and ∂ kJ an k is unambiguous when we allow the kinetic coefficients to vary in space. 61 Next, we move the total divergence of the anomalous current (i.e. −∂ k J an k ) to the left hand side of Eq. (69) and define the total current as the sum of the first moment v Fǧk /3 (cf. Eq. 68) and J an k : Here the total spin Hall angle and spin swapping coefficient are given by the following expressions, Equation (74) is the result announced in Eq. (2) of Sec. II. Interestingly, the term describing the spin Hall effect in Eq. (74), i.e. the term proportional to θ, has exactly the same form as the spin Hall term obtained in superconductors with intrinsic spin-orbit coupling 30 . This means that in systems with both extrinsic and intrinsic SOC the coefficient describing the charge-spin coupling is simply the sum of the two contributions.
As we have already discussed in the previous subsection, the first (second) term in Eq. (75) corresponds to the skew-scattering ( side-jump) contribution to the spin Hall angle θ. One half of the side-jump contribution comes from the SOC correction to the anisotropic part of the GFǧ (1) k while another half appears from the anomalous current.
Similarly, the spin swapping coefficient κ also receives contributions from two independent scattering mechanisms. The first term in Eq. (76) exactly reproduces the swapping coefficient identified by Lifshitz and Dyakonov 33 . In addition we found another contribution given by the second term in Eq. (76) which is proportional to (p F l) −1 and relies on the nonlocality of the selfenergy. This term gives rise to a "non-local" component to the spin current swapping effect, which is formally similar to the side-jump component of the spin Hall effect. In fact the one half of the "non-local" contribution to κ comes from the "normal" and another half from the anomalous currents. Note also that the side jump contribution to θ and the "non-local" contribution to κ scale in exactly the same way with respect to the impurity concentration -both are proportional to n im and thus inversely proportional to the Drude conductivity. In a homogeneous system where gradient of kinetic coefficients vanished ∂ k θ(r) = ∂ k κ(r) = 0, the divergence of the current (or the source of diffusion) takes the usual form ∂ kJk = −∂ k (Dǧ∂ kǧ ). Hence, a flow (current) induced by SOC typically occurs at material boundaries and/or interface of two materials where ∂ k θ(r) and ∂ k κ(r) are finite. Although our quasiclassical theory cannot describe boundary effects which occurs at the scales smaller than the mean free path, it provides an unambiguous definition of the generalized matrix current whose conservation define a boundary condition for the kinetic equations. The same philosophy has been used extensively to model spin-charge conversion also in the normal state 55,56,62 . In the next section we compile all the results above and obtain the generalized Usadel equation.

V. DISCUSSION
Substitution of Eq. (68), (70) and (71) into Eq. (69) leads to the main result of our paper, the generalize Usadel equation: (77) Recallǧ =ǧ(r, t 1 , t 2 ) and the generalized currentJ k is defined in Eq. (74) and the SOC induced torque on the right hand side is given by In order to discern the physics behind these expressions, it is helpful to first study Eq. (77) in a (normal state) metal. The advantage of using quasiclassical equation is that it describes both superconducting and normal state in a coherent manner. Indeed, the normal state diffusion equation can be readily obtained by setting the retarded and advanced GFs to +τ 3 δ(t 1 −t 2 ), −τ 3 δ(t 1 −t 2 ) respectively, and the time t 1 = t 2 = t in the Keldysh GF. Then, from Eq. (77), we can obtain the well-known spin diffusion equation after multiplying it by the vector σ and taking the trace 63 where S a = −(π/4)N F Tr σ a τ 3ǧ K is the a-component of the non-equilibrium spin density. It measures the deviation from the equilibrium spin-density. The equilibrium spin density can induce, for example, from a static Zeeman field. Similarly, one can obtain the diffusion equation for the charge density, n = −(π/4)N F Trǧ K n by taking the trace over spin σ.
Notice that the termŤ on the right hand side of Eq. (77) does not contribute to the normal state spin diffusion equation, Eq. (79). In other words, the right hand side of Eq. (79) which describes spin torque contains only the well-known Elliott-Yafet spin relaxation term. Mathematically, the reason for the vanishingŤ contribution is due to the fact that this term contains products of two derivatives of GFs. In order to obtain the Keldysh component ofŤ , one necessarily needs to differentiate at least one retarded or advanced GF which are constant in space. We will see that this is different in the superconducting state whereŤ = 0.
Next, we discuss the equation of motion for the currents. By taking the trace over spin in Eq. (74), one obtains the charge current, which in the normal state reads: The spin current is obtained by multiplying Eq. (74) with σ b and taking the trace: The first term on the right hand side of these two equations are the charge and spin diffusion currents respectively. The terms proportional to θ couple charge and spin degrees of freedom and they describe the spin Hall effect and inverse spin Hall effect. They stem from the anti-commutator in Eq. (74). The third term on the right hand side of Eq. (81), leads to the swap of the spin and direction indexes of the spin-current 64 . Thus, the generalized Usadel equation leads, in the normal state, to the known effects regarding spin-charge coupling.
Having established the known results in the normal phase, let us now discuss Eq. (77) in the superconducting phase which is characterized by non-trivial excitation spectrum, i.e. the retarded and advanced Green functions are non-trivial functions of the time difference t 1 − t 2 and of the space coordinate. More importantly, their matrix structure in the Nambu space (c.f. Eq. 13 ) carries finite anomalous (off-diagonal) components which leads to, among many other phenomena, the equilibrium supercurrent and the equilibrium magnetoelectric effect 31 which we will discuss below. In the normal metallic phase, the charge and spin current is related to different components of the gradients of charge and spin densities as in Eqs. (80-81). However, this is not generally the case in superconductor due to the the non-trivial dependence of retarded and advance GFs on energy. For example, while the second (spin Hall) term on the right hand side of Eq. (74) can indeed be written as a spin or charge (spectral) density, the third (swapping) term proportional to κ cannot, due to the productǧ∂ kǧ .
We start by inspecting the expression for the 8 × 8 matrix current described in Eq. (74). The charge current is obtained by taking trace in Eq. (74) after multiplying it with the Nambu matrix τ 3 . In the absence of SOC, only the first term in Eq. (74) is non-vanishing. Notice that even in equilibrium, this term can result in a finite current if the anomalous GFf and its time-reversal conjugatef c are different. For instance, this can happen in a bulk superconductor with a finite phase gradient which leads to the well-known supercurrent. Similarly, there exists equilibrium magnetoelectric effects in superconductor with SOC, i.e. spin to charge conversion in the absence of spin-injection (pumping) field. For example, a static Zeeman field can polarize the condensate and creates triplet correlations. These correlations enter the second (SH) term in Eq. (74) and contribute to the charge current. Reciprocally, a charge supercurrent may induce a spin current in a superconductor. These effects have been studied in Ref. 31 within the first Born approximation.
Interestingly, the expression for the current, Eq. (74), has exactly the same form as in the Born approximation studied in Ref. 31. However, unlike Ref. 31, the analysis in the present work is valid to all orders in the scalar elastic impurity potential and includes resonant skew scattering and side-jump mechanism. This results in the renormalization of the kinetic coefficients (spin Hall angle and swapping coefficient) discussed in Eqs. (75)-(76).
More importantly, we found a hitherto unknown term in the Usadel equation:Ť in the right hand side of Eq. (77) and defined in Eq. (78). SinceŤ is a commutator with the spin Pauli matrix, it vanishes under trace in spin-space. In other words, it does not modify the (spectral) charge diffusion equation and only enters the spin diffusion equation as a spin-orbit torque whose magnitude is characterized by the spin Hall angle and swap current coefficient. As mentioned above, this torque vanishes in the normal state [cf. Eq. (81)] hence it describes a spin torque unique for the superconducting state.
We now illustrate with an example the consequences of this new term by considering a superconductor with a single magnetic vortex. The latter is described by a spatial dependent order parameter ∆ = |∆(r)|e inϕ , where r is the radial component of the position vector r, ϕ = tan −1 (y/x) is the polar angle and n is the vorticity or the topological charge. The axis of the vortex is in z-direction. Let us assume that the superconductor is subject to a homogeneous spin-splitting (Zeeman) field −h a σ a where h a is the unit vector pointing along the direction of the Zeeman field. In the absence of SOC, the normal and anomalous components of the GF, Eq. (13), have the following form: Note the ϕ dependence enters only in the anomalous GFf (r, ϕ). Next, we substitute Eq. (82) and (83) into Eq. (77) to compute the spin torque to linear order in SOC. The result reads: Here the function F (r) is a function of the radial components and can be computed by solving the Usadel equation, but its explicit form is not relevant for this discussion. Due to the Nambu matrix τ 3 in the trace of Eq. (84), F (r) contains products of the anomalous GF and its time-conjugated with derivatives with respect to r and ϕ. Therefore, this spin torque only appears in the superconducting state. For the example consider here the amplitude of the torque decays away from the vortex core, i.e. with increasing r. From Eq. (84), we conclude that the spin-torque generated byŤ is proportional to the vector product between the angular momentum of the condensate n z (where n is the vorticity of the vortex) and the triplet vector h. Let us emphasis all the spin-charge conversion (or spinorbitronic) effects we discussed in this article occurs at leading order in spatial non-uniformity (∂ 2 i ) of the Usadel equation. The characteristic scale of the spin torquě T is determined by the spin Hall angle θ and swap current coefficient κ defined in the normal metallic state. For example,Ť can be observed in a thin layer of superconducting proximitized Pt or Ta (θ ∼ 0.01 − 0.1 from Ref. 19), which are the typical metals used for measuring the spin Hall effect. Notice also that the spin Hall angle in superconducting Nb, although small, has been recently quantified in the experiment of Ref .65

VI. SUMMARY AND OUTLOOK
We have systematically extended the Usadel equation to incorporate spin-orbit coupling disorder in diffusive superconducting systems. In addition to the spin Hall effect and the swap-current effect describe by the current op-eratorJ k , we identified a non-linear spin-orbit coupling induced torqueŤ that has no counterpart in the normal metallic state. Interestingly, the torque is parameterized by the same spin Hall angle θ and swap current coefficient κ that one would obtain in the normal metallic state. Note that our generalization of Usadel equation accounts for spatially varying kinetic coefficients. By imposing suitable boundary conditions on the Usadel equation, it can readily describe both equilibrium and out-of equilibrium phenomena related to spin-charge coupling in diffusive hybrid-structures made by superconductor, ferromagnetic and normal metal. In this appendix, we detailed the derivation of the self-energyΣ 1a , as shown in Fig. 3a of the main text. It can be diagrammatically read off from Fig. 3a as follow, Here t so (r) = −iλ 2 ajk σ a ∂ j t 0 (r)∂ k is the disorder SOC induced by the impurity potential t 0 (r). For shortrange randomly distributed impurities, the correlation between the (scalar) disorder potential t 0 (r 1 ) and the matrix-valued T-matrix,Ť (0) (r 2 ) is given by the fol- . Substitute the impurity average into Eq. (A1), we arrive at the following: where the T-matrix (without SOC) is given in Eq. (33): Next, we substitute Eq. (A3) into the self-energy to arrive at the following equation: Here the derivative always acts on its immediate neighbour.Σ o (r 1 , r 2 ) andΣ e (r 1 , r 2 ) are characterized by two scattering rates on the Fermi level: ω 1 (r 1 ) = 2πn im N F Re t * F (r 1 )B(r 1 ) and ω 2 (r 1 ) = 2πn im N F Im t * F (r 1 )B(r 1 ) where B(r 1 ) = λ 2 p 2 F t 0 (r 1 ) is the SOC scattering vertex. In deriving the above, we used Re t F (r 1 ) = Re t * F (r 1 ) and Im t F (r 1 ) = −Im t * F (r 1 ). To proceed further, we shift all quantities to the center of mass coordinate r = (r 1 + r 2 )/2 and relative coordinate s = r 1 − r 2 : Here r ± = r ± s/2 and ∂ k ( ∂ s k ) are the spatial derivative along direction k on the variable r (s). Using these equations, the self-energy can be written down as follow: Next, we Fourier transform the relative coordinate of the self-energy (s) into momentum p: In order to do so, we express all the Green functions in the Wigner coordinates: Let us begin by evaluatingΣ o (r, p). We substitute the first equation in Eq. (A13) into Eq. (A10) and apply the Fourier transform in Eq. (A12) to arrive at the following: The intricate momentum summation above can be done in the quasiclassical limit. For p −1 F to be much smaller than the typical variation of spectral weights and observables, the Green functions will be peaked at the Fermi level so their momentum summation can be simplify as follow: Hereǧ (r) is the isotropic (zeroth-moment) of the Green function whileǧ k (r) is the (first-moment) of the Green function along spatial direction k. Since the relevant self-energy that enters the collision integral has its momentum parked at the Fermi level, we setΣ o (r, p = p F n) =Σ o (r, n) and arrive at the equation quoted in the main text: Note thatΣ o (r, p) is now an algebraic equation expressed in terms of the zeorth and first moment of the quasiclassical Green function, this is a tremendous simplification compare to the previous equation whereΣ o (r, p) is a functional of the Green function.Σ e (r, p) can be derived following the same procedure: we substitute Eqs. (A13) and into Eq. (A11) and apply the Fourier transform Eq. (A12) to arrive at the following: Next, we use Eq. (A15) to perform the momentum integration and set p = p F n to arrive at the result in the main text.

Appendix B: Collision integral of the Usadel Equation
In this appendix, we explain in detail how we arrive at the result I Noteǧ =ǧ(n, r) is the Eilenberger GF and readers should not be confused with the zeroth moment (Usadel) GF g =ǧ(r). To linear order in SOC, there are two class of self-energy diagrams shown in Fig. 3 of the main text which are denoted asΣ 1a andΣ 1b . The corresponding contribution toǏ 0 are calculated below.

Contribution ofΣ 1a to the collision integral I
(1) 0 We first consider the self-energy shown in Fig. 3a which we defined asΣ 1a . From Appendix. A and the main text, Σ 1a =Σ o +Σ e reads as followš Σ o (n, r) = ω 1 ajk 2 Σ e (n, r) = −i ajk ω 2 2 1 3 n j (σ aǧ kǧ −ǧǧ k σ a ) − n j 2p F (σ a (i∂ kǧ )ǧ +ǧ (i∂ kǧ )σ a ) + i∂ k 6p F (σ aǧ jǧ +ǧǧ j σ a ) − ajk ∂ j ω 2 12p F σ aǧ kǧ +ǧǧ k σ a (B3) In the following, we substitute each of the above contributions into the collision integral of the Usadel equation: Let us begin withΣ o (n, r). We substitute Eq. (B2) into Eq. (B4) and arrive at the following result: The collision integral of the Usadel equation in Ref. 31 is obtained by considering only the first line above. However, there is typo in the supplementary material of Ref. 31 where we skipped the first term inǏ (o) . We go beyond Ref. 31 and include spatial variation of ω 1 . Note that the last two term in the first line of Eq. (B5) arised from the anticommutator in Eq. (B4). They contribute to the same order of magnitude as the third and fourth term in the first line of Eq. (B5). Therefore, as emphasized in the main text, when the system is subjected to disorder SOC, it is crucial to retain the Poisson bracket in the collision integral of the kinetic equation.
Recall that the constitutive relation obtained in the main text is given by the following: To proceed further, we substitute Eq. (B6) into (B5). To linear in SOC, this amounts to replacingǧ k = −lǧ∂ kǧ in Eq. (B5) and the resulting equation reads: Importantly, the space derivative (i.e. ∂ k ) in the last term applies to the quantities inside the bracket as chain-rules, i.e. it is a total divergence. Recall that ω and l depend on the spatially coordinate r. In deriving the above, we used the property ǧ , ∂ kǧ = 0 which follows from the normalization conditionǧ 2 = 1 we found in Eq. (59). Next, we use D = v F l/3 to simplify the coefficients above to arrive at the following equation: akj 2ω 1 τ p F l σ a ,ǧ∂ kǧ ∂ jǧ + 2 3 ω 1 τ i ∂ kǧ ∂ jǧ , σ a ] − ∂ k Dω 1 τ 2p F l akj ∂ jǧ , σ a (B8) The first two terms correspond to part of the spin-orbit torqueŤ while the last term corresponds to the anticommutator part ofJ an k , c.f. Eq. (71) Next, we follow the same logic above to derive the collision integral from Eq. (B3). First, we substitute Eq. (B3) into Eq. (B4) and replaceǧ k = −lǧ∂ kǧ , the resulting equation reads: I (e) = ajk ω 2 − l 2 18 σ a , ∂ jǧ ∂ kǧ − l 2 18 ∂ jǧ [σ a ,ǧ]∂ kǧ + l 6p F i σ a , ∂ jǧ ∂ kǧ − i∂ k 6p F l ǧ∂ jǧ , σ a − ajk l 6p F i ǧ∂ jǧ , σ a ∂ k ω 2 (B9) In the above, the last two terms can be combined into a total divergence. Next, we use D = v F l/3 to simplify the coefficients above and arrive at the following equation: I (e) = D 4 akj 2 3 ω 2 τ σ a ,ǧ∂ kǧ ∂ jǧ − ω 2 τ p F l i ∂ kǧ ∂ jǧ , σ a ] −∂ k akj Dω 2 τ 2p F l i [σ a ,ǧ∂ jǧ ] − ajk Dω 2 τ 6 ∂ jǧ [σ a ,ǧ]∂ kǧ (B10) Let us now combine Eq. (B7) and Eq. (B10) and arrive at the following: The last term stems from the second term inǏ (e) . As we will show in next section, as a consequence of the optical theorem, this term is canceled exactly by the contribution coming from the self-energyΣ 1b shown in Fig. 3b of the main text.
2. Contribution ofΣ 1b to the collision integral I (1) 0 The expression forΣ 1b is derived in the main text and given by Eqs. (53)-(54) which we repeat here: As stated in the main text,Σ 1b does not carry the external momentum. As a result it does not contribute to the collision integralǏ k of the constitutive relation as it would lead to terms of the order of 3 1 and 2 1 2 which are irrelevant in our approximation. Secondly, whenΣ 1b is substituted into the collision integral of Eq. B1, the anticommutator part, which is proportional to ∂ p Σ(p, r), vanishes and we are left with the following simple resultš I 1b = −i Σ 1b ,ǧ . (B14) Next, using the normalization condition (i.e.ǧ 2 = 1) we showed from Eq. (59), the collision integral can be greatly simplified. First, note that the last term in Eq. (B12) does not contribute toǏ 1b because ǧK +Kǧ,ǧ =ǧKǧ +Kǧ 2 −ǧ 2K −ǧKǧ = 0.
Note also that because ǧKǧ,ǧ = −[K,ǧ] the second term in Eq. B12 reduces to the form similar to the first term. Hence, by collecting these results, the collision integral can be simplified as followš To linear order in SOC, we substituteǧ k = −lǧ∂ kǧ to evaluate K ,ǧ and the resulting collision integral reads: This has exactly the same matrix structure as the last term in Eq. (B11). To see how the coefficients add up, we invoke the optical theorem of Eq. (35) which follows from unitarity of the scattering S-matrix: This brings our collision integral to the following form I 1b = ajk l 2 ω 2 18 ∂ jǧ [σ a ,ǧ]∂ kǧ = ajk Dω 2 τ 6 ∂ jǧ [σ a ,ǧ]∂ kǧ , where we used the definition ω 2 = 2πn im N F Imt * F λ 2 p 2 F t 0 .
Finally, by adding the derivedǏ 1b to the collision integral fromΣ 1a =Σ o +Σ e , Eq. (B11), we arrive at the result quoted in Eq. (71) of the main text:Ǐ