Superconductor spintronics: Modeling spin and charge accumulation in out-of-equilibrium NS junctions subjected to Zeeman magnetic fields

We study the spin and charge accumulation in junctions between a superconductor and a ferromagnet or a normal metal in the presence of a Zeeman magnetic field, when the junction is taken out of equilibrium by applying a voltage bias. We write down the most general form for the spin and charge current in such junctions, taking into account all spin-resolved possible tunneling processes. We make use of these forms to calculate the spin accumulation in NS junctions subjected to a DC bias, and to an AC bias, sinusoidal or rectangular. We observe that in the limit of negligeable changes on the superconducting gap, the NS dynamical conductance is insensitive to spin imbalance. Therefore to probe the spin accumulation in the superconductor, one needs to separate the injection and detection point, i. e. the electrical spin detection must be non-local. We address also the effect of the spin accumulation induced in the normal leads by driving a spin current and its effects on the detection of the spin accumulation in the superconductor. Finally, we investigate the out-of-equilibrium spin susceptibility of the SC, and we show that it deviates drastically from it's equilibrium value.


I. INTRODUCTION
Chargeless spin accumulation has recently been realized experimentally, 1-3 enabling an estimate of the spinrelaxation time in mesoscopic superconductors. This is remarkably long, of the order of nanoseconds, significantly longer than the charge-relaxation time [4][5][6][7][8][9][10] . The measurement of this spin relaxation time using the injection of spin from a ferromagnet (FM) into a superconductor (SC) has recently been presented yielding an estimation of the order of 10ns. 1, 11 The SC was subject to a Zeeman field which splits its density of states (DOS) for the up and down quasiparticles. 12,13 To test that this spin relaxation time value is accurate, a new experiment has been developed in which a normal-superconducting junction is subject to both a constant (DC) and an time dependent (AC) bias. 14 The ferromagnetic lead is replaced by a normal one, as in Ref. 1 it was shown that spin accumulation was mostly due to the Zeeman split rather than to the ferromagnetic character of the lead. In the new experiments the SC is much cleaner, with a DOS much closer to the BCS form. The time-averaged spin chemical potential difference resulting from the spin accumulation is measured as a function of the applied chemical potential, the intensity of the AC signal, and as a function of frequency. 14 Here we present the model which allows us to fit the experimental data in Ref. 14 , extract the value of the spinrelaxation time, and explain the experimental dependence of the accumulated spin with respect to various parameters. We model the injection between the SC and the FM/metallic lead using Fermi's golden rule (FGR) and taking into account all possible tunneling processes when an electron from the lead is injected into the SC (the electron can enter the SC as a quasiparticle with the same spin, or as a quasihole of opposite spin). [15][16][17] Our model generalizes the theory presented in Ref. 1 which was an incomplete approximation, as some of the terms have not been taken into account, and which was valid only for small values of the spin imbalance. We subsequently write down and solve self-consistently the time-dependent semi-classical equations of motion for the spin accumulation in the SC. The relation between the spin accumulation in the SC and what is seen in a non-local measurement of this accumulation performed using a FM lead 14 is obtained by imposing the condition that the electric current through the non-local probe is zero. The bias voltage to which the probe is subjected to satisfy this condition is the measured non-local voltage that we can relate theoretically with the applied voltage.
A crucial observation to be made is that a frequency dependence for the time-averaged non-local measured voltage is conditioned by the existence of non-linearities in the system (such as the non-linearity of the conductance of the detection junction, and the non-linear dependence of the spin accumulation in the SC on the difference of chemical potentials between the spin up and down quasiparticles). Such non-linearities are automatic in a SC for not too small values of the spin accumulation due to the strong non-linear form of the BCS DOS. For a fully linear system the measured spin accumulation is independent of frequency, and the spin-relaxation time cannot be obtained from frequency-domain measurements.
The paper is organized as follows. In Sec. II we present the theory of spin injection into a SC, and we write down the relation between the injection electric and spin current as a function of the applied voltage. In Sec. III we present the experimental setup, and we write down the semi-classical equations of motion for the spin. In Sec. IV we solve these equations for a DC bias, and in Sec. V we study sinusoidal and rectangular AC biases. We conclude in section VI.

