Controlling spin supercurrents via nonequilibrium spin injection

We propose a mechanism whereby spin supercurrents can be manipulated in superconductor/ferromagnet proximity systems via nonequilibrium spin injection. We find that if a spin supercurrent exists in equilibrium, a nonequilibrium spin accumulation will exert a torque on the spins transported by this current. This interaction causes a new spin supercurrent contribution to manifest out of equilibrium, which is proportional to and polarized perpendicularly to both the injected spins and the equilibrium spin current. This is interesting for several reasons: as a fundamental physical effect; due to possible applications as a way to control spin supercurrents; and timeliness in light of recent experiments on spin injection in proximitized superconductors.

There are four relevant species of quasiparticles in the systems that we will consider: namely electrons and holes, which each have two distinct spin projections. These have the densities where Ψ † σ and Ψ σ are standard creation and annihilation operators. For comparison, the propagators are defined as [S1-S3]: where the subscripts σ and σ denote possible spin projections. Combining these definitions, we see that the quasiparticle densities are directly related to the equal-coordinate propagators: These expressions can be used to calculate the spin-resolved density of electrons and holes, respectively. Note that holes carry both opposite charge and opposite spin compared to electrons [S26]. The charge and spin accumulations are then found by multiplying each quasiparticle density with their respective charges or spins, and summing up their contributions: ρ e e 1 2 [+n e↑ + n e↓ − n h↑ − n h↓ ], (S10) where we use the convention that e is the electron charge (e < 0). The prefactors 1/2 are required to prevent double-counting, and can be explained as follows. If we add one physical electron to the system, then the charge of the system increases by e. However, the number of electrons increases by one, and the number of holes decreases by one, meaning that the difference between electrons and holes increases by two. Thus, when the charge density ρ e is described in terms of both electrons and holes, we need an extra factor 1/2 to get the right physical charge. The same logic applies to the spin density ρ z . We can rewrite the results in terms of the propagators above, and recognize the remaining sum as a trace over spins: There is nothing special about the spin-z axis, so it is straightforward to generalize this result to arbitrary spin-projections: where σ = (σ 1 , σ 2 , σ 3 ) is the Pauli vector. From the definition of the Keldysh propagator above, we can also use the identity AB * = B † A † to show that G K * σσ = −G K σσ . This means that G K σσ is imaginary, which makes ρ e , ρ s ∼ iG K manifestly real. For later convenience, we will therefore write this out explicitly:

B. Quasiparticle currents
Now that we know the charge and spin accumulations, the next step is to find the corresponding currents. To derive these, we go back to the quasiparticle densities defined in Eq. (S1): To rigorously derive expressions for the charge and spin currents, we will use the definitions above to look for quasiparticle continuity equations on the form where j τσ is the particle-and spin-resolved current density we are interested in, while q τσ represents possible source terms. We start by differentiating the densities with respect to time: We can rewrite the above using the Heisenberg equation of motion for the field operators. Note that any contributions to the continuity equation arising from non-derivative terms in the Hamiltonian-such as a superconducting gap or an exchange field-can be incorporated into the source term q. Thus, for the purposes of deriving current equations, it is sufficient to consider only derivative terms. Whether or not the currents we derive are conserved currents can be checked at the end of the derivation, by substituting the Usadel equation into the final quasiclassical current equations [S4, S5]. If we for simplicity disregard gauge fields for now, the equations reduce to: We then substitute these back into the equations for ∂ t n τσ : Thanks to cancellation of cross-terms, these can be factorized: Comparing this to Eq. (S20), we conclude that: As a mathematical trick, let us now use different coordinates Ψ σ = Ψ σ (r, t) and Ψ † σ = Ψ † σ (r , t ) for the field operators, where we let r → r and t → t in the end. In this case, the differential operators acting on the field operators can be factored out of the expectation value without ambiguity: We are now ready to define the charge and spin current densities. In correspondence with Eq. (S10), we define these as: Substituting Eq. (S31) into Eqs. (S33) and (S34), comparing the results to the Keldysh propagator in Eq. (S7), and recognizing the results as traces in spin space, we conclude: Generalizing to all spin projections, we obtain the final results: We wish to point out that these currents are manifestly real. From the definition of the Keldysh propagator, we see that: But which set of coordinates we chose to call (r, t) and (r , t ) was arbitrary, and should not affect the physical results, since we are considering the limit r , t → r, t anyway. This means that we can interchange the coordinates (r, t) and (r , t ) on the righthand side of the equation, as long as we do this consistently for every factor simultaneously. The coordinate interchange leads to a sign flip in (∇ − ∇) which cancels the minus sign inside the brackets, and makes the two sides of the equation equal. This lets us conclude that (∇ − ∇)G K * = (∇ − ∇)G K , which in turn implies that the charge and spin currents are real. For later convenience, we can therefore rewrite the above as

