Tunneling Conductance and Spin Transport in Clean Ferromagnet-Ferromagnet-Superconductor Heterostructures

# Tunneling Conductance and Spin Transport in Clean Ferromagnet-Ferromagnet-Superconductor Heterostructures

Chien-Te Wu    Oriol T. Valls Also at Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455 School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455    Klaus Halterman Michelson Lab, Physics Division, Naval Air Warfare Center, China Lake, California 93555
July 17, 2019
###### Abstract

We present a transfer matrix approach that combines the Blonder-Tinkham-Klapwijk (BTK) formalism and self-consistent solutions to the Bogolibuov-de Gennes (BdG) equations and use it to study the tunneling conductance and spin transport in ferromagnet ()-superconductor () trilayers () as functions of bias voltage. The self-consistency ensures that the spin and charge conservation laws are properly satisfied. We consider forward and angularly averaged conductances over a broad range of the strength of the exchange fields and thicknesses, as the relative in-plane magnetization angle, , between the two ferromagnets varies. The -dependence of the self-consistent conductance curves in the trilayers can differ substantially from that obtained via a non-self-consistent approach. The zero bias forward conductance peak exhibits, as varies, resonance effects intricately associated with particular combinations of the geometrical and material parameters. We find, when the magnetizations are non-collinear, signatures of the anomalous Andreev reflections in the subgap regions of the angularly averaged conductances. When is half-metallic, the angularly averaged subgap conductance chiefly arises from anomalous Andreev reflection. The in-plane components of the spin current are strongly bias dependent, while the out-of-plane spin current component is only weakly dependent upon voltage. The components of the spin current aligned with the local exchange field of one of the F layers are conserved in that layer and in the S region, while they oscillate in the other layer. We compute the spin transfer torques, in connection with the oscillatory behavior of spin currents, and verify that the spin continuity equation is strictly obeyed in our method.

###### pacs:
74.45.+c,74.78.Fk,75.75.-c

## I Introduction

Over the last two decades, significant progress in fabrication techniques has allowed the development of spintronics devices, such as spin valves,igor () that utilize both charge and spin degrees of freedom. Traditional spin valves consist of magnetic materials only. There is another important type of spintronics devices, involving ferromagnet (F)-superconductor (S) heterostructures. These heterostructures have also received much attention because of the fundamental physics related to the interplay between ferromagnetic and superconducting order. Their potential applications in spintronics include magnetic memory technology where information storage is accomplished via control of the magnetic moment bit. It is then crucial to have precise control over the magnetization direction. Spin transfer torque (STT) is one effect that affords such control. The generation of spin-polarized supercurrents may be used to obtain a superconducting STT acting on the magnetization of a ferromagnet. This effect may be utilized in high density nanotechnologies that require magnetic tunnel junctions. Thus, the dissipationless nature of the supercurrent flow offers a promising avenue in terms of low energy nanoscale manipulation of superconducting and magnetic orderings.

Although ferromagnetism and -wave superconductivity seem incompatible because of the inherently opposite natures of their order parameter spin configurations, superconductivity can still be induced in the F layers of F-S layered structures by the superconducting proximity effects.Buzdin2005 () In essence, the superconducting proximity effects describe the leakage of superconductivity into a non-superconducting normal (N) or magnetic metal, as well as its depletion in S near the interface. However, proximity effects in F-S systems are very different from those in N-S structures due to the inherent exchange field in the F materials. As a consequence of this exchange field, the Cooper pair acquires a non-zero center-of-mass momentumdemler (); Buzdin2005 (); Halterman2001 (); Halterman2002 () and the overall Cooper pair wavefunction oscillates spatially in the F regions. Owing to this oscillatory nature, many new physical phenomena emerge in F-S heterostructures such as oscillations of the superconducting transition temperature, , with the thickness of the F layers. Buzdin1990 (); Buzdin2005 (); demler (); Halterman2004 ()