II. FM-SC JUNCTION
In this section we introduce the model to describe the FM-SC junction and we compute the tunneling currents flowing through it.

A. Theoretical model
The total Hamiltonian for the junction is : H = H F + H S + H T . The Hamiltonian for the FM lead is with q = 2 q 2 /2m and µ F the chemical potential in the FM. Actually, the "ferromagnetic" character of the lead will be encoded in the spin-dependent tunnel amplitude . The BCS Hamiltonian H S is given by with p = 2 p 2 /2m, ∆ is the superconducting energy gap and µ B H is the Zeeman energy. The tunneling Hamiltonian is The first term of H T describes the transfer of electrons with spin σ and momentum q into the superconductor with a spin dependent amplitude T σ p,q . We assume that there is no spin flip at the interface but that the momentum is not conserved, p → q. It is more convenient to express the operators for electrons in the SC in terms of quasiparticle operators using the Josephson's definition of the Bogoliubov-Valatin transformation 18,19 where the γ † e(h),p are creation operators of electronlike (holelike) excitations. Since a quasiparticle has probability u 2 p to be an electron and probability v 2 k to be a hole, the quasiparticle charge is given by q p = u 2 p − v 2 p . The corresponding charge carried by the condensate is 2v 2 p = 1 − q p . 15 With the proper choice of the coherence factors u p and v p , the BCS Hamiltonian can be written as where , and E p,σ = 2 p + ∆ 2 − σµ B H the excitation energy. 12,13 Introducingσ = −σ, the tunneling Hamiltonian becomes B. Tunnel current We now calculate the particle tunneling current flowing through the junction by using Fermi's Golden rule. 15,20 The allowed tunneling processes and the corresponding probabilities are described in Table I. For example, for a γ † e,pσ c q,σ process (an electron c q,σ is annihilated in the SC to create a quasiparticle γ † e,pσ in the SC), the total probability is a product of the tunnel probability given by the tunnel Hamiltonian, |T σ p,q | 2 u 2 p , the probability to have a filled state in the FM to tunnel from, f ( q ), and the probability to have an empty state in the SC to tunnel into, 1 − f pσ (E p,σ ).
The tunneling average current through the junction for a given spin can be written as where we assumed T σ p,q = T σ q,p . To simplify this formula we note that for each state with p + > k F , energy E p + ,σ and u p + , there exists another state p − < k F with the same energy E p − ,σ = E p + ,σ (note however that p + = − p − ). This implies for the coherence factors that u 2 p ± = v 2 p ∓ . Moreover we can reasonably assume that |T p + ,q | = |T p − ,q |. By separating the sum over p into two sums over p ± , and by noting that only the coherence factors depend on p ± (with u p ± = v p ∓ , u 2 p ± + v 2 p ± = 1) we obtain where we assumed that there is no charge accumulation, meaning that the Fermi-Dirac distribution in the superconductor is identical for momenta greater and lesser than the Fermi momentum k F . In the following the sign + of the momentum p will be omitted for brevity. We can now convert the momentum summation into an energy integral, using For the energy range we consider, it is reasonable to assume that the FM DOS and the tunneling probabilities are roughly independent of energy. Performing the momentum-energy conversion for the q in the FM, and subsequently the resulting energy integral, we obtain Here ρ F is the total density of states of the ferromagnetic lead (integrated over the volume of the FM). The conversion of the summation over the momentum p in the SC into an energy integral is more tricky. This is because the two first terms of the above expression correspond to the injection of an electron as an electron-like excitation of energy E p,σ . The last two terms correspond to the conversion of an electron into a hole-like excitation at energy −E p,σ . The SC densities of states are different for the two processes, due to the presence of the Zeeman field: ρ S (E p,σ ) = ρ S (−E p,σ ). Converting the momentum summation over p into an energy integral thus leads to  where E σ = E − σµ B H, ρ 0 is the DOS of the superconductor at the Fermi energy, and is the usual normalized BCS density of states, with θ the Heaviside step function. In real materials this DOS needs to be modified to round-off singularities and to take into account the existence of states inside the gap. 21 Such states can be induced for instance by disorder (static impurities, magnetic impurities, etc.). Several theoretical models which describe the presence of such gap states can be used in order to fit properly the experimentally measured DOS. The most commonly used is the Dynes model for the DOS, 22 in which a phenomenological parameter δ is added to the energy as an imaginary part corresponding to a quasiparticle lifetime We can also use the Fulde DOS which also takes into account the presence of spin orbit effects, disorder, etc., in conjunction with an applied magnetic field. 23 It is however very difficult to properly fit the experimentally measured DOS using this models. Sometimes is most useful to use in the numerical calculations directly the DOS measured experimentally. The total spin and electric tunneling current can be written as I e = e σ I σ and I s = /2 σ σI σ where Here µ s refers to the shift of the chemical potential due to the spin imbalance.
It is important to note the difference between the above formulas obtained taking into account all possible tunneling processes and those obtained from the semiconductor model for a superconductor. 17 In the semiconductor model an electron injected into the superconductor enters as a quasiparticle with the same spin, and consequently access solely one spin density of states. This gives rise to terms in the tunneling current containing solely products between samespin densities of state in the electrode and in the SC. On the other hand, more generally, a quasiparticle injected into a superconductor can enter either as a quasiparticle of the same spin or as a quasihole of opposite spin. Thus the injection has access to both spin densities of states, which gives rise to terms in the tunneling currents containing both products of same-spin DOS in the electrode and in the SC, as well as of opposite-spin DOS, as it can be clearly seen in the above equation. The two models are of course different only if the DOS of the spin up and down electrons are different, as it is the case here due to the presence of an applied Zeeman magnetic field; in this case the semiconductor model does not suffice to capture the physics of the injection into the SC.