C. Quasiclassical and diffusive limits
To derive equations we can use together with the Usadel equation, we now follow the standard prescription for taking the quasiclassical and diffusive limits [S1-S3, S6]. The net change to the Keldysh propagator and its derivative are then: 1 2m where v p/m is interpreted as the quasiparticle velocity, N 0 is the density of states at the Fermi level, and · · · F refers to the average over the Fermi surface. From the derivation of the Usadel equation, we also know that in the diffusive limit the Fermi-surface averages can be written [S6, S7] where∇ is a gauge-covariant derivative including the electromagnetic vector potential and spin-orbit interactions [S6, S7, S20],ǧ s is the isotropic propagator, and D is the diffusion constant. We drop the subscripts on the isotropic propagatorsǧ s , and substitute the above into the accumulations and currents: where we have reintroduced the matrix currentǏ D(ǧ∇ǧ). Note that these equations only depend on the "electronic" part of the propagators in Nambu space, which in reality contains information about both the electrons and holes in the system.
All these results can be written as integrals over only positive energies using the symmetries of the Nambu-space matriceŝ In other words, the negative-energy contributions can be recast in terms of the lower-right blocks; and since take the real part of the results, the complex conjugations are irrelevant. The remaining structure can be recognized as a trace over Nambu space, yielding the final quasiclassical transport equations Note thatσÎ K should be interpreted as an outer product between two vectors, which results in a rank-2 tensor. This is because a general description of spin transport requires both a direction of transport ∼Î K and a spin orientation ∼σ.

D. Higher-order gauge contributions
The equations of motion for the field operators also include first-order derivative terms in systems with electromagnetic [S6, S9] or spin-orbit [S7, S10, S11, S20] gauge fields. If we ignore all other terms in the Hamiltonian, these derivative terms give the following Heisenberg equations: where we implicitly sum over the spin index σ . Going through the same kind of derivations as without the gauge fields, we find that we basically just have to make the following replacement in the results right before taking the quasiclassical limit: Note that the gauge fields also affects charge and spin transport in a different way, since they also appear as covariant deriva-

A. Supercurrents vs. resistive currents
As shown in previous sections, the total spin current J s can in the quasiclassical limit be calculated as an energy integral, where the spectral spin current j s Re Tr[τ 3σÎ K ]/8 and the matrix currentÎ Dǧ∇ǧ. If we substitute the parametrization g K =ĝ Rĥ −ĥĝ A into the definition of the matrix current, we find that its Keldysh component can be expanded aŝ The terms on the first line may be finite even for a homogeneous distribution functionĥ, and produces spin currents even in equilibrium. Furthermore, they are sensitive to the phase-winding of the superconducting condensate viaĝ R ∇ĝ R andĝ A ∇ĝ A . We therefore identify this as a supercurrent contribution. The terms on the second line, however, are proportional to ∇ĥ. This current contribution both requires an inhomogeneous distribution function, and is insensitive to the phase-winding of the superconducting condensate, and has to be a resistive current.
In this work, we are primarily interested in generating a spin supercurrent from a nonequilibrium spin accumulation. We therefore limit our attention to systems with a positionindependent distribution functionĥ that has an excited spin mode. Since we assume ∇ĥ = 0, the second line of Eq. (S59) disappears, and only the supercurrent contribution remains: As for the distribution function, it can be parametrized aŝ where h s points along the net quantization axis of the accumulated spins, and the magnitudes of the modes above are (S63) Note that the energy mode h 0 and spin mode h s are odd and even functions of energy, respectively. We have parametrized the spin mode in terms of a spin voltage V s (V ↑ − V ↓ )/2, where V σ are the effective potentials experienced by spinσ quasiparticles [S12-S14]. The spin mode h s is related to the spin accumulation in Eq. (S52) by an energy integral where N( ) is the density of states [S13].