It is of fundamental importance that superconducting proximity effects are governed by Andreev reflection,Andreev () which is a process of electron-to-hole conversion at N-S or F-S interfaces, and it involves the creation or annihilation of a Cooper pair. Therefore, consideration of Andreev reflection is central when studying the transport properties of N-Sbtk (); tanaka () or F-S systems.beenakker (); zv1 (); zv2 () Of particular interestbtk (); tanaka (); beenakker (); zv1 (); zv2 () is the behavior of the tunneling conductance in the subgap region, where hybrid systems can carry a supercurrent due to Andreev reflection. In conventional Andreev reflection, the reflected hole has opposite spin to the incident particle. Accordingly, the exchange field in the F materials that causes the splitting of spin bands has a significant effect on the tunneling conductance in the subgap region. Most important, the qualitative behavior of the conductance peak in the zero bias limit is strongly influenced by the degree of conduction electron spin polarization in the F materials.beenakker (); zv1 (); zv2 (); mazin () Experimentally, this concept has been applied to quantify the spin polarization. raychaudhuri (); upad (); chalsani (); soulen (); hacohen ()

An intriguing phenomenon in F-S structures is the induction of triplet pairing correlations.berg86 (); Bergeret2007 (); Wang2010 (); Hubler2012 (); hv2p () These correlations are very important when studying transport phenomena such as those found in SFS Josephson junctions.Keizer2006 (); robinson (); khaire () In contrast to the short proximity lengthHalterman2002 () of singlet Cooper pair condensates into F materials, the triplet pairing correlations are compatible with the exchange fields and hence largely immune to the pair breaking effect produced by the latter. However, for such correlations to be induced F-S structures must possess a spin-flip mechanism. Examples include a spin-dependent scattering potential at the F-S interface Halterman2009 (); Eschrig2008 () and the introduction of another magnetic layer with a misoriented magnetic moment such as superconducting spin valves.Halterman2007 () The pairing state of induced triplet correlations is at variance with the effects of conventional Andreev reflection, responsible for the generation of singlet Cooper pairs. Thus, recent studieslinder2009 (); visani (); niu (); ji (); feng () on the tunneling conductance propose the existence of anomalous Andreev reflection, that is, a reflected hole with the same spin as the incident particle can be Andreev reflected under the same circumstances as the generation of triplet pairing correlations becomes possible. In this view, triplet proximity effects are correlated with the process of this anomalous Andreev reflection. This will be confirmed and discussed in this work.

Another important geometry for a superconducting spin valve consists of a conventional spin valve with a superconductor layer on top: a trilayer. By applying an external magnetic field, or switching via STT, one is able to control the relative orientation of the intrinsic magnetic moments and investigate the dependencealejandro (); leksin (); zdravkov () of physical properties such as on the misorientation angle between the two magnetic layers. Due to the proximity effects, is often found to be minimized when the magnetizations are approximately perpendicular to each other,cko () reflecting the presence of long range triplet correlations, induced in trilayers. Their existence has been verified both theoreticallycko () and experimentally.alejandro () The non-monotonic behavior of as a function of has also been shown to be quantitativelyalejandro () related to the long range triplet correlations, with excellent agreement between theory and experiment.

Motivated by these important findings, we will investigate here, in a fully self-consistent manner, the dependence of the tunneling conductance and other transport quantities of these trilayers. Non-self-consistent theoretical studies of tunneling conductance have been performed on trilayers in previous work.ji (); cheng () However, as we shall see in Sec. II, only self-consistent methods guarantee that conservation laws are not violated and (see Sec. III) only then can one correctly predict the proximity effects on the angular dependence of transport properties. The spin-polarized tunneling conductance of F-S bilayers only, was studied in Refs. zv1, ; zv2, ; kashiwaya, ; yamashita, . Also, in traditional spin valves e.g. - layered structures, the spin-polarized current generated in the layer can transfer angular momentum to the layer when their magnetic moments are not parallel to each otherigor () via the effect of STT.berger (); myers () As a result, the spin current is not a conserved quantity and one needs a general law that relates local spin current to local STT.linder2009 () The transport properties of structures, in particular the dependence on applied bias of the spin-transfer torque and the spin-polarized tunneling conductance have been previously studied.romeo (); bozovic (); linder2009 ()