A. Experimental setup
The setup that we consider is presented in Fig. 1. A voltage bias V is applied between a lead and the SC at point A (injection junction). The spin accumulation S(t) in the SC is measured at point B (detection junction) by applying a voltage difference V m between a ferromagnetic lead and the SC and imposing the condition that the electrical current flowing through the detection junction is zero.

B. Semiclassical equations of motion
We assume that the time dependent spin accumulation S(t) in the superconductor satisfies a simple diffusion equation of motion where τ s is the spin relaxation time in the SC. As described in the previous section, I i s , the spin current in the injection junction, is a function of the applied voltage V between the SC and the injection lead. In much of our analysis, as The experimental setup -schematics (left), and actual setup (right) with a scale bar of 1µm. Connection between the two figures: S is the superconducting finger, J3 is the detection junction with F corresponding to F 2. J1 and J2 are two injection junctions between N 1/N 2 (corresponding to F1) and S, situated at different distances from J3. 14 well as in the experiments presented in Ref. 14, the injection lead is considered to be non-FM (i.e. T i ↓ = T i ↑ = T i ) and the injected spin current can be written as The spin accumulation in the superconductor can be written as This relates the accumulated spin to the chemical potential difference between the up and down spins in the SC. The equation of motion can be formally integrated Note that I i s is also a function of µ s . The equations (17), (18), and (19) form a self-consistent system of integral equations which can be solved numerically to determine µ S as a function of V for all times t. An exception to this can be made for very small values of µ s when we can neglect the dependence of I i s on µ s , and we calculate the accumulated spin and correspondingly of µ s directly from Eq. (19), as done in Ref. 1.
When the applied voltage is constant in time the system will end up in a stationary state in which the injected spin current will exactly compensate the spin relaxation in the superconductor such that dS(t)/dt = 0. Imposing this condition in the equations of motion yields I i s = S/τ s , where now both S and I i s are time-independent. Equations (17) and (18) form a self-consistent system of equations which can be solved numerically to determine the stationary value of µ S as a function of V . To determine the measured voltage V m as a function of µ s , and correspondingly on the value of the applied voltage, V , we need to impose the extra condition that in the detection junction the total electrical current I d e = 0, where I d e , as described in the previous section, is given by Note that the currents through the injector and the detector satisfy the same equations but with different parameters.
Here, for the FM detector the transmission coefficients are spin-dependent T d If this is not the case, V m = 0 for all accumulated µ s . Conversely, the measured V m value depends strongly on the polarization of the detector , with a typical experimental P d ≈ 0.1 for the Cobalt lead.