B. Expansion in Pauli matrices
Once we substitute Eq. (S61) into Eq. (S60), there are a few subtleties to be careful about. To handle these, without yet introducing all the details of the singlet/triplet-decomposition, we first expandĝ R ∇ĝ R directly in terms of Pauli matrices: The first four terms parametrizes a general block-diagonal matrix, while the last termˆ represents off-block-diagonal parts. Since the distributionĥ can always be chosen to be block-diagonal,ˆ does not contribute to the trace in Eq. (S60). The other coefficients are found by taking appropriate traces: We parametrizeĝ R † ∇ĝ R † using coefficients α, β, γ, δ that are defined in the same manner as above.
We will now argue that the parameter δ is identically zero. By differentiating the normalization condition (ĝ R ) 2 = 1, one can show that the retarded propagator anticommutes with its gradient, {ĝ R , ∇ĝ R } = 0. This identity can be rewritten Let us now trace both sides of the equation, and use the cyclic rule Tr[ÂB] = Tr[BÂ] on the right-hand side, Sinceσ 0τ0 is an identity matrix, we see from Eq. (S65) that: In other words, δ = 0 is always satisfied, as any other conclusion would violate the normalization condition (ĝ R ) 2 = 1. Next, to clarify another subtlety, we need to derive some trace identities. By explicitly writing out the matrix products and usingσ = diag(σ, σ * ), one can show that Products of spin matrices in general satisfy (a · σ)(b · σ) = (a · b) +i(a × b) · σ; multiplying by σ and taking the trace, we find the associated trace rule Tr[(a · σ)(b · σ)σ] = +2i(a × b). However, if we complex-conjugate before taking the trace, we uncover another identity Tr[(a · σ * )(b · σ * )σ * ] = −2i(a × b). A geometric motivation for the sign difference is that if the basis σ = (σ 1 , σ 2 , σ 3 ) defines a right-handed coordinate system, then σ * = (σ 1 , −σ 2 , σ 3 ) has to define a left-handed one-and this inverts the right-hand rule that cross-products usually satisfy.
With the aid of the results above, we see that This is the subtle trap alluded to above: due to the way we defineσ = diag(σ, σ * ), the generalization of the Pauli crossproduct identity to matrices in Nambu space requires an extra factorτ 3 in the trace to produce a nonzero result. We now substitute Eqs. (S61) and (S64) into Eq. (S60). With the identities above, we see that the only contributions are: By multiplying Eq. (S66) by appropriate Pauli matrices, taking traces, and using Tr[Â † ] = Tr[Â] * , one can show that α = −α * . This makes α−α real and α+α imaginary, so both contributions are compatible with the normalization condition. We could also use this information to eliminate the underlined coefficients, but this would make it harder to see how mixed singlet/tripletterms cancel later in the derivation. Interestingly, all spin supercurrent contributions depend on the same coefficient α, and do not couple to the other traces ofĝ R ∇ĝ R . The physically observable spin supercurrent is found by integrating the spectral current over all positive and negative energies. We also know that h 0 and h s are odd and even functions of energy, respectively. We can therefore let α (+ ) → ∓α (− ) = ∓α * (+ ) in the spectral current without changing the total spin supercurrent: This form of the result will be useful later, as it makes it clearer which parts of the non-underlined and underlined coefficients cancel for symmetry reasons. Conveniently, this also makes the h 0 and h s contributions take very similar forms.
We now proceed with an expansion of the propagators in terms of physically meaningful components. Following the same kind of parametrization as Ref. S15, we can writê Here, f s represents the spin-singlet pair amplitude, while f t is the spin-triplet amplitude. On the other hand, we can interpret g s and g t as the spin-independent and spin-dependent parts of the density of states, respectively [S15]. In our notation, this means that the density of states for particles with spinprojection p is given by N = N 0 Re[g s + g t · p]. In equilibrium, the spin accumulation is found by integrating h 0 g t over energies, giving another interpretation of g t . Outside of equilibrium, we of course get another kind of spin accumulation due to a nonzero spin mode h s , which we are interested in here. Using Eq. (S75) and the identity σ 2 σσ 2 = −σ * , we find that the diagonal components ofĝ R ∇ĝ R in Nambu space are where the subscripts [· · · ] i, j are matrix indices in Nambu space.
Using the identity (a · σ)(b · σ) = (a · b) + i(a × b) · σ and its conjugate (a · σ * )(b · σ * ) = (a · b) − i(a × b) · σ * , we can sort the above into spin-independent and spin-dependent terms, Since we defineσ = diag(σ, σ * ), Eq. (S65) tells us that the coefficient α that we require can be expressed as Together with the expansion ofĝ R ∇ĝ R above, and standard trace identities for Pauli matrices, we then obtain Let us now calculate the corresponding coefficient α from the matrixĝ R † ∇ĝ R † . Taking the complex-transpose of Eq. (S75), we see thatĝ R † changed as follows compared toĝ R : (S81) Other than these transformations, the parametrization is clearly identical, and the derivation of α becomes identical as well. If we in the end results also choose to let → − , corresponding to a combination of complex-conjugation and tilde-conjugation, the net transformation rules become We can therefore simply perform the changes above to Eq. (S79) to get the corresponding equations forα * : We are now ready to calculate the spectral spin supercurrent in terms of the singlet/triplet-decomposition. Adding up Eqs. (S79) and (S83), we see that all mixed singlet/triplet terms drop out, and we are left with only the cross-product terms: Substituting this into Eq. (S74), we immediately see that: We have shown earlier in the derivation that both contributions are compatible with the normalization condition. The fact that they did not cancel during the last simplification above, shows that both contributions are compatible with the energy symmetries of h 0 and h s . Finally, we know that the contents of the brackets g t × ∇g t − f t × ∇f t can be nonzero, since this is the source of equilibrium spin currents.
The final result shows that if one in equilibrium has a spin supercurrent j eq s , then a nonequilibrium spin mode h s gives rise to a new component j neq s ∼ j eq s × h s . This can intuitively be interpreted as the injected spins h s exerting some kind of torque on the spins transported by the equilibrium current j eq s , thus producing a component j neq s that is spin-polarized in a direction perpendicular to both. This analogy is not perfect: it leaves out the Im and Re operations in Eq. (S86), and the fact that the cross-product relation is between spectral currents and accumulations. However, the intuition provided by this picture is sufficient to explain the results in the main manuscript.