Here, we consider charge transport and both spin current and spin-transfer torque in trilayers. In previous theoretical work, such as that mentioned above, when computing tunneling conductance of N-S and F-S structures, using methods based on the Blonder-Tinkham-Klapwijk (BTK) procedurebtk (); tanaka (); zv1 (); zv2 (); bozovic (); linder2007 (); linder2009 () and quasi-classical approximations,Eschrig2010 () the superconducting pair amplitude was assumed to be a step function: a constant in S, dropping abruptly to zero at the N-S or F-S interface and then vanishing in the non-superconducting region. This assumption neglects proximity effects. Only qualitative predictions on the behavior of the tunneling conductance can be reliably made. Still, results exhibit many interesting features especially in F-S systems.zv1 (); zv2 () However, to fully account for the proximity effects, in the transport properties, one must use a self-consistent pair potential. This is because that reveals realistic information regarding the leakage and depletion of superconductivity. Also, as we shall discuss below, self-consistent solutions guarantee that conservation laws are satisfied. In Ref. bv, , the tunneling conductance of F-S bilayers was extracted via self-consistent solutions of Bogoliubov-de Gennes (BdG) equations.degennes () However, the numerical methods used there required awkward fitting procedures that led to appreciable uncertainties and precluded their application to trilayers. The findings indicated that the self-consistent tunneling conductances for the bilayer are quantitatively different from those computed in a non-self-consistent framework, thus demonstrating the importance of properly accounting for proximity effects in that situation. Here we report on a powerful self-consistent approach and use it to compute the tunneling conductance of trilayers. It is based on the BTK method, incorporated into a transfer matrix procedure similar to that usedstrinati () in Josephson junction calculations and simple F-S junctions within a Hubbard modelting (). As we shall demonstrate, this approach not only has the advantage of being more numerically efficient but also can be used to compute spin transport quantities. Thus, we are able to address many important points regarding both charge and spin transport in trilayers, including the spin currents and spin-transfer torque, the proximity effects on the tunneling conductance, and the correlation between the anomalous Andreev reflection and the triplet correlations.

This paper is organized as follows: we present our self-consistent approach, and its application to compute the tunneling conductance, the spin-transfer torques, the spin current, and the proper way to ensure that conservation laws are satisfied, in Sec. II. In Sec. III we present the results. In Subsec. III.1, we briefly compare the results of F-S bilayers obtained in our self-consistent approach with non-self-consistent ones. The rest of Subsec. III.2 includes our results for trilayers, that is, the main results of this work. The dependence on the tunneling conductance of trilayers on the angle is extensively discussed as a function of geometrical and material parameters. Results for the effect of the anomalous Andreev reflection, the spin-transfer torque, and the spin current are also presented. We conclude with a recapitulation in Sec. IV.

## Ii Methods

### ii.1 Description of the system

The geometry of our system is depicted in Fig. 1. We denote the outer ferromagnet as and the middle layer as . We choose our coordinate system so that the interfaces are parallel to the plane, and infinite in extent, while the system has a finite width in the direction.

