Tunneling Conductance and Spin Transport in Clean Ferromagnet-Ferromagnet-Superconductor Heterostructures
aa r X i v : . [ c ond - m a t . s up r- c on ] S e p Tunneling Conductance and Spin Transport in Clean Ferromagnet-Ferromagnet-SuperconductorHeterostructures
Chien-Te Wu, ∗ Oriol T. Valls, † and Klaus Halterman ‡ School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455 Michelson Lab, Physics Division, Naval Air Warfare Center, China Lake, California 93555 (Dated: November 6, 2018)We present a transfer matrix approach that combines the Blonder-Tinkham-Klapwijk (BTK) formalism andself-consistent solutions to the Bogolibuov-de Gennes (BdG) equations and use it to study the tunneling conduc-tance and spin transport in ferromagnet (F)-superconductor (S) trilayers (F F S) as functions of bias voltage.The self-consistency ensures that the spin and charge conservation laws are properly satisfied. We considerforward and angularly averaged conductances over a broad range of the strength of the exchange fields andF 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 obtainedvia a non-self-consistent approach. The zero bias forward conductance peak exhibits, as φ varies, resonanceeffects 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 re-gions of the angularly averaged conductances. When F is half-metallic, the angularly averaged subgap con-ductance chiefly arises from anomalous Andreev reflection. The in-plane components of the spin current arestrongly 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 conservedin that layer and in the S region, while they oscillate in the other layer. We compute the spin transfer torques, inconnection with the oscillatory behavior of spin currents, and verify that the spin continuity equation is strictlyobeyed in our method. PACS numbers: 74.45.+c,74.78.Fk,75.75.-c
I. INTRODUCTION
Over the last two decades, significant progress in fabrica-tion techniques has allowed the development of spintronicsdevices, such as spin valves, that utilize both charge and spindegrees of freedom. Traditional spin valves consist of mag-netic materials only. There is another important type of spin-tronics devices, involving ferromagnet (F)-superconductor (S)heterostructures. These heterostructures have also receivedmuch attention because of the fundamental physics related tothe interplay between ferromagnetic and superconducting or-der. Their potential applications in spintronics include mag-netic memory technology where information storage is ac-complished via control of the magnetic moment bit. It is thencrucial to have precise control over the magnetization direc-tion. Spin transfer torque (STT) is one effect that affords suchcontrol. The generation of spin-polarized supercurrents maybe used to obtain a superconducting STT acting on the magne-tization of a ferromagnet. This effect may be utilized in highdensity nanotechnologies that require magnetic tunnel junc-tions. Thus, the dissipationless nature of the supercurrent flowoffers a promising avenue in terms of low energy nanoscalemanipulation of superconducting and magnetic orderings.Although ferromagnetism and s -wave superconductivityseem incompatible because of the inherently opposite naturesof their order parameter spin configurations, superconductiv-ity can still be induced in the F layers of F-S layered struc-tures by the superconducting proximity effects. In essence,the superconducting proximity effects describe the leakageof superconductivity into a non-superconducting normal (N)or magnetic metal, as well as its depletion in S near the in- terface. However, proximity effects in F-S systems are verydifferent from those in N-S structures due to the inherent ex-change field in the F materials. As a consequence of this ex-change field, the Cooper pair acquires a non-zero center-of-mass momentum and the overall Cooper pair wavefunctionoscillates spatially in the F regions. Owing to this oscilla-tory nature, many new physical phenomena emerge in F-S het-erostructures such as oscillations of the superconducting tran-sition temperature, T c , with the thickness of the F layers. It is of fundamental importance that superconducting prox-imity effects are governed by Andreev reflection, which isa process of electron-to-hole conversion at N-S or F-S inter-faces, and it involves the creation or annihilation of a Cooperpair. Therefore, consideration of Andreev reflection is cen-tral when studying the transport properties of N-S or F-Ssystems. Of particular interest is the behavior of thetunneling conductance in the subgap region, where hybrid sys-tems can carry a supercurrent due to Andreev reflection. Inconventional Andreev reflection, the reflected hole has oppo-site spin to the incident particle. Accordingly, the exchangefield in the F materials that causes the splitting of spin bandshas a significant effect on the tunneling conductance in thesubgap region. Most important, the qualitative behavior of theconductance peak in the zero bias limit is strongly influencedby the degree of conduction electron spin polarization in the Fmaterials.
Experimentally, this concept has been appliedto quantify the spin polarization.
An intriguing phenomenon in F-S structures is the induc-tion of triplet pairing correlations.
These correlations arevery important when studying transport phenomena such asthose found in SFS Josephson junctions.
In contrast tothe short proximity length of singlet Cooper pair conden-sates into F materials, the m = ± and the introduction of another magnetic layer with a misori-ented magnetic moment such as F SF superconducting spinvalves. The pairing state of m = ± conventional Andreevreflection, responsible for the generation of singlet Cooperpairs. Thus, recent studies on the tunneling conductancepropose the existence of anomalous
Andreev reflection, thatis, a reflected hole with the same spin as the incident parti-cle can be Andreev reflected under the same circumstances asthe generation of m = ± F S trilayer. By applying an ex-ternal magnetic field, or switching via STT, one is able tocontrol the relative orientation of the intrinsic magnetic mo-ments and investigate the dependence of physical prop-erties such as T c on the misorientation angle φ between thetwo magnetic layers. Due to the proximity effects, T c is of-ten found to be minimized when the magnetizations are ap-proximately perpendicular to each other, reflecting the pres-ence of long range triplet correlations, induced in F F S tri-layers. Their existence has been verified both theoretically and experimentally. The non-monotonic behavior of T c asa function of φ has also been shown to be quantitatively related to the long range triplet correlations, with excellentagreement between theory and experiment.Motivated by these important findings, we will investigatehere, in a fully self-consistent manner, the φ dependence of thetunneling conductance and other transport quantities of theseF F S trilayers. Non-self-consistent theoretical studies of tun-neling conductance have been performed on F F S trilayers inprevious work.
However, as we shall see in Sec. II, onlyself-consistent methods guarantee that conservation laws arenot violated and (see Sec. III) only then can one correctlypredict the proximity effects on the angular dependence oftransport properties. The spin-polarized tunneling conduc-tance of F-S bilayers only, was studied in Refs. 12, 13, 41, and42. Also, in traditional spin valves e.g. F -F layered struc-tures, the spin-polarized current generated in the F layer cantransfer angular momentum to the F layer when their mag-netic moments are not parallel to each other via the effect ofSTT. As a result, the spin current is not a conserved quan-tity and one needs a general law that relates local spin cur-rent to local STT. The transport properties of F SF struc-tures, in particular the dependence on applied bias of the spin-transfer torque and the spin-polarized tunneling conductancehave been previously studied. Here, we consider charge transport and both spin cur- rent and spin-transfer torque in F F S trilayers. In pre-vious theoretical work, such as that mentioned above,when computing tunneling conductance of N-S and F-Sstructures, using methods based on the Blonder-Tinkham-Klapwijk (BTK) procedure and quasi-classicalapproximations, the superconducting pair amplitude wasassumed to be a step function: a constant in S, droppingabruptly to zero at the N-S or F-S interface and then van-ishing in the non-superconducting region. This assumptionneglects proximity effects. Only qualitative predictions onthe behavior of the tunneling conductance can be reliablymade. Still, results exhibit many interesting features espe-cially in F-S systems. However, to fully account for theproximity effects, in the transport properties, one must usea self-consistent pair potential. This is because that revealsrealistic information regarding the leakage and depletion ofsuperconductivity. Also, as we shall discuss below, self-consistent solutions guarantee that conservation laws are sat-isfied. In Ref. 49, the tunneling conductance of F-S bilay-ers was extracted via self-consistent solutions of Bogoliubov-de Gennes (BdG) equations. However, the numerical meth-ods used there required awkward fitting procedures that ledto appreciable uncertainties and precluded their application totrilayers. The findings indicated that the self-consistent tun-neling conductances for the bilayer are quantitatively differ-ent from those computed in a non-self-consistent framework,thus demonstrating the importance of properly accounting forproximity effects in that situation. Here we report on a pow-erful self-consistent approach and use it to compute the tun-neling conductance of F F S trilayers. It is based on the BTKmethod, incorporated into a transfer matrix procedure similarto that used in Josephson junction calculations and simpleF-S junctions within a Hubbard model . As we shall demon-strate, this approach not only has the advantage of being morenumerically efficient but also can be used to compute spintransport quantities. Thus, we are able to address many impor-tant points regarding both charge and spin transport in F F Strilayers, including the spin currents and spin-transfer torque,the proximity effects on the tunneling conductance, and thecorrelation between the anomalous Andreev reflection and thetriplet correlations.This paper is organized as follows: we present our self-consistent approach, and its application to compute the tun-neling conductance, the spin-transfer torques, the spin current,and the proper way to ensure that conservation laws are sat-isfied, in Sec. II. In Sec. III we present the results. In Sub-sec. III A, we briefly compare the results of F-S bilayers ob-tained in our self-consistent approach with non-self-consistentones. The rest of Subsec. III B includes our results for trilay-ers, that is, the main results of this work. The dependenceon the tunneling conductance of F F S trilayers on the angle φ is extensively discussed as a function of geometrical andmaterial parameters. Results for the effect of the anomalousAndreev reflection, the spin-transfer torque, and the spin cur-rent are also presented. We conclude with a recapitulation inSec. IV. (cid:20) (cid:21) \][ ¥ (cid:3) FIG. 1. (Color online) Schematic of the F F S trilayer. The exchangefield, h , denoted by a black solid arrow, is along the + z direction inthe outer magnetic layer (F ) while within the inner magnetic layer(F ), h is oriented at an angle φ in the x − z plane. The outer magneticlayer and the superconducting layer are connected to electrodes thatare biased with a finite voltage V . II. METHODSA. Description of the system
The geometry of our system is depicted in Fig. 1. We de-note the outer ferromagnet as F and the middle layer as F .We choose our coordinate system so that the interfaces areparallel to the x − z plane, and infinite in extent, while thesystem has a finite width d = d F + d F + d S in the y direction.The Hamiltonian appropriate to our system is, H e f f = Z d r X α ψ † α ( r ) H ψ α ( r ) + X α, β (cid:16) i σ y (cid:17) αβ ∆ ( r ) ψ † α ( r ) ψ † β ( r ) + H . c . − X α, β ψ † α ( r ) ( h · σ ) αβ ψ β ( r ) , (1)where H is the single-particle Hamiltonian, h is a Stonerexchange field that characterizes the magnetism, and σ arePauli matrices. The superconducting pair potential ∆ ( r ) ≡ g ( r ) (cid:10) ψ ↑ ( r ) ψ ↓ ( r ) (cid:11) is the product of pairing constant, g ( r ),in the singlet channel, and the pair amplitude. We be-gin by writing down the BdG equations, which we willsolve self-consistently for our F F S trilayers. By per-forming the generalized Bogoliubov transformation , ψ σ = P n (cid:16) u n σ γ n + η σ v ∗ n σ γ † n (cid:17) , where σ = ( ↑ , ↓ ) and η σ ≡ −
1) forspin-down (up), the Hamiltonian [Eq. (1)] can be diagonal-ized. We can then for our geometry rewrite Eq. (1) as a quasi-one-dimensional eigensystem: H − h z − h x ∆ − h x H + h z ∆ ∆ − ( H − h z ) − h x ∆ − h x − ( H + h z ) u n ↑ u n ↓ v n ↑ v n ↓ = ǫ n u n ↑ u n ↓ v n ↑ v n ↓ , (2)where the u n σ and v n σ are respectively the quasiparticle andquasihole amplitudes with spin σ . The exchange field van-ishes in the S region, while in F it is directed along z , h = h ˆ z ≡ h , and in F it can rotate in the x − z plane, h = h (sin φ ˆ x + cos φ ˆ z ) ≡ h . The single-particle Hamil-tonian now reads H = − (1 / m )( d / dy ) + ǫ ⊥ − E F ( y ),where ǫ ⊥ ≡ k ⊥ / m denotes the transverse kinetic energy inthe x − z plane. Also, E F ( y ) = E FS ≡ k FS / m in the super-conducting region and E F ( y ) = E FM ≡ k FM / m in the fer-romagnetic layers. Throughout this paper, we assume naturalunits ~ = k B = E FS .To take into account the more realistic situation where the Fmaterials can in general have different bandwidths than the Slayer, we define (as in Ref. 49) a mismatch parameter Λ via E FM ≡ Λ E FS .We are aiming here to solve the problem in a fully self con-sistent manner. The self-consistent pair potential ∆ ( y ) can beexpressed in terms of the quasi-particle and quasi-hole wave-functions. Accordingly, ∆ ( y ) = g ( y )2 X n ′ (cid:2) u n ↑ ( y ) v ∗ n ↓ ( y ) + u n ↓ ( y ) v ∗ n ↑ ( y ) (cid:3) tanh (cid:18) ǫ n T (cid:19) , (3)where the primed sum is over all eigenstates with energies ǫ n smaller than a characteristic Debye energy, and g ( y ) is thesuperconducting coupling constant in the S region and van-ishes elsewhere. We obtain the self-consistent pair potentialby solving Eqs. (2) and (3) following the iterative numericalprocedures discussed in previous work. B. Application of the BTK method
The BTK formalism is a procedure to extract the trans-mitted and reflected amplitudes, and hence the conductance,from solutions to the BdG equations. This is accomplished bywriting down the appropriate eigenfunctions in different re-gions. In this subsection, we review the relevant aspects ofthe formalism for the non-self-consistent case (a step func-tion pair potential) with the objective of establishing notationand methodology to describe, in the next subsection, the pro-cedure to be used in the self-consistent case.Consider first a spin-up quasi-particle with energy ǫ , inci-dent into the left side labeled “F ”, in Fig. 1). Since the ex-change fields in the F and F layers can be non-collinear, itfollows from Eq. (2) that the spin-up (-down) quasi-particlewavefunction 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 F layer is a linear combination ofthe original incident spin-up quasi-particle wavefunctions andvarious types of reflected wavefunctions, namely reflectedspin-up and spin-down quasi-particle and quasi-hole wave-functions (via both ordinary and Andreev reflections). We usea single column vector notation to represent these combina-tions, Ψ F , ↑ ≡ e ik + ↑ y + b ↑ e − ik + ↑ y b ↓ e − ik + ↓ y a ↑ e ik −↑ y a ↓ e ik −↓ y . (4)If the incident particle has spin down, the correspondingwavefunction in F is Ψ F , ↓ ≡ b ↑ e − ik + ↑ y e ik + ↓ y + b ↓ e − ik + ↓ y a ↑ e ik −↑ y a ↓ e ik −↓ y . (5)In these expressions k ± σ are quasi-particle ( + ) and quasi-hole( − ) wavevectors in the longitudinal direction y , and satisfy therelation, k ± σ m = h Λ (1 − η σ h m ) ± ǫ − k ⊥ i / , (6)where m = m =
2, used later.As mentioned above, all energies are in units of E FS and,in addition, we measure all momenta in units of k FS . Inthis simple case, one can easily distinguish the physicalmeaning of each individual wavefunction. For instance inEq. (4), a ↓ (0 , , , T e ik −↓ y is the reflected spin-down quasi-hole wavefunction. The quasi-hole wavefunctions are the timereversed solutions of the BdG equations and carry a positivesign in the exponent for a left-going wavefunction. The rele-vant angles can be easily found in terms of wavevector com-ponents. Thus, e.g., the incident angle θ i (for spin-up) at theF − F interface is θ i = tan − (cid:16) k ⊥ / k + ↑ (cid:17) , and the Andreev re-flected angle θ − r ↓ for reflected spin-down quasi-hole wavefunc-tion is θ − r ↓ = tan − (cid:16) k ⊥ / k −↓ (cid:17) . The conservation of transversemomentum leads to many important features when oneevaluates the angularly averaged tunneling conductance, aswe will see below. For the intermediate layer F , the eigen-function in general contains both left- and right-moving planewaves, that is, Ψ F ≡ c f + ↑ e ik + ↑ y + c f + ↑ e − ik + ↑ y + c g + ↑ e ik + ↓ y + c g + ↑ e − ik + ↓ y c f + ↓ e ik + ↑ y + c f + ↓ e − ik + ↑ y + c g + ↓ e ik + ↓ y + c g + ↓ e − ik + ↓ y c f −↑ e ik −↑ y + c f −↑ e − ik −↑ y + c g −↑ e ik −↓ y + c g −↑ e − ik −↓ y c f −↓ e ik −↑ y + c f −↓ e − ik −↑ y + c g −↓ e ik −↓ y + c g −↓ e − ik −↓ y , (7)where k ±↑ and k ±↓ are defined in Eq. (6). The ± indices aredefined as previously, and the up and down arrows refer to F .The eigenspinors f and g that correspond to spin parallel orantiparallel to h respectively, are given, for 0 ≤ φ ≤ π/
2, by the expression, f + ↑ f + ↓ ! = N − cos φ sin φ ! = f −↑ − f −↓ ! ; g + ↑ g + ↓ ! = N − sin φ + cos φ ! = − g −↑ g −↓ ! (8)with the normalization constant N = p / + cos φ . Thesespinors reduce to those for pure spin-up and spin-down quasi-particles and holes when φ =
0, corresponding to a uniformmagnetization along z . One can also easily see that the partic-ular wavefunction of Eq. (7), c (cid:16) f + ↑ , f + ↓ , , (cid:17) T e ik + ↑ y denotes aquasi-particle with spin parallel to the exchange field in F .When π/ < φ ≤ π , these eigenspinors read f + ↑ f + ↓ ! = N sin φ − cos φ ! = − f −↑ f −↓ ! ; g + ↑ g + ↓ ! = N − + cos φ sin φ ! = g −↑ − g −↓ ! (9)with N = p / − cos φ .In this subsection where we are still assuming a non-self-consistent stepwise potential equal to ∆ throughout the Sregion and to zero elsewhere, we have the superconduct-ing coherence factors, √ u = (cid:20)(cid:18) ǫ + q ǫ − ∆ (cid:19) /ǫ (cid:21) / and √ v = (cid:20)(cid:18) ǫ − q ǫ − ∆ (cid:19) /ǫ (cid:21) / . In this case the right-goingeigenfunctions on the S side can be written as, Ψ S ≡ t u e ik + y + t v e − ik − y t u e ik + y + t v e − ik − y t v e ik + y + t u e − ik − y t v e ik + y + t u e − ik − y , (10)where, k ± = (cid:20) ± q ǫ − ∆ − k ⊥ (cid:21) / are quasi-particle (+) andquasi-hole (-) wavevectors in the S region. By using continu-ity of the four-component wavefunctions and their first deriva-tives at both interfaces, one can obtain all sixteen unknowncoefficients in the above expressions for the wavefunctionsby solving a set of linear equations of the form M F x F ,σ = M F x F at the F − F interface and ˜ M F x F = M S x S at theF − S interface simultaneously, where x TF , ↑ = (cid:0) , b ↑ , , b ↓ , , a ↑ , , a ↓ (cid:1) (11a) x TF , ↓ = (cid:0) , b ↑ , , b ↓ , , a ↑ , , a ↓ (cid:1) (11b) x TF = ( c , c , c , c , c , c , c , c ) (11c) x TS = ( t , , t , , t , , t , , (11d)and M F , M F , ˜ M F , and M S are appropriate 8 × a σ and b σ which are used to compute the conductance,as discussed in the next two subsections. C. Transfer matrix self consistent method
The non-self-consistent step potential assumption is largelyunrealistic. Proximity effects lead to a complicated oscillatorybehavior of the superconducting order parameter in the F lay-ers and to the generation of triplet pairs as discussedin Sec. I. The concomitant depletion of the pair amplitudesnear the F-S interface means that unless the superconductor isthick enough, the pair amplitude does not saturate to its bulkvalue even deep inside the S regions. Furthermore, as we shallemphasize below, lack of self consistency may lead to viola-tion of charge conservation: hence, while non-self-consistentapproximations might be sometimes adequate for equilibriumcalculations, their use must be eschewed for transport. There-fore, one should generally use a self-consistent pair potentialthat is allowed to spatially vary, as required by Eq. (3), andhence results in a minimum in the free energy of the system.We begin by extending the BTK formalism to the spatiallyvarying self-consistent pair potential obtained as explainedbelow Eq. (3). Although the self-consistent solutions of theBdG equations reveal that the pair amplitudes are non-zero inthe non-superconducting regions due to the proximity effects,the pair potential vanishes in these regions since g ( y ) ≡ and F regions. To deal withthe spatially varying pair potential on the S side, we divideit into many very thin layers with microscopic thicknesses oforder k − FS . We treat each layer as a very thin superconduc-tor with a constant pair potential, ∆ i , as obtained from theself-consistent procedure. We are then able to write the eigen-functions of each superconducting layer corresponding to thatvalue of the pair potential. For example, in the i -th layer, theeigenfunction should contain all left and right going solutions,and it reads: Ψ S i ≡ t i u i e ik + i y + ¯ t i u i e − ik + i y + t i v i e − ik − i y + ¯ t i v i e ik − i y t i u i e ik + i y + ¯ t i u i e − ik + i y + t i v i e − ik − i y + ¯ t i v i e ik − i y t i v i e ik + i y + ¯ t i v i e − ik + i y + t i u i e − ik − i y + ¯ t i u i e ik − i y t i v i e ik + i y + ¯ t i v i e − ik + i y + t i u i e − ik − i y + ¯ t i v i e ik − i y , (12)where, k ± i = (cid:20) ± q ǫ − ∆ i − k ⊥ (cid:21) / , and ∆ i represents thestrength of the normalized self consistent pair potential in the i -th superconducting layer. The superconducting coherencefactors u i and v i depend on ∆ i in the standard way. All thecoefficients in Eq. (12) are unknown, and remain to be de-termined. However, in the outermost S layer (rightmost inour convention) the eigenfunctions are of a form identical toEq. (10) but with different locally constant pair potential.We see then that the price one has to pay for including theproximity effects is the need to compute a very large num-ber of coefficients. To do so, we adopt here a transfer matrixmethod to solve for these unknowns. If one considers the in-terface between the i -th and the ( i + M i x i = M i + x i + , where, for a generic i , x Ti = ( t i , ¯ t i , t i , ¯ t i , t i , ¯ t i , t i , ¯ t i ) , (13)and the matrices, ˜ M i and M i + , can be written as discussedin connection with Eq. (11). The coefficients in the ( i + i -th layeras x i + = M − i + ˜ M i x i . In the same way, for the interface be-tween the ( i − i -th layer, we can write x i = M − i ˜ M i − x i − . From the above relations, one canwrite down the relation between x i + and x i − , i.e. x i + = M − i + ˜ M i M − i ˜ M i − x i − . By iteration of this procedure, onecan “transfer” the coefficients layer by layer and eventuallyrelate the coefficients of the rightmost layer, x n , to those ofthe leftmost layer in S and then on to the inner ferromagneticlayer F : x n = M − n ˜ M n − M − n − · · · ˜ M M − ˜ M F x F (14)By solving Eq. (14) together with M F x F = M F x F , weobtain all the coefficients in the F region, where the wave-function is formally still described by the expressions givenin Eqs. (4) and (5). Of course, all coefficients involved, in-cluding the energy dependent a σ and b σ values from which(see below) the conductance is extracted, are quite differentfrom those in a non-self-consistent calculation. These differ-ences will be reflected in our results. One can also provethat, when the pair potential in S is a constant (non-self-consistent), then M i + = ˜ M i and therefore Eq. (14) becomes x n = x = M − ˜ M F x F . This is formally identical to that wehave seen in our discussion of the non-self-consistent formal-ism.This efficient technique, besides allowing us to determineall the reflected and transmitted amplitudes in the outermostlayers, permits us to perform a consistency check by recom-puting the self-consistent solutions to the BdG equations (theeigenfunctions). Once we have determined the amplitudes x F , x F , and x n , we can use them to find the amplitudes inany intermediate layer by “transferring” back the solutions.For example, the coefficients x n − can be found by using x n = M − n ˜ M n − x n − if we know the coefficient x n for therightmost layer. Knowledge of these coefficients in every re-gion yields again the self-consistent wavefunctions of the sys-tem. These of course should be the same as the eigenfunc-tions found in the original procedure. Although the numericalcomputations involved in this consistency check are rather in-tensive, it is worthwhile to perform them: we have verifiedthat, by plugging these solutions into Eq. (3) and consider-ing all possible solutions with all possible incident angles tothe BdG equations, the output pair potential obtained fromthe transport calculation is the same as the input pair poten-tial obtained by direct diagonalization. This would obviouslynot have been the case if the initial pair potential had not beenfully self consistent to begin with. The reflected and transmit-ted amplitudes calculated from the self-consistent solutionsare in general very different from the non-self-consistent onesand lead to different quantitative behavior of the tunnelingconductance, as we shall discuss in section III. D. Charge conservation
We discuss now the important issue of the charge conserva-tion laws. In transport calculations, it is fundamental to assurethat they are not violated . From the Heisenberg equation ∂∂ t h ρ ( r ) i = i Dh H e f f , ρ ( r ) iE . (15)By computing the above commutator, we arrive at the follow-ing continuity condition ∂∂ t h ρ ( r ) i + ∇ · j = − e Im h ∆ ( r ) D ψ †↑ ( r ) ψ †↓ ( r ) Ei . (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 simplyan expression for the divergence of the current. In our quasione-dimensional system, and in terms of our wavefunctions,the conservation law can be rewritten as: ∂ j y ( y ) ∂ y = e Im ∆ ( y ) X n h u ∗ n ↑ v n ↓ + u ∗ n ↓ v n ↑ i tanh (cid:18) ǫ n T (cid:19) (17)When the system is in equilibrium the self-consistency con-dition on the pair potential causes the right hand side ofEqs. (16) or (17) to vanish. This would not necessarilybe the case if a non-self-consistent solution were used. It was shown that charge conservation is only guaranteedwhen self consistency is adhered to in microscopic Josephsonjunctions. Current-voltage calculations for N-S heterostruc-tures show that self-consistency is crucial to properly accountfor all of the Andreev scattering channels arising when thecurrent is constant throughout the system. While non-self-consistent solutions are less computationally demanding, theirvalidity when calculating transport quantities in the nonequi-librium regime is always suspect.In the problem we are considering, there exists a finite volt-age bias V between the two leads of the system (see Fig. 1).This finite bias leads to a non-equilibrium quasi-particle dis-tribution and results of course in a net current. Still, chargeconservation must hold. To see how this works in this non-equilibrium case we first write down the net quasi-particlecharge density in the T → | k k · · · i caused by the fi-nite bias V . Thus, this excited state contains all single particlestates | k j i ( j = , , · · · ) with energies less than eV . For sim-plicity, let us first consider the contribution by a single-particlestate. We use | k i to characterize this single particle state withan incident wavevector k = k ⊥ + k ˆ y and energy ǫ k . The chargedensity associated with it is written as ρ = − e X σ D k (cid:12)(cid:12)(cid:12) ψ † σ ψ σ (cid:12)(cid:12)(cid:12) k E (18) = − e X n σ (cid:16) | u n σ | D k (cid:12)(cid:12)(cid:12) γ † n γ n (cid:12)(cid:12)(cid:12) k E + | v n σ | D k (cid:12)(cid:12)(cid:12) γ n γ † n (cid:12)(cid:12)(cid:12) k E(cid:17) = − e X n σ (cid:16) | u n σ | D k (cid:12)(cid:12)(cid:12) γ † n γ n (cid:12)(cid:12)(cid:12) k E + | v n σ | D k (cid:12)(cid:12)(cid:12) − γ † n γ n (cid:12)(cid:12)(cid:12) k E(cid:17) = − e X n σ | v n σ | − e X n σ (cid:16) | u n σ | − | v n σ | (cid:17) δ n k = − e X n σ | v n σ | − e X σ (cid:16) | u k σ | − | v k σ | (cid:17) The first term represents the ground state charge density. For ageneric excited state, | k k · · · i , that can contain many single-particle states, one need to sum over all single-particle statesfor the charge density such that ρ = − e X n σ | v n σ | − e X ǫ k < eV X σ (cid:16) | u k σ | − | v k σ | (cid:17) . (19) The quasi-particle current density from this generic excitedstate can also be computed, j y = − e m X ǫ k < eV X σ * − i ψ † σ ∂∂ y ψ σ + i ∂∂ y ψ † σ ! ψ σ + k (20) = − em Im X n σ v n σ ∂ v ∗ n σ ∂ y + X ǫ k < eV X σ u ∗ k σ ∂ u k σ ∂ y + v ∗ k σ ∂ v k σ ∂ y ! = − em Im X ǫ k < eV X σ u ∗ k σ ∂ u k σ ∂ y + v ∗ k σ ∂ v k σ ∂ y ! , where h ... i k is a shorthand notation of h k | ... | k i . The first termin the second line vanishes because it represents the net cur-rent for the system in the ground state with a real pair poten-tial. The right hand side of the continuity equation, Eq. (17),becomes − e Im h ∆ P ǫ k < eV (cid:16) u ∗ k ↑ v k ↓ + v k ↑ u ∗ k ↓ (cid:17)i and is responsi-ble for the interchange between the quasi-particle current den-sity and the supercurrent density . We have numerically veri-fied that by properly including these terms, all of our numeri-cal results for the current density are constant throughout thewhole system. E. Extraction of the conductance
We are now in a position to compute the differential tun-neling conductances. We begin by discussing the extractionof the conductance from the BTK theory. As we mentionedin the previous subsection, the finite bias V and the resultingnon-equilibrium distribution leads to an electric current flow-ing in the junction. In the BTK theory, this current can beevaluated from the following expression, I ( V ) = Z G ( ǫ ) (cid:2) f ( ǫ − eV ) − f ( ǫ ) (cid:3) d ǫ, (21)where f is the Fermi function. The energy dependent tun-neling conductance, G ( ǫ ) = ∂ I /∂ V | V = ǫ in the low- T limit, isgiven as: G ( ǫ, θ i ) = X σ P σ G σ ( ǫ, θ i ) (22) = X σ P σ + k −↑ k + σ | a ↑ | + k −↓ k + σ | a ↓ | − k + ↑ k + σ | b ↑ | − k + ↓ k + σ | b ↓ | , where we have used, as is customary, natural units of conduc-tance ( e / h ). In the above expression the different k compo-nents are as explained in subsection II B (see e.g. Eq. (6)) andthe a σ and b σ are as defined in Eqns. (4) and (5). These coef-ficients, which are of course energy dependent, are calculatedusing the self-consistent transfer matrix technique of subsec-tion II C. Therefore, even though Eq. (22) is formally the samein the self-consistent and non-self-consistent cases, the resultsfor the reflection amplitudes or probabilities involved, | a ↑ | , | a ↓ | , | b ↑ | , and | b ↓ | are different in these two schemes. Theangle θ i is the incident angle, discussed in terms of k compo-nents below Eq. (6). The weight factor P σ ≡ (1 − h η σ ) / layer by using Eq. (20) and we havebeen able to verify that the resulting current density is identi-cal to the terms inside the bracket in the expression of G ( ǫ ),Eq. (22). In other words, in the low- T limit the continuum-limit version of Eq. (20) is equivalent to Eq. (21).The conductance results Eq. (22) also depend on the in-cident angle of electrons, θ i . Experimentally, one can mea-sure the forward conductance, θ i =
0, via point contacts or, inmost other experimental conditions, an angular average. Con-sequently, it is worthwhile to compute the angularly averagedconductance by using the following definitions, h G σ ( ǫ ) i = R θ c σ d θ i cos θ i G σ ( ǫ, θ i ) R θ c σ d θ i cos θ i , (23)and h G i = X σ P σ h G σ i , (24)where the critical angle θ c σ is in general different for spin-up and spin-down bands. This critical angle arises from theconservation of transverse momentum and the correspondingSnell law: q(cid:16) k + σ + k ⊥ (cid:17) sin θ i = q(cid:16) k + σ ′ + k ⊥ (cid:17) sin θ + r σ ′ = q(cid:16) k − σ ′ + k ⊥ (cid:17) sin θ − r σ ′ = sin θ S , (25)where we continue to measure wavevectors in units of k FS .The angles θ ± r σ satisfy tan − (cid:16) k ⊥ / k ± σ (cid:17) , and the σ and σ ′ areeach ↑ or ↓ . The last equality in Eq. (25) represents the caseof the transmitted wave in S, and θ S is the transmitted angle.Although the self-consistent pair potential varies in S and sodo the quasi-particle (hole) wavevectors, we here need onlyconsider the transmitted angle θ S in the rightmost layer: thisfollows in the same way as the usual Snell’s law in a layeredsystem, as given in elementary textbooks. From Eq. (25), onecan determine the critical angles for different channels. Con-sider, e.g., a spin-up electron incident from F without anyFermi wavevector mismatch, i.e. Λ =
1. Since we are onlyconcerned with the case that the bias of tunneling junctionsis of the order of superconducting gap and therefore muchsmaller than the Fermi energy, the approximate magnitude ofthe incident wavevector is √ + h , the Andreev approxima-tion. We substitute this and similar expressions into Eq. (25)and, with the help of Eq. (6), we obtain p + h sin θ i = p − h sin θ − r ↓ = sin θ S . (26)One can straightforwardly verify that, when the relation θ i > sin − [((1 − h ) / (1 + h )) / ] is satisfied for the incident angle,the conventional Andreev reflection becomes an evanescentwave . In this case, the conventional Andreev reflection does not contribute to the angular averaging. On the other hand, ifthe energy ǫ of the incident electron is less than the saturatedvalue of the superconducting pair amplitude in S, all the con-tribution to the conductance from the transmitted waves in Salso vanishes because k ± acquires an imaginary part. How-ever, even the condition that ǫ is greater than the saturatedsuperconducting amplitude does not guarantee that the contri-bution from the transmitted waves to the conductance is non-vanishing. One still needs to consider the transmitted criticalangle sin − [1 / (1 + h ) / ]. We define the critical angle θ c σ to bethe largest one among all the reflected and transmitted criticalangles. It is obvious that the critical angles θ c σ are differentfor spin-up and spin-down bands when h , F. Spin transport
We consider now the spin-transfer torque and the spin cur-rent. As the charge carriers that flow through our system,along the y direction in our convention, are spin polarized,the STT provides an additional probe of the spin degree offreedom. Unlike the charge current, that must be a constantthroughout the system, the spin current density is generallynot a conserved quantity in the ferromagnet regions as we willdemonstrate below. The discussion in Sec. II D on how theBTK formalism deals with the charge current can be extendedto compute these spin dependent transport quantities. Weneed here the continuity equation for the local magnetization m ≡ − µ B P σ D ψ † σ σ ψ σ E , where µ B is the Bohr magneton. Byusing the Heisenberg equation ∂∂ t h m ( r ) i = i Dh H e f f , m ( r ) iE we obtain the relation: ∂∂ t h m i i + ∂∂ y S i = τ i , i = x , y , z (27)where τ is the spin-transfer torque, τ ≡ m × h , and the spincurrent density S i is given by S i ≡ i µ B m X σ * ψ † σ σ i ∂ψ σ ∂ y − ∂ψ † σ ∂ y σ i ψ σ + . (28)The spin current density reduces from a tensor form to a vectorbecause of the quasi-one-dimensional nature of our geometry.From Eq. (27), we can see that S is a local physical quantityand τ is responsible for the change of local magnetizationsdue to the flow of spin-polarized currents. As we shall see inSec. III, the conservation law (with the source torque term) forthe spin density is fundamental and one has to check it is notviolated when studying these transport quantities.In the low- T limit and with the presence of a fi-nite bias, the non-equilibrium local magnetizations m i ≡ P ǫ k < eV P σ − µ B h ψ † σ σ i ψ σ i k in Eq. (27) reads m x = − µ B X n (cid:16) − v n ↑ v ∗ n ↓ − v n ↓ v ∗ n ↑ (cid:17) (29a) + X ǫ k < eV (cid:16) u ∗ k ↑ u k ↓ + v k ↑ v ∗ k ↓ + u ∗ k ↓ u k ↑ + v k ↓ v ∗ k ↑ (cid:17) m y = − µ B i X n (cid:16) v n ↑ v ∗ n ↓ − v n ↓ v ∗ n ↑ (cid:17) (29b) − i X ǫ k < eV (cid:16) u ∗ k ↑ u k ↓ + v k ↑ v ∗ k ↓ − u ∗ k ↓ u k ↑ − v k ↓ v ∗ k ↑ (cid:17) m z = − µ B X n (cid:16) | v n ↑ | − | v n ↓ | (cid:17) (29c) + X ǫ k < eV (cid:16) | u k ↑ | − | v k ↑ | − | u k ↓ | + | v k ↓ | (cid:17) , where the first summations in the expressions for m i denotethe ground state local magnetizations. The second summa-tions appear as a consequence of the finite bias between elec-trodes. The expressions for the corresponding spin currents, S i ≡ i µ B m X ǫ k < eV X σ * ψ † σ σ i ∂ψ σ ∂ y − ∂ψ † σ ∂ y σ i ψ σ + k , (30)becomes S x = − µ B m Im X n − v n ↑ ∂ v ∗ n ↓ ∂ y − v n ↓ ∂ v ∗ n ↑ ∂ y ! (31a) + X ǫ k < eV u ∗ k ↑ ∂ u k ↓ ∂ y + v k ↑ ∂ v ∗ k ↓ ∂ y + u ∗ k ↓ ∂ u k ↑ ∂ y + v k ↓ ∂ v ∗ k ↑ ∂ y ! S y = µ B m Re X n − v n ↑ ∂ v ∗ n ↓ ∂ y + v n ↓ ∂ v ∗ n ↑ ∂ y ! (31b) + X ǫ k < eV u ∗ k ↑ ∂ u k ↓ ∂ y + v k ↑ ∂ v ∗ k ↓ ∂ y − u ∗ k ↓ ∂ u k ↑ ∂ y − v k ↓ ∂ v ∗ k ↑ ∂ y ! S z = − µ B m Im X n v n ↑ ∂ v ∗ n ↑ ∂ y − v n ↓ ∂ v ∗ n ↓ ∂ y ! (31c) + X ǫ k < eV u ∗ k ↑ ∂ u k ↑ ∂ y − v k ↑ ∂ v ∗ k ↑ ∂ y − u ∗ k ↓ ∂ u k ↓ ∂ y + v k ↓ ∂ v ∗ k ↓ ∂ y ! . The first summations in Eq. (31) represent the static spin cur-rent densities when there is no bias. The static spin currentdoes not need to vanish, since a static spin-transfer torquemay exist near the boundary of two magnets with misalignedexchange fields. The finite bias leads to a non-equilibriumquasi-particle distribution for the system and results in non-static spin current densities that are represented by the secondsummation in Eq. 31. Obviously, the spin-transfer torque hasto vanish in the superconductor where the exchange field iszero. It is conventional to normalize m to − µ B ( N ↑ + N ↓ ),where the number densities N ↑ = k FS (1 + h m ) / / (6 π ) and N ↓ = k FS (1 − h m ) / / (6 π ). Following this convention, we nor-malize τ to − µ B ( N ↑ + N ↓ ) E FS and S to − µ B ( N ↑ + N ↓ ) E FS / k FS . h=0.866 G E h=0.7 G ( E = ) Λ FIG. 2. (Color online) Bias dependence of the results for the forwardconductance, G , in thick F-S bilayers (see text). The values of h are indicated. In both main panels the solid and dashed curves show G , in units of e / h for non-self-consistent and self-consistent results,respectively. The bias E is in units of the S bulk gap ∆ . In the toppanel the (red) lower curves are for a mismatch parameter Λ = . Λ = .
5, and the (blue) higher curvesfor
Λ =
1. In the bottom panel, the (purple) top curves are for
Λ = .
41, the (blue) curves are as in the top panel, and the (black) lowerones for
Λ = .
71. The inset (see text) shows G ( E =
1) vs. Λ in theself consistent calculation (dots) and the non-self-consistent result(line). III. RESULTS
The forward scattering conductances G are computed byconsidering a particle incident with an angle θ i (cid:27) E ≡ eV is in units of the zerotemperature gap, ∆ , in bulk S material and e / h is used as thenatural unit of conductance. When the F and F regions aremade of same F material, i.e., h = h and k ±↑ , ( ↓ = k ±↑ , ( ↓ )2 ,we will use h (not to be confused with Planck’s constant) and k ±↑ , ( ↓ ) to denote their exchange fields and wavevectors. Thisis the case we will mostly study. All results are for the low- T limit. All of the lengths are measured in unit of k − FS anddenoted by capital letters, e.g. D S denotes k FS d S . A. Bilayers
We begin with a brief discussion of self-consistent resultsfor the tunneling conductance in F-S bilayers, contrastingthem with non-self-consistent results. We assume that the Slayer is very thick so that the pair amplitude saturates to itsbulk value deep inside the S region. In this subsection, the di-mensionless superconducting coherence length Ξ is taken tobe 50 and the thicknesses D F and D S of the F and S layers areboth 15 Ξ . By computing the pair amplitudes via the directdiagonalization method, we have verified that they indeedsaturate to their bulk value with this large ratio of D S to Ξ .As discussed in Sec. I, the replacement of non-magneticmetals with ferromagnets in a bilayer leads to strong suppres-sion of the Andreev reflection in the subgap region. The de-crease of the zero bias conductance (ZBC) strongly dependson the magnitude of the exchange field in F. This dependenceis used to measure the degree of spin-polarization of magneticmaterials experimentally. However, in early theoreticalwork, it was shown that to accurately determine the degreeof spin-polarization, one has to consider the Fermi wavevec-tor mismatch (FWM), Λ , as well as the interfacial barriers.The ZBC peak is very sensitive to both spin-polarization andFWM and the dependence cannot be characterized by a singleparameter.We display in Fig. 2 forward conductance vs. bias resultsfor both the self-consistent and non-self-consistent calcula-tions, at two different values of the exchange fields and sev-eral FWM values. One sees at once that the self-consistentresults approach the non-self-consistent ones in the zero biaslimit, while deviating the most for energies near the supercon-ducting gap. The ZBC decreases with increasing h and withdecreasing Λ . Also, larger h indeed leads to a conspicuousreduction in the subgap conductance and so does the intro-duction of FWM. One can conclude that the behavior of theZBC can not be characterized by only one parameter, either h or Λ . Instead, one should expand the fitting parameter spaceto determine the degree of spin polarization.In the non-self-consistent framework, the conductance atthe superconducting gap ( E = Λ at a given h . However, earlier work predicted that thisconclusion is invalid in self-consistent approach, and that theconductance at the superconducting gap varies monotonicallywith increasing Λ . Here we verify this via our self-consistenttransfer matrix method. The inset in the bottom panel of Fig. 2clearly shows this dependence on Λ . Figure 2 also shows thatthe self-consistent results (dashed curves) on subgap conduc-tances are in general lower than those obtained in the non-self-consistent framework (solid curves) for a strong exchangefield. On the other hand, in the high bias limit, the self-consistent results become similar to the non-self-consistentones. This is simply because the particle does not experiencemuch of a difference between a step-like pair potential and asmooth pair potential when it is incident with high enough en-ergy. Finally, clear cusps appear at the superconducting gapvalue in some cases, e.g., the forward scattering conductancecurve at h = .
866 and
Λ =
1. This is consistent with what isfound in previous work for thick bilayers. G E φ =0 ° φ =90 ° φ =180 ° FIG. 3. (Color online) Comparison between the self-consistent andnon-self-consistent forward scattering conductances of F F S trilay-ers. The solid and the dashed lines are for non-self-consistent andself-consistent results respectively. The (red) curves, highest at thecritical bias (CB) are for φ = ◦ . The (blue) curves, lowest at CB,are for φ = ◦ . We have D F = D F =
12, and D S =
180 (seetext).
B. Trilayers
We now discuss our results for F F S trilayers of finitewidths. First, we discuss the dependence of the tunneling con-ductances on the angle φ between h and h (see below Eq. (2)and Fig. 1). An important reason for considering trilayers withfinite widths is the strong dependence of the superconductingtransition temperatures T c on the angle φ due to proximityeffects and induced long-range triplet correlations. Fieldinduced switching effects also make these structures attrac-tive candidates for memory elements. The non-monotonic be-havior of T c ( φ ) with its minimum being near φ = ◦ , wasextensively discussed in Ref. 39. This angular dependencehas been shown to be related to the induced triplet pairingcorrelations . The superconducting transition temperaturesare also predicted to be positively correlated with the singletpair amplitudes deep inside the S regions . Therefore, it is ofparticular importance to consider systems of finite size to takeinto view the whole picture of proximity effects on the angu-lar dependence of the tunneling conductance. For the resultsshown in this subsection, we assume the absence of FWM( Λ =
1. Forward Scattering
As a typical example of our results, we show in Fig. 3 re-sults for the φ dependence of the forward scattering conduc-tances. The exchange field we use here for both F layers is h = .
3, and the thicknesses of the F and F layers corre-spond to D F =
10 and D F =
12 respectively, while the Slayer has width D S = = . Ξ . Results obtained via thenon-self-consistent approach are plotted for comparison. Inthe non-self-consistent framework where the single parame-0 G h=0.5 φ =0 °φ =30 °φ =60 °φ =90 °φ =120 °φ =150 °φ =180 ° h=0.45 h=0.6 h=0.35 E h=0.4 h=0.425 FIG. 4. (Color online) Forward scattering conductance of F F S trilayers for several angles φ as indicated in the legend. The top panels are for D F = D F =
12, and D S =
180 and the bottom panels for D F = D F =
18, and D S = h is indicated.For the left panels, the conductances at CB decrease with increasing φ . For the other panels, the ZBC (see text) decreases as φ increases. ter ∆ describes the stepwise pair potential, one sees in Fig. 3that for all values of the angle φ the conductance curves dropwhen the bias is at ∆ , corresponding to E = φ non-monotonically, with φ = ◦ corresponding to the largest and φ = ◦ to the smallest bias values. Since the CB depends onthe strength of the superconducting gap deep inside the S re-gions, the non-monotonicity of the CB in Fig. 3 is correlatedwith the non-monotonicity of T c . The CB never reaches unity,in these trilayers, due to their finite size. Accordingly, thisfeature of the correct self-consistent results implies that onecannot adequately determine the angular dependence of theforward conductance in the non-self-consistent framework.This feature also provides experimentalists with another wayto measure the strength of the superconducting gap for differ-ent angles in these trilayers by determining the CB in a setof conductance curves. The remaining results shown in thissection are all computed self-consistently.In Fig. 4, we present more results for the dependence ofthe forward scattering conductances on φ . In the top panelsthe thicknesses of each layer and the coherence length are thesame as Fig. 3. In the bottom panels we increase the thick-ness of the inner magnetic layer to D F =
18 while D F , D S ,and Ξ remain unchanged. For each row of Fig. 4, results forthree different exchange fields are plotted. In the top left panel( h = .
5) we see that the angular dependence of the CB (or themagnitude of the saturated pair amplitudes) is monotonic with φ . Although this monotonicity is not common, we have ver-ified that it is consistent with the theoretical results for T c ( φ )for the same particular case. The more usual non-monotonic dependence is found in all other panels, as discussed in theprevious paragraph. In every case, we have also checked thatthe magnitude of the CB reflects the magnitude of the self-consistent pair amplitudes deep inside the superconductor.For the ZBC, we see that the degree of its angular depen-dence is very sensitive to h . In the top left panel, with h = . φ . On the other hand, theZBC in the top right panel, h = .
6, drops by almost a fac-tor of two as φ varies from the relative parallel (P) orienta-tion, φ = ◦ , to the antiparallel (AP) orientation, φ = ◦ .This is a consequence of interference between the spin-up andspin-down wavefunctions under the influence of the rotatedexchange field in the middle layer. In the top left panel, wesee that the conductance at CB decreases with increasing an-gle. In other words, the zero bias conductance peak (ZBCP)becomes more prominent as φ is increased. However, for thetop middle panel, h = .
45, the development of the ZBCP isless noticeable when the angle φ is increased. In the top rightpanel, h = .
6, the ZBCP evolves into a zero bias conductancedip (ZBCD) as φ varies from φ = ◦ to φ = ◦ , with a clearfinite bias conductance peak (FBCP) appearing just below theCB. This behavior is reminiscent of that which occurs whena barrier, or mismatch, are present. In the bottom panels of thisfigure, corresponding to a larger value of D F one can observesimilar features. For example, a slight change from h = . h = . φ > ◦ . Thelocation of the FBCP also moves closer to the CB value when φ increases. That these features of the ZBC depend on boththe strength of exchange field (reflected in k ±↑ and k ±↓ ) and thethickness of the F layer indicates that the ZBC shows thecharacteristics of a resonance scattering phenomenon as in an1 G ( E = , φ = ° ) h D F2 =18 D F2 =12 FIG. 5. (Color online) Resonance effects in the forward scatteringconductance at zero bias for trilayers at φ = ◦ . In the top panel,the trilayers have same thicknesses as in the top panels of Fig. 4, andin the bottom panel, they are as in the bottom panels of Fig. 4.The(blue) dots are the results from our computations and the (red) curvesfrom Eq. (33). elementary quantum mechanical barrier. The main differenceis that the scattering problem here involves the intricate inter-ference between quasi-particle and quasi-hole spinors.When the bias is high enough, the tunneling conductanceapproaches its normal state value. Thus, one can extractthe magnetoresistance from the conductance at E =
2. Weonly discuss here the magnetoresistance’s qualitative behav-ior. One can define a measure of the magnetoresistance as, M G ( E , φ ) ≡ G ( E , φ = ◦ ) − G ( E , φ ) G ( E , φ = ◦ ) . (32)For all results shown in the panels of Fig. 4, the conductanceat E = φ , i.e., it is a monotonicfunction of φ , the standard behavior for conventional, non-superconducting, spin-valves. Furthermore, one can also seethat M G ( E = , φ = ◦ ) increases with exchange field.Therefore, the behavior of the magnetoresistance at large biasis as one would expect in the present self-consistent BTKframework. However, the behavior of M G ( E = , φ = ◦ )that is associated with the behavior of the ZBC is generally anon-monotonic function of h .We next investigate the high sensitivity of the ZBC to h byexamining its resonances for two different F widths arrangedin an AP magnetic configuration ( φ = ◦ ). To do so, weperformed an analytic calculation of the ZBC in the non-self-consistent framework in situations where (as discussed in con- nection with Fig. 3) the results nearly coincide with those ofself-consistent calculations. We find that the ZBC at φ = ◦ , G ( E = , φ = ◦ ) ≡ G ZB , for a given h and D F is: G ZB = k ↑ k ↓ A + (cid:0) h − h − h k ↑ k ↓ (cid:1) cos (cid:2) (cid:0) k ↑ − k ↓ (cid:1) D F (cid:3) . (33)The expression for A in Eq. (33) is: a sin (cid:2)(cid:0) k ↑ + k ↓ (cid:1) D F (cid:3) + a (cid:2) cos (cid:0) k ↑ D F (cid:1) − cos (cid:0) k ↓ D F (cid:1)(cid:3) + a , (34)where a = h (1 − k ↑ k ↓ ) , a = h , and a = h + ( − + h − k ↑ k ↓ ) . Here we have omitted the ± indices for the quasi-particle and quasi-hole wavevectors, since we are in the zerobias limit. In Fig. 5, we plot Eq. (33) as a function of h for D F =
12 (top panel) and 18 (bottom panel). In this zero biaslimit, the (blue) circles (self-consistent numerical results) areon top of the (red) curves (analytic results). As the thicknessof the intermediate layer increases, the number of resonancemaxima and minima increases. Therefore, the resonance be-havior of the ZBC is more sensitive to h for larger D F , as wehave seen in Fig. 4. For a given D F , the ZBC drops consid-erably as φ varies from φ = ◦ to φ = ◦ when h is near theminimum of the resonance curve (rightmost panels of Fig. 4).On the other hand, when h is near the resonance maximum(leftmost panels of Fig. 4), the ZBC is a very weak functionof φ provided that h is not too strong. By examining the de-nominator of Eq. (33), we find that the terms involved in A areless important than the last term. This is because the wave-length ( k ↑ − k ↓ ) − associated with that term is the dominantcharacteristic wavelength in the theory of proximity effectsin F-S structures. In both panels of Fig. 5, we see that theZBC for φ = ◦ vanishes in the half-metallic limit. To showthis analytically, one can use the conservation of probabilitycurrents and write down the relation, valid when the bias issmaller than the superconducting gap: k −↑ k + σ | a ↑ | + k −↓ k + σ | a ↓ | + k + ↑ k + σ | b ↑ | + k + ↓ k + σ | b ↓ | = . (35)By combining Eq. (35) with Eq. (22), it becomes clear that thesubgap conductances arise largely from Andreev reflection. Inthe half-metallic limit, conventional Andreev reflection is for-bidden due to the absence of an opposite-spin band: this leadsto zero ZBC at φ = ◦ . Same-spin Andreev reflection (seediscussion in the paragraph above Eq. (4)) is not allowed incollinear magnetic configurations. Equation (35) also reflectsanother important feature of the ZBC: the contributions to G atzero bias from the spin-up and down channels are identical ex-cept for the weight factor P σ : one can prove analytically thatthe sum of first two terms (related to Andreev reflection) inEq. (35) is spin-independent. As a result, the sum of last twoterms, related to ordinary reflection, is also spin-independent,and so is the ZBC.We briefly consider here one example where the two F ma-terials in the trilayers have different field strengths. In thisexample all the thicknesses and the coherence length are as inthe top panels of Fig. 4. In Fig. 6, we plot the forward scatter-ing conductance for several φ at h = . h = .
1. One2 G E φ =0 °φ =30 °φ =60 °φ =90 °φ =120 °φ =150 °φ =180 ° FIG. 6. (Color online) Forward scattering conductance of a F F Strilayer with differing magnetic materials corresponding to exchangefields of h = . h = .
1. Various magnetic orientations, φ , areconsidered as shown. Geometry and other parameters are as in thetop panels of Fig. 4. can quickly identify that the ZBC here is a non-monotonicfunction of φ with it maximum at the orthogonal relative mag-netization angle, φ = ◦ . In contrast, results at equal ex-change field strengths usually demonstrate monotonic behav-ior, as previously shown. However, many features are still thesame, such as the formation of a FBCP when φ > ◦ . For φ = ◦ and φ = ◦ , the conductance curves are not mono-tonically decreasing, as was the case at h = h . There, when h = h and φ < ◦ , we always see monotonically decreasingbehavior because the scattering effect due to misoriented mag-netizations is not as great as at φ > ◦ . Also, when h , h ,we have to include in our considerations another scattering ef-fect that comes from the mismatch between k ±↑ , ( ↓ and k ±↑ , ( ↓ .Specifically, when φ = ◦ , the only important scattering effectis that due to mismatch from h , h and it leads to suppres-sion of the ZBC at φ = ◦ . However, we see that the scat-tering due to the misoriented magnetic configuration ( φ , ◦ )compensates the effect of mismatch from h , h and ZBCis maximized when φ = ◦ . Qualitatively, one can examineEqs. (8) and (9) and verify that the spinor at φ = ◦ is com-posed of both pure spin-up and spin-down spinors with equalweight, apart from phase factors. As a result, the scatteringeffect due to mismatch from k ±↑ , ( ↓ and k ±↑ , ( ↓ is reduced. Wealso verified that, when the strength of h is increased towards h , the locations for the maximum of the ZBC( φ ) curves grad-ually move from φ = ◦ at h = . φ = ◦ at h = .
2. Angularly averaged conductance
We now present results for the angularly averaged conduc-tance, h G i as defined in Eq. (23). The details of the angularaveraging are explained under Eq. (25). The angularly aver-aged conductance is relevant to a much wider range of experi-mental results than the forward conductance, which is relevantstrictly only for some point contact experiments. This is par-ticularly true if one recalls that the critical angle θ c σ and the weight factor for angular averaging in Eq. (23) used in thiswork can be modified based on a real experimental set-up oron the geometry of the junction.In Fig. 7, we present results for h G i at D F =
12 (left pan-els) and D F =
18 (right panels). All curves are obtained with D F =
10 and D S = = . Ξ at the values of h indi-cated in each panel. Results are plotted over the entire rangeof φ values. The CB values obtained for h G i are again non-monotonic functions of φ and the non-monotonicity matchesthat of the saturated pair amplitudes, for the reasons previ-ously given. The CB values for h G i in these cases are thesame as those for the forward scattering conductance. Onecan also see that the resonance phenomenon is washed out inthe angularly averaged conductance. For example, the reso-nance curve in the top panel of Fig. 5 tells us that h ≈ . h ≈ . G . However, inthe top left panel of Fig. 7, the ZBC is no longer a weak func-tion of φ and it gradually decreases when φ is increased. Nearthe resonance minimum, h = .
6, bottom left panel of Fig. 7,we can see a trace of the appearance of the FBCP when φ isabove 90 ◦ . This FBCP in h G i is not as prominent as that in theforward scattering G , due to the averaging.The magnetoresistance measure M G ( E = , φ ) is larger for h G i than for the forward scattering conductance. For example, M G ( E = , φ = ◦ ) in the forward scattering conductancefor h = . D F =
12 is half of that in h G i . As for the zerobias magnetoresistance M G ( E = , φ = ◦ ) in h G i , it is ofabout the same order as M G ( E = , φ = ◦ ) and it does notdepend on where it is located in the resonance curve, Fig. 5(recall that M G ( E = , φ = ◦ ) for the forward scatteringconductance almost vanishes at the resonance maximum).In the right panels of Fig. 7, we plot results for a larger D F with values of h = .
35 (near a resonance maximum)and h = .
725 (near a resonance minimum). They share verysimilar features with the thinner D F case in the left panels.However, for h = . φ shrink to almost or less than unity and they are justbarely higher than the conductance at E = h G i by dividing the contribution from all angles into two ranges:the range above and the range below the conventional An-dreev critical angles θ Ac [see discussion below Eq. (26)]. Con-sider e.g., the case of spin-up incident quasi-particles. When θ Ac ≡ sin − [ √ (1 − h ) / (1 + h )] < θ i < sin − [ √ / (1 + h )],the conventional Andreev reflected waves become evanescentwhile the transmitted waves are still traveling waves abovethe CB. When θ i > sin − [ √ / (1 + h )], both the conventionalAndreev reflected waves and the transmitted waves becomeevanescent. Here, θ c ↑ = sin − [ √ / (1 + h )] is the upper limitin Eq. (23).The case of spin-down incident quasi-particles is trivial, be-cause the dimensionless incident momentum is √ − h whichis less than both the conventional Andreev reflected wavevec-3 E D F2 =12 and h=0.6 〈 G 〉 D F2 =12 and h=0.3 φ =0 °φ =30 °φ =60 °φ =90 °φ =120 °φ =150 °φ =180 ° D F2 =18 and h=0.725D F2 =18 and h=0.35 FIG. 7. (Color online) Bias dependence of the angularly averaged conductance of F F S trilayers for several angles φ (see legend). In the leftpanels, D F = D S = = . Ξ , as in the top panels of Fig. 4. In the right panels, D F =
10 as in the bottom panels of Fig. 4. In allcases, the ZBC decreases with increasing φ . tor, √ + h , and the transmitted wavevector, (unity in our con-ventions). Therefore, all the reflected and transmitted wavesabove the CB are traveling waves. As a result, we shouldconsider all possible incident angles and the upper limit ofEq. (23) is π/
2. Let us therefore focus on the nontrivial spin-up component of h G i . In Fig. 8 we separately plot the contri-butions to h G ↑ i from angles in the range above θ Ac (top panels)and below (bottom panels) for the field values and geometryin the left panels of Fig. 7, in particular D F =
12. Thesecontributions we will denote as h G ↑ ( E ) i above and h G ↑ ( E ) i below respectively. The h G ↑ ( E ) i below contributions, in the bottompanels of Fig. 8 are, for both h = . h = .
6, similarto the result for their total forward scattering counterpart (seeFig. 3 and the top right panel of Fig. 4). Of course, the angu-lar averaging leads to a smearing of the pronounced featuresoriginally in the forward scattering G . Qualitatively, the simi-larity comes from the propagating nature of all possible wavesexcept the transmitted waves below the CB when θ i < θ Ac .Therefore, the forward scattering G is just a special examplewith the incident angle perpendicular to the interface.In the subgap region, the contribution to h G ↑ ( E ) i above isvanishingly small although small humps appear when theexchange fields in the two F layers are non-collinear, i.e., φ , , π . These small humps are generated by the process of anomalous, equal-spin Andreev reflection. This processis possible in trilayers because, in a non-collinear magneticconfiguration, a spin up quasiparticle can Andreev reflect asa spin-up hole. This can be seen from the matrix form of theBdG equations, Eq. (2). The occurrence of anomalous An-dreev reflection leads to some important physics which weshall discuss in the next sub-subsection. One can see fromFig. 8, that when the exchange fields are strictly parallel oranti-parallel to each other, anomalous Andreev reflection doesnot arise.Above θ Ac , the conventional Andreev-reflected wave isevanescent and it does not contribute to h G ↑ i . When the bias isabove the saturated pair amplitude, contributions to h G ↑ i fromthe upper range are provided by both the transmitted wavesand by anomalous Andreev reflected waves. Recall that or-dinary transmitted waves are propagating when E is greaterthan the saturated pair amplitudes. We also see that h G ↑ i above decreases with increasing φ . At φ = ◦ , h G ↑ i is vanish-ingly small due to the effect of a large mismatch from theanti-parallel exchange field. Note also that the contributionfrom above θ Ac is less in the h = . h = .
6. Thisis mainly due to a smaller fraction of states at h = . θ Ac . On the other hand, the contri-bution from below θ Ac is larger in the h = . C on t r i bu ti on s t o 〈 G ↑ 〉 E below θ cA and h=0.3 0.1 0.2 0.3 0.4 above θ cA and h=0.3 φ =0 °φ =30 °φ =60 °φ =90 °φ =120 °φ =150 °φ =180 ° below θ cA and h=0.6above θ cA and h=0.6 FIG. 8. (Color online) Contributions (see text) to the spin-up angularly averaged conductance, h G ↑ i , from angular ranges above (top panels)and below (bottom panels) the Andreev critical angle θ Ac . Several values of φ are considered, as indicated. The top panel results at φ = ◦ arevanishingly small. The geometric and exchange field values are as in the left panels of Fig. 7. For the top panels, the plotted values at E = φ . For the bottom panels, their values at both E = E = φ . crease of h G ↑ i above and the decrease of h G ↑ i below from h = . h = .
3. Anomalous Andreev reflection
As we have seen, equal-spin (anomalous) Andreev reflec-tion (ESAR) can be generated when the magnetic configu-ration is non-collinear. We have previously shown that con-ventional Andreev reflection is forbidden when θ i > θ Ac = sin − ( √ (1 − h ) / (1 + h )). Thus, θ Ac vanishes in the half-metallic limit. In that case, conventional Andreev reflec-tion is not allowed for any incident angle θ i and the sub-gap h G ↑ i arises only from ESAR. For this reason, in thissub-subsection we present results for a trilayer structure thatconsists of one half-metal ( h =
1) and a much weaker( h = .
1) ferromagnet. The weaker ferromagnet servesthe purpose of generating ESAR. A somewhat similar ex-ample that has been extensively discussed in the literatureis that of half metal-superconductor bilayers with spin-flipinterface.
There the spin-flip interface plays the samerole as the weaker ferromagnet here. Another interestingphenomenon also related to ESAR is the induction of triplet pairing correlations in F-S structures.
To induce thistype of triplet pairing, F-S systems must be in a non-collinearmagnetic configuration such as F F S or F SF trilayers with φ , , π . Hence, the mechanism behind induced triplet pair-ing correlations is also responsible for ESAR and these twophenomena are closely related.In Fig. 9, we plot the h G i of this particular system for sev-eral φ . The geometrical parameters are again D F = D F =
12, and D S = h G i = h G ↑ i becausethe weight factor P ↓ = φ = ◦ and φ = ◦ the CB value is about 0.65 and, below the CB(in the subgap region), h G i vanishes because the conventionalAndreev reflection is completely suppressed and ESAR is notallowed in the collinear cases. For φ = ◦ and φ = ◦ , theCB is near 0.4 and 0.5 respectively and all of the subgap h G i is due to ESAR. The CB values for φ = ◦ , 90 ◦ , and 120 ◦ are0.15, 0.12, and 0.15. For these three angles, a FBCP clearlyforms, arising from the ESAR in the subgap region.To examine the conductance in the subgap region, which isin this case due only to ESAR, we choose the φ = ◦ angleand plot, in Fig. 10, the contributions to G (for this case G and h G i are very similar) from the reflected spin-up particle andthe reflected spin-up hole wavefunctions. The spin-down par-ticle and spin-down hole wavefunctions are evanescent and5 〈 G 〉 E φ =0 °φ =30 °φ =60 °φ =90 °φ =120 °φ =150 °φ =180 ° FIG. 9. (Color online) The angularly averaged conductance of F F Strilayers with exchange field h = h = . φ . See text for discussion. C on t r i b . t o G ( E , φ = ° ) E Contrib. from reflected spin-up particleContrib. from reflected spin-up holeTotal contrib. to G(E, φ =150 ° ) FIG. 10. (Color online) Contributions to G ( E , φ = ◦ ), computedfor the parameter values used in Fig. 9, from the spin-up quasiparticleand spin-up quasihole ESAR (see text for discussion). The total G isalso shown. do not contribute to the conductance. Thus, Eq. (22) reads G = + ( k −↑ / k + ↑ ) | a ↑ | − | b ↑ | . The quantities plotted are thesecond ((green) curve) and third ((red) curve, highest at theorigin) terms in this expression. The value of G is also plot-ted. One sees that the reflected ESAR amplitudes decay veryquickly above the CB. However, these reflected amplitudesare quite appreciable in the subgap region. In other words,the supercurrent in the subgap region contains signatures ofthe triplet correlations. This confirms the simple picture thatabove the CB the current flowing throughout the junction isgoverned by the transport of quasiparticles. However, belowthe CB it is dominated by ESAR.
4. Spin current densities and spin-transfer torques
Finally, we now report on spin-dependent transport quan-tities, including the spin current, the spin-transfer torques,and their connections to the local magnetization at finite bias. -0.4-0.2 0 0.2 0.4-40 -30 -20 -10 0 10 20 S z ( - ) Y δ m z ( - ) -3-2-1 0 1-40 -30 -20 -10 0 10 S y ( - ) -0.4 0 0.4 0.8 1.2-40 -30 -20 -10 0 δ m y ( - ) -0.4-0.2 0 0.2 0.4-40 -30 -20 -10 0 10 20 S x ( - ) E=0E=0.4E=0.8E=1.2E=1.6E=1.96 δ m x ( - ) FIG. 11. (Color online) The components of the spin current density, S x , S y , and S z , calculated from Eq. (31) are plotted vs. Y ≡ k F y forseveral values of the bias E ≡ eV (main panels). We have φ = ◦ , h = . D F = D F = D S = = Ξ . The F -S inter-face is at Y = -F interface at Y = −
30. Vertical linesat these interfaces in the top and bottom panels help locate the dif-ferent regions. Only the central portion of the Y range is included(see text). The ranges included depend on the component. The in-sets show the change in each component of the local magnetization, δ m ( E ) ≡ m ( E ) − m (0), also as a function of Y . The values of E areas in the main plot, the ranges included may be different. An objective here is to demonstrate that the conservation lawEq. (27) which in the steady state is simply: ∂∂ y S i = τ i , i = x , y , z , (36)is satisfied in our self consistent calculations for F F S trilay-ers. We consider these spin dependent quantities in a trilayerwith h = . φ = ◦ . Thus, the internal field in the outer electrodeF is along the z axis, while that in F is along x . The thick-nesses are D F = D F =
30, and D S = = Ξ .A set of results is shown in Fig. 11. There, in the three mainplots, we display the three components of the spin current den-sity, computed from Eqs. (31) and normalized as explainedbelow that equation. They are plotted as functions of the di-mensionless position Y ≡ k F y for several values of the bias V , E ≡ eV . The F -S interface and the F -F interface are locatedat Y = Y = −
30, respectively. For clarity, only the rangeof Y corresponding to the ”central” region near the interfacesis included in these plots: the shape of the curves deeper intoS or F can be easily inferred by extrapolation. From thesemain panels, one sees that the current is spin-polarized in the x -direction (the direction of the exchange field in F ) to theright of F -F interface, including the S region. Furthermore, S x is found to be a constant except in the F region, where itexhibits oscillatory behavior. This indicates the existence of anon-vanishing, oscillating spin-transfer torque in the F layer,as we will verify below. We also see that S x vanishes whenthe bias is less than the superconducting gap in bulk S ( E < S x with V is simi-lar to that of the ordinary charge current in an N-S tunnelingjunction with a very strong barrier where there is no currentuntil V > ∆ . This phenomenon is very different from whatoccurs in ordinary spin valves (F -F ), where the spin currentis not blocked below any finite characteristic bias.The S y component, along the normal to the layers, is shownin the middle main panel of Fig. 11. It depends extremelyweakly on the bias E . It is very small except near the interfacebetween the two ferromagnets but there it is about an order ofmagnitude larger than the other two components. Hence onlya somewhat smaller Y range is shown. Unlike the S x and S z components, S y does not vanish even when there is no bias ap-plied to the trilayer (the (red) curve in this panel). From theseobservations, one can infer that S y is largely derived from itsstatic part with only a very small contribution from the effectof finite bias. The emergence of a static spin current is dueto the leakage of the local magnetization m z into the F layerand of m x into the F layer. This explains why the static spincurrent S y is mostly localized near the F -F interface. The S z component (lower panel) is constant in the F region, as onewould expect. It oscillates in the F region, and vanishes in theS layer. As opposed to the S x component, S z is non-vanishing,although very small, when E < E >
1. The oscillatory behavior of S z , again, is relatedto the local spin-transfer torque as we will verify below.We can summarize the behavior of the spin current vector,in this φ = ◦ configuration, as follows: when E >
1, the spincurrent, which is initially (at the left side) spin-polarized in the + z direction, is twisted to the x direction under the action ofthe spin torques discussed below, as it passes through the sec-ond magnet, which therefore acts as a spin filter. The currentremains then with its spin polarization in the + x direction as itflows through the superconductor. Thus in this range of E thetrilayer switches the polarization of the spin current. On theother hand, when E <
1, the small z -direction spin-polarizedcurrent tunneling into the superconductor is gradually con- -0.2-0.15-0.1-0.05 0 0.05-40 -30 -20 -10 0 10 20 τ z ( - ) Y -0.15-0.1-0.05 0 0.05-40 -30 -20 -10 0 -15-10-5 0 5-50 -40 -30 -20 -10 0 10 τ y ( - ) -15-10-5 0 5-50 -40 -30 -20 -10 -0.15-0.1-0.05 0 0.05 0.1 0.15 0.2-50 -40 -30 -20 -10 0 10 τ x ( - ) E=0E=0.4E=0.8E=1.2E=1.6E=1.96 -0.1-0.05 0 0.05 0.1 0.15-50 -40 -30 -20 -10
FIG. 12. (Color online) The components of the spin-transfer torque τ ≡ m × h plotted vs. Y for several bias values. All parametersand geometry are as in Fig. 11. Vertical lines, denoting interfaces,are in the top and bottom panels. The insets show (for bias E = . verted into supercurrent and becomes spin-unpolarized.In the insets of the three panels of Fig. 11, we illustratethe behavior with bias of the corresponding component of thelocal magnetization as it is carried into S. Specifically, weplot the components of the vector difference between the localmagnetization with and without bias, δ m ( E ) ≡ m ( E ) − m ( E = Y . The range of Y is chosen to display thesalient aspects of the behavior of this quantity, and it is not thesame as in the main plots, nor is it the same for each compo-nent. The bias values are the same as in the main plots, how-ever. The magnetizations are computed from Eqs. (29) andnormalized in the usual way, as discussed below Eqs. (31).7In these units, and at h = . m in the magnetic layers is about 0.15. This scaleshould be kept in mind.The behavior of the x component is nontrivial in the F andS regions, and the corresponding Y range is included in thetop panel inset. When the applied bias is below the bulk Sgap value, δ m x ( E ) penetrates into the S layer with a decaylength ∼ Ξ =
50. This decay length is much longer than thatfound for the static magnetization, m (0) . When the bias isabove the gap, δ m x ( E ) penetrates even more deeply into theS layer, with a clearly very different behavior than for E <
1. This long-range propagation is of course consistent withthe behavior of S x , as S x , the spin current polarized in the x direction, appears only when E >
1. The magnitude of δ m y ismuch smaller than that of δ m x or δ m z . It peaks near the F -F interface and that range of Y is emphasized in the middle inset.Its overall scale monotonically increases with increasing bias.It damps away from the interface in an oscillatory manner. Asto δ m z , which can conveniently be plotted in the same Y range,it decays with a very short decay length and oscillates in F .The overall damped oscillatory behavior of δ m y and δ m z inthe F region reflects the precession, as a function of position,of the spin density around the local exchange field that pointstoward the + x direction. This phenomenon is well known inspin-valves. The oscillation periods for δ m y , δ m z , S x , and S z are very similar and of the order of 1 / ( hk FS ).Next, we investigate the spin-transfer torque, τ ≡ m × h .This quantity, computed from the normalized values of h and m , is plotted as a function of position in Fig. 12 for the samesystem as in Fig. 11. Results are shown for each of its threecomponents in the main panels of the figure. One sees that atzero bias, E =
0, both τ x and τ z vanish identically. In the F layer, τ x increases in magnitude with increasing E . It vanishesin F and in S. The behavior of τ z is, as one would expect, theconverse: it vanishes in F and S, and its magnitude increasesin F . The oscillatory behavior of τ x and τ z is consistent, as weshall see below, with the results for S x and S y . The componentnormal to the layers, τ y , is nonvanishing only near the F -F interface, although its peak there attains a rather high value,nearly two orders of magnitude larger than the peak value ofthe other components. It is independent of bias, consistentwith the behavior of S y .In the insets, we verify, for each component, that Eq. (36)is satisfied, that is, that our self-consistent methods strictlypreserve the conservation laws in this nontrivial case. (Wehave already mentioned that we have verified that the chargeor particle current are independent of y ). We specifically con-sider the bias value E = . x component ofthe spin-transfer torque, τ x (blue dashes), taken from the cor-responding main plot, and the derivative of the spin current, ∂ S x /∂ Y (blue circles), obtained by numerically differentiatingthe corresponding result in the top panel of Fig. 11. Clearly,the curves are in perfect agreement. (One can easily checkthat with the normalizations and units chosen there should beno numerical factor between the two quantities). In the secondpanel, the same procedure is performed for the y component,although in this case, because of the very weak dependence of both S y and τ y on bias, the value of the latter is hardly rel-evant. Nevertheless, despite the evident difficulty in comput-ing the numerical derivative of the very sharply peaked S y , theagreement is excellent. For τ z , its vanishing in the F regionis in agreement with the constant spin current in that layer.The conservation law Eq. (36) is verified in the inset for thiscomponent, again at bias E = .
6. Just as for the x compo-nent, the dots and the line are on top of each other. Thus theconservation law for each component is shown to be perfectlyobeyed.The results of this sub-subsection can be summarized as fol-lows: the finite bias leads to spin currents. As opposed to theordinary charge currents, these spin currents are generally notconserved locally because of the presence of the spin-transfertorques which act as source terms and are responsible for thechange of spin-density. But a self-consistent calculation must still contain exactly the correct amount of non-conservation,that is, Eq.(36) must be satisfied. It is therefore of fundamen-tal importance to verify that it is, as we have. IV. CONCLUSIONS
In summary, we have investigated important transport prop-erties of F F S trilayers, including tunneling conductancesand spin transport. To properly take into account the prox-imity effects that lead to a spatially varying pair potential,we have incorporated a transfer matrix method into the BTKformalism. This allows us to use self-consistent solutions ofthe BdG equations. This technique also enables us to com-pute spin transport quantities including spin transfer torqueand spin currents. We have shown that in F-S bilayers theself-consistent calculations lead to conductances at the super-conducting gap that increase with the Fermi wavevector mis-match whereas non-self-consistent ones predict they are in-sensitive to this parameter. In F F S trilayers, we have foundthat the critical bias CB (where tunneling conductance curvesdrop) for different relative magnetization angles, φ , dependson the strength of the superconducting order parameter nearthe interface. The angular dependence of the critical bias re-flects that of the transition temperatures T c , which are usuallynonmonotonic functions of φ . For forward scattering in theseF F S trilayers, we found that the dependence of the zero biasconductance peak (ZBCP) on φ is related to both the strengthof the exchange fields and the thickness of the F layers. Thisremarkable behavior can be explained via quantum interfer-ence effects. At the resonance minimum, the ZBCP drops sig-nificantly and monotonically from φ = ◦ to φ = ◦ . On theother hand, the φ dependence of the ZBCP is very weak whenit is at its resonance maximum. For asymmetric cases where h , h , we found that the ZBCP is a nonmonotonic func-tion of φ with its value at φ = π/ h G i , and found that features of resonance effectsare then somewhat washed out due to the averaging. However,by studying h G i in the subgap regions, we found that anoma-lous (equal spin) Andreev reflection (ESAR) arises when φ corresponds to noncollinear orientations. The emergence of8ESAR is correlated with the well-known induced triplet pair-ing correlations in proximity coupled F-S structures. Whenthe outer magnet is a half metal, the h G i signatures arisechiefly from the process of ESAR. We have also studied thebias dependence of the spin currents and spin transfer torquesand their general behavior in F F S trilayers with φ = ◦ (the exchange fields in F and F point toward the z and x di-rections, respectively). The spin current components are ingeneral non-conserved quantities. The S z component, parallelto the local exchange field in the F layer, does not change inthe F region but shows damped oscillatory behavior in the F layer and eventually vanishes in the S region. However, S x isa constant throughout the F and S regions and oscillates in F layers. We found that S y (the component normal to the layers)depends very weakly on the bias, and thus its spatial depen-dence arises largely from a static effect. The bias dependenceof S x in the S region is very similar to that of the tunnelingcharge current in normal/superconductor systems with highbarriers: S x vanishes in the subgap regions and arises right above the gap. The behavior of m is consistent with that of S .We found that m x , parallel to the local exchange fields in F ,spreads out over the S regions when the bias is larger than thesuperconducting gap. We have also investigated the bias de-pendence of the spin transfer torques, and we have carefullyverified that the appropriate continuity equation for the spincurrent is strictly obeyed in our self-consistent approach. Ourmethod can be extended to include the effects of interfacialscattering and wavevector mismatch. It can also be used forfurther study of the intricate phenomena associated with spintransport in these systems. ACKNOWLEDGMENTS
Portions of this work were supported by IARPA grant No.N66001-12-1-2023. CTW thanks the University of Minnesotafor a Dissertation Fellowship. The authors thank I. Krivorotov(Irvine) for helpful discussions. ∗ [email protected] † [email protected]; Also at Minnesota Supercomputer Institute,University of Minnesota, Minneapolis, Minnesota 55455 ‡ [email protected] I. ˇZuti´c, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. , 323(2004). A. I. Buzdin, Rev. Mod. Phys. , 935 (2005). E. A. Demler, G. B. Arnold, and M. R. Beasley, Phys. Rev. B ,15174 (1997). K. Halterman and O. T. Valls, Phys. Rev. B , 014509 (2001). K. Halterman and O. T. Valls, Phys. Rev. B , 224516 (2002). Buzdin, A. I., and M. Y. Kuprianov, Pisma Zh. Eksp. Teor. Phys. , 1089-1091 [JETP Lett. , 487-491 (1990)]. K. Halterman and O. T. Valls, Phys. Rev. B , 104516 (2004). A. F. Andreev, Sov. Phys. JETP , 1228 (1964). G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Phys. Rev. B ,4515 (1982). S. Kashiwaya, Y. Tanaka, M. Koyanagi, and K. Kajimura, Phys.Rev. B , 2667 (1996). M. J. M. de Jong and C. W. J. Beenakker, Phys. Rev. Lett. , 1657(1995). I. ˇZuti´c and O. T. Valls, Phys. Rev. B , 6320 (1999). I. ˇZuti´c and O. T. Valls, Phys. Rev. B , 1555 (2000). I. I. Mazin, Phys. Rev. Lett. , 1427 (1999). R. J. Soulen Jr., J. M. Byers, M. S. Osofsky, B. Nadgorny, T. Am-brose, S. F. Cheng, P. R. Broussard, C. T. Tanaka, J. Nowak, J. S.Moodera, A. Barry, and J. M. D. Coey, Science , (5386):85-88(1998). S. K. Upadhyay, A. Palanisami, R. N. Louie, and R. A. Buhrman,Phys. Rev. Lett. , 3247 (1998). P. Raychaudhuri, A. P. Mackenzie, J. W. Reiner, and M. R.Beasley, Phys. Rev. B , 020411 (2003). P. Chalsani, S. K. Upadhyay, O. Ozatay, and R. A. Buhrman, Phys.Rev. B , 094417 (2007). S. Hacohen-Gourgy, B. Almog, and G. Deutscher, Appl. Phys.Lett. 92, 152502 (2008). F.S. Bergeret, A.F Volkov, and K.B. Efetov, Phys. Rev. Lett. ,3140 (2001); Phys. Rev. B , 064513 (2003); Rev. Mod. Phys. , 1321-1373 (2005). F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Appl. Phys. A ,599 (2007) J. Wang, M. Singh, M. Tian, N. Kumar, B. Liu, C. Shi, J. K. Jain,N. Samarth, T. E. Mallouk, and M. H. W. Chan, Nat. Phys. , 389(2010). F. H¨ubler, M. J. Wolf, T. Scherer, D. Wang, D. Beckmann, and H.v. L¨ohneysen, Phys. Rev. Lett. , 087004 (2012). K. Halterman, O. T. Valls, and P. H. Barsic, Phys. Rev. B ,174511 (2008). R. S. Keizer, S. T. B. Goennenwein, T. M. Klapwijk, G. Miao, G.Xiao, and A. Gupta, Nature , 825 (2006). J. W. A. Robinson, J. D. S. Witt, and M. G. Blamire, Science ,59 (2010) T. S. Khaire, M. A. Khasawneh, W. P. Pratt, Jr., and N. O. Birge,Phys. Rev. Lett. , 137002 (2010). K. Halterman and O. T. Valls, Phys. Rev. B , 104502 (2009). M. Eschrig and T. L¨ofwander, Nature Physics , 138 (2008). K. Halterman, P. H. Barsic, and O. T. Valls, Phys. Rev. Lett. ,127002 (2007). J. Linder, T. Yokoyama, and A. Sudbø, Phys. Rev. B , 224504(2009). C. Visani, Z. Sefrioui, J. Tornos, C. Leon, J. Briatico, M. Bibes,A. Barth´el´emy, J. Santamar´ıa, and Javier E. Villegas, Nature Phys. , 539 (2012). Z. P. Niu, Europhys. Lett. Y.-Q. Ji, Z.-P. Niu, C.-D. Feng, and D.-Y. Xing, Chinese Phys.Lett. , 691 (2008) C. D. Feng, Z. M. Zheng, R. Shen, B. Wang, and D. Y. Xing, Phys.Rev. B , 224510 (2010) A. A. Jara, C. Safranski, I. N. Krivorotov, C.-T. Wu. A. N. Malmi-Kakkada, O. T. Valls, and K. Halterman, Phys. Rev. B , 184502(2014). P.V. Leksin, N. N. Garif’yanov, I. A. Garifullin, Ya.V. Fominov,J. Schumann, Y. Krupskaya, V. Kataev, O. G. Schmidt, and B.B¨uchner, Phys. Rev. Lett. , 057005 (2012). V. I. Zdravkov, J. Kehrle, G. Obermeier, D. Lenk, H.-A. Krug vonNidda, C. M¨uller, M. Yu. Kupriyanov, A. S. Sidorenko, S. Horn,R. Tidecks , and L. R. Tagirov, Phys. Rev. B , 144507 (2013). C.-T. Wu, K. Halterman, and O. T. Valls, Phys. Rev. B , 014523(2012). Q. Cheng and B. Jin, Physica C: Superconductivity , 29(2012). S. Kashiwaya, Y. Tanaka, N. Yoshida, and M. R. Beasley, Phys.Rev. B , 3572 (1999). T. Yamashita, S. Takahashi, H. Imamura, and S. Maekawa, Phys.Rev. B , 172509 (2002). L. Berger, Phys. Rev. B , 9353 (1996). E. B. Myers, D. C. Ralph, J. A. Katine, R. N. Louie, and R. A.Buhrman, Science , 867 (1999). M. Boz˘ovi´c and Z. Radovi´c, Phys. Rev. B , 134524 (2002). F. Romeo and R. Citro, Phys. Rev. B , 024531 (2011). J. Linder and A. Sudbø, Phys. Rev. B , 134509 (2007). R. Grein, T. L¨ofwander, G. Metalidis, and M. Eschrig, Phys. Rev.B , 094508 (2008). P. H. Barsic and O. T. Valls Phys. Rev. B , 014502 (2009). P. G. de Gennes,
Superconductivity of Metals and Alloys (Addison-Wesley, Reading, MA, 1989). A. Spuntarelli, P. Pieri, and G. C. Strinati, Physics Reports ,111 (2010). J.-X. Zhu, C.S. Ting, Phys. Rev. B , 1456 (2000). G. Baym and L.P. Kadanoff, Phys. Rev. , 287 (1961). With a current present, one cannot adopt the expedient of choosingreal wavefunctions. P.F. Bagwell, Phys. Rev. B , 6841 (1993). F. Sols and J. Ferrer, Phys. Rev. B , 15913 (1994). J. Sanchez-Canizares and F. Sols, Phys. Rev. B , 531 (1997). Y.V. Fominov, A.A. Golubov, T.Y. Karminskaya, M.Y.Kupriyanov, R.G. Deminov, and L.R. Tagirov, JETP Letters , 308 (2010). S. Oh, D. Youm, and M. R. Beasley, Appl. Phys. Lett. , 2376(1997). M. Eschrig, J. Kopu, J. C. Cuevas, and G. Schr¨on, Phys. Rev.Lett. , 137003 (2003). J. Linder, M. Cuocco, and Asle Sudbø, Phys. Rev. B , 174526(2010). D.C. Ralph and M.D. Stiles, J. Magn, Magn. Mater.320