First principles supercurrent calculation in realistic magnetic Josephson junctions
Hervé Ness, Ivan A. Sadovskyy, Andrey E. Antipov, Mark van Schilfgaarde, Roman M. Lutchyn
FFirst principles supercurrent calculation in realistic magnetic Josephson junctions
Hervé Ness, ∗ Ivan A. Sadovskyy, Andrey E. Antipov, Mark van Schilfgaarde,
1, 3 and Roman M. Lutchyn Department of Physics, King’s College London, Strand Campus, London WC2R 2LS, UK Microsoft Quantum, Microsoft Station Q, University of California, Santa Barbara, California 93106, USA National Renewable Energy Laboratory, Golden, Colorado 80401, USA
We investigate the transport properties of magnetic Josephson junctions. In order to capture re-alistic material band structure effects, we develop a numerical method combining density functionaltheory and Bogoliubov-de Gennes model. We demonstrate the capabilities of this method by study-ing Nb/Ni/Nb junctions in the clean limit. The supercurrent through the junctions is calculated asa function of the ferromagnetic Ni thickness, magnetization, and crystal orientation. We identifytwo generic mechanisms for the supercurrent decay with ferromagnet thickness: (i) large exchangesplitting may gap out minority or majority carriers leading to the suppression of Andreev reflectionin the junction, (ii) loss of synchronization between different modes due to the significant dispersionof the quasiparticle velocity with the transverse momentum. Our results are in good agreementwith recent experimental studies of Nb/Ni/Nb junctions. The present approach opens a path formaterial composition optimization in magnetic Josephson junctions and superconducting magneticspin valves.
Keywords: magnetic Josephson junction, π -junction, JMRAM, ab initio calculations, critical current,scattering theory. I. INTRODUCTION
Coherent quantum tunneling of Cooper pairs througha thin barrier is one of the first examples of macroscopicquantum coherent phenomena. Predicted by Josephsonmore than 50 years ago [1], it has important applicationsin quantum circuits used in metrology, quantum sensingand quantum information processing [2].Most of the previous studies focused on conventionalJosephson junctions (JJs) consisting of two s -wave su-perconductors (S) that are connected by an insulating(I) or a normal (N) region [3, 4]. The flow of supercur-rent through a JJ depends on the superconducting phasedifference φ between two superconductors and, in gen-eral, is characterized by the current-phase relationship J ( φ ) (CPR). In conventional JJs CPR should be peri-odic with π , I ( φ ) = I ( φ + 2 π ) which follows from theBCS theory [5]. This result is a manifestation of a e charge of Cooper pairs, and is used in metrology to mea-sure electron charge. Time-reversal symmetry requiresthat I ( φ ) = − I ( − φ ) which imposes a constraint that thesupercurrent should be zero for φ = πn where n is integer.In general, CPR can be expanded in Fourier harmonics, I ( φ ) = (cid:80) n I n sin( nφ ) .In many cases, however, CPR is well approximated bythe first harmonic I ( φ ) ≈ I c sin( φ ) with I c being the max-imum supercurrent that can flow through the junction,i.e. the critical current. At a microscopic level, the super-current through a short SNS junction is determined bybound states forming in the constriction due to Andreevreflection at the NS interfaces. In the Andreev reflectionprocess an incident electron-like quasiparticle with spin ↑ gets reflected at the NS interface as a hole-like quasi- ∗ [email protected] particle with spin ↓ and a Cooper pair is emitted into thecondensate. When time-reversal symmetry is not broken,electrons and holes propagate with the same velocity inthe normal region. In this case no phase shift accumu-lates between this pair of quasiparticles along their tra-jectories in the N region, and the sign of I c is fixed. When I c > we refer to this case as -junction.In a magnetic Josephson junction (MJJ), exchangesplitting breaks time-reversal symmetry and leads toan interesting interplay of superconductivity and mag-netism [4, 7–10]. In superconductor-ferromagnet-superconductor (SFS) junctions [Fig. 1(a)] the corre-lated quasi-particles and quasi-holes forming Andreevbound states propagate through the junction under theexchange field of the ferromagnet (F). In many ferromag-nets, such as Fe or Ni, the exchange splitting is large (ofthe order of eV) and significantly perturbs the band struc-ture of a metal and, consequently, significantly modifiesFermi velocities of minority (spin ↓ ) and majority (spin ↑ )carriers. Strong time-reversal symmetry breaking leadsto the appearance of characteristic superconducting cor-relations with oscillatory dependence determined by thedifference in wave numbers, k ↑ F − k ↓ F [11, 12]. This ef-fect opens a possibility for the supercurrent reversal asa function of the thickness of the ferromagnetic region,the so-called Josephson π -junction [13, 14]. The correla-tion between the phase shift of the supercurrent and themagnetization provides a possibility for realizing mag-netic spin valves, see Fig. 1(b), which may have promisingnovel applications for cryogenic superconducting digitaltechnologies [6, 15, 16]. Understanding the microscopicphysics of MJJs is of great scientific interest as well astechnological importance. – π transitions in SFS junc-tions have been extensively studied experimentally sincethe early 2000’s and have been observed in different ma-terial systems [14, 17–23]. While qualitatively these ob-servations are consistent with the previous phenomeno- a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n S u p e r c o d u c t o r F e rr o m a g n e t S u p e r c o d u c t o r S u p e r c o d u c t o r F e rr o m a g n e t N o r m a l m e t a l F e rr o m a g n e t S u p e r c o d u c t o r a bNi(111)Nb(110) Nb(110)c zx y w Figure 1: (a) Schematic view of the SFS junction. (b) SFNFSjunction. Arrows indicate possible magnetization of ferromag-nets. The supercurrent through the spin-valve JJ depends onthe relative magnetization of the ferromagnets, which gov-erns the properties of the Josephson magnetic random-accessmemory (JMRAM) [6]. (c) Ball-and-stick representation ofthe Nb(110)/Ni(111)/Nb(110) junction with 5 layers of Ni.Nb atoms are light grey, Ni atoms are blue. The top andbottom atomic planes of Nb(110) are repeated periodicallyin the z -direction to create the semi-infinite Nb leads of thejunction through which the current flows. Periodic boundaryconditions are used in the xy -plane. The corresponding re-ciprocal space defines two-dimensional k (cid:107) = ( k x , k y ) vectors,i.e. the transverse modes, used in the calculations. logical theories [24–40] the microscopic understanding ofthe current-phase relationship is still lacking: the originof the supercurrent suppression with junction thicknessand the appearance of different phase offsets remains un-clear. The supercurrent suppression is often associatedwith presence of disorder in the ferromagnetic region.However, significant supercurrent suppression can alsoappear in relatively clean metals (e.g., Ni) whose meanfree path is larger than the junction thickness. The dif-ferent phase offsets may originate from the NS interfaceproperties or from inhomogeneous magnetization in theferromagnet.In this paper we develop a microscopic theory for thesupercurrent in realistic MJJs. We use a combinationof density functional theory (DFT) and Bogoliubov-deGennes (BdG) model to investigate the – π transitionin realistic material stacks of Nb/Ni/Nb junctions in theclean limit. We also study the suppression of the crit-ical current in the MJJs and identify two microscopicmechanisms for its suppression, both a consequence ofthe asymmetry of the majority and minority carriers inthe F region.First, there is an asymmetry in the structure of the k E V ex E ( i ) F E ( ii ) F k F k F ba S F S w k k
Hole ( , )Electron (+ , )
Figure 2: (a) Simplified band structure of a ferromagnet,with majority and minority bands split by V ex . Fermi level E (i) F corresponds to small V ex . k ↑ F and k ↓ F are the Fermi mo-menta for majority and minority carriers, respectively. E (ii) F corresponds to large exchange splitting, where the minorityband is pushed above E (ii) F . Thus, the propagation of minorityquasiparticles is characterized by an imaginary momentum k ↓ F and is suppressed. (b) Supercurrent in SFS junction is carriedby Andreev bound states localized in the junction. Solid redline represents quasi-classical trajectory corresponding to anAndreev bound state. The spectrum of Andreev states de-pends on the relative superconducting phase difference acrossthe junction as well as the phase, δϕ = | k ↑ F − k ↓ F | w/ , accu-mulated due to the difference of Fermi momenta for majorityand minority carriers. Note that in scenario (ii) the propa-gation of minority carriers is suppressed leading to an overallexponential decay of the supercurrent with w . This is to becontrasted with the normal transport through the junction. Fermi surface, see Fig. 2(a). For certain bands and mo-menta, the Fermi surface present in one spin channel maybe absent in the other. This is typical in ferromagneticmaterials like Fe, Co, and Ni because the bandwidth for d -electrons is relatively small and is often comparable tothe exchange splitting. As a result, the wave number ofone of the constituent quasi-particles forming Andreevbound states in MJJ becomes imaginary and the super-current becomes suppressed. We label this scenario asmechanism (i). In a second mechanism (ii), we showthere is a dephasing of a harmonic signal originatingfrom different Fourier components of the supercurrentdue the Fermi velocity dispersion, see Fig. 3(a). Boththese mechanisms lead to an exponential suppression ofthe supercurrent, which was previously believed to occurdue to the presence of disorder in the magnetic region.We show that band structure effects are important andmay be even dominant in many cases.The paper is organized as follows. We describe ourfindings and summarize our main results in Sec. II.In Sec. III we discuss the numerical method we devel-oped to perform first principles supercurrent calculations.Detailed discussion of the main results is presented inSec. IV. We conclude with Sec. V. Some technical detailsare relegated to Appendices A–C. II. QUALITATIVE DISCUSSION AND MAINRESULTS
In this section we describe basic concepts for the su-percurrent flow in MJJs and summarize our results. Ourmain qualitative conclusions are supported by micro-scopic calculations for Nb/Ni/Nb MJJs. Ni appears tohave fairly long mean free path l MFP ≈ Å (see estima-tions in Appendix A), which is comparable or larger thanthe typical thickness of the ferromagnet used in recentexperiments [21, 22]. Therefore, the motion of quasipar-ticles in Nb/Ni/Nb junction is quasi-ballistic, and ourmethod is applicable to this system. Most of this paperfocuses on clean Nb/Ni/Nb junctions.First it is illuminating to consider a toy model, a one-dimensional SFS junction, and calculate the supercurrentin such a system for different Fermi energies, see Fig. 2(a).Using the results of Ref. 41, one finds that both majorityand minority spin bands are both occupied in the limit V ex /E F (cid:28) [i.e. scenario (i) in Fig. 2(a)] the supercur-rent does not decay with ferromagnet thickness at zerotemperature, I ( φ ) = 2 e ∆ (cid:126) cos δϕ sin φ , < φ < π − δϕ, − sin δϕ cos φ , π − δϕ < φ < π + 2 δϕ, − cos δϕ sin φ , π + 2 δϕ < φ < π. (1)Here perfect interface transparency T is assumed, T ≈ .The phase offset δϕ originates from the Fermi momen-tum difference of a quasi-particle and a quasi-hole form-ing Andreev bound state in the junction, see Fig. 2(b),and is given by δϕ = | k ↑ F − k ↓ F | w/ ≈ V ex w/v F . The -and π -junction regimes can be clearly identified at δϕ = 0 and δϕ = π/ , respectively. At the intermediate values < δϕ < π/ , the CPR is anharmonic which is a genericfeature at – π transition as shown below. In the lowtransparency regime, T (cid:28) , one would expect qualita-tively similar results with the maximal critical currentbeing suppressed I c ∼ ( e ∆ / (cid:126) ) T but still independent ofthe ferromagnet thickness, w .In the case of large exchange splitting V ex /E F (cid:29) [scenario (ii) in Fig. 2(a)] the minority band may becomeunoccupied. Given that minority carriers are gapped outand their propagation through the junction is suppressed,the supercurrent decays exponentially with w [41], I ( φ ) ≈ e ∆ (cid:126) exp( − κw ) (cid:104) − E F V z sin ( kw ) (cid:105) sin φ. (2)Here κ = (cid:112) m ∗ ( V ex − E F ) / (cid:126) and k = (cid:112) m ∗ ( V ex + E F ) / (cid:126) with m ∗ being effective electron mass. One may notice the drastic difference between normal-state and super-conducting transport in this case — the former is weaklyaffected (because majority and minority contributionsare additive) whereas the supercurrent is strongly sup-pressed. In this case, the measurement of normal-statejunction resistance does not necessarily predict the mag-nitude of the supercurrent through the junction.We now generalize above results for the realistic three-dimensional (3D) geometry and material composition ofthe MJJ. In the clean limit, the supercurrent I ( φ ) in theshort-junction limit ( w much smaller than the coherencelength of the superconductor) is obtained from the spec-trum of the Andreev bound states ε ν ( φ, k (cid:107) ) localized inthe junction [3] which now also depends on the parallelmomentum k (cid:107) . The supercurrent density J ( φ ) per junc-tion area A is J ( φ ) = I ( φ ) /A . For the junction withperiodic atomic structure in xy -plane the supercurrentdensity at zero temperature is given by J ( φ ) = − e (cid:126) (cid:90) BZ d k (cid:107) (2 π ) (cid:88) ν> ∂ε ν ( φ, k (cid:107) ) ∂φ , (3)where the k (cid:107) integration is performed over the Brillouinzone (BZ) of the corresponding surface supercell of area, A , and the sum is carried over positive quasiparticle en-ergies, ε ν ( φ, k (cid:107) ) > . Note that we use spin-resolved ε ν and therefore omit spin prefactor 2 in Eq. (3). Thederivative is periodic in φ and can be represented as aFourier series − e (cid:126) ∂ε ν ( φ, k (cid:107) ) ∂φ = (cid:88) n (cid:62) I nν ( k (cid:107) ) sin (cid:2) nφ + δϕ nν ( k (cid:107) ) (cid:3) , so that Eq. (3) can be written as J ( φ ) = (cid:90) BZ d k (cid:107) (2 π ) (cid:88) ν> (cid:88) n (cid:62) I nν sin (cid:2) nφ + δϕ nν ( k (cid:107) ) (cid:3) . As we will show below, away from – π transition the su-percurrent is dominated by the first ( n = 1 ) harmonic.Therefore, we focus henceforth on the first harmonic con-tribution and drop n index in the following discussion.Next, one may notice that the supercurrent amplitudes I ν ( k (cid:107) ) and phase offsets δϕ ν ( k (cid:107) ) depend on the parallelmomentum k (cid:107) . One may include this dependence anddefine an effective Fermi energy E F ( k (cid:107) ) that counts theenergy corresponding to k (cid:107) in each band of the ferromag-net from the bottom of the band. Depending on E F ( k (cid:107) ) and V ex , either scenario (i) or (ii) of Fig. 2(a) may berealized.Thus, it is important to compare the exchange splittingwith the bandwidth of d -character states in transitionmetals to make sure that a perturbation theory in V ex isjustified. Assuming it is the case, one can estimate thephase offset δϕ ν ( k (cid:107) ) by expanding in exchange splittingto find δϕ ν ( k (cid:107) ) ≈ V ex /v z ( k (cid:107) ) . (4) w [Å]32101234 J [ A / m ] w [Å]0.00.51.01.5 c d eba G / A [ m m ] G / AG / AG / AG / A J [ A / m ] J ( ) J sin( ) J sin(2 ) Figure 3: (a) First Fourier harmonic J of the supercurrent density [Eq. (3), black circles] and its fit [Eq. (14), solid blackline] as a function of Ni layer thickness, w , calculated for the Nb(110)/Ni(111)/Nb(110) junctions shown in Fig. 1(c). Greensemitransparent curves correspond to j fit1 ( k (cid:107) ) for individual k (cid:107) . (b) Normal state conductance per unit of area as a function of w for majority ( G ↑ ) and minority ( G ↓ ) spins [Eq. (10)] as well as G = G ↑ + G ↓ . G does not depend on w and is approximatedby the single value (cid:104) G (cid:105) . (c)–(e) Supercurrent as a function of phase difference φ for (c) 4 atomic layers of Ni (strong -junctionregime), (d) 8 layers (intermediate regime), and (e) 13 layers (strong π -junction regime). In the - and π -junction regimes, the J component dominates. In the intermediate regime higher-order terms may prevail. In general, the dependence of v z on k (cid:107) is compli-cated, especially in spd -transition metals. The combi-nation of complicated amplitude and phase offset de-pendence on k (cid:107) leads to a non-trivial supercurrent de-pendence on the ferromagnet thickness w . As shownin Fig. 3(a), the critical current decays with w for theNb(110)/Ni(111)/Nb(110) junctions. We analyze thedetails of the decay and perform an exponential fit inSec. IV. At small thicknesses, below 30 Å, this decayoriginates from the evanescent modes corresponding togapped out minority or majority carriers which cannotpropagate through the junction, see Table I. This is themechanism (i) noted at the end of Sec. I. At larger thick-nesses (i.e. w (cid:38) Å) a loss of synchronization betweendifferent modes due to the dispersion of v z ( k (cid:107) ) becomesimportant. This second mechanism (ii) has been previ-ously discussed in the literature [7, 13] under assumptionsof a single spherical Fermi surface and a small uniformexchange splitting V ex in the magnetic region. Withinthese assumptions, one finds that critical current shoulddecay algebraically with the thickness w of a magnetic layer [13]. However, as we show below most of theseassumptions do not apply to realistic SFS junctions in-volving transition metals. Thus, in order to understandCPR in realistic MJJs, one needs to use accurate ab ini-tio methods, which capture the physical effects describedabove.To make a connection between the decay seen in Fig. 3and the mechanisms responsible for it, Fig. 4(a) presentsthe energy band structure of bulk Ni. It is probably thehighest fidelity band structure available: it very closelyreproduces ARPES data in both majority and minorityspin bands, with exchange splitting V ex = 0 . eV [42], andshould be an excellent predictor of the real Fermi surfaceand velocities in Ni.Here we use this potential to analyze the bulk Ni Fermisurface [Fig. 4(b)] and Fermi velocities (Table I). Count-ing from the bottom s -band, bands of Ni d characterare bands 2 to 6. These bands are nearly full: onlymajority band 6 ↑ and minority bands ↓ – ↓ cross theFermi level E F . Only band 6 has both majority andminority carriers at the Fermi surface. Bands 6 ↑ and Majority ( ↑ ) Minority ( ↓ ) v min F v max F (cid:104) v F (cid:105) (cid:112) (cid:104) v F (cid:105) ρ ( E F ) ∆ E m ∗ /m e κ v min F v max F (cid:104) v F (cid:105) (cid:112) (cid:104) v F (cid:105) ρ ( E F ) Band [ ] [ ] [ ] [ ] [eV − Å − ] [eV] [Å] [ ] [ ] [ ] [ ] [eV − Å − ]2 − . − . − . − . v F for the majority and minority Fermisurfaces in bulk Ni. ρ ( E F ) is the density of states at the Fermi level, ∆ E is the energy of the valence band maximum relative to E F for the majority bands not crossing E F , m ∗ ( m e ) is the effective (bare) electron mass, and κ estimates the decay exponentof the evanescent mode for a given band at E F . a b E [ e V ] Γ Γ L 6 356 Γ L X − − MajorityMinorityAverage ab initio
ARPES
Figure 4: (a) Electronic band structure of bulk Ni (solid lines)calculated from first principles including many-body effects,see Ref. 42. It is the highest fidelity available and is very closeto ARPES data (diamonds) in both majority (red) and minor-ity (green) spin bands, with exchange splitting V ex = 0 . eV.(b) Majority (solid line) and minority (dashed line) Fermisurfaces of bulk Ni in the k (cid:107) plane (with k z = 0 ) correspond-ing to the (111) plane direction used for the stacking of theNb/Ni/Nb junctions, as discussed in the text. Axes corre-spond to two Γ -L lines: the Fermi surface has a three-foldsymmetry in the entire plane. Majority band ↑ is depictedas solid red-blue line with color interpolating between redand blue depending on the Fermi velocity (cid:126) − ∂E/∂k , whichranges between × m/s (red) and × m/s (blue). Mi-nority bands are shown by dashed lines: ↓ (blue), ↓ (cyan), ↓ (green), and ↓ (red). ↓ have roughly the same shape and the energy split-ting is approximately constant and equal to V ex (see theright panel of Fig. 4 in Ref. 42). Beyond this, however,the correspondence between the Ni band structure anda simple parabolic band structure deviate substantially,in two ways that critically affect the analysis. First,the Fermi velocity, v F = (cid:126) − ( ∂E/∂k ) | k = k F , is not fixedeven for a single pocket: it varies in band ↑ by a fac-tor of 2 [see Fig. 4 and Table I]. Accordingly the split-ting k ↑ F − k ↓ F between bands ↑ and ↓ varies by factor oftwo as expected from the twofold variation in v F . Sec- ond, bands ↓ – ↓ have no majority counterpart at E F ,indicating that the wave number of bands ↑ – ↑ is com-plex. This is the origin for the exponential decay in sce-nario (ii) in Fig. 2(a) as noted above: a large portionof Andreev levels are carried by Cooper pairs made ofsingle-particle wave functions with one or both of k ↑ F and k ↓ F having an imaginary component. The magnitude of Im k depends on the particular mode and k (cid:107) leading toa distribution of decay exponents. The slowest decay ineach of these evanescent modes can be estimated fromthe distance ∆ E of the closest approach to E F and theeffective mass m ∗ /m e , using (cid:126) k / m ∗ = ∆ E and de-cay κ = 2 π/ Im( k min ) , see Table I. ( m ∗ /m e is found to behighly anisotropic, so only the effective transport mass m ∗ = 3 [1 /m ∗ + 1 /m ∗ + 1 /m ∗ ] − is shown.) κ is only arough measure of the evanescent mode decay for a givenband. We discuss in detail the distribution of decay ex-ponents and phase offsets in Sec. IV.Let us now focus on the mechanism (ii) for the super-current decay, i.e. loss of synchronization between differ-ent transverse modes. This mechanism is well-known indiffusive systems where quasiparticle trajectory is ran-dom and thus the phase offset δϕ accumulated alongsuch a trajectory also gets randomized. Thus, upon av-eraging Eq. (3) over different disorder realizations, oneends up with exponentially decaying critical current [7].Previously, such a suppression of the supercurrent withjunction thickness, w , was often associated with impurityscattering in the ferromagnet. Here we demonstrate thatthis dephasing mechanism can also appear in clean sys-tems (where quasiparticle trajectory is well defined) dueto the dispersion of the velocity v z with an in-plane mo-mentum k (cid:107) . Specifically, we find that, in the Nb/Ni/Nbjunctions, this mechanism becomes relevant for junctionsthicker than Å, see Fig. 3 and Sec. IV. In SFS junctionsthe combination of disorder in the ferromagnet, interfacescattering as well band-structure-induced dephasing ul-timately determines the magnitude of the supercurrent.However, we believe that band structure effects providean upper bound on the magnitude of the critical currentas a function of w .In Sec. IV, we present numerical results which supportthe qualitative discussion presented above. III. METHOD
We develop a numerical method to perform realis-tic simulations of MJJs using a combination of first-principles DFT and BdG calculations. The former isused to obtain the normal-state properties (e.g., bandstructure, Fermi velocities, magnetization) and to calcu-late the normal scattering matrices through the inhomo-geneous 3D realistic junctions. As a next step we takesuperconductivity into account and calculate supercur-rent through the stack assuming the short junction limit.
A. Normal transport: ab initio description
To calculate the normal scattering matrix we usethe
Questaal package for electronic structure calcu-lations based on the linear muffin-tin orbital (LMTO)method [43]. It calculates the full non-linear, i.e. nonequilibrium, transport properties of an infinite system de-scribing a central (C) region cladded by two semi-infiniteleft (L) and right (R) leads [44, 45], as represented below: L (cid:122) (cid:125)(cid:124) (cid:123) . . . | L | L | C (cid:122) (cid:125)(cid:124) (cid:123) PL | PL | . . . | PL L − | R (cid:122) (cid:125)(cid:124) (cid:123) R | R | . . .
The LCR system [Fig. 1(c)] can be partitioned into aninfinite stack of principal layers (PLs) which interactonly with their nearest-neighbors. This is possible be-cause the screened LMTO structure constants are short-ranged [46]. In the present case the C region consists ofthe ferromagnet, plus two layers of Nb at the LC andCR interfaces respectively. This is the range over whichthe perturbation from C significantly modifies the po-tential in the L or R region. To construct the Nb/Ni/Nbstack, coincident site lattices for Ni and Nb must be found(details of how this was accomplished are given in Ap-pendix B). Planes of coincident site lattices are stackedto form the Nb/Ni/Nb structures. Figure 5(a) shows theNb and Ni planes we used, which are denoted here as‘surface supercells.’The electronic current flows along the z direction, per-pendicular to the PLs lying in the xy -plane (transversedirection), see Fig. 1(c). Periodic boundary conditionsare used within each PL. The corresponding reciprocalspace defines the two-dimensional (2D) k (cid:107) vectors, i.e.the transverse modes, used in the calculations. The k (cid:107) mesh is discretized and integrals over k (cid:107) are performednumerically.The electronic structure of the C region can be sepa-rated from L and R regions through self-energies, Σ L and Package website is questaal.org. Σ R , that modify the Hamiltonian of C region. They aremost easily calculated if the potential of each PL in theL or R region are identical all through the bulk region.This is the reason for adding a few Nb layers folded intothe C region. Thus the periodically repeating unit cellsin the L and R regions can be safely assumed to havethe potential of the bulk crystal. To construct the self-energies, the potentials of the PL in an infinite stack areneeded. These potentials are functions only of the PL intheir own region, and may be calculated in several ways. Σ L and Σ R are obtained from ‘surface’ Green’s function(a fictitious system which consists of a semi-infinite stackof PL, each with the same potential). Note that the po-tential of the C region is calculated self-consistently. Thisis important, as the local moments of Ni are small at theboundary layers, and build up gradually, see Fig. 5(b).With the potential in hand, the normal-state transportcan be calculated using scattering formalism [45, 47]. Forthis, knowing the retarded Green’s function, G r , of thejunction is sufficient. However, this is not the case forthe Josephson current: the individual eigenfunctions arerequired. Within a Green’s function framework, G r mustbe organized by normal modes which correspond to theeigenstates of the L and R leads for a prescribed energy E . In the PL representation, the Hamiltonian has beendiscretized into the linear combinations of the LMTO ba-sis functions, and the normal modes are represented aseigenvectors of these basis functions. The Schrödingerequation becomes a difference equation in the PL basisfunctions [48]. Eigenvectors are calculated by solving aquadratic eigenvalue problem [48, 49], whose eigenvaluescorrespond to exp( ± ik z,n a ) , where a is the thickness ofthe PL. The wave number k z,n of the normal mode n canbe complex, but to correspond to a propagating mode k z,n must be real. By solving the equation as a functionof the energy E , one gets all the eigenvalues and eigenvec-tors which provide the information needed to constructthe self-energies Σ L and Σ R (for each k (cid:107) and each spin σ ).Note that, in the mode basis, the imaginary part of theself-energies is proportional to the (band) velocity of themodes, and is diagonal for non-degenerate modes [49, 50].The retarded Green’s function G r of the C region (con-nected to the L and R leads) can be written as a matrixin the normal mode basis, G r σ ( E, k (cid:107) ) = (cid:110) [ G r C ,σ ( E, k (cid:107) )] − − Σ L ,σ ( E, k (cid:107) ) − Σ R ,σ ( E, k (cid:107) ) (cid:111) − , (5)where G r C is the Green’s function of the isolated C region.In this basis, G r is decomposed into four blocks, G r σ ( E, k (cid:107) ) = (cid:20) G r LL ,σ ( E, k (cid:107) ) G r LR ,σ ( E, k (cid:107) ) G r RL ,σ ( E, k (cid:107) ) G r RR ,σ ( E, k (cid:107) ) (cid:21) , (6)upon projecting onto the propagating modes of the L andR leads. These four quantities and the mode velocitiescompletely determine the normal state transport proper-ties of the junctions.The transmission matrices are defined by the off-diagonal parts of Eq. (6). More specifically, the transmis-sion coefficients [ t LR ,σ ] nm , connecting L and R regions,are given by [47] [ t LR ,σ ] nm ( E, k (cid:107) ) = i (cid:113) | [ v L ,σ ] n ( E, k (cid:107) ) |× [ G r LR ,σ ] nm ( E, k (cid:107) ) (cid:113) | [ v R ,σ ] m ( E, k (cid:107) ) | , (7)where [ v L ,σ ] n and [ v R ,σ ] m are the velocity matrix elementsfor propagating modes n and m in the L and R leads, re-spectively. The transmission matrix t RL can be obtainedfrom Eq. (7) by replacing L ↔ R . The reflection coeffi-cients are given by the diagonal blocks of Eq. (6). Forinstance, on the L side [ r LL ,σ ] nn (cid:48) ( E, k (cid:107) ) = i (cid:113) | [ v L ,σ ] n ( E, k (cid:107) ) |× [ G r LL ,σ ] nn (cid:48) ( E, k (cid:107) ) (cid:113) | [ v L ,σ ] n (cid:48) ( E, k (cid:107) ) | − δ nn (cid:48) . (8)The reflection matrix r RR on the R side is obtained fromEq. (8) by replacing L ↔ R .We define the normal scattering matrix as S = r LL , ↑ t LR , ↑ r LL , ↓ t LR , ↓ t RL , ↑ r RR , ↑ t RL , ↓ r RR , ↓ , (9)where we omit E - and k (cid:107) -dependence for brevity.The linear-response normal conductance, G σ , per spinis given by G σ A = e h (cid:90) BZ d k (cid:107) (2 π ) (cid:88) n,m (cid:12)(cid:12) [ t LR ,σ ] nm ( E F , k (cid:107) ) (cid:12)(cid:12) . (10)It is calculated at the Fermi energy, E = E F . The totalconductance is given by G = G ↑ + G ↓ . B. Superconducting transport: scattering matrixapproach
Equation (9) is the normal-state scattering matrix formetal-ferromagnet-metal (NFN) structure taking into ac-count reflection at both NF interfaces. To account for su-perconductivity, we introduce a step-like superconduct-ing pairing potential ∆ = 3 . meV and use the An-dreev approximation to account for electron-hole scatter-ing processes [3, 10]. This approach combines the detailsof the atomic structure of Nb/Ni/Nb and superconduc-tivity within the mean-field approximation.The direct contact between S and F layers, as shownin Fig. 1(a), would lead to an interaction between them.The back-action of the ferromagnet on the superconduc-tor (i.e. the inverse proximity effect) results in a spatialdependence of the pairing potential near the SF inter-faces [51, 52]. In the case of a clean SFS junction model, the self-consistent BdG calculations have been discussedin Refs. [53–55]. In typical experimental systems the su-perconductor is disordered, so disorder effect on pairingpotential needs to be considered, see, e.g., Ref. [56]. Fur-thermore, in recent experiments S and F layers are sepa-rated by an intermediate spacer layer, which significantlyreduces this inverse proximity effect. Thus, to under-stand inverse proximity effect in realistic SFS devices onewould need to include both of the abovementioned ingre-dients in the model as well as to take into account theinhomogeneous magnetization in the ferromagnet (dis-cussed in the next section) which is outside the scope ofthis paper. For the sake of clarity, we focus here on arealistic band structure in the ferromagnet and its effecton the supercurrent in the SFS structures. Henceforth,we also neglect the orbital effects of the fringe magneticfield, created by the ferromagnet.Since ∆ (cid:28) E F , spin-resolved Andreev reflection at SN L and N R S interfaces is described by r A ( φ ) = e iφ/ e iφ/ e − iφ/ e − iφ/ , (11)where φ is the phase difference between left and rightsuperconducting leads and is the identity matrix. Inthe short junction limit, the main contribution to thesupercurrent comes from Andreev bound states localizedin the junction having energy ε ν . The energy spectrumof Andreev states can be obtained using the followingequation [57], α ( ε ) (cid:20) r ∗ A ( φ ) r A ( φ ) 0 (cid:21)(cid:20) S ( E F + ε, k (cid:107) ) 00 S ∗ ( E F − ε, k (cid:107) ) (cid:21) Ψ in = Ψ in , (12)where α ( ε ) = (cid:112) − ε / ∆ + iε/ ∆ . The vector Ψ in =[ ψ L → e ↑ , ψ L → e ↓ , ψ R ← e ↑ , ψ R ← e ↓ , ψ L → h ↑ , ψ L → h ↓ , ψ R ← h ↑ , ψ R ← h ↓ ] T corre-sponds to the electron- and hole-like (e/h) waves in N L and N R regions incident on the F region from the left ( → )and from the right ( ← ).Simulations show that S is weakly-dependent on E in the range [ E F − ∆ , E F + ∆] . Therefore, we expand S ( E, k (cid:107) ) in E − E F and keep only the leading term, i.e. S ( E, k (cid:107) ) ≈ S ( E F , k (cid:107) ) = S ( k (cid:107) ) . Using this approxima-tion, one can simplify the quantization condition (12)and reduce it to the matrix eigenvalue problem (see de-tails in Ref. 58). This approach allows one to reliablycalculate the Andreev bound states spectrum, ε ν ( φ, k (cid:107) ) .The zero-temperature supercurrent, J , through the junc-tion is given by Eq. (3). IV. RESULTS
It is illuminating to compare our numerical simula-tions for the supercurrent with the experimental mea-surements involving quasi-ballistic MJJs. As previouslydiscussed, we believe that the Nb/Ni/Nb junctions repre-sent a good model system for which experimental data isreadily available [16, 21, 22]. The best-performing stacksconsist of Nb(110)/Cu/Ni(111)/Cu/Nb(110). Cu spacerlayers seem to be essential to get strong supercurrent,likely because it prevents intermixing of the Ni and Nb.Our model junction simulates this geometry, via super-cells in the plane normal to the stack to account for thelattice mismatch (see Fig. 1 and Appendix B), thoughwe do not include the Cu layers. We anticipate thatCu spacers will mainly affect transmission matrix ele-ments rather than the dependence of the supercurrenton ferromagnet thickness, which is the main focus of thiswork. Furthermore, as discussed before, the Cu spacerswill suppress the direct interaction between the ferromag-net and the superconductor and reduce inverse proximityeffect justifying Andreev approximation for the boundaryconditions, see Eq. (11). Therefore, we consider only thesimplified Nb/Ni/Nb stack and vary the number of lay-ers (atomic planes) of Ni. Additionally, we also considereffect of different crystallographic orientation of the Niplanes and investigate Nb/Ni(110)/Nb junctions in Ap-pendix C.To make a Nb/Ni superlattice, the unit cells of theNb and Ni regions in the plane normal to the interfacemust be coincident. This is complicated by the severelattice constant mismatch, and also the incompatibilityof the (110) and (111) atomic planes. It is necessary toconstruct superlattices with Nb(110) and Ni(111) bothrotated to the z axis, and with lattice vectors in the planecoincident. A supercell with nearly coincident vectorswas found (see Appendix B for details). By applyinga small shear strain to the Ni, the lattice vectors aremade exactly coincident. Figure 5(a) shows the Nb(110)surface supercell and the Ni(111) surface supercell withequal lattice vectors used to match the Nb/Ni interfaces.Each atomic plane of Ni(111) contains 14 atoms and eachatomic plane of Nb(110) consists of 10 atoms. The atomicstructure of the Nb(110)/Ni(111)/Nb(110) for 5 layers(atomic planes) of Ni is shown in Fig. 1(c).Next, we performed self-consistent DFT calculationswithin the local density approximation (LDA) in order toobtain the relaxed structure and corresponding electronicstructure. For the smallest structures we performed aconstrained optimization. Only the atoms in the planesclosest to the Nb/Ni interfaces are allowed to relax tofacilitate stacking of arbitrarily large cells. The Nb/Niinterplanar spacing has also been optimized to minimizethe total energy. Once the structure is determined, one can determinethe normal-state thermodynamic and transport proper-ties of the junction, e.g., calculate the magnetization pro-file and spin-resolved conductance through the junctionas a function of Ni thickness. For transport calculations, Description of the Nb/Ni/Nb supercells optimization at the
Questaal website.
Ni(111) layer M [ B ] a bNb(110) Figure 5: (a) Top view of the Nb(110) and Ni(111) sur-face supercells used to build the Nb/Ni interfaces and theNb(110)/Ni(111)/Nb(110) stacks shown in Fig. 1(c). The sur-face supercell is defined from two 2D vectors a = [10 . , Åand a = [6 . , . Å with periodic boundary condition inthe 2D plane. The corresponding reciprocal space defines the2D k (cid:107) vectors used in the calculations. Each atomic planeof Ni (Nb) contains 14 (10) atoms of Ni (Nb). (b) Mag-netic moment profile of Nb(110)/Ni(111)/Nb(110) junctionsfor different thickness of Ni (from 3 to 9 layers). The valueof the moment is an averaged over the moments of the 14 Niatoms in each atomic plane. Note the magnetic dead layerat the Nb/Ni interfaces and that all moments vanish for theshortest junction made of 3 layers of Ni. we use a layer transport technique [44] which employsthe atomic spheres approximation (ASA). Careful checkswere made of ASA band structures of elemental Nb andNi, and also superlattices, to confirm that they are verysimilar to the full potential LMTO DFT-LDA ones.We find that the magnetic properties of Ni are sensi-tive to their local environment, indicative of the itinerantferromagnetism. As shown in Fig. 5(b), the magnetiza-tion profile is non-uniform in the junction with averagedmagnetic moments per atom being suppressed near theNb interface. For thickness larger then 4 layers, one re-covers the bulk value of ∼ . µ B in the middle layers,away from the Nb/Ni interfaces. The averaged momentdrops down towards the edges and becomes considerablyreduced down to ∼ . µ B at the interface with Nb. Thestrong reduction of magnetism is exemplified for the shortjunction with 3 layers of Ni where the moments on theNi atoms have completely vanished. Such a non-uniformmagnetic moment dependence in Nb/Ni/Nb junctions af-fect superconducting properties of the SFS junctions in anon-trivial way. For example, Nb/Ni/Nb junctions thin-ner than 4 layers of Ni behave as essentially SNS junc-tions.It is well known that the LDA tends to overestimatelocal moments M in itinerant magnets [59] because spinfluctuations reduce the average moment [60], and under-estimate M when local moments are very large [42]. ForNi, LDA yields M in good agreement with the experi-ment, but this is likely an artifact of an accidental can-cellation of errors. Most important for transport is theexchange splitting V ex , which the LDA predicts to be0.6 eV, about twice larger than the experimental value of0.3 eV [61]. It is possible to reproduce both M and V ex atthe same time, but a high-level theory is needed to sur- w [Å]10 e J c A / G e J c A / Ge | J | A / Ge | J fit1 | A / G Figure 6: Comparison of critical current density, J c =max φ | J ( φ ) | (blue squares), absolute value of first Fouriercomponent for supercurrent, | J | [see Eq. (13), red circles],and its fitting | J fit1 | [Eq. (14), black curve]. All these quanti-ties are ‘normalized’ by normal-state conductance, G . mount both kinds of errors inherent the LDA [42]. Thehigh cost and poor scaling of such a theory is not practi-cal for these junctions, so we elect to stay within the LDAand scale the self-consistently calculated V ex . This wasthe approach Karlsson and Aryasetiawan used to calcu-late the spin wave spectra in Ni [62]. Scaling of V ex can beaccomplished using different approaches, e.g., by addingsome effective magnetic field to simulate the effect of spinfluctuations. Since Ni is a simple case with a nearly linearrelation between M and V ex , the band structure hardlydepends on the details in which the LDA potential ismodified. Here we first perform fully self-consistent cal-culations. Then, to construct the potential for transportproperties, we rescale the spin component of the densityby a constant factor, which we denote as M/M . Thisenables parametric studies of transport as a function of V ex . M/M = 0 . yields the observed V ex = 0 . eV, andwe use this scaling unless stated otherwise.The conductance per unit of area, G/A , is shown inFig. 3(b). It is weakly dependent on the thickness, w , ofthe magnetic layer, as expected for a metallic system inthe absence of disorder.We now focus on superconducting properties. The de-pendence of the supercurrent on the phase difference φ for 5, 8, and 11 Ni layers is shown in Figs. 3(c)–3(e). Onecan present current-phase relation, J ( φ ) in a form of aFourier series, J ( φ ) = (cid:88) n (cid:62) J n sin( nφ ) . (13)In the -junction mode, the first term in the Fourier seriesdominates with J > [see solid black line in Fig. 3(c)].In π -junction case [Fig. 3(e)], the supercurrent is alsomostly defined by the first harmonic but with J < .Close to the – π transition J dies out, so that the be-havior is governed by higher Fourier harmonics, e.g. for8 Ni layers supercurrent has mostly second harmonic, J ,shown by dashed black line in Fig. 3(d). Figure 6 shows the critical current density, J c =max φ | J ( φ ) | , normalized by normal-state conductance, eJ c A/G ∆ , as a function of w (blue squares). Far fromthe – π transitions, the critical current coincides with theabsolute value of the first Fourier harmonic, e | J | A/G ∆ (red circles). Since J contains a sign of the current andhas better numerical stability than J c , we use this quan-tity for the analysis. We exclude very thin junctions (3Ni layers or less) from the analysis since the magneticproperties are suppressed there. The J dependence on w can be fit by the followingexpression, J fit1 ( w ) = Θ J exp( − w/ξ J ) cos (cid:2) π ( w + δ J ) /λ J (cid:3) , (14)where Θ J = 3 .
20 A /µ m , ξ J = 41 . Å, λ J = 23 . Å, and δ J = − . Å are fitting parameters. We interpret ξ J asa decay length, λ J as the ‘half-period’ of the oscillationin J as a function of w , and δ J as a measure of thesuppressed magnetization in Ni layers near the Ni/Nbboundaries. J fit1 ( w ) accurately fits the discrete points J as shown in Figs. 3(a), 6, 7, and 10 by black solid lineand black circles, accordingly.In order to gain insight into the evolution of J with w ,let us resolve contributions from different k (cid:107) . For this,we rewrite Eq. (3) as J ( φ ) = A (cid:90) BZ d k (cid:107) (2 π ) j ( φ, k (cid:107) ) , (15a) j ( φ, k (cid:107) ) = − e (cid:126) A (cid:88) ν> ∂ε ν ( φ, k (cid:107) ) ∂φ . (15b)Here A = 76 . Å is the area of the surface supercellshown in Fig. 5(a). Similar to Eq. (13), we denote thefirst φ -harmonic of j ( φ, k (cid:107) ) as j ( k (cid:107) ) . The evolution ofthe first Fourier harmonic of the supercurrent J and j ( k (cid:107) ) as a function of w are shown in Figure 7. Calcu-lations were performed for two different sets of k (cid:107) with290 discrete k (cid:107) -points (full black circles) and 4142 k (cid:107) -points (empty black circles). One can see that both setsgive the same result for J , establishing that the k (cid:107) in-tegration is well converged. In Fig. 7, ‘errorbars’ denotethe standard deviation in j ( k (cid:107) ) with respect to the k (cid:107) -summed average, J ( φ ) . Individual j ( k (cid:107) ) are shown bysemi-transparent horizontal dashes. The important ob-servation is that while J decays with w , the dispersionin j ( k (cid:107) ) does not change significantly. For the first data point in Fig. 6 corresponding to 3 layers ofNi the magnetization is completely suppressed [see Fig. 5(b)]and ratio eJ c A/G ∆ ≈ . is significantly higher than the one forthicker Ni regions with non zero magnetic moments. We comparethis value with the result for the short disordered SNS junction.Combination of analytical energy spectrum [57] with Dorokhovdistribution of channel transmissions [63, 64] leads to a ratio of2.1 (horizontal dashed line in Fig. 6). We attribute this differenceto the fact that there is no interfacial disorder in our model. w [Å]32101234 J [ A / m ] J (290 k -points) J (4142 k -points) J fit1 Figure 7: First Fourier harmonic, J , of the supercurrent as a function of junction thickness, w . Full black circles correspondto J calculated with 290 k (cid:107) -points, empty circles (mostly superposed onto the full black circles) correspond to 4142 k (cid:107) -points.Solid black line is J fit1 ( w ) [ J fit given by Eq. (14)]. Green semitransparent dashes show J fit1 ( w ) contributions [Eq. (15b)] for theindividual k (cid:107) -points. The vertical ‘errorbars’ correspond to the standard deviation of j ( k (cid:107) ) with respect to J . The standarddeviation is the same for both sets of k (cid:107) -point indicating that these results are independent of the chosen discretization. w = 6.1 Å4 L -junction24.2 Å13 L 0-junction48.3 Å25 L -junction70.4 Å36 L 0-junction94.5 Å48 L -junction118.6 Å60 L k x k y j ( k || ) / m a x k || | j ( k || ) | Figure 8: Colorplot of j ( k (cid:107) ) with k (cid:107) = ( k x , k y ) . Each panel corresponds to the local extrema of J fit1 ( w ) shown in Fig. 7.For 4 and 13 layers all the k (cid:107) channels contribute with the same sign. For 48 layers and larger, different k (cid:107) channels losesynchronization and contribute to the total supercurrent with different signs. Figure 8 shows colorplots of j ( k (cid:107) ) corresponding tolocal extrema of J fit1 ( w ) (4, 13, 25, 36, 48, and 60 layerslabeled in Fig. 7). For small w , one can see that mostof all the k (cid:107) contributions to J have the same sign, i.e.positive in -junction regime and negative in π -junctionregime. In this regime the decay is predominantly due toevanescent modes decaying into the junction. For larger w , the dephasing mechanism becomes important sincethe phase offset spread grows with w . One can observethe apparition of contributions of the opposite sign for w (cid:38) Å. This dephasing mechanism is mainly due tothe variation of the Fermi velocity with k (cid:107) , and becomesmore important with increasing w .In order to study the distribution of the phase offsetsand decay exponents for different modes, we fit the indi-vidual j ( k (cid:107) ) using an expression analogous to Eq. (14).The set of the resulting fitted curves j fit1 ( w, k (cid:107) ) for 4142 k (cid:107) -points are shown by the green semitransparent curvesin Fig. 3(a). Here to minimize the numerical ‘noise,’ j ( w, k (cid:107) ) curves are smoothed over the 2D k (cid:107) space usingGaussian filter with σ k (cid:107) = 0 . Å − which is of the orderof the Fermi wave vector in Nb. Thus, each data point j ( w, k (cid:107) ) approximately corresponds to a transverse con-ducting channel. This fitting procedure works reason-ably well, e.g., the relationship in Eq. (15a) holds if onereplaces J ( φ ) by J fit1 ( w ) and j ( φ, k (cid:107) ) by j fit1 ( w, k (cid:107) ) for w (cid:38) Å.The distribution of the fitting parameters for j fit1 ( w, k (cid:107) ) is shown in Fig. 9. Histograms for decaylengths and half-periods reveal a complicated picture de-scribing different contributions to the supercurrent in realmaterials. First of all, in Fig. 9(a) one can see the dis-tribution of the half-periods, λ j , which is similar to aGaussian distribution with a mean value (cid:104) λ j (cid:105) = 23 . Å1 F r e q u e n c y [ a . u . ] b j [Å]15 20 25 30 35 a j [Å] Figure 9: Analysis of the fitting parameters of j fit1 ( w ) for indi-vidual k (cid:107) , shown in Fig. 3(a) by semitransparent green lines.Distribution of (a) half-periods, λ j , and (b) decay lengths, ξ j . and standard deviation . Å. The mean value (cid:104) λ j (cid:105) isvery close to λ J [see text after Eq. (14)] while the spreadin λ j leads to dephasing and is responsible for the expo-nential decay of J at large w . Indeed, it is well-knownthat the average of an oscillatory function with respectto a random fluctuating phase (described by a Gaussiandistribution) results in an exponentially decaying func-tion.In addition to the dephasing mechanism, the decay ofthe supercurrent originates from the evanescent modes.The histogram for decay lengths, ξ j is shown in Fig. 9(b).Here small ξ j corresponds to fast-decaying j fit1 ( w, k (cid:107) ) ,large ξ j is responsible for non-decaying modes (i.e. modeswith the decay exponents larger than the junction thick-ness). The right-skewed distribution of the decay expo-nents has a mean value of (cid:104) ξ j (cid:105) = 108 Å which is muchlarger than ξ J in Eq. (14). The shoulder at small ξ j presumably corresponds to the evanescent mode decaycomprising of d bands 2 and 3, see Table I, whereas thetail at large ξ j originates predominantly from the band6. Overall, one can see that a fit with a single decay ex-ponent, discussed in Eq. (14), is quite oversimplified fora Nb/Ni/Nb junction considered here.We now turn to the discussion of the effect of exchangesplitting energy on the supercurrent in MJJs. So farwe have used M/M = 0 . , which yields the experi-mentally observed V ex = 0 . eV. It is interesting to in-vestigate how a ferromagnet with a different V ex (butotherwise the same band structure as Ni) would affectthe w -dependence of J . In Fig. 10 we show the re-sults for parametric variations in M/M . One can seein Fig. 10(a) that the half-period, λ J , and the decaylength, ξ J , strongly depend on M . Here black pointscorrespond to M/M = 0 . , and the self-consistent cal-culations with no rescaling correspond to M/M = 1 . Inorder to understand how λ J and ξ J depend on M , we per-form the fitting procedure Eq. (14) for different magneticmoments and plot, in Fig. 10(b), the supercurrent den-sity as a function of the rescaled thickness, ( w + δ J ) / λ J with ˜ λ J = ( M /M ) 11 . Å. Remarkably J , as a functionof the rescaled thickness, collapses to the same universalcurve. The inset demonstrates that the fitted half-period, λ J , is proportional to /M showing that the oscillation w [Å]2101234 J [ A / m ] M / M = 0.125 M / M = 0.25 M / M = 0.33 M / M = 0.5 M / M = 1 w + J )/2 J J [ A / m ] J [Å]050100 J [ Å ] ab Figure 10: J for different rescaling of the magnetic moments M/M as a function of (a) the thickness, w , and (b) therescaled thickness, ( w + δ J ) / λ J , where δ J is the fitting pa-rameter in Eq. (14) and ˜ λ J = ( M /M ) 11 . Å. J values areshown by circles; corresponding fittings J fit1 ( w ) are shown bylines of the same color. (Inset) Half-period, λ J , fitted usingEq. (14) versus ˜ λ J . period scales linearly with the inverse of V ex in this pa-rameter range. Finally, we have also performed calculations forNb/Ni/Nb junctions built from stacking the Ni atomicplanes in the (110) orientation instead of the (111) ori-entation. The results are given in Appendix C. Qualita-tively, the same physics hold for both stacks built from(111) and (110) Ni planes. However, our calculationsshow that the actual value for the period of oscillationand for the current decay depend crucially on the detailsof the electronic structure of the junctions, such as therelative crystal orientation.
V. CONCLUSION
In this paper we have studied the magnetic and super-conducting properties of SFS junctions using the com-bination of density functional theory and Bogoliubov-deGennes model. We believe that this method allows one Deviations from linear regime become significant for
M/M (cid:38) . . For the clarity of the data, we do not show the results for M/M > in Fig. 10(b). w . We identified two generic mechanisms for the decay ofthe supercurrent with w : (i) exchange-splitting inducedgap opening for minority or majority carriers and (ii) de-phasing between different modes due to the significantquasiparticle velocity dispersion with the transverse mo-mentum. It was previously believed that disorder in theferromagnet is mainly responsible for the supercurrentdecay in SFS junctions. In the present work we haveshown that band structure effects also contribute to thecritical current suppression and thus provide an upperbound for the supercurrent in ideal (i.e. disorder-free)structure.We identified that the Nb/Ni/Nb junction is a suit-able system for comparison with the simulations becauseof the long mean free path in Ni relative to the junctionthickness and the quasi-ballistic nature of quasiparticlepropagation in the ferromagnet. We have found goodagreement with published experimental data for the half-period of the critical current oscillations: λ J ≈ Å [seeEq. (14) and text after it] versus ≈ Å in experiment,Ref. 21. We have also found that the critical currentdecays exponentially with the ferromagnet thickness w .This is to be contrasted with previously assumed alge-braic decay based on results for the clean SFS junctionsusing simple parabolic-like band structure. We believethat in measured Nb/Ni/Nb junctions with w (cid:46) Å themechanism (i) is likely to be responsible for the supercur-rent decay. This finding is crucial for material and geom-etry optimization of MJJs and superconducting magneticspin valves.Understanding the interplay of band structure effectsand disorder in MJJs is an interesting open problem. Webelieve that interfacial disorder due to, for example, sur-face roughness will mix different k (cid:107) modes and will leadto a larger spread of half-periods. This, in turn, will fur-ther enhance the dephasing mechanism (ii) of the super-current decay discussed here. Strong disorder in the bulk(i.e. mean free path much smaller than junction thickness w ) would lead to the diffusive motion of quasiparticles inthe ferromagnet which is a significant departure from thequasi-ballistic junction limit considered here. We thinkthat bulk disorder would induce even more dephasing be-tween different modes because phase offsets in this casewill depend on different random trajectories of minorityand majority carriers. We, therefore, believe that bulkdisorder will lead to even stronger decay of the supercur-rent with junction thickness, w . ACKNOWLEDGMENTS
HN and MvS acknowledge financial support from Mi-crosoft Station Q via a sponsor agreement between KCLand Microsoft Quantum. HN is much indebted to Dimi-tar Pashov for stimulating discussions about developing the
Questaal package. The authors express their grati-tude to Mason Thomas for the organizational help anddiscussions at the early stages of the project. The au-thors acknowledge stimulating discussions with NormanBirge, Anna Herr, Tom Ambrose, Nick Rizzo, and DonMiller.
Appendix A: Mean free path estimate
In this section we provide an estimate for the mean freepath l MFP in the Ni ferromagnet. Our approach is similarto Ref. 65. We use the Kubo formula for conductivity σ xx = e (cid:88) n,σ τ nσ (cid:104) v nσ (cid:105) ρ n ( E F ) , where index n labels the Ni bands, σ is the spin projec-tion, and τ nσ is the scattering time in each band assumedto be momentum independent. Mean square velocity, (cid:104) v nσ (cid:105) , and density of states at the Fermi level, ρ n ( E F ) ,are obtained from our ab initio model calculations andare provided in Table I.We use the available experimental data [66, 67] forthick Ni samples w (cid:38) Å: the low temperature linearresistance / ( σ ↑ xx + σ ↓ xx ) ≈
33 nΩ · m , and bulk spin scat-tering asymmetry, β F = ( σ ↑ xx − σ ↓ xx ) / ( σ ↑ xx + σ ↓ xx ) = 0 . .Only a single ↑ band of the majority spin is occupied.Thus, one can readily estimate mean free path for major-ity carriers l ↑ MFP = (cid:104) v ↑ (cid:105) τ ↑ ≈ Å. For the minority chan-nel we assume that the conductance at large thicknessesis dominated by the most mobile band, i.e. ↓ [68]. Thisgives an estimate for the mean free path l ↓ MFP ≈ Å,similar to its exchange-split partner ↑ . These estimatesare consistent with the available ARPES [69] and com-putational data [65]. As mentioned in the main text, themean free path for ↑ and ↓ carriers exceeds the typicaljunction thicknesses measured experimentally. Appendix B: Nb(110)/Ni(111)/Nb(110) stacks
Here we describe how coincident site lattices are con-structed for the Nb(110)/Ni(111) interface. Supercells ofboth Ni and Nb are separately constructed in the follow-ing way, to make a coincident site lattice.The primitive unit cells ( fcc in the Ni case, with a 0 Klattice constant 3.515 Å, and bcc in the Nb case with a0 K lattice constant 3.295 Å) are rotated. For Ni, therotation is compactly described in terms of three Eu-ler angles: rotation about z by π/ , rotation about x (cid:48) by arccos( (cid:112) / , rotation about z (cid:48)(cid:48) by arccos(4 / √ .Axes xyz are shown in Fig. 1(c); x (cid:48) y (cid:48) z (cid:48) and x (cid:48)(cid:48) y (cid:48)(cid:48) z (cid:48)(cid:48) areaxes after first and second rotation, respectively. Fromthe first two rotations the [1¯11] axis becomes the new z (cid:48)(cid:48) axis. The last rotation is needed to make it approx-imately coincident with Nb. The Nb is rotated about (1 , − , by − π/ , about z (cid:48) by π/ , and about z (cid:48) again3 ca b Ni(111)Nb(110)Nb(110) w [Å]210123 J [ A / m ] Nb/Ni(111)/NbNb/Ni(110)/Nb
Figure 11: (a) Top view of the Nb(110)/Ni(110) in-terface supercell. Ni (Nb) atoms are shown in lightblue (grey). A repetition ( × ) of the surface super-cell is shown. The corresponding Nb(110) and Ni(110)cells contain 2 atoms each. (b) Transverse view of theNb(110)/Ni(110)/Nb(110) junction with 5 layers of Ni.(c) J component for the Nb(110)/Ni(110)/Nb(110) and theNb(110)/Ni(111)/Nb(110) junctions versus the thickness w .Both currents have a decaying oscillatory behavior, with ahalf period of oscillation of λ J ≈ Å for the Ni(110) caseand of λ J ≈ Å for the Ni(111) case. by arccos(5 / √ . After rotations both Ni (1¯11) and Nb (¯1¯10) planes are normal to z , which is the propagationdirection.Next, superlattices must be constructed. They are gen-erated by scaling the primitive lattice vectors by the fol-lowing integer multiples, for Ni and Nb respectively: − , − − − . The Ni supercell contains 14 atoms, all in a single planeperpendicular to z ; the Nb supercell consists of 2 planeswith 10 atoms per plane [Fig. 5(a)]. The lattice vectorstransverse to z are nearly coincident, but not identicallyso. To render them coincident, we opt to shear the Ni bythe following linear transformation, (cid:20) .
997 00 .
050 1 . (cid:21) . This shear is a measure of the remaining mismatch ofthe undistorted Ni and Nb lattices. It is close enough tounity to have a minor effect on the Ni band structure.Along this axis, Nb planes form an ABAB . . . stack-ing pattern; the Ni are stacked ABCABC . . .
Ni/Nb in-terfaces are formed by stacking varying numbers of Ni planes on a Nb substrate, and relaxing the supercells tominimize the total energy. This was done for a few smallsuperlattices with 3, 4, and 5 Ni planes, restricting therelaxation to the two Ni and two Nb frontier planes. Inthis way we can build up structures of arbitrarily manyNi planes, sandwiching unrelaxed planes between frontierplanes. Thus the entire set of structures has three fami-lies of interfaces: those with integer numbers n , n + 1 or n + 2 of Ni planes.Figure 3(b) shows that normal conductance, G , expe-riences ∼ . deviations from its mean (cid:104) G (cid:105) followingthe pattern of these three families. Overall, both normalconductance and supercurrent (see Fig. 7) are reasonablysmooth functions of the thickness w , which suggests thatthe discreteness of the lattice and the details of latticerelaxation plays a minor role. Appendix C: Nb(110)/Ni(110)/Nb(110) junctions
In this section, we present results for the Nb/Ni/Nbjunctions built from stacking the Ni atomic planes withthe (110) orientation instead of the (111) orientation. Al-though the (110) plane orientation is not the most stableconfiguration for the Nb/Ni interface and most experi-mental studies focus on (111) orientation, it is interestingto ascertain whether there is an effect of crystallographicorientation on the supercurrent. Indeed, the supercur-rent decay and half-period of oscillations with w do de-pend on crystal orientation.We study the Nb(110)/Ni(110)/Nb(110) trilayershown in Fig. 11(a). The inter-plane distance betweenthe Nb and Ni planes at the interfaces is taken to be theaverage of the Nb and Ni inter-plane distances. For thesake of simplicity, we do not perform atomic relaxationsfor this system. The lattice parameter of Nb (3.295 Å)has been increased to match the lattice parameter of Ni(3.515 Å) to simplify the construction of periodic super-cell. This has a slight effect on the Nb band structure.Calculations of the trilayers were performed self-consistently, and for transport we reconstruct the po-tential by rescaling the magnetic moments by 0.5, as inthe Ni(111) case. The – π and π – transitions occuraround w ≈ Å and 41 Å [shown in red in Fig. 11(c)],yielding a half-period for the oscillations of λ J ≈ Å.This period is significantly larger than λ J ≈ Å forthe Nb(110)/Ni(111)/Nb(110) case [shown in black inFig. 11(c)]. The difference of ≈ Å between the twohalf periods represents roughly 5 inter-plane distancesfor Ni(110) [3 inter-plane distances for Ni(111)], and isnot induced by the unrelaxed atomic structures at theNb(110)/Ni(110) interfaces. It is noteworthy that thecalculated period for the Ni(111) trilayers better coin-cides with available experimental data [21] than for theNi(110) case.4 [1] B. Josephson, Possible new effects in superconductivetunnelling, Phys. Lett. , 251 (1962).[2] P. A. Warburton, The josephson effect: 50 years of sci-ence and technology, Phys. Educ. , 669 (2011).[3] C. W. J. Beenakker, Three “universal” mesoscopicJosephson effects, in Transport phenomena in mesoscopicsystems , Vol. 109, edited by H. Fukuyama and T. Ando(Springer, Berlin, Heidelberg, 1992) pp. 235–253.[4] A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, Thecurrent-phase relation in Josephson junctions, Rev. Mod.Phys. , 411 (2004).[5] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Micro-scopic theory of superconductivity, Phys. Rev. , 162(1957).[6] I. M. Dayton, T. Sage, E. C. Gingrich, M. G. Loving,T. F. Ambrose, N. P. Siwak, S. Keebaugh, C. Kirby, D. L.Miller, A. Y. Herr, Q. P. Herr, and O. Naaman, Experi-mental demonstration of a Josephson magnetic memorycell with a programmable π -junction, IEEE Magn. Lett. , 1 (2018).[7] A. I. Buzdin, Proximity effects in superconductor-ferromagnet heterostructures, Rev. Mod. Phys. , 935(2005).[8] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Oddtriplet superconductivity and related phenomena insuperconductor-ferromagnet structures, Rev. Mod. Phys. , 1321 (2005).[9] M. G. Blamire and J. W. A. Robinson, The interfacebetween superconductivity and magnetism: understand-ing and device prospects, J. Condens. Matter Phys. ,453201 (2014).[10] M. Eschrig, Spin-polarized supercurrents for spintronics:a review of current progress, Rep. Prog. Phys , 104501(2015).[11] P. Fulde and R. A. Ferrell, Superconductivity in a strongspin-exchange field, Phys. Rev. , A550 (1964).[12] A. I. Larkin and Y. N. Ovchinnikov, Inhomogeneous stateof superconductors, Sov. Phys. JETP , 762 (1965).[13] A. L. Buzdin, L. N. Bulaevskil, and S. V. Panyukov,Critical-current oscillations as a function of the exchangefield and thickness of the ferromagnetic metal (F) in anS-F-S Josephson junction, JETP Lett. , 178 (1982).[14] V. V. Ryazanov, V. A. Oboznov, A. Y. Rusanov, A. V.Veretennikov, A. A. Golubov, and J. Aarts, Coupling oftwo superconductors through a ferromagnet: Evidencefor a π junction, Phys. Rev. Lett. , 2427 (2001).[15] C. Bell, G. Burnell, C. W. Leung, E. J. Tarte, D.-J.Kang, and M. G. Blamire, Controllable Josephson cur-rent through a pseudospin-valve structure, Appl. Phys.Lett , 1153 (2004).[16] E. C. Gingrich, B. M. Niedzielski, J. A. Glick, Y. Wang,D. L. Miller, R. Loloee, W. P. P. Jr, and N. O. Birge,Controllable 0- π Josephson junctions containing a ferro-magnetic spin valve, Nat. Phys. , 564 (2016).[17] T. Kontos, M. Aprili, J. Lesueur, F. Genêt, B. Stephani-dis, and R. Boursier, Josephson junction through a thinferromagnetic layer: Negative coupling, Phys. Rev. Lett. , 137007 (2002).[18] H. Sellier, C. Baraduc, F. m. c. Lefloch, and R. Calem-czuk, Temperature-induced crossover between and π states in S/F/S junctions, Phys. Rev. B , 054531 (2003).[19] J. W. A. Robinson, S. Piano, G. Burnell, C. Bell, andM. G. Blamire, Critical current oscillations in strongferromagnetic π junctions, Phys. Rev. Lett. , 177003(2006).[20] T. S. Khaire, W. P. Pratt, and N. O. Birge, Critical cur-rent behavior in Josephson junctions with the weak fer-romagnet PdNi, Phys. Rev. B , 094523 (2009).[21] B. Baek, M. L. Schneider, M. R. Pufall, and W. H. Rip-pard, Phase offsets in the critical-current oscillations ofJosephson junctions based on Ni and Ni-(Ni Fe ) x Nb y barriers, Phys. Rev. Applied , 064013 (2017).[22] B. Baek, M. L. Schneider, M. R. Pufall, and W. H. Rip-pard, Anomalous supercurrent modulation in Josephsonjunctions with Ni-based barriers, IEEE Trans. Appl. Su-percond , 1 (2018).[23] V. Aguilar, D. Korucu, J. A. Glick, R. Loloee, W. P.Pratt, and N. O. Birge, Spin-polarized triplet supercur-rent in Josephson junctions with perpendicular ferromag-netic layers, Phys. Rev. B , 024518 (2020).[24] W. L. McMillan, Theory of superconductor—normal-metal interfaces, Phys. Rev. , 559 (1968).[25] T. Wolfram, Tomasch oscillations in the density of statesof superconducting films, Phys. Rev. , 481 (1968).[26] I. O. Kulik, Macroscopic quantization and the proximityeffect in S-N-S junctions, J. Exp. Theor. Phys. , 944(1970).[27] J. Demers and A. Griffin, Scattering and tunneling ofelectronic excitations in the intermediate state of super-conductors, Can. J. Phys. , 285 (1971).[28] A. Griffin and J. Demers, Tunneling in the normal-metal-insulator-superconductor geometry using the Bogoliubovequations of motion, Phys. Rev. B , 2202 (1971).[29] O. Entin-Wohlman, Effect of a barrier at thesuperconducting-normal metal interface, J. Low Temp.Phys , 777 (1977).[30] G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Tran-sition from metallic to tunneling regimes in supercon-ducting microconstrictions: Excess current, charge im-balance, and supercurrent conversion, Phys. Rev. B ,4515 (1982).[31] A. Furusaki and M. Tsukada, Dc Josephson effect andAndreev reflection, Solid State Commun , 299 (1991).[32] A. Furusaki, H. Takayanagi, and M. Tsukada, Joseph-son effect of the superconducting quantum point contact,Phys. Rev. B , 10563 (1992).[33] A. Furusaki, DC Josephson effect in dirty SNS junctions:Numerical study, Physica B , 214 (1994).[34] M. J. M. de Jong and C. W. J. Beenakker, Andreev re-flection in ferromagnet-superconductor junctions, Phys.Rev. Lett. , 1657 (1995).[35] Y. Tanaka and S. Kashiwaya, Theory of Joseph-son effect in superconductor-ferromagnetic-insulator-superconductor junction, Physica C , 357 (1997).[36] I. Žutić and O. T. Valls, Spin-polarized tunneling in ferro-magnet/unconventional superconductor junctions, Phys.Rev. B , 6320 (1999).[37] Z. Radović, N. Lazarides, and N. Flytzanis, Josephson ef-fect in double-barrier superconductor-ferromagnet junc-tions, Phys. Rev. B , 014501 (2003).[38] J. Cayssol and G. Montambaux, Incomplete an- dreev reflection in a clean superconductor-ferromagnet-superconductor junction, Phys. Rev. B , 012507(2005).[39] F. Konschelle, J. Cayssol, and A. I. Buzdin, Nonsinu-soidal current-phase relation in strongly ferromagneticand moderately disordered sfs junctions, Phys. Rev. B , 134505 (2008).[40] A. F. Tzortzakakis and N. Flytzanis, Josephson junctionswith spin-orbit and spin-flip interactions , Master’s thesis,University of Crete (2019).[41] M. Cheng and R. M. Lutchyn, Josephson cur-rent through a superconductor/semiconductor-nanowire/superconductor junction: Effects of strongspin-orbit coupling and Zeeman splitting, Phys. Rev. B , 134522 (2012).[42] L. Sponza, P. Pisanti, A. Vishina, D. Pashov, C. Weber,M. van Schilfgaarde, S. Acharya, J. Vidal, and G. Kotliar,Self-energies in itinerant magnets: A focus on Fe and Ni,Phys. Rev. B , 041112 (2017).[43] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson,K. D. Belashchenko, A. Chantis, F. Jamet, and M. vanSchilfgaarde, Questaal: A package of electronic structuremethods based on the linear muffin-tin orbital technique,Comput. Phys. Commun , 107065 (2020).[44] S. V. Faleev, F. Léonard, D. A. Stewart, and M. vanSchilfgaarde, Ab initio tight-binding LMTO method fornonequilibrium electron transport in nanosystems, Phys.Rev. B , 195422 (2005).[45] Y. Meir and N. S. Wingreen, Landauer formula for thecurrent through an interacting electron region, Phys.Rev. Lett. , 2512 (1992).[46] O. K. Andersen and O. Jepsen, Explicit, first-principlestight-binding theory, Phys. Rev. Lett. , 2571 (1984).[47] D. S. Fisher and P. A. Lee, Relation between conductivityand transmission matrix, Phys. Rev. B , 6851 (1981).[48] A.-B. Chen, Y.-M. Lai-Hsu, and W. Chen, Difference-equation approach to the electronic structures of sur-faces, interfaces, and superlattices, Phys. Rev. B , 923(1989).[49] Y. Fujimoto and K. Hirose, First-principles treatmentsof electron transport properties for nanoscale junctions,Phys. Rev. B , 195315 (2003).[50] M. Wimmer, Quantum transport in nanostructures:From computational concepts to spintronics in grapheneand magnetic tunnel junctions , Ph.D. thesis, UniversitätRegensburg (2008).[51] K. Halterman and O. T. Valls, Proximity effects atferromagnet-superconductor interfaces, Phys. Rev. B ,014509 (2001).[52] K. Halterman and O. T. Valls, Proximity effectsand characteristic lengths in ferromagnet-superconductorstructures, Phys. Rev. B , 224516 (2002).[53] K. Halterman and O. T. Valls, Layered ferromagnet-superconductor structures: The π state and proximity effects, Phys. Rev. B , 014517 (2004).[54] K. Halterman, P. H. Barsic, and O. T. Valls, Oddtriplet pairing in clean superconductor/ferromagnet het-erostructures, Phys. Rev. Lett. , 127002 (2007).[55] K. Halterman, O. T. Valls, and C.-T. Wu, Charge andspin currents in ferromagnetic Josephson junctions, Phys.Rev. B , 174516 (2015).[56] V. O. Yagovtsev, N. G. Pugach, and M. Eschrig,The inverse proximity effect in strong ferromag-net–superconductor structures, Supercond. Sci. Technol. , 025003 (2021).[57] C. W. J. Beenakker, Universal limit of critical-currentfluctuations in mesoscopic josephson junctions, Phys.Rev. Lett. , 3836 (1991).[58] B. van Heck, S. Mi, and A. R. Akhmerov, Single fermionmanipulation via superconducting phase differences inmultiterminal Josephson junctions, Phys. Rev. B ,155450 (2014).[59] A. Aguayo, I. I. Mazin, and D. J. Singh, Why Ni Al is anitinerant ferromagnet but Ni Ga is not, Phys. Rev. Lett. , 147201 (2004).[60] T. Moriya, Spin fluctuations in itinerant electron mag-netism (Springer-Verlag, Berlin, 1985).[61] F. J. Himpsel, J. A. Knapp, and D. E. Eastman, Exper-imental energy-band dispersions and exchange splittingfor ni, Phys. Rev. B , 2919 (1979).[62] K. Karlsson and F. Aryasetiawan, A many-body ap-proach to spin-wave excitations in itinerant magnetic sys-tems, J. Phys. Condens. Matter , 7617 (2000).[63] O. N. Dorokhov, On the coexistence of localized and ex-tended electronic states in the metallic phase, Solid StateCommun , 381 (1984).[64] P. A. Mello, P. Pereyra, and N. Kumar, Macroscopic ap-proach to multichannel disordered conductors, Ann Phys , 290 (1988).[65] D. Gall, Electron mean free path in elemental metals, J.Appl. Phys , 085101 (2016).[66] C. E. Moreau, I. C. Moraru, N. O. Birge, and W. P.Pratt, Measurement of spin diffusion length in sputteredNi films using a special exchange-biased spin valve geom-etry, Appl. Phys. Lett , 012101 (2007).[67] J. Bass, Giant magnetoresistance: Experiment, in Hand-book of spin transport and magnetism , edited by I. Zuticand E. Y. Tsymbal (CRC Press, 2011) pp. 69–94.[68] E. Y. Tsymbal, D. Pettifor, and S. Maekawa, Giant mag-netoresistance: Theory, in
Handbook of spin transportand magnetism , edited by I. Zutic and E. Y. Tsymbal(CRC Press, 2011) pp. 95–114.[69] D. Y. Petrovykh, K. N. Altmann, H. Höchst, M. Laub-scher, S. Maat, G. J. Mankey, and F. J. Himpsel, Spin-dependent band structure, Fermi surface, and carrier life-time of permalloy, Appl. Phys. Lett73