The Hamiltonian appropriate to our system is,

 Heff = ∫d3r{∑αψ†α(r)H0ψα(r) (1) + −

where is the single-particle Hamiltonian, is a Stoner exchange field that characterizes the magnetism, and are Pauli matrices. The superconducting pair potential is the product of pairing constant, , in the singlet channel, and the pair amplitude. We begin by writing down the BdG equations, which we will solve self-consistently for our FFS trilayers. By performing the generalized Bogoliubov transformationdegennes (), , where and for spin-down (up), the Hamiltonian [Eq. (1)] can be diagonalized. We can then for our geometry rewritecko () Eq. (1) as a quasi-one-dimensional eigensystem:

 ⎛⎜ ⎜ ⎜ ⎜⎝H0−hz−hx0Δ−hxH0+hzΔ00Δ−(H0−hz)−hxΔ0−hx−(H0+hz)⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜⎝un↑un↓vn↑vn↓⎞⎟ ⎟ ⎟ ⎟⎠ =ϵn⎛⎜ ⎜ ⎜ ⎜⎝un↑un↓vn↑vn↓⎞⎟ ⎟ ⎟ ⎟⎠, (2)

where the and are respectively the quasiparticle and quasihole amplitudes with spin . The exchange field vanishes in the S region, while in it is directed along , , and in it can rotate in the plane, . The single-particle Hamiltonian now readscko () , where denotes the transverse kinetic energy in the plane. Also, in the superconducting region and in the ferromagnetic layers. Throughout this paper, we assume natural units and measure all energies in units of . To take into account the more realistic situation where the F materials can in general have different bandwidths than the S layer, we define (as in Ref. bv, ) a mismatch parameter via .

We are aiming here to solve the problem in a fully self consistent manner. The self-consistent pair potential can be expressed in terms of the quasi-particle and quasi-hole wavefunctions. Accordingly,

 Δ(y)=g(y)2∑n′[un↑(y)v∗n↓(y)+un↓(y)v∗n↑(y)]tanh(ϵn2T), (3)

where the primed sum is over all eigenstates with energies smaller than a characteristic Debye energy, and is the superconducting coupling constant in the region and vanishes elsewhere. We obtain the self-consistent pair potential by solving Eqs. (2) and (3) following the iterative numerical procedures discussed in previous work.hv2p (); cko ()

### ii.2 Application of the BTK method

The BTK formalism is a procedure to extract the transmitted and reflected amplitudes, and hence the conductance, from solutions to the BdG equations. This is accomplished by writing down the appropriate eigenfunctions in different regions. In this subsection, we review the relevant aspects of the formalismbtk () for the non-self-consistent case (a step function pair potential) with the objective of establishing notation and methodology to describe, in the next subsection, the procedure to be used in the self-consistent case.

Consider first a spin-up quasi-particle with energy , incident into the left side labeled “”, in Fig. 1). Since the exchange fields in the and layers can be non-collinear, it follows from Eq. (2) that the spin-up (-down) quasi-particle wavefunction is not just coupled to the spin-down (-up) quasi-hole wavefunction, as is the case of F-S bilayers. Indeed, the wavefunction in the layer is a linear combination of the original incident spin-up quasi-particle wavefunctions and various types of reflected wavefunctions, namely reflected spin-up and spin-down quasi-particle and quasi-hole wavefunctions (via both ordinary and Andreev reflections). We use a single column vector notation to represent these combinations,

 ΨF1,↑≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝eik+↑1y+b↑e−ik+↑1yb↓e−ik+↓1ya↑eik−↑1ya↓eik−↓1y⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (4)

If the incident particle has spin down, the corresponding wavefunction in is

 ΨF1,↓≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝b↑e−ik+↑1yeik+↓1y+b↓e−ik+↓1ya↑eik−↑1ya↓eik−↓1y⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (5)

In these expressions are quasi-particle and quasi-hole wavevectors in the longitudinal direction , and satisfy the relation,

 k±σm=[Λ(1−ησhm)±ϵ−k2⊥]1/2, (6)

where (as used above) or , used later. As mentioned above, all energies are in units of and, in addition, we measure all momenta in units of . In this simple case, one can easily distinguish the physical meaning of each individual wavefunction. For instance in Eq. (4), is the reflected spin-down quasi-hole wavefunction. The quasi-hole wavefunctions are the time reversed solutions of the BdG equations and carry a positive sign in the exponent for a left-going wavefunction. The relevant angles can be easily found in terms of wavevector components. Thus, e.g., the incident angle (for spin-up) at the interface is , and the Andreev reflected angle for reflected spin-down quasi-hole wavefunction is . The conservation of transverse momentum leads to many important featuresbv (); zv1 () when one evaluates the angularly averaged tunneling conductance, as we will see below. For the intermediate layer , the eigenfunction in general contains both left- and right-moving plane waves, that is,

 ΨF2≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝c1f+↑eik+↑2y+c2f+↑e−ik+↑2y+c3g+↑eik+↓2y+c4g+↑e−ik+↓2yc1f+↓eik+↑2y+c2f+↓e−ik+↑2y+c3g+↓eik+↓2y+c4g+↓e−ik+↓2yc5f−↑eik−↑2y+c6f−↑e−ik−↑2y+c7g−↑eik−↓2y+c8g−↑e−ik−↓2yc5f−↓eik−↑2y+c6f−↓e−ik−↑2y+c7g−↓eik−↓2y+c8g−↓e−ik−↓2y⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (7)

where and are defined in Eq. (6). The indices are defined as previously, and the up and down arrows refer to . The eigenspinors and that correspond to spin parallel or antiparallel to respectively, are given, for , by the expression,

 (f+↑f+↓)=1N(11−cosϕsinϕ)=(f−↑−f−↓);(g+↑g+↓)=1N(−sinϕ1+cosϕ1)=(−g−↑g−↓) (8)

with the normalization constant . These spinors reduce to those for pure spin-up and spin-down quasi-particles and holes when , corresponding to a uniform magnetization along . One can also easily see that the particular wavefunction of Eq. (7), denotes a quasi-particle with spin parallel to the exchange field in . When , these eigenspinors read

 (f+↑f+↓)=1N(sinϕ1−cosϕ1)=(−f−↑f−↓);(g+↑g+↓)=1N(1−1+cosϕsinϕ)=(g−↑−g−↓) (9)

with .

In this subsection where we are still assuming a non-self-consistent stepwise potential equal to throughout the S region and to zero elsewhere, we have the superconducting coherence factors, and . In this case the right-going eigenfunctions on the S side can be written as,

 ΨS≡⎛⎜ ⎜ ⎜ ⎜ ⎜⎝t1u0eik+y+t4v0e−ik−yt2u0eik+y+t3v0e−ik−yt2v0eik+y+t3u0e−ik−yt1v0eik+y+t4u0e−ik−y⎞⎟ ⎟ ⎟ ⎟ ⎟⎠, (10)

where, are quasi-particle (+) and quasi-hole (-) wavevectors in the S region. By using continuity of the four-component wavefunctions and their first derivatives at both interfaces, one can obtain all sixteen unknown coefficients in the above expressions for the wavefunctions by solving a set of linear equations of the form at the interface and at the interface simultaneously, where

 xTF1,↑=(1,b↑,0,b↓,0,a↑,0,a↓) (11a) xTF1,↓=(0,b↑,1,b↓,0,a↑,0,a↓) (11b) xTF2=(c1,c2,c3,c4,c5,c6,c7,c8) (11c) xTS=(t1,0,t2,0,t3,0,t4,0), (11d)

and , , , and are appropriate matrices, which are straightforward to write down. Use of these coefficients gives us all the reflected and transmitted amplitudes and which are used to compute the conductance, as discussed in the next two subsections.

### ii.3 Transfer matrix self consistent method

The non-self-consistent step potential assumption is largely unrealistic. Proximity effects lead to a complicated oscillatory behavior of the superconducting order parameter in the F layers and to the generationBuzdin2005 (); Keizer2006 (); Bergeret2007 (); Wang2010 (); visani (); Hubler2012 (); Halterman2007 (); hv2p () of triplet pairs as discussed in Sec. I. The concomitant depletion of the pair amplitudes near the F-S interface means that unless the superconductor is thick enough, the pair amplitude does not saturate to its bulk value even deep inside the S regions. Furthermore, as we shall emphasize below, lack of self consistency may lead to violation of charge conservation: hence, while non-self-consistent approximations might be sometimes adequate for equilibrium calculations, their use must be eschewed for transport. Therefore, one should generally use a self-consistent pair potential that is allowed to spatially vary, as required by Eq. (3), and hence results in a minimum in the free energy of the system.

We begin by extending the BTK formalism to the spatially varying self-consistent pair potential obtained as explained below Eq. (3). Although the self-consistent solutions of the BdG equations reveal that the pair amplitudes are non-zero in the non-superconducting regions due to the proximity effects, the pair potential vanishes in these regions since there. Therefore, one can still use Eqs. (4) and (5), with  (7), for the wavefunctions in the and regions. To deal with the spatially varying pair potential on the S side, we divide it into many very thin layers with microscopic thicknesses of order . We treat each layer as a very thin superconductor with a constant pair potential, , as obtained from the self-consistent procedure. We are then able to write the eigenfunctions of each superconducting layer corresponding to that value of the pair potential. For example, in the -th layer, the eigenfunction should contain all left and right going solutions, and it reads:

 ΨSi≡⎛⎜ ⎜ ⎜ ⎜ ⎜⎝t1iuieik+iy+¯t1iuie−ik+iy+t4ivie−ik−iy+¯t4ivieik−iyt2iuieik+iy+¯t2iuie−ik+iy+t3ivie−ik−iy+¯t3ivieik−iyt2ivieik+iy+¯t2ivie−ik+iy+t3iuie−ik−iy+¯t3iuieik−iyt1ivieik+iy+¯t1ivie−ik+iy+t4iuie−ik−iy+¯t4ivieik−iy⎞⎟ ⎟ ⎟ ⎟ ⎟⎠, (12)