IV. RESULTS FOR AN APPLIED DC VOLTAGE
In Fig. 2 we present a typical dependence of V m on V . Similar to Ref. 1 this follows qualitatively the form of I s (V ) which exhibits the same main features as the BCS DOS 1 (e.g. two coherence peaks, a quasi-null value at small V 's and a saturation at large V 's). To understand the difference between our approach and previous approximations we will separately take into account the effects of two important factors, the non-linearity and the self-consistence. The non-linearity comes into place when B and µ s are large, and the above current formulas cannot be Taylor expanded in these parameters. The self-consistence expresses the back-reaction of the accumulated µ s on the injection current; when µ s is small with respect to the applied voltage, we can neglect it in Eq. (17) and solve the equations of motion in a non-self-consistent matter.

The self-consistence :
We begin by analyzing the effect of the self-consistence. In Fig. 2 we present the dependence of V m on V obtained both self-consistently (full line) and non-self-consistently (dashed line). It appears that the main difference is quantitative, i.e. the self-consistence is introducing an overall correcting factor which does not depend strongly on V . To test this assumption we write down the self-consistent and non-self-consistent solutions of the equations of motion in the linear limit (small µ s ). In this regime we can perform a Taylor expansion of Eq. (18) in µ S , S = ρ 0 ρ(µ B H)µ nsc s , which together with the condition I i s = S/τ s yields We can also solve the equations of motion self-consistently by making a Taylor expansion of Eq. (17) in µ s with g ns = (2πe 2 / )ρ F |T i | 2 ρ 0 is the normalized conductance of the injection junction. Noting that S = 2 ρ 0 ρ(µ B H)µ sc s , we find and µ sc s /µ nsc s = 1/(1 − 2τ s g ns /e 2 ρ 0 ). In our numerical calculations we will set g ns /e 2 ρ 0 = 1, so that τ s is measured in units of e 2 ρ 0 /g ns . Indeed, it seems that in the linear limit, the non-self-consistent and self-consistent approaches differ by a simple numerical factor, which converges to 1 when τ s 1. This observations has been checked numerically in Fig. 3 where we have plotted the accumulated µ s calculated using the non-self-consistent and the self-consistent approach (with a correcting factor of 1/(1 − 2τ s ) taken into account). Indeed we see that the two give the same result in the linear (small V ) regime.
It appears thus that the effect of solving the equations of motion in a self-consistent or non-self-consistent manner is mainly quantitative (an overall numerical factor), which is however very important if we are interested in extracting the value of the spin relaxation time from a fit of the experimental data for an applied DC voltage. A factor of 10 difference arising in the model from the self-consistency can give thus a factor of 10 error in the estimation of the spin-relaxation time.
The non-linearity of the detection : The second important factor to be understood is the non-linearity, in particular the non-linearity of the detection. We note that the chemical potential describing the spin accumulation µ s depends qualitatively similar on the applied V as the measured V m . The differences come from the non-linearity of the detection junction. To understand this, we have plotted in Fig. 4  of the detection junction, thus for V ≤ ∆/2, corresponding to a small µ s , we have a linear dependence of V m on µ s , V m ∝ P d µ s as expected, while for V > ∆/2 the linearity does not hold (for reminder: ∆ = 0.22meV and P d = 2%). We should also note that, as we will show in the next section, the non-linearities in the system are a crucial ingredient in observing a frequency dependence of the measured non-local signal, i.e V m depends on the frequency of the applied AC voltage only because of the non-linear form of the SC DOS.

V. RESULTS FOR AN APPLIED AC VOLTAGE
A. Time-dependent behavior