III. NUMERICAL MODEL
As summarized in the main article, we perform the numerical calculations using the Usadel formalism [S1-S3, S13, S16]. This is formulated in terms of the 8×8 quasiclassical propagatoř which satisfies the identitiesĝ K =ĝ Rĥ −ĥĝ A andĝ A =τ 3ĝ R †τ 3 . Together, these identities show that we have to determine two 4 × 4 matrices to knowǧ: the retarded propagatorĝ R , which determines the spectral properties of a material; and the distribution functionĥ, which describes the occupation numbers of the states in the material. Both of these matrices can be functions of position r and quasiparticle energy .
In general, the distribution functionĥ follows from solving a kinetic equation that can be derived directly from the full 8 × 8 Usadel equation. We present a complete derivation of a kinetic equation and relevant boundary conditions in Ref. S12, which is valid for quite general superconducting structures. The result is formulated as an explicit and linear differential equation, which can be easily and efficiently implemented in a numerical Usadel solver. Related derivations can be found in Refs. [S2, S3, S13, S17-S19]. However, instead of solving the kinetic equation explicitly, we have made two simplifying assumptions about the distribution functionĥ. The first is that it is roughly constant throughout the superconductor, which is reasonable as long as the superconductor is not too thick compared to its spin relaxation length. The second is that the distribution function can be modelled using a spin voltage, which we justify in the Discussion in the main article. These assumptions imply that we can treatĥ as a constant parameter, which simplifies our model system from a 2D to 1D geometry, thus making it much more feasible to attack numerically.
In the numerical simulations, we also approximate the effect of inelastic scattering using a Dynes parameter: → + 0.01i∆ 0 . The ferromagnetic insulators in our model system were treated as spin-active interfaces, i.e. boundary conditions that account for spin-dependent phase shifts for quasiparticles in the superconductor that are reflected from the insulating interface [S21-S24]. This boundary condition takes the form where the plus and minus signs describe boundary conditions at the left and right interfaces of the superconductor, respectively. Here, L is the length of the superconductor along the x-axis, G is the bulk conductance of the superconductor in its non-superconducting state, G ϕ describes the effect of spindependent phase shifts when quasiparticles are reflected from a magnetic interface, and m is a unit vector that describes the magnetization direction at the same interface. In addition to the equations above, we need to solve a selfconsistency equation for the order parameter ∆(x) [S20], Here, f K s is the singlet part of the anomalous part of the Keldysh propagatorĝ K . In practice, it can be evaluated at all energies from the calculatedĝ K at positive energies using the identities f K s (x, − ) = 1 4 Tr (+iσ 2 ) (τ 1 + iτ 2 )ĝ K (x, + ) * .

(S92)
These follow from parametrizingĝ K in a similar manner as Eq. (S75), then invoking the definition of tilde-conjugatioñ f s (x, ) f * s (x, − ), and finally using standard trace identities. The actual numerical implementation was done using the Riccati parametrization of the retarded propagator [S20, S25], where γ is a 2 × 2 complex matrix, and N (1 − γγ) −1 . This matrix structure is defined in a way that automatically satisfies the normalization condition (ĝ R ) 2 = 1, and accounts for the particle-hole symmetries of the propagator, thus reducing the number of independent variables one has to solve for numerically. This parametrization also has the additional benefits that the Riccati parameter γ is single-valued and has a bounded norm, thus resulting in a more stable numerical solution procedure. How to reformulate the equations forĝ R in terms of γ is described in e.g. Ref. S20. In practice, we alternate between calculating γ from Eqs. (S88) and (S89), and updating ∆(x) using Eqs. (S90)-(S92), until they converge to a satisfactory degree. The simulation code itself is publicly available from github.com/jabirali/ .