where, , and represents the strength of the normalized self consistent pair potential in the -th superconducting layer. The superconducting coherence factors and depend on in the standard way. All the coefficients in Eq. (12) are unknown, and remain to be determined. However, in the outermost S layer (rightmost in our convention) the eigenfunctions are of a form identical to Eq. (10) but with different locally constant pair potential.

We see then that the price one has to pay for including the proximity effects is the need to compute a very large number of coefficients. To do so, we adopt here a transfer matrix method to solve for these unknowns.strinati () If one considers the interface between the -th and the -th layer, we have the linear relation , where, for a generic ,

 xTi=(t1i,¯t1i,t2i,¯t2i,t3i,¯t3i,t4i,¯t4i), (13)

and the matrices, and , can be written as discussed in connection with Eq. (11). The coefficients in the -th layer can be obtained in terms of those in the -th layer as . In the same way, for the interface between the -th layer and the -th layer, we can write . From the above relations, one can write down the relation between and , i.e. . By iteration of this procedure, one can “transfer” the coefficients layer by layer and eventually relate the coefficients of the rightmost layer, , to those of the leftmost layer in S and then on to the inner ferromagnetic layer :

 xn=M−1n~Mn−1M−1n−1⋯~M1M−11~MF2xF2 (14)

By solving Eq. (14) together with , we obtain all the coefficients in the region, where the wavefunction is formally still described by the expressions given in Eqs. (4) and (5). Of course, all coefficients involved, including the energy dependent and values from which (see below) the conductance is extracted, are quite different from those in a non-self-consistent calculation. These differences will be reflected in our results. One can also prove that, when the pair potential in S is a constant (non-self-consistent), then and therefore Eq. (14) becomes . This is formally identical to that we have seen in our discussion of the non-self-consistent formalism.

This efficient technique, besides allowing us to determine all the reflected and transmitted amplitudes in the outermost layers, permits us to perform a consistency check by recomputing the self-consistent solutions to the BdG equations (the eigenfunctions). Once we have determined the amplitudes , , and , we can use them to find the amplitudes in any intermediate layer by “transferring” back the solutions. For example, the coefficients can be found by using if we know the coefficient for the rightmost layer. Knowledge of these coefficients in every region yields again the self-consistent wavefunctions of the system. These of course should be the same as the eigenfunctions found in the original procedure. Although the numerical computations involved in this consistency check are rather intensive, it is worthwhile to perform them: we have verified that, by plugging these solutions into Eq. (3) and considering all possible solutions with all possible incident angles to the BdG equations, the output pair potential obtained from the transport calculation is the same as the input pair potential obtained by direct diagonalization. This would obviously not have been the case if the initial pair potential had not been fully self consistent to begin with. The reflected and transmitted amplitudes calculated from the self-consistent solutions are in general very different from the non-self-consistent ones and lead to different quantitative behavior of the tunneling conductance, as we shall discuss in section III.

### ii.4 Charge conservation

We discuss now the important issue of the charge conservation laws. In transport calculations, it is fundamental to assure that they are not violated baym (). From the Heisenberg equation

 (15)

By computing the above commutator, we arrive at the following continuity condition

 ∂∂t⟨ρ(r)⟩+∇⋅j=−4eIm[Δ(r)⟨ψ†↑(r)ψ†↓(r)⟩]. (16)

In the steady state, which is all that we are considering here, the first term on the left is omitted. Eqn. (16) is then simply an expression for the divergence of the current. In our quasi one-dimensional system, and in terms of our wavefunctions, the conservation law can be rewritten as:

 ∂jy(y)∂y=2eIm{Δ(y)∑n[u∗n↑vn↓+u∗n↓vn↑]tanh(ϵn2T)} (17)

When the system is in equilibrium the self-consistency condition on the pair potential causes the right hand side of Eqs. (16) or (17) to vanish. This would not necessarily be the case if a non-self-consistentimaginary () solution were used.bagwell () It was shown that charge conservation is only guaranteed when self consistency is adhered to in microscopic Josephson junctions.sols2 () Current-voltage calculations for N-S heterostructures show that self-consistency is crucial to properly account for all of the Andreev scattering channels arising when the current is constant throughout the system.sols () While non-self-consistent solutions are less computationally demanding, their validity when calculating transport quantities in the nonequilibrium regime is always suspect.

In the problem we are considering, there exists a finite voltage bias between the two leads of the system (see Fig. 1). This finite bias leads to a non-equilibrium quasi-particle distribution and results of course in a net current. Still, charge conservation must hold. To see how this works in this non-equilibrium case we first write down the net quasi-particle charge density in the limit (the case we consider here) by considering the excited state caused by the finite bias . Thus, this excited state contains all single particle states () with energies less than . For simplicity, let us first consider the contribution by a single-particle state. We use to characterize this single particle state with an incident wavevector and energy . The charge density associated with it is written as

 ρ =−e∑σ⟨k∣∣ψ†σψσ∣∣k⟩ (18) =−e∑nσ(|unσ|2⟨k∣∣γ†nγn∣∣k⟩+|vnσ|2⟨k∣∣1−γ†nγn∣∣k⟩) =−e∑nσ|vnσ|2−e∑nσ(|unσ|2−|vnσ|2)δnk =−e∑nσ|vnσ|2−e∑σ(|ukσ|2−|vkσ|2)

The first term represents the ground state charge density. For a generic excited state, , that can contain many single-particle states, one need to sum over all single-particle states for the charge density such that

 ρ=−e∑nσ|vnσ|2−e∑ϵk

The quasi-particle current density from this generic excited state can also be computed,

 jy (20) =−emIm⎡⎣∑nσvnσ∂v∗nσ∂y+∑ϵk

where is a shorthand notation of . The first term in the second line vanishes because it represents the net current for the system in the ground state with a real pair potential. The right hand side of the continuity equation, Eq. (17), becomes and is responsible for the interchange between the quasi-particle current density and the supercurrent densitybtk (). We have numerically verified that by properly including these terms, all of our numerical results for the current density are constant throughout the whole system.

### ii.5 Extraction of the conductance

We are now in a position to compute the differential tunneling conductances. We begin by discussing the extraction of the conductance from the BTK theory. As we mentioned in the previous subsection, the finite bias and the resulting non-equilibrium distribution leads to an electric current flowing in the junction. In the BTK theory, this current can be evaluated from the followingbtk () expression,

 I(V)=∫G(ϵ)[f(ϵ−eV)−f(ϵ)]dϵ, (21)

where is the Fermi function. The energy dependent tunneling conductance, in the low- limit, is given as:

 G(ϵ,θi)=∑σPσGσ(ϵ,θi) (22) =∑σPσ(1+k−↑1k+σ1|a↑|2+k−↓1k+σ1|a↓|2−k+↑1k+σ1|b↑|2−k+↓1k+σ1|b↓|2),

where we have used, as is customary, natural units of conductance . In the above expression the different components are as explained in subsection II.2 (see e.g. Eq. (6)) and the and are as defined in Eqns. (4) and (5). These coefficients, which are of course energy dependent, are calculated using the self-consistent transfer matrix technique of subsection II.3. Therefore, even though Eq. (22) is formally the same in the self-consistent and non-self-consistent cases, the results for the reflection amplitudes or probabilities involved, , , , and are different in these two schemes. The angle is the incident angle, discussed in terms of components below Eq. (6). The weight factor accounts for the number of available states for spin-up and spin-down bands in the outer electrode. The tunneling conductance can also be interpreted as the transmission coefficient for electrical current. The method enables us also to compute the current density directly from the wavefunctions, Eqs. (4) and (5), in the layer by using Eq. (20) and we have been able to verify that the resulting current density is identical to the terms inside the bracket in the expression of , Eq. (22). In other words, in the low- limit the continuum-limit version of Eq. (20) is equivalent to Eq. (21).

The conductance results Eq. (22) also depend on the incident angle of electrons, . Experimentally, one can measure the forward conductance, , via point contacts or, in most other experimental conditions, an angular average. Consequently, it is worthwhile to compute the angularly averaged conductance by using the following definitions,

 ⟨Gσ(ϵ)⟩=∫θcσ0dθicosθiGσ(ϵ,θi)∫θcσ0dθicosθi, (23)

and

 ⟨G⟩=∑σPσ⟨Gσ⟩, (24)