Numerical results
We now apply a time-dependent (AC) sinusoidal voltage of frequency ω and amplitude V rf on the injector, such that V (t) = V + V rf cos(ωt). The expressions for the spin currents do not change (see Eq. (17)). The only change comes from the dependence on time of the spin imbalance. We can solve numerically the self-consistent integral equations of motion (Eqs. (16)- (20)) to obtain the time-dependent V m (t), µ s (t) and accumulated spin S(t), as a function of applied V for various values of frequencies and AC amplitudes.
We begin by plotting the accumulated spin S(t) as a function of time. The time-dependence of the spin accumulation can be understood easily by thinking that the superconductor acts as a capacitor (its charge could be viewed as the spin imbalance). Indeed, for large frequencies the capacitor is loading (the spins are accumulating up to a maximal value) but its decreasing never happens because the spin relaxation time is larger than the period of the oscillations. On the contrary, for small frequencies voltage the spins can relax because of the large period of the oscillations. In Fig. 5 we plot S(t) for three values of frequency, ω = 0.8/τ s , ω = 0.2/τ s and ω = 0.04/τ s . All other parameters are the same as in the previous section. We note that the average of the oscillations is independent of frequency, while their amplitude is not. The larger the frequency, the more the behavior of S(t) approaches that of a charging capacitor with smaller and smaller oscillations around the saturation value. It would thus seem that experimentally one cannot see a frequency dependence for the time-averaged spin accumulation. However, in an actual experiment one does not measure S but V m , which can exhibit a strong non-linear behavior with S. We should then expect that if the amplitude of the oscillations in S(t) depends on frequency the time average of V m depends on frequency via rectification effects. In Fig. 6 we plot the time dependence of V m for three different frequencies and we see that indeed both the amplitude of the oscillations and the time average depend on frequency.

Taylor expansion
To understand the above numerical results we study a few limiting cases that can be solved analytically. For a small applied AC voltage (V rf V ), with V (t) = V + V rf cos(ωt), we can use a Taylor expansion, and the spin current can be expressed as Injecting Eq. (24) in Eq. (19) gives us an expansion for the spin imbalance in powers of V rf : S(t) = S 0 (t) + S 1 (t) + S 2 (t) + .... We focus on the first two terms in the expansion, but the next orders can be studied in a similar fashion. For times much more larger than τ s we obtain (see Appendix A).
We see from the above formula that the average accumulated spin is independent of frequency, consistent with the numerical analysis in the previous section (see Fig. 5). However, the amplitude of the oscillations does depend on frequency. By having a quick look at the above expression we can directly point out two regimes. Thus, for ωτ s 1 we should expect a constant behavior for the dependence of the amplitude of the oscillations with the frequency, as τs 1+τ 2 s ω 2 −→ τ s . The opposite regime where ωτ s 1 gives us a 1/ω dependence. This behavior is confirmed by a numerical analysis of the dependence of the amplitude of oscillations with the frequency where we indeed note a crossover behavior for a frequency ω ≈ 1/τ s . However this analysis in valid only when V rf is very small, and in the regime in which the system is well described by the non-self-consistent calculation. Also, since the timescales involved are very short, it is much harder to have experimentally access to the amplitude of the oscillations than to the time averages, and in what follows we will focus rather on time-averaged quantities than on time-dependent ones. We expect that the non-linearity will give rise to a frequency dependence even when averaging over time, allowing us to detect directly this spin relaxation time in the frequency domain. Experimentally, one uses 0.5 − 50M Hz frequencies since we expect a spin relaxation time of the order of nanoseconds.

B. Time-averaged quantities
In this section we study the dependence of the average measured V m as a function of the applied voltage for different AC amplitudes and frequencies. We begin by plotting V m and its derivative dV m /dV as a function of V for different frequencies at fixed AC amplitude (see Fig. 7). All frequencies are given in units of 1/τ s . Subsequently, in Fig. 8   i.e a flattening of the V m dependence on V , with an eventual extra peak arising at V = 0, a doubling of the peaks in the dV m /dV dependence on V , whose position depend quasi-linearly on V rf , and a saturation of V m and dV m /dV with increasing the frequency, are qualitatively similar to what is measured in Ref. 14. The frequency at which the saturation occurs is directly related to the inverse of the spin-relaxation time. Thus our theory seems to describe very well the frequency and amplitude dependence of the measured accumulation. We have checked that while the measured V m and µ s do depend on the frequency (because of the non-linearities in the system), the accumulated spin S does not, as described also in the previous section. We have calculated the accumulated spin S (see Fig. 11) for different values of V and frequency and we have seen that S is indeed unaffected by the frequency.