where the critical angle is in general different for spin-up and spin-down bands. This critical angle arises from the conservation of transverse momentum and the corresponding Snell law:

 √(k+σ12+k2⊥)sinθi=√(k+σ′12+k2⊥)sinθ+rσ′ (25) =√(k−σ′12+k2⊥)sinθ−rσ′=sinθS,

where we continue to measure wavevectors in units of . The angles satisfy , and the and are each or . The last equality in Eq. (25) represents the case of the transmitted wave in S, and is the transmitted angle. Although the self-consistent pair potential varies in S and so do the quasi-particle (hole) wavevectors, we here need only consider the transmitted angle in the rightmost layer: this follows in the same way as the usual Snell’s law in a layered system, as given in elementary textbooks. From Eq. (25), one can determine the critical angles for different channels. Consider, e.g., a spin-up electron incident from without any Fermi wavevector mismatch, i.e. . Since we are only concerned with the case that the bias of tunneling junctions is of the order of superconducting gap and therefore much smaller than the Fermi energy, the approximate magnitude of the incident wavevector is , the Andreev approximation. We substitute this and similar expressions into Eq. (25) and, with the help of Eq. (6), we obtain

 √1+h1sinθi=√1−h1sinθ−r↓=sinθS. (26)

One can straightforwardly verify that, when the relation is satisfied for the incident angle, the conventional Andreev reflection becomes an evanescent wavezv2 (). In this case, the conventional Andreev reflection does not contribute to the angular averaging. On the other hand, if the energy of the incident electron is less than the saturated value of the superconducting pair amplitude in S, all the contribution to the conductance from the transmitted waves in S also vanishes because acquires an imaginary part. However, even the condition that is greater than the saturated superconducting amplitude does not guarantee that the contribution from the transmitted waves to the conductance is nonvanishing. One still needs to consider the transmitted critical angle . We define the critical angle to be the largest one among all the reflected and transmitted critical angles. It is obvious that the critical angles are different for spin-up and spin-down bands when .

### ii.6 Spin transport

We consider now the spin-transfer torque and the spin current. As the charge carriers that flow through our system, along the direction in our convention, are spin polarized, the STT provides an additional probe of the spin degree of freedom. Unlike the charge current, that must be a constant throughout the system, the spin current density is generally not a conserved quantity in the ferromagnet regions as we will demonstrate below. The discussion in Sec. II.4 on how the BTK formalism deals with the charge current can be extended to compute these spin dependent transport quantities. We need here the continuity equation for the local magnetization , where is the Bohr magneton. By using the Heisenberg equation we obtain the relation:

 ∂∂t⟨mi⟩+∂∂ySi=τi,i=x,y,z (27)

where is the spin-transfer torque, , and the spin current density is given by

 Si≡iμB2m∑σ⟨ψ†σσi∂ψσ∂y−∂ψ†σ∂yσiψσ⟩. (28)

The spin current density reduces from a tensor form to a vector because of the quasi-one-dimensional nature of our geometry. From Eq. (27), we can see that is a local physical quantity and is responsible for the change of local magnetizations due to the flow of spin-polarized currents. As we shall see in Sec. III, the conservation law (with the source torque term) for the spin density is fundamental and one has to check it is not violated when studying these transport quantities.

In the low- limit and with the presence of a finite bias, the non-equilibrium local magnetizations in Eq. (27) reads

 mx= −μB[∑n(−vn↑v∗n↓−vn↓v∗n↑) (29a) +∑ϵk

where the first summations in the expressions for denote the ground state local magnetizations. The second summations appear as a consequence of the finite bias between electrodes. The expressions for the corresponding spin currents,

 Si≡iμB2m∑ϵk

becomes

 Sx= −μBmIm[∑n(−vn↑∂v∗n↓∂y−vn↓∂v∗n↑∂y) (31a) +∑ϵk

The first summations in Eq. (31) represent the static spin current densities when there is no bias. The static spin current does not need to vanish, since a static spin-transfer torque may exist near the boundary of two magnets with misaligned exchange fields. The finite bias leads to a non-equilibrium quasi-particle distribution for the system and results in non-static spin current densities that are represented by the second summation in Eq. 31. Obviously, the spin-transfer torque has to vanish in the superconductor where the exchange field is zero. It is conventional to normalize toHalterman2007 () , where the number densities and . Following this convention, we normalize to and to .