C. Different type of AC voltage: rectangular pulses
To get an analytical understanding of the numerical results presented in the previous section we consider also a different type of AC signal, for example a chain of rectangular pulses. In this case we can also calculate analytically the form of the spin imbalance, if we make the assumption that the self-consistent effects are negligible. The pulse has the following shape with T = 2π/ω the period of the signal, V rf its amplitude, and the width of the pulse. The difference of spin imbalance between the stationary regime (V rf = 0) and the time-dependent one (V rf = 0) can be calculated exactly using Eq. (19): Performing the integral over time (see Appendix B) leads to The average of Eq. (28) can be written as δS = 1/T T 0 δS(t)dt and leads to If the pulse constitutes a constant fraction of the period (ω is constant), the accumulated average S should be independent of frequency. We have checked that is indeed the case by a numerical analysis. Also this is consistent with our previous observations for the sinusoidal signal. An interesting observation that we make is that, as it can be seen from Eq. (29) the dependence of δS on V is given generically by I i s (V + V rf ) − I i S (V ). In Fig. 12 we plot the dependence of the excess accumulated spin δS as a function of V obtained numerically for a specific value of V rf and of frequency. We also sketch the behavior of showing that indeed, to first approximation, the behavior of δS follows qualitatively the behavior of I i s (V + V rf ) − I i S (V ). A simple generalization can be made to understand qualitatively the behavior of δS with V for the sinusoidal signal. To first approximation a sinusoidal signal is equivalent to a superposition of two V rf and −V rf pulses, with an ω = 1/2. We would then expect an overall dependence of δS qualitatively similar with Fig. 13 we plot the dependence of the excess accumulated spin as a function of V obtained numerically for a specific value of V rf = 0.2meV and ω = 0.04/τ s for a sinusoidal signal, and we also sketch the behavior of , as a function of V for a sinusoidal signal with V rf = 0.2meV and the frequency ω = 0.04/τs. Note the qualitatively similar behavior of the two curves (same as in the previous figure, a constant has been introduced to uniformize the two curves).

Differences between accumulation and relaxation times
In this section we consider the possibility that the time for accumulation (loading) and relaxation (unloading) are different. Such phenomenon could be detected by applying a time dependent voltage with the following shape (see This is because if A 1 > 0 and A 2 < 0 the first step corresponds to the loading of the superconductor and is controlled by τ 1 , and the second one to the unloading and is controlled by τ 2 . A similar calculation as before can be performed leading to the following form for the average spin imbalance By setting 1 = 2 = , A 1 = A 2 = A and τ 1 = τ 2 = τ we effectively restore the previous situation with τ s being replaced by 2τ = τ 1 + τ 2 . We can see that in this limit the dependence is still linear with frequency. One important observation to make is that the difference between τ 1 and τ 2 can be measured directly by applying an AC voltage with A 1 = −A 2 and 1 = 2 . If the two times are different, the average excess accumulation will be non-zero, which is not the case when τ 1 = τ 2 .

VI. CONCLUSIONS
We have calculated the accumulated spin in a junction between a superconductor and a normal or ferromagnetic contact, taken out of equilibrium, and in the presence of a Zeeman magnetic field. We have studied the effects of DC and AC applied voltages, for different amplitudes and frequencies. We have found that for an applied DC voltage, the dependence of the measured spin accumulation V m on the applied voltage V follows roughly the form of the DOS, i.e. two coherence peaks, a reduction at small biases, and a saturation for large applied voltages. When an AC component is added to the applied voltage one observes a flattening of the measured spin voltage with the applied V , with an eventual extra peak arising at V = 0, a doubling of the peaks in the dV m /dV dependence on V , whose position depends quasi-linearly with the amplitude of the AC voltage, and a saturation of V m and dV m /dV with increasing the frequency. The frequency at which the saturation occurs is directly related to the inverse of the spin-relaxation time. Our theoretical predictions show a very good qualitative agreement with the experimental measurements in Ref. 14. A detailed comparison between our theoretical analysis and these experimental measurements allows the extraction of the spin relaxation time in the SC, in particular from the dependence of the measured accumulation as a function of frequency. We find that this value is of the order of ns.
Note the dependence on frequency of the prefactor τs 1+τ 2 s ω 2 , corresponding to a frequency dependence for the amplitude of the oscillations in S(t).
Appendix B -Analytical calculation of S(t) for a rectangular pulse.