# Supercurrent decay in ballistic magnetic Josephson junctions

Hervé Ness  Department of Physics, King’s College London, Strand Campus,
London WC2R 2LS, UK  Ivan A. Sadovskyy  Microsoft Quantum, Microsoft Station
Q, University of California, Santa Barbara, California 93106, USA  Andrey E.
Antipov  Microsoft Quantum, Microsoft Station Q, University of California,
Santa Barbara, California 93106, USA  Mark van Schilfgaarde  Department of
Physics, King’s College London, Strand Campus, London WC2R 2LS, UK  National
Renewable Energy Laboratory, Golden, Colorado 80401, USA  Roman M. Lutchyn
Microsoft Quantum, Microsoft Station Q, University of California, Santa
Barbara, California 93106, USA

###### Abstract

We investigate transport properties of ballistic magnetic Josephson junctions
and establish that suppression of supercurrent is an intrinsic property of the
junctions, even in absence of disorder. By studying the role of ferromagnet
thickness, magnetization, and crystal orientation we show how the supercurrent
decays exponentially with thickness and identify two mechanisms responsible
for the effect: (i) large exchange splitting may gap out minority or majority
carriers leading to the suppression of Andreev reflection in the junction,
(ii) loss of synchronization between different modes due to the significant
dispersion of the quasiparticle velocity with the transverse momentum. Our
results for Nb/Ni/Nb junctions are in good agreement with recent experimental
studies. Our approach combines density functional theory and Bogoliubov-de
Gennes model and opens a path for material composition optimization in
magnetic Josephson junctions and superconducting magnetic spin valves.

Magnetic Josephson junction, $\pi$-junction, JMRAM, ab initio calculations,
scattering theory.

##  1 Introduction

Coherent quantum tunneling of Cooper pairs through a thin barrier is one of
the first examples of macroscopic quantum coherent phenomena. Predicted by
Josephson more than 50 years ago [1], it has important applications in quantum
circuits used in metrology, quantum sensing and quantum information processing
[2].

Most of the previous studies focused on conventional Josephson junctions (JJs)
consisting of two $s$-wave superconductors (S) that are connected by an
insulating (I) or a normal (N) region [3, 4]. The flow of supercurrent through
a JJ depends on the superconducting phase difference $\phi$ between two
superconductors and, in general, is characterized by the current-phase
relationship $J(\phi)$ (CPR). In conventional JJs CPR should be periodic with
$2\pi$, $I(\phi)=I(\phi+2\pi)$ which follows from the BCS theory [5]. This
result is a manifestation of a $2e$ charge of Cooper pairs, and is used in
metrology to measure electron charge. Time-reversal symmetry requires that
$I(\phi)=-I(-\phi)$ which imposes a constraint that the supercurrent should be
zero for $\phi=\pi n$ where $n$ is integer. In general, CPR can be expanded in
Fourier harmonics, $I(\phi)=\sum_{n}I_{n}\sin(n\phi)$.

In many cases, however, CPR is well approximated by the first harmonic
$I(\phi)\approx I_{\mathrm{c}}\sin(\phi)$ with $I_{\mathrm{c}}$ being the
maximum supercurrent that can flow through the junction, i.e. the critical
current. At a microscopic level, the supercurrent through a short SNS junction
is determined by bound states forming in the constriction due to Andreev
reflection at the NS interfaces. In the Andreev reflection process an incident
electron-like quasiparticle with spin $\uparrow$ gets reflected at the NS
interface as a hole-like quasiparticle with spin $\downarrow$ and a Cooper
pair is emitted into the condensate. When time-reversal symmetry is not
broken, electrons and holes propagate with the same velocity in the normal
region. In this case no phase shift accumulates between this pair of
quasiparticles along their trajectories in the N region, and the sign of
$I_{\mathrm{c}}$ is fixed. When $I_{\mathrm{c}}>0$ we refer to this case as
$0$-junction.

Figure 1:  1 Schematic view of the SFS junction. 1 SFNFS junction. Arrows
indicate possible magnetization of ferromagnets. The supercurrent through the
spin-valve JJ depends on the relative magnetization of the ferromagnets, which
governs the properties of the Josephson magnetic random-access memory (JMRAM)
[6]. 1 Ball-and-stick representation of the Nb(110)/Ni(111)/Nb(110) junction
with 5 layers of Ni. Nb atoms are light grey, Ni atoms are blue. The top and
bottom atomic planes of Nb(110) are repeated periodically in the $z$-direction
to create the semi-infinite Nb leads of the junction through which the current
flows. Periodic boundary conditions are used in the $xy$-plane. The
corresponding reciprocal space defines two-dimensional
$\mathbf{k}_{\parallel}=(k_{x},k_{y})$ vectors, i.e. the transverse modes,
used in the calculations.

In a magnetic Josephson junction (MJJ), exchange splitting breaks time-
reversal symmetry and leads to an interesting interplay of superconductivity
and magnetism [4, 7, 8, 9, 10]. In superconductor-ferromagnet-superconductor
(SFS) junctions [Fig. 1] the correlated quasi-particles and quasi-holes
forming Andreev bound states propagate through the junction under the exchange
field of the ferromagnet (F). In many ferromagnets, such as Fe or Ni, the
exchange splitting is large (of the order of eV) and significantly perturbs
the band structure of a metal and, consequently, significantly modifies Fermi
velocities of minority (spin $\downarrow$) and majority (spin $\uparrow$)
carriers. Strong time-reversal symmetry breaking leads to the appearance of
characteristic superconducting correlations with oscillatory dependence
determined by the difference in wave numbers, $k_{\mathrm{\scriptscriptstyle
F}}^{\uparrow}-k_{\mathrm{\scriptscriptstyle F}}^{\downarrow}$ [11, 12]. This
effect opens a possibility for the supercurrent reversal as a function of the
thickness of the ferromagnetic region, the so-called Josephson $\pi$-junction
[13, 14]. The correlation between the phase shift of the supercurrent and the
magnetization provides a possibility for realizing magnetic spin valves, see
Fig. 1, which may have promising novel applications for cryogenic
superconducting digital technologies [15, 16, 6]. Understanding the
microscopic physics of MJJs is of great scientific interest as well as
technological importance. $0$-$\pi$ transitions in SFS junctions have been
extensively studied experimentally since the early 2000’s and have been
observed in different material systems [14, 17, 18, 19, 20, 21, 22, 23, 24,
25]. While qualitatively these observations are consistent with the previous
phenomenological theories [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42] the roles of the microscopic band structure arising from the
atomic lattice on the supercurrent suppression with junction thickness remain
unclear. Our primary goal is to address these essential points. The
supercurrent suppression that is exponential in the junction length is often
associated with presence of disorder in the ferromagnetic region [43, 7].
However, significant supercurrent suppression can also appear in relatively
clean metals (e.g., Ni) whose mean free path is larger than the junction
thickness [44, 22, 23].

Here we study the suppression of the critical current in the MJJs shown in
Fig. 1 and identify two microscopic mechanisms for its suppression, both a
consequence of the band structure asymmetry of the majority and minority
carriers in the F region. First, there is an asymmetry in the structure of the
Fermi surface, see Fig. 2. For certain bands and momenta, the Fermi surface
present in one spin channel may be absent in the other. This is typical in
ferromagnetic materials like Fe, Co, and Ni because the bandwidth for
$d$-electrons is relatively small and is often comparable to the exchange
splitting. As a result, the wave number of one of the constituent quasi-
particles forming Andreev bound states in MJJ becomes imaginary and the
supercurrent becomes suppressed. We label this scenario as mechanism (i). In a
second mechanism (ii), we show there is a dephasing of a harmonic signal
originating from different Fourier components of the supercurrent due the
Fermi velocity dispersion, see Fig. 3. Both these mechanisms lead to an
exponential suppression of the supercurrent, which was previously believed to
occur due to the presence of disorder in the magnetic region. We show that
band structure effects are important and may be even dominant in many cases.

In order to capture realistic band structure, we develop a microscopic theory
for the supercurrent in realistic MJJs. We use a combination of density
functional theory (DFT) and Bogoliubov-de Gennes (BdG) model to investigate
the $0$-$\pi$ transition in realistic material stacks of Nb/Ni/Nb junctions in
the clean limit. This method allows one to predict and explain key properties
of MJJs such as the period and decay of the critical current oscillations with
the ferromagnet thickness.

Figure 2:  2 Simplified band structure of a ferromagnet, with majority and
minority bands split by $V_{\mathrm{ex}}$. $k_{\mathrm{\scriptscriptstyle
F}}^{\uparrow}$ and $k_{\mathrm{\scriptscriptstyle F}}^{\downarrow}$ are the
Fermi momenta for majority and minority carriers, respectively. Fermi level
$E_{\mathrm{\scriptscriptstyle F}}^{\mathrm{(i)}}$ corresponds to large
$V_{\mathrm{ex}}$, where the minority band is pushed above
$E_{\mathrm{\scriptscriptstyle F}}^{\mathrm{(i)}}$.
$E_{\mathrm{\scriptscriptstyle F}}^{\mathrm{(ii)}}$ corresponds to small
$V_{\mathrm{ex}}$. Thus, the propagation of minority quasiparticles is
characterized by an imaginary momentum $k_{\mathrm{\scriptscriptstyle
F}}^{\downarrow}$ and is suppressed. 2 Supercurrent in SFS junction is carried
by Andreev bound states localized in the junction. Solid red line represents
quasi-classical trajectory corresponding to an Andreev bound state. The
spectrum of Andreev states depends on the relative superconducting phase
difference across the junction as well as the phase,
$\delta\varphi=|k_{\mathrm{\scriptscriptstyle
F}}^{\uparrow}-k_{\mathrm{\scriptscriptstyle F}}^{\downarrow}|w$, accumulated
due to the difference of Fermi momenta for majority and minority carriers.
Note that in scenario (i) the propagation of minority carriers is suppressed
leading to an overall exponential decay of the supercurrent with $w$. This is
to be contrasted with the normal transport through the junction.

The paper is organized as follows. We describe our findings and summarize our
main results in Sec. 2. In Sec. 3 we discuss the numerical method we developed
to perform first principles supercurrent calculations. Detailed discussion of
the main results is presented in Sec. 4. We conclude with Sec. 5. Some
technical points are relegated to Appendices A–E.

Figure 3:  3 First Fourier harmonic $J_{1}$ of the supercurrent density [Eq.
(3), black circles] and its fit [Eq. (14), solid black line] as a function of
Ni layer thickness, $w$, calculated for the Nb(110)/Ni(111)/Nb(110) junctions
shown in Fig. 1. Green semitransparent curves correspond to
$j_{1}^{\mathrm{fit}}(\mathbf{k}_{\parallel})$ for individual
$\mathbf{k}_{\parallel}$. 3 Normal state conductance per unit of area as a
function of $w$ for majority ($G_{\uparrow}$) and minority ($G_{\downarrow}$)
spins [Eq. (10)] as well as $G=G_{\uparrow}+G_{\downarrow}$. $G$ does not
depend on $w$ and is approximated by the single value $\langle G\rangle$. 3–3
Supercurrent as a function of phase difference $\phi$ for 3 4 atomic layers of
Ni (strong $0$-junction regime), 3 8 layers (intermediate regime), and 3 13
layers (strong $\pi$-junction regime). In the $0$\- and $\pi$-junction
regimes, the $J_{1}$ component dominates. In the intermediate regime higher-
order terms may prevail.

##  2 Qualitative discussion and main results

In this section we describe basic concepts for the supercurrent flow in MJJs
and summarize our results. Our main qualitative conclusions are supported by
microscopic calculations for Nb/Ni/Nb MJJs. Ni appears to have fairly long
mean free path $l_{\mathrm{\scriptscriptstyle MFP}}\approx 60$ Å (see
estimations in Appendix A), which is comparable or larger than the typical
thickness of the ferromagnet used in recent experiments [22, 23]. Therefore,
the motion of quasiparticles in Nb/Ni/Nb junction is quasi-ballistic, and our
method is applicable to this system. Most of this paper focuses on clean
Nb/Ni/Nb junctions.

First it is illuminating to consider a toy model, a one-dimensional SFS
junction, and calculate the supercurrent in such a system for different Fermi
energies, see Fig. 2. Using the results of Ref. 45, one finds that both
majority and minority spin bands are both occupied in the limit
$V_{\mathrm{ex}}/E_{\mathrm{\scriptscriptstyle F}}\ll 1$ [i.e. scenario (ii)
in Fig. 2] the supercurrent does not decay with ferromagnet thickness at zero
temperature,

$\displaystyle
I(\phi)=\frac{2e\Delta}{\hbar}\begin{cases}\cos\delta\varphi\sin\cfrac{\phi}{2},&0<\phi<\pi-2\delta\varphi,\\\
-\sin\delta\varphi\cos\cfrac{\phi}{2},&\pi-2\delta\varphi<\phi<\pi+2\delta\varphi,\\\
-\cos\delta\varphi\sin\cfrac{\phi}{2},&\pi+2\delta\varphi<\phi<2\pi.\end{cases}$
(1)

Here perfect interface transparency $\mathcal{T}$ is assumed,
$\mathcal{T}\approx 1$. The phase offset $\delta\varphi$ originates from the
Fermi momentum difference of a quasi-particle and a quasi-hole forming Andreev
bound state in the junction, see Fig. 2, and is given by
$\delta\varphi=|k_{\mathrm{\scriptscriptstyle
F}}^{\uparrow}-k_{\mathrm{\scriptscriptstyle F}}^{\downarrow}|w\approx
V_{\mathrm{ex}}w/\hbar v_{\mathrm{\scriptscriptstyle F}}$. The $0$\- and
$\pi$-junction regimes can be clearly identified at $\delta\varphi=0$ and
$\delta\varphi=\pi/2$, respectively. At the intermediate values
$0<\delta\varphi<\pi/2$, the CPR is anharmonic which is a generic feature at
$0$-$\pi$ transition as shown below. In the low transparency regime,
$\mathcal{T}\ll 1$, one would expect qualitatively similar results with the
maximal critical current being suppressed
$I_{\mathrm{c}}\sim(e\Delta/\hbar)\mathcal{T}$ but still independent of the
ferromagnet thickness, $w$.

In the case of large exchange splitting
$V_{\mathrm{ex}}/E_{\mathrm{\scriptscriptstyle F}}\gg 1$ [scenario (i) in Fig.
2] the minority band may become unoccupied. Given that minority carriers are
gapped out and their propagation through the junction is suppressed, the
supercurrent decays exponentially with $w$ [45],

$I(\phi)\approx\frac{2e\Delta}{\hbar}\exp(-\kappa
w)\Bigl{[}1-\frac{E_{\mathrm{\scriptscriptstyle
F}}}{8V_{\mathrm{ex}}}\sin^{2}(kw)\Bigr{]}\sin\phi.$ (2)

Here $\kappa=\sqrt{2m^{*}(V_{\mathrm{ex}}\\!-\\!E_{\mathrm{\scriptscriptstyle
F}})}/\hbar$ and
$k=\sqrt{2m^{*}(V_{\mathrm{ex}}\\!+\\!E_{\mathrm{\scriptscriptstyle
F}})}/\hbar$ with $m^{*}$ being effective electron mass. One may notice the
drastic difference between normal-state and superconducting transport in this
case — the former is weakly affected (because majority and minority
contributions are additive) whereas the supercurrent is strongly suppressed.
In this case, the measurement of normal-state junction resistance does not
necessarily predict the magnitude of the supercurrent through the junction.

We now generalize above results for the realistic three-dimensional (3D)
geometry and material composition of the MJJ. In the clean limit, the
supercurrent $I(\phi)$ in the short-junction limit ($w$ much smaller than the
coherence length of the superconductor) is obtained from the spectrum of the
Andreev bound states $\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})$
localized in the junction [3] which now also depends on the parallel momentum
$\mathbf{k}_{\parallel}$. The supercurrent density $J(\phi)$ per junction area
$A$ is $J(\phi)=I(\phi)/A$. For the junction with periodic atomic structure in
$xy$-plane the supercurrent density at zero temperature is given by

$J(\phi)=-\frac{e}{\hbar}\int\limits_{\mathrm{\scriptscriptstyle
BZ}}\\!\frac{d\mathbf{k}_{\parallel}}{(2\pi)^{2}}\,\sum_{\nu>0}\frac{\partial\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})}{\partial\phi},$
(3)

where the $\mathbf{k}_{\parallel}$ integration is performed over the Brillouin
zone (BZ) of the corresponding surface supercell of area, $A$, and the sum is
carried over positive quasiparticle energies,
$\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})>0$. Note that we use spin-
resolved $\varepsilon_{\nu}$ and therefore omit spin prefactor 2 in Eq. (3).
The derivative is periodic in $\phi$ and can be represented as a Fourier
series

$-\frac{e}{\hbar}\,\frac{\partial\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})}{\partial\phi}=\sum_{n\geqslant
1}I_{n\nu}(\mathbf{k}_{\parallel})\sin\bigl{[}n\phi+\delta\varphi_{n\nu}(\mathbf{k}_{\parallel})\bigr{]},$

so that Eq. (3) can be written as

$J(\phi)=\int\limits_{\mathrm{\scriptscriptstyle
BZ}}\\!\frac{d\mathbf{k}_{\parallel}}{(2\pi)^{2}}\,\\!\sum_{\nu>0}\,\sum_{n\geqslant
1}I_{n\nu}\sin\bigl{[}n\phi+\delta\varphi_{n\nu}(\mathbf{k}_{\parallel})\bigr{]}.$

As we will show below, away from $0$-$\pi$ transition the supercurrent is
dominated by the first ($n=1$) harmonic. Therefore, we focus henceforth on the
first harmonic contribution and drop $n$ index in the following discussion.
Next, one may notice that the supercurrent amplitudes
$I_{\nu}(\mathbf{k}_{\parallel})$ and phase offsets
$\delta\varphi_{\nu}(\mathbf{k}_{\parallel})$ depend on the parallel momentum
$\mathbf{k}_{\parallel}$. One may include this dependence and define an
effective Fermi energy $E_{\mathrm{\scriptscriptstyle
F}}(\mathbf{k}_{\parallel})$ that counts the energy corresponding to
$\mathbf{k}_{\parallel}$ in each band of the ferromagnet from the bottom of
the band. Depending on $E_{\mathrm{\scriptscriptstyle
F}}(\mathbf{k}_{\parallel})$ and $V_{\mathrm{ex}}$, either scenario (i) or
(ii) of Fig. 2 may be realized.

Thus, it is important to compare the exchange splitting with the bandwidth of
$d$-character states in transition metals to make sure that a perturbation
theory in $V_{\mathrm{ex}}$ is justified. Assuming it is the case, one can
estimate the phase offset $\delta\varphi_{\nu}(\mathbf{k}_{\parallel})$ by
expanding in exchange splitting to find

$\delta\varphi_{\nu}(\mathbf{k}_{\parallel})\approx V_{\mathrm{ex}}w/\hbar
v_{z\nu}(\mathbf{k}_{\parallel}).$ (4)

In general, the dependence of $v_{z\nu}$ on $\mathbf{k}_{\parallel}$ is
complicated, especially in $spd$-transition metals. The combination of
complicated amplitude and phase offset dependence on $\mathbf{k}_{\parallel}$
leads to a non-trivial supercurrent dependence on the ferromagnet thickness
$w$. As shown in Fig. 3, the critical current decays with $w$ for the
Nb(110)/Ni(111)/Nb(110) junctions. We analyze the details of the decay and
perform an exponential fit in Sec. 4. At small thicknesses, below 30 Å, this
decay originates from the evanescent modes corresponding to gapped out
minority or majority carriers which cannot propagate through the junction, see
Table 1. This is the mechanism (i) discussed above. At larger thicknesses
(i.e. $w\gtrsim 50$ Å) a loss of synchronization between different modes due
to the dispersion of $v_{z\nu}(\mathbf{k}_{\parallel})$ becomes important.
This second mechanism (ii) has been previously discussed in the literature
[13, 7] under assumptions of a single spherical Fermi surface and a small
uniform exchange splitting $V_{\mathrm{ex}}$ in the magnetic region. Within
these assumptions, one finds that critical current should decay algebraically
with the thickness $w$ of a magnetic layer [13]. However, as we show below
most of these assumptions do not apply to realistic SFS junctions involving
transition metals. Thus, in order to understand CPR in realistic MJJs, one
needs to use accurate ab initio methods, which capture the physical effects
described above.

Figure 4:  4 Electronic band structure of bulk Ni (solid lines) calculated
from first principles including many-body effects, see Ref. 46. It is the
highest fidelity available and is very close to ARPES data (diamonds) in both
majority (red) and minority (green) spin bands, with exchange splitting
$V_{\mathrm{ex}}=0.3$ eV. 4 Majority (solid line) and minority (dashed line)
Fermi surfaces of bulk Ni in the $\mathbf{k}_{\parallel}$ plane (with
$k_{z}=0$) corresponding to the (111) plane direction used for the stacking of
the Nb/Ni/Nb junctions, as discussed in the text. Axes correspond to two
$\Gamma$-L lines: the Fermi surface has a three-fold symmetry in the entire
plane. Majority band $6^{\uparrow}$ is depicted as solid red-blue line with
color interpolating between red and blue depending on the Fermi velocity
$\hbar^{-1}\partial E/\partial k$, which ranges between $3\times 10^{5}$ m/s
(red) and $6\times 10^{5}$ m/s (blue). Minority bands are shown by dashed
lines: $3^{\downarrow}$ (blue), $4^{\downarrow}$ (cyan), $5^{\downarrow}$
(green), and $6^{\downarrow}$ (red).

To make a connection between the decay seen in Fig. 3 and the mechanisms
responsible for it, Fig. 4 presents the energy band structure of bulk Ni. It
is probably the highest fidelity band structure available: it very closely
reproduces ARPES data in both majority and minority spin bands, with exchange
splitting $V_{\mathrm{ex}}=0.3$ eV [46], and should be an excellent predictor
of the real Fermi surface and velocities in Ni.

| Majority ($\uparrow$) | Minority ($\downarrow$)
---|---|---
| $v_{\mathrm{\scriptscriptstyle F}}^{\mathrm{min}}$ | $v_{\mathrm{\scriptscriptstyle F}}^{\mathrm{max}}$ | $\langle v_{\mathrm{\scriptscriptstyle F}}\rangle$ | $\sqrt{\langle v_{\mathrm{\scriptscriptstyle F}}^{2}\rangle}$ | $\rho(E_{\mathrm{\scriptscriptstyle F}})$ | $\Delta{E}$ | $m^{*}\\!/m_{\mathrm{e}}$ |  $\kappa$ | $v_{\mathrm{\scriptscriptstyle F}}^{\mathrm{min}}$ | $v_{\mathrm{\scriptscriptstyle F}}^{\mathrm{max}}$ | $\langle v_{\mathrm{\scriptscriptstyle F}}\rangle$ | $\sqrt{\langle v_{\mathrm{\scriptscriptstyle F}}^{2}\rangle}$ | $\rho(E_{\mathrm{\scriptscriptstyle F}})$
Band | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [eV-1Å-3] | [eV] |  | [Å] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [$10^{5}\frac{\mathrm{m}}{\mathrm{s}}$] | [eV-1Å-3]
2 |  |  |  |  |  | $-1.59$ | 0.68 | 12 |  |  |  |  |
3 |  |  |  |  |  | $-0.22$ | 3.90 | 13 | 2.7 | 3.6 | 3.3 | 3.6 | 0.004
4 |  |  |  |  |  | $-0.11$ | 0.61 | 47 | 0.8 | 1.8 | 1.1 | 1.2 | 0.015
5 |  |  |  |  |  | $-0.10$ | 2.95 | 22 | 0.6 | 3.0 | 1.5 | 1.8 | 0.173
6 | 3.5 | 6.0 | 4.6 | 5.2 | 0.029 |  |  |  | 0.6 | 3.0 | 2.3 | 2.6 | 0.045
Table 1:  Minimum, maximum, average, and root mean square of the Fermi
velocities $v_{\mathrm{\scriptscriptstyle F}}$ for the majority and minority
Fermi surfaces in bulk Ni. $\rho(E_{\mathrm{\scriptscriptstyle F}})$ is the
density of states at the Fermi level, $\Delta{E}$ is the energy of the valence
band maximum relative to $E_{\mathrm{\scriptscriptstyle F}}$ for the majority
bands not crossing $E_{\mathrm{\scriptscriptstyle F}}$, $m^{*}$
($m_{\mathrm{e}}$) is the effective (bare) electron mass, and $\kappa$
estimates the decay exponent of the evanescent mode for a given band at
$E_{\mathrm{\scriptscriptstyle F}}$.

Here we use this potential to analyze the bulk Ni Fermi surface [Fig. 4] and
Fermi velocities (Table 1). Counting from the bottom $s$-band, bands of Ni $d$
character are bands 2 to 6. These bands are nearly full: only majority band 6↑
and minority bands $3^{\downarrow}$–$6^{\downarrow}$ cross the Fermi level
$E_{\mathrm{\scriptscriptstyle F}}$. Only band 6 has both majority and
minority carriers at the Fermi surface. Bands 6↑ and 6↓ have roughly the same
shape and the energy splitting is approximately constant and equal to
$V_{\mathrm{ex}}$ (see the right panel of Fig. 4 in Ref. 46). Beyond this,
however, the correspondence between the Ni band structure and a simple
parabolic band structure deviate substantially, in two ways that critically
affect the analysis. First, the Fermi velocity, $v_{\mathrm{\scriptscriptstyle
F}}={\hbar}^{-1}(\partial E/\partial k)|_{k=k_{\mathrm{\scriptscriptstyle
F}}}$, is not fixed even for a single pocket: it varies in band $6^{\uparrow}$
by a factor of 2 [see Fig. 4 and Table 1]. Accordingly the splitting
$k^{\uparrow}_{\mathrm{\scriptscriptstyle
F}}-k^{\downarrow}_{\mathrm{\scriptscriptstyle F}}$ between bands
$6^{\uparrow}$ and $6^{\downarrow}$ varies by factor of two as expected from
the twofold variation in $v_{\mathrm{\scriptscriptstyle F}}$. Second, bands
$3^{\downarrow}$–$5^{\downarrow}$ have no majority counterpart at
$E_{\mathrm{\scriptscriptstyle F}}$, indicating that the wave number of bands
$3^{\uparrow}$–$5^{\uparrow}$ is complex. This is the origin for the
exponential decay in scenario (i) in Fig. 2 as noted above: a large portion of
Andreev levels are carried by Cooper pairs made of single-particle wave
functions with one or both of $k_{\mathrm{\scriptscriptstyle F}}^{\uparrow}$
and $k_{\mathrm{\scriptscriptstyle F}}^{\downarrow}$ having an imaginary
component. The magnitude of $\mathrm{Im}\,k$ depends on the particular mode
and $\mathbf{k}_{\parallel}$ leading to a distribution of decay exponents. The
slowest decay in each of these evanescent modes can be estimated from the
distance $\Delta{E}$ of the closest approach to $E_{\mathrm{\scriptscriptstyle
F}}$ and the effective mass $m^{*}/m_{e}$, using
$\hbar^{2}k_{\mathrm{min}}^{2}/2m^{*}=\Delta{E}$ and decay
$\kappa=2\pi/\mathrm{Im}(k_{\mathrm{min}})$, see Table 1. ($m^{*}/m_{e}$ is
found to be highly anisotropic, so only the effective transport mass
$m^{*}=3\,[1/m_{1}^{*}+1/m_{2}^{*}+1/m_{3}^{*}]^{-1}$ is shown.) $\kappa$ is
only a rough measure of the evanescent mode decay for a given band. We discuss
in detail the distribution of decay exponents and phase offsets in Sec. 4.

Let us now focus on the mechanism (ii) for the supercurrent decay, i.e. loss
of synchronization between different transverse modes. This mechanism is well-
known in diffusive systems where quasiparticle trajectory is random and thus
the phase offset $\delta\varphi$ accumulated along such a trajectory also gets
randomized. Thus, upon averaging Eq. (3) over different disorder realizations,
one ends up with exponentially decaying critical current [7]. Previously, such
a suppression of the supercurrent with junction thickness, $w$, was often
associated with impurity scattering in the ferromagnet. Here we demonstrate
that this dephasing mechanism can also appear in clean systems (where
quasiparticle trajectory is well defined) due to the dispersion of the
velocity $v_{z\nu}$ with an in-plane momentum $\mathbf{k}_{\parallel}$.
Specifically, we find that, in the Nb/Ni/Nb junctions, this mechanism becomes
relevant for junctions thicker than $50$ Å, see Fig. 3 and Sec. 4. In SFS
junctions the combination of disorder in the ferromagnet, interface scattering
as well band-structure-induced dephasing ultimately determines the magnitude
of the supercurrent. However, we believe that band structure effects provide
an upper bound on the magnitude of the critical current as a function of $w$.

In Sec. 4, we present numerical results which support the qualitative
discussion presented above.

##  3 Method

We develop a numerical method to perform realistic simulations of MJJs using a
combination of first-principles DFT and BdG calculations. The former is used
to obtain the normal-state properties (e.g., band structure, Fermi velocities,
magnetization) and to calculate the normal scattering matrices through the
inhomogeneous 3D realistic junctions. As a next step we take superconductivity
into account and calculate supercurrent through the stack assuming the short
junction limit.

###  3.1 Normal transport: ab initio description

To calculate the normal scattering matrix we use the Questaal package for
electronic structure calculations based on the linear muffin-tin orbital
(LMTO) method [47]. It calculates the full non-linear, i.e. non equilibrium,
transport properties of an infinite system describing a central (C) region
cladded by two semi-infinite left (L) and right (R) leads [48, 49], as
represented below:

$\overbrace{\ldots\,|\,\mathcal{L}\,|\,\mathcal{L}\,}^{\textstyle\mathrm{L}}\,|\,\overbrace{\,\mathrm{PL}_{0}\,|\,\mathrm{PL}_{1}\,|\,\ldots\,|\,\mathrm{PL}_{L-1}\,}^{\textstyle\mathrm{C}}\,|\,\overbrace{\,\mathcal{R}\,|\,\mathcal{R}\,|\,\ldots}^{\textstyle\mathrm{R}}$

The LCR system [Fig. 1] can be partitioned into an infinite stack of principal
layers (PLs) which interact only with their nearest-neighbors. This is
possible because the screened LMTO structure constants are short-ranged [50].
In the present case the C region consists of the ferromagnet, plus two layers
of Nb at the LC and CR interfaces respectively. This is the range over which
the perturbation from C significantly modifies the potential in the L or R
region. To construct the Nb/Ni/Nb stack, coincident site lattices for Ni and
Nb must be found (details of how this was accomplished are given in Appendix
B). Planes of coincident site lattices are stacked to form the Nb/Ni/Nb
structures. Figure 5 shows the Nb and Ni planes we used, which are denoted
here as ‘surface supercells.’

The electronic current flows along the $z$ direction, perpendicular to the PLs
lying in the $xy$-plane (transverse direction), see Fig. 1. Periodic boundary
conditions are used within each PL. The corresponding reciprocal space defines
the two-dimensional (2D) $\mathbf{k}_{\parallel}$ vectors, i.e. the transverse
modes, used in the calculations. The $\mathbf{k}_{\parallel}$ mesh is
discretized and integrals over $\mathbf{k}_{\parallel}$ are performed
numerically.

The electronic structure of the C region can be separated from L and R regions
through self-energies, $\Sigma_{\mathrm{\scriptscriptstyle L}}$ and
$\Sigma_{\mathrm{\scriptscriptstyle R}}$, that modify the Hamiltonian of C
region. They are most easily calculated if the potential of each PL in the L
or R region are identical all through the bulk region. This is the reason for
adding a few Nb layers folded into the C region. Thus the periodically
repeating unit cells in the L and R regions can be safely assumed to have the
potential of the bulk crystal. To construct the self-energies, the potentials
of the PL in an infinite stack are needed. These potentials are functions only
of the PL in their own region, and may be calculated in several ways.
$\Sigma_{\mathrm{\scriptscriptstyle L}}$ and
$\Sigma_{\mathrm{\scriptscriptstyle R}}$ are obtained from ‘surface’ Green’s
function (a fictitious system which consists of a semi-infinite stack of PL,
each with the same potential). Note that the potential of the C region is
calculated self-consistently. This is important, as the local moments of Ni
are small at the boundary layers, and build up gradually, see Fig. 5.

With the potential in hand, the normal-state transport can be calculated using
scattering formalism [49, 51]. For this, knowing the retarded Green’s
function, $\mathcal{G}^{\mathrm{r}}$, of the junction is sufficient. However,
this is not the case for the Josephson current: the individual eigenfunctions
are required. Within a Green’s function framework, $\mathcal{G}^{\mathrm{r}}$
must be organized by normal modes which correspond to the eigenstates of the L
and R leads for a prescribed energy $E$. In the PL representation, the
Hamiltonian has been discretized into the linear combinations of the LMTO
basis functions, and the normal modes are represented as eigenvectors of these
basis functions. The Schrödinger equation becomes a difference equation in the
PL basis functions [52]. Eigenvectors are calculated by solving a quadratic
eigenvalue problem [52, 53], whose eigenvalues correspond to $\exp(\pm
ik_{z,n}a)$, where $a$ is the thickness of the PL. The wave number $k_{z,n}$
of the normal mode $n$ can be complex, but to correspond to a propagating mode
$k_{z,n}$ must be real. By solving the equation as a function of the energy
$E$, one gets all the eigenvalues and eigenvectors which provide the
information needed to construct the self-energies
$\Sigma_{\mathrm{\scriptscriptstyle L}}$ and
$\Sigma_{\mathrm{\scriptscriptstyle R}}$ (for each $\mathbf{k}_{\parallel}$
and each spin $\sigma$). Note that, in the mode basis, the imaginary part of
the self-energies is proportional to the (band) velocity of the modes, and is
diagonal for non-degenerate modes [53, 54].

The retarded Green’s function $\mathcal{G}^{\mathrm{r}}$ of the C region
(connected to the L and R leads) can be written as a matrix in the normal mode
basis,

$\mathcal{G}^{\mathrm{r}}_{\sigma}(E,\mathbf{k}_{\parallel})=\Bigl{\\{}[\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
C}},\sigma}(E,\mathbf{k}_{\parallel})]^{-1}\\\
-\Sigma_{{\mathrm{\scriptscriptstyle
L}},\sigma}(E,\mathbf{k}_{\parallel})-\Sigma_{{\mathrm{\scriptscriptstyle
R}},\sigma}(E,\mathbf{k}_{\parallel})\Bigr{\\}}^{-1},$ (5)

where $\mathcal{G}^{\mathrm{r}}_{\mathrm{\scriptscriptstyle C}}$ is the
Green’s function of the isolated C region. In this basis,
$\mathcal{G}^{\mathrm{r}}$ is decomposed into four blocks,

$\mathcal{G}^{\mathrm{r}}_{\sigma}(E,\mathbf{k}_{\parallel})=\left[\begin{array}[]{cc}\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
LL}},\sigma}(E,\mathbf{k}_{\parallel})&\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
LR}},\sigma}(E,\mathbf{k}_{\parallel})\\\
\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
RL}},\sigma}(E,\mathbf{k}_{\parallel})&\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
RR}},\sigma}(E,\mathbf{k}_{\parallel})\end{array}\right]\\!,$ (6)

upon projecting onto the propagating modes of the L and R leads. These four
quantities and the mode velocities completely determine the normal state
transport properties of the junctions.

The transmission matrices are defined by the off-diagonal parts of Eq. (6).
More specifically, the transmission coefficients
$[t_{{\mathrm{\scriptscriptstyle LR}},\sigma}]_{nm}$, connecting L and R
regions, are given by [51]

$[t_{{\mathrm{\scriptscriptstyle
LR}},\sigma}]_{nm}(E,\mathbf{k}_{\parallel})=i\sqrt{|[v_{{\mathrm{\scriptscriptstyle
L}},\sigma}]_{n}(E,\mathbf{k}_{\parallel})|}\\\
\times[\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
LR}},\sigma}]_{nm}(E,\mathbf{k}_{\parallel})\,\sqrt{|[v_{{\mathrm{\scriptscriptstyle
R}},\sigma}]_{m}(E,\mathbf{k}_{\parallel})|},$ (7)

where $[v_{{\mathrm{\scriptscriptstyle L}},\sigma}]_{n}$ and
$[v_{{\mathrm{\scriptscriptstyle R}},\sigma}]_{m}$ are the velocity matrix
elements for propagating modes $n$ and $m$ in the L and R leads, respectively.
The transmission matrix $t_{\mathrm{\scriptscriptstyle RL}}$ can be obtained
from Eq. (7) by replacing $\mathrm{L}\leftrightarrow\mathrm{R}$. The
reflection coefficients are given by the diagonal blocks of Eq. (6). For
instance, on the L side

$[r_{{\mathrm{\scriptscriptstyle
LL}},\sigma}]_{nn^{\prime}}(E,\mathbf{k}_{\parallel})=i\sqrt{|[v_{{\mathrm{\scriptscriptstyle
L}},\sigma}]_{n}(E,\mathbf{k}_{\parallel})|}\\\
\times[\mathcal{G}^{\mathrm{r}}_{{\mathrm{\scriptscriptstyle
LL}},\sigma}]_{nn^{\prime}}(E,\mathbf{k}_{\parallel})\,\sqrt{|[v_{{\mathrm{\scriptscriptstyle
L}},\sigma}]_{n^{\prime}}(E,\mathbf{k}_{\parallel})|}-\delta_{nn^{\prime}}.$
(8)

The reflection matrix $r_{\mathrm{\scriptscriptstyle RR}}$ on the R side is
obtained from Eq. (8) by replacing $\mathrm{L}\leftrightarrow\mathrm{R}$.

We define the normal scattering matrix as

$S=\left[\begin{array}[]{cc|cc}r_{{\mathrm{\scriptscriptstyle
LL}},\uparrow}&0&t_{{\mathrm{\scriptscriptstyle LR}},\uparrow}&0\\\
0&r_{{\mathrm{\scriptscriptstyle
LL}},\downarrow}&0&t_{{\mathrm{\scriptscriptstyle LR}},\downarrow}\\\
\hline\cr t_{{\mathrm{\scriptscriptstyle
RL}},\uparrow}&0&r_{{\mathrm{\scriptscriptstyle RR}},\uparrow}&0\\\
0&t_{{\mathrm{\scriptscriptstyle
RL}},\downarrow}&0&r_{{\mathrm{\scriptscriptstyle
RR}},\downarrow}\end{array}\right]\\!,$ (9)

where we omit $E$\- and $\mathbf{k}_{\parallel}$-dependence for brevity.

The linear-response normal conductance, $G_{\sigma}$, per spin is given by

$\frac{G_{\sigma}}{A}=\frac{e^{2}}{h}\int\limits_{\mathrm{\scriptscriptstyle
BZ}}\\!\frac{d\mathbf{k}_{\parallel}}{(2\pi)^{2}}\,\sum_{n,m}\bigl{|}[t_{{\mathrm{\scriptscriptstyle
LR}},\sigma}]_{nm}(E_{\mathrm{\scriptscriptstyle
F}},\mathbf{k}_{\parallel})\bigr{|}^{2}.$ (10)

It is calculated at the Fermi energy, $E=E_{\mathrm{\scriptscriptstyle F}}$.
The total conductance is given by $G=G_{\uparrow}+G_{\downarrow}$. Figure 3
shows that $G$, $G_{\uparrow}$, and $G_{\downarrow}$ are almost independent of
the junction thickness, $w$.

###  3.2 Superconducting transport: scattering matrix approach

Equation (9) is the normal-state scattering matrix for metal-ferromagnet-metal
(NFN) structure taking into account reflection at both NF interfaces. To
account for superconductivity, we introduce a step-like superconducting
pairing potential $\Delta=3.1$ meV and use the Andreev approximation to
account for electron-hole scattering processes [3, 10]. This approach combines
the details of the atomic structure of Nb/Ni/Nb and superconductivity within
the mean-field approximation.

The direct contact between S and F layers, as shown in Fig. 1, would lead to
an interaction between them. The back-action of the ferromagnet on the
superconductor (i.e. the inverse proximity effect) results in a spatial
dependence of the pairing potential near the SF interfaces [55, 56, 57]. In
the case of a clean SFS junction model, the self-consistent BdG calculations
have been discussed in Refs. 58, 59, 60, 61, 62. In typical experimental
systems the superconductor is disordered, so disorder effect on pairing
potential needs to be considered, see, e.g., Ref. 63. Furthermore, in recent
experiments S and F layers are separated by an intermediate spacer layer,
which significantly reduces this inverse proximity effect. Thus, to understand
inverse proximity effect in realistic SFS devices one would need to include
both of the abovementioned ingredients in the model as well as to take into
account the inhomogeneous magnetization in the ferromagnet (discussed in the
next section) which is outside the scope of this paper. For the sake of
clarity, we focus here on a realistic band structure in the ferromagnet and
its effect on the supercurrent in the SFS structures. Henceforth, we also
neglect the orbital effects of the fringe magnetic field, created by the
ferromagnet.

Since $\Delta\ll E_{\mathrm{\scriptscriptstyle F}}$, spin-resolved Andreev
reflection at SNL and NRS interfaces is described by

$r_{\mathrm{\scriptscriptstyle
A}}(\phi)=\left[\begin{array}[]{cc|cc}0&\mathbf{1}\,e^{i\phi/2}&0&0\\\
\mathbf{1}\,e^{i\phi/2}&0&0&0\\\ \hline\cr 0&0&0&\mathbf{1}\,e^{-i\phi/2}\\\
0&0&\mathbf{1}\,e^{-i\phi/2}&0\end{array}\right]\\!,$ (11)

where $\phi$ is the phase difference between left and right superconducting
leads and $\mathbf{1}$ is the identity matrix. In the short junction limit,
the main contribution to the supercurrent comes from Andreev bound states
localized in the junction having energy $\varepsilon_{\nu}$. The energy
spectrum of Andreev states can be obtained using the following equation [64],

$\alpha(\varepsilon)\\!\left[\\!\begin{array}[]{cc}0&\\!\\!r_{\mathrm{\scriptscriptstyle
A}}^{*}(\phi)\\\ r_{\mathrm{\scriptscriptstyle
A}}(\phi)\\!\\!&0\end{array}\\!\right]\\!\\!\left[\\!\begin{array}[]{cc}S(E_{\mathrm{\scriptscriptstyle
F}}\\!+\\!\varepsilon,\mathbf{k}_{\parallel})\\!\\!\\!&0\\\
0&\\!\\!\\!S^{*}(E_{\mathrm{\scriptscriptstyle
F}}\\!-\\!\varepsilon,\mathbf{k}_{\parallel})\end{array}\\!\right]\\!\Psi^{\mathrm{in}}\\!=\\!\Psi^{\mathrm{in}},$
(12)

where
$\alpha(\varepsilon)=\sqrt{1-\varepsilon^{2}/\Delta^{2}}+i\varepsilon/\Delta$.
The vector
$\Psi^{\mathrm{in}}=[\psi_{\mathrm{e\uparrow}}^{{\mathrm{\scriptscriptstyle
L}}\rightarrow},$ $\psi_{\mathrm{e\downarrow}}^{{\mathrm{\scriptscriptstyle
L}}\rightarrow},$ $\psi_{\mathrm{e\uparrow}}^{{\mathrm{\scriptscriptstyle
R}}\leftarrow},$ $\psi_{\mathrm{e\downarrow}}^{{\mathrm{\scriptscriptstyle
R}}\leftarrow},$ $\psi_{\mathrm{h\uparrow}}^{{\mathrm{\scriptscriptstyle
L}}\rightarrow},$ $\psi_{\mathrm{h\downarrow}}^{{\mathrm{\scriptscriptstyle
L}}\rightarrow},$ $\psi_{\mathrm{h\uparrow}}^{{\mathrm{\scriptscriptstyle
R}}\leftarrow},$ $\psi_{\mathrm{h\downarrow}}^{{\mathrm{\scriptscriptstyle
R}}\leftarrow}]^{\mathrm{T}}$ corresponds to the electron- and hole-like (e/h)
waves in NL and NR regions incident on the F region from the left
($\rightarrow$) and from the right ($\leftarrow$).

Simulations show that $S$ is weakly-dependent on $E$ in the range
$[E_{\mathrm{\scriptscriptstyle F}}-\Delta,E_{\mathrm{\scriptscriptstyle
F}}+\Delta]$. Therefore, we expand $S(E,\mathbf{k}_{\parallel})$ in
$E-E_{\mathrm{\scriptscriptstyle F}}$ and keep only the leading term, i.e.
$S(E,\mathbf{k}_{\parallel})\approx S(E_{\mathrm{\scriptscriptstyle
F}},\mathbf{k}_{\parallel})=S(\mathbf{k}_{\parallel})$. Using this
approximation, one can simplify the quantization condition (12) and reduce it
to the matrix eigenvalue problem (see details in Ref. 65). This approach
allows one to reliably calculate the Andreev bound states spectrum,
$\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})$. The zero-temperature
supercurrent, $J$, through the junction is given by Eq. (3). Figures 3–3 show
$J$ as a function of a phase difference, $\phi$. Figure 3 shows the first
Fourier harmonic of the supercurrent as a function of the junction thickness,
$w$.

##  4 Results

It is illuminating to compare our numerical simulations for the supercurrent
with the experimental measurements involving quasi-ballistic MJJs. As
previously discussed, we believe that the Nb/Ni/Nb junctions represent a good
model system for which experimental data is readily available [16, 22, 23].
The best-performing stacks consist of Nb(110)/Cu/Ni(111)/Cu/Nb(110). Cu spacer
layers 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 supercells in the plane normal to the stack to account for the
lattice mismatch (see Fig. 1 and Appendix B), though we do not include the Cu
layers. We anticipate that Cu spacers will mainly affect transmission matrix
elements rather than the dependence of the supercurrent on ferromagnet
thickness, which is the main focus of this work. Furthermore, as discussed
before, the Cu spacers will suppress the direct interaction between the
ferromagnet and the superconductor and reduce inverse proximity effect
justifying Andreev approximation for the boundary conditions, see Eq. (11).
Therefore, we consider only the simplified Nb/Ni/Nb stack and vary the number
of layers (atomic planes) of Ni. Additionally, we also consider effect of
different crystallographic orientation of the Ni planes and investigate
Nb/Ni(110)/Nb junctions in Appendix C.

To make a Nb/Ni superlattice, the unit cells of the Nb and Ni regions in the
plane normal to the interface must be coincident. This is complicated by the
severe lattice constant mismatch, and also the incompatibility of the (110)
and (111) atomic planes. It is necessary to construct superlattices with
Nb(110) and Ni(111) both rotated to the $z$ axis, and with lattice vectors in
the plane coincident. A supercell with nearly coincident vectors was found
(see Appendix B for details). By applying a small shear strain to the Ni, the
lattice vectors are made exactly coincident. Figure 5 shows the Nb(110)
surface supercell and the Ni(111) surface supercell with equal lattice vectors
used to match the Nb/Ni interfaces. Each atomic plane of Ni(111) contains 14
atoms and each atomic plane of Nb(110) consists of 10 atoms. The atomic
structure of the Nb(110)/Ni(111)/Nb(110) for 5 layers (atomic planes) of Ni is
shown in Fig. 1.

Figure 5:  5 Top view of the Nb(110) and Ni(111) surface supercells used to
build the Nb/Ni interfaces and the Nb(110)/Ni(111)/Nb(110) stacks shown in
Fig. 1. The surface supercell is defined from two 2D vectors
$\mathbf{a}_{1}=[10.8,0]$ Å and $\mathbf{a}_{2}=[6.03,7.11]$ Å with periodic
boundary condition in the 2D plane. The corresponding reciprocal space defines
the 2D $\mathbf{k}_{\parallel}$ vectors used in the calculations. Each atomic
plane of Ni (Nb) contains 14 (10) atoms of Ni (Nb). 5 Magnetic moment profile
of Nb(110)/Ni(111)/Nb(110) junctions for different thickness of Ni (from 3 to
9 layers). The value of the moment is an averaged over the moments of the 14
Ni atoms in each atomic plane. Note the magnetic dead layer at the Nb/Ni
interfaces and that all moments vanish for the shortest junction made of 3
layers of Ni.

Next, we performed self-consistent DFT calculations within the local density
approximation (LDA) in order to obtain the relaxed structure and corresponding
electronic structure. For the smallest structures we performed a constrained
optimization. Only the atoms in the planes closest to the Nb/Ni interfaces are
allowed to relax to facilitate stacking of arbitrarily large cells. The Nb/Ni
interplanar spacing has also been optimized to minimize the total energy, see
Appendix E for more details.

Once the structure is determined, one can determine the normal-state
thermodynamic and transport properties of the junction, e.g., calculate the
magnetization profile and spin-resolved conductance through the junction as a
function of Ni thickness. For transport calculations, we use a layer transport
technique [48] which employs the atomic spheres approximation (ASA). Careful
checks were made of ASA band structures of elemental Nb and Ni, and also
superlattices, to confirm that they are very similar to the full potential
LMTO DFT-LDA ones.

We find that the magnetic properties of Ni are sensitive to their local
environment, indicative of the itinerant ferromagnetism. As shown in Fig. 5,
the magnetization profile is non-uniform in the junction with averaged
magnetic moments per atom being suppressed near the Nb interface. For
thickness larger than 4 layers, one recovers the bulk value of $\sim
0.6\mu_{\mathrm{\scriptscriptstyle B}}$ in the middle layers, away from the
Nb/Ni interfaces. The averaged moment drops down towards the edges and becomes
considerably reduced down to $\sim 0.1\mu_{\mathrm{\scriptscriptstyle B}}$ at
the interface with Nb. The strong reduction of magnetism is exemplified for
the short junction with 3 layers of Ni where the moments on the Ni atoms have
completely vanished. Such a non-uniform magnetic moment dependence in Nb/Ni/Nb
junctions affect superconducting properties of the SFS junctions in a non-
trivial way. For example, Nb/Ni/Nb junctions thinner than 4 layers of Ni
behave as essentially SNS junctions.

It is well known that the LDA tends to overestimate local moments $M$ in
itinerant magnets [66] because spin fluctuations reduce the average moment
[67], and underestimate $M$ when local moments are very large [46]. For Ni,
LDA yields $M$ in good agreement with the experiment, but this is likely an
artifact of an accidental cancellation of errors. Most important for transport
is the exchange splitting $V_{\mathrm{ex}}$, which the LDA predicts to be 0.6
eV, about twice larger than the experimental value of 0.3 eV [68]. It is
possible to reproduce both $M$ and $V_{\mathrm{ex}}$ at the same time, but a
high-level theory, potentially including spin-orbit coupling, is needed to
surmount both kinds of errors inherent the LDA [46, 69]. The high cost and
poor scaling of such a theory is not practical for these junctions, so we
elect to stay within the LDA and scale the self-consistently calculated
$V_{\mathrm{ex}}$. This was the approach Karlsson and Aryasetiawan used to
calculate the spin wave spectra in Ni [70]. Scaling of $V_{\mathrm{ex}}$ can
be accomplished using different approaches, e.g., by adding some effective
magnetic field to simulate the effect of spin fluctuations. Since Ni is a
simple case with a nearly linear relation between $M$ and $V_{\mathrm{ex}}$,
the band structure hardly depends on the details in which the LDA potential is
modified. Here we first perform fully self-consistent calculations. Then, to
construct the potential for transport properties, we rescale the spin
component of the density by a constant factor, which we denote as $M/M_{0}$.
This enables parametric studies of transport as a function of
$V_{\mathrm{ex}}$. $M/M_{0}=0.5$ yields the observed $V_{\mathrm{ex}}=0.3$ eV,
and we use this scaling unless stated otherwise.

The conductance per unit of area, $G/A$, is shown in Fig. 3. It is weakly
dependent on the thickness, $w$, of the magnetic layer, as expected for a
metallic system in the absence of disorder.

Figure 6:  Comparison of critical current density,
$J_{\mathrm{c}}=\max_{\phi}|J(\phi)|$ (blue squares), absolute value of first
Fourier component for supercurrent, $|J_{1}|$ [see Eq. (13), red circles], and
its fitting $|J_{1}^{\mathrm{fit}}|$ [Eq. (14), black curve]. All these
quantities are ‘normalized’ by normal-state conductance, $G$.

We now focus on superconducting properties. The dependence of the supercurrent
on the phase difference $\phi$ for 5, 8, and 11 Ni layers is shown in Figs.
3–3. One can present current-phase relation, $J(\phi)$ in a form of a Fourier
series,

$J(\phi)=\sum\limits_{n\geqslant 1}J_{n}\sin(n\phi).$ (13)

In the $0$-junction mode, the first term in the Fourier series dominates with
$J_{1}>0$ [see solid black line in Fig. 3]. In $\pi$-junction case [Fig. 3],
the supercurrent is also mostly defined by the first harmonic but with
$J_{1}<0$. Close to the $0$-$\pi$ transition $J_{1}$ dies out, so that the
behavior is governed by higher Fourier harmonics [71], e.g. for 8 Ni layers
supercurrent has mostly second harmonic, $J_{2}$, shown by dashed black line
in Fig. 3.

Figure 6 shows the critical current density,
$J_{\mathrm{c}}=\max_{\phi}|J(\phi)|$, normalized by normal-state conductance,
$eJ_{\mathrm{c}}A/G\Delta$, as a function of $w$ (blue squares). Far from the
$0$-$\pi$ transitions, the critical current coincides with the absolute value
of the first Fourier harmonic, $e|J_{1}|A/G\Delta$ (red circles). Since
$J_{1}$ contains a sign of the current and has better numerical stability than
$J_{\mathrm{c}}$, we use this quantity for the analysis. We exclude very thin
junctions (3 Ni layers or less) from the analysis since the magnetic
properties are suppressed there.111For the first data point in Fig. 6
corresponding to 3 layers of Ni the magnetization is completely suppressed
[see Fig. 5] and ratio $eJ_{\mathrm{c}}A/G\Delta\approx 2.2$ is significantly
higher than the one for thicker Ni regions with non zero magnetic moments. We
compare this value with the result for the short disordered SNS junction.
Combination of analytical energy spectrum [64] with Dorokhov distribution of
channel transmissions [72, 73] leads to a ratio of 2.1 (horizontal dashed line
in Fig. 6). We attribute this difference to the fact that there is no
interfacial disorder in our model.

The $J_{1}$ dependence on $w$ can be fit by the following expression,

$J_{1}^{\mathrm{fit}}(w)=\Theta_{J}\exp(-w/\xi_{J})\cos\bigl{[}\pi(w+\delta_{J})/\lambda_{J}\bigr{]},$
(14)

where $\Theta_{J}=3.20\,\mathrm{A}/\mu\mathrm{m}^{2}$, $\xi_{J}=41.1$ Å,
$\lambda_{J}=23.2$ Å, and $\delta_{J}=-2.75$ Å are fitting parameters. We
interpret $\xi_{J}$ as a decay length, $\lambda_{J}$ as the ‘half-period’ of
the oscillation in $J_{1}$ as a function of $w$, and $\delta_{J}$ as a measure
of the suppressed magnetization in Ni layers near the Ni/Nb boundaries.
$J_{1}^{\mathrm{fit}}(w)$ accurately fits the discrete points $J_{1}$ as shown
in Figs. 3, 6, 7, and 10 by black solid line and black circles, accordingly.

Figure 7:  First Fourier harmonic, $J_{1}$, of the supercurrent as a function
of junction thickness, $w$. Full black circles correspond to $J_{1}$
calculated with 290 $\mathbf{k}_{\parallel}$-points, empty circles (mostly
superposed onto the full black circles) correspond to 4142
$\mathbf{k}_{\parallel}$-points. Solid black line is $J_{1}^{\mathrm{fit}}(w)$
[$J_{1}$ fit given by Eq. (14)]. Green semitransparent dashes show
$J_{1}^{\mathrm{fit}}(w)$ contributions [Eq. (15b)] for the individual
$\mathbf{k}_{\parallel}$-points. The vertical ‘errorbars’ correspond to the
standard deviation of $j_{1}(\mathbf{k}_{\parallel})$ with respect to $J_{1}$.
The standard deviation is the same for both sets of
$\mathbf{k}_{\parallel}$-point indicating that these results are independent
of the chosen discretization.  Figure 8:  Colorplot of
$j_{1}(\mathbf{k}_{\parallel})$ with $\mathbf{k}_{\parallel}=(k_{x},k_{y})$.
Each panel corresponds to the local extrema of $J_{1}^{\mathrm{fit}}(w)$ shown
in Fig. 7. For 4 and 13 layers all the $\mathbf{k}_{\parallel}$ channels
contribute with the same sign. For 48 layers and larger, different
$\mathbf{k}_{\parallel}$ channels lose synchronization and contribute to the
total supercurrent with different signs.

In order to gain insight into the evolution of $J_{1}$ with $w$, let us
resolve contributions from different $\mathbf{k}_{\parallel}$. For this, we
rewrite Eq. (3) as

$\displaystyle J(\phi)$
$\displaystyle=A\int\limits_{\mathrm{\scriptscriptstyle
BZ}}\frac{d\mathbf{k}_{\parallel}}{(2\pi)^{2}}\,j(\phi,\mathbf{k}_{\parallel}),$
(15a) $\displaystyle j(\phi,\mathbf{k}_{\parallel})$
$\displaystyle=-\frac{e}{\hbar}\,\frac{1}{A}\,\sum_{\nu>0}\frac{\partial\varepsilon_{\nu}(\phi,\mathbf{k}_{\parallel})}{\partial\phi}.$
(15b)

Here $A=76.75$ Å2 is the area of the surface supercell shown in Fig. 5.
Similar to Eq. (13), we denote the first $\phi$-harmonic of
$j(\phi,\mathbf{k}_{\parallel})$ as $j_{1}(\mathbf{k}_{\parallel})$. The
evolution of the first Fourier harmonic of the supercurrent $J_{1}$ and
$j_{1}(\mathbf{k}_{\parallel})$ as a function of $w$ are shown in Figure 7.
Calculations were performed for two different sets of $\mathbf{k}_{\parallel}$
with 290 discrete $\mathbf{k}_{\parallel}$-points (full black circles) and
4142 $\mathbf{k}_{\parallel}$-points (empty black circles). One can see that
both sets give the same result for $J_{1}$, establishing that the
$\mathbf{k}_{\parallel}$ integration is well converged. In Fig. 7, ‘errorbars’
denote the standard deviation in $j_{1}(\mathbf{k}_{\parallel})$ with respect
to the $\mathbf{k}_{\parallel}$-summed average, $J_{1}(\phi)$. Individual
$j_{1}(\mathbf{k}_{\parallel})$ are shown by semi-transparent horizontal
dashes. The important observation is that while $J_{1}$ decays with $w$, the
dispersion in $j_{1}(\mathbf{k}_{\parallel})$ does not change significantly.

Figure 8 shows colorplots of $j_{1}(\mathbf{k}_{\parallel})$ corresponding to
local extrema of $J_{1}^{\mathrm{fit}}(w)$ (4, 13, 25, 36, 48, and 60 layers
labeled in Fig. 7). For small $w$, one can see that most of all the
$\mathbf{k}_{\parallel}$ contributions to $J_{1}$ have the same sign, i.e.
positive in $0$-junction regime and negative in $\pi$-junction regime. In this
regime the decay is predominantly due to evanescent modes decaying into the
junction. For larger $w$, the dephasing mechanism becomes important since the
phase offset spread grows with $w$. One can observe the apparition of
contributions of the opposite sign for $w\gtrsim 50$ Å. This dephasing
mechanism is mainly due to the variation of the Fermi velocity with
$\mathbf{k}_{\parallel}$, and becomes more important with increasing $w$.

In order to study the distribution of the phase offsets and decay exponents
for different modes, we fit the individual $j_{1}(\mathbf{k}_{\parallel})$
using an expression analogous to Eq. (14). The set of the resulting fitted
curves $j_{1}^{\mathrm{fit}}(w,\mathbf{k}_{\parallel})$ for 4142
$\mathbf{k}_{\parallel}$-points are shown by the green semitransparent curves
in Fig. 3. Here to minimize the numerical ‘noise,’
$j_{1}(w,\mathbf{k}_{\parallel})$ curves are smoothed over the 2D
$\mathbf{k}_{\parallel}$ space using Gaussian filter with
$\sigma_{\mathbf{k}_{\parallel}}=0.01$ Å-1 which is of the order of the Fermi
wave vector in Nb. Thus, each data point $j_{1}(w,\mathbf{k}_{\parallel})$
approximately corresponds to a transverse conducting channel. This fitting
procedure works reasonably well, e.g., the relationship in Eq. (15a) holds if
one replaces $J(\phi)$ by $J_{1}^{\mathrm{fit}}(w)$ and
$j(\phi,\mathbf{k}_{\parallel})$ by
$j_{1}^{\mathrm{fit}}(w,\mathbf{k}_{\parallel})$ for $w\gtrsim 10$ Å.

The distribution of the fitting parameters for
$j_{1}^{\mathrm{fit}}(w,\mathbf{k}_{\parallel})$ is shown in Fig. 9.
Histograms for decay lengths and half-periods reveal a complicated picture
describing different contributions to the supercurrent in real materials.
First of all, in Fig. 9 one can see the distribution of the half-periods,
$\lambda_{j}$, which is similar to a Gaussian distribution with a mean value
$\langle\lambda_{j}\rangle=23.2$ Å and standard deviation $2.8$ Å. The mean
value $\langle\lambda_{j}\rangle$ is very close to $\lambda_{J}$ [see text
after Eq. (14)] while the spread in $\lambda_{j}$ leads to dephasing and is
responsible for the exponential decay of $J_{1}$ at large $w$. Indeed, it is
well-known that the average of an oscillatory function with respect to a
random fluctuating phase (described by a Gaussian distribution) results in an
exponentially decaying function.

In addition to the dephasing mechanism, the decay of the supercurrent
originates from the evanescent modes. The histogram for decay lengths,
$\xi_{j}$ is shown in Fig. 9. Here small $\xi_{j}$ corresponds to fast-
decaying $j_{1}^{\mathrm{fit}}(w,\mathbf{k}_{\parallel})$, large $\xi_{j}$ is
responsible for non-decaying modes (i.e. modes with the decay exponents larger
than the junction thickness). The right-skewed distribution of the decay
exponents has a mean value of $\langle\xi_{j}\rangle=108$ Å which is much
larger than $\xi_{J}$ in Eq. (14). The shoulder at small $\xi_{j}$ presumably
corresponds to the evanescent mode decay comprising of $d$ bands 2 and 3, see
Table 1, whereas the tail at large $\xi_{j}$ originates predominantly from the
band 6. Overall, one can see that a fit with a single decay exponent,
discussed in Eq. (14), is quite oversimplified for a Nb/Ni/Nb junction
considered here.

Figure 9:  Analysis of the fitting parameters of $j_{1}^{\mathrm{fit}}(w)$ for
individual $\mathbf{k}_{\parallel}$, shown in Fig. 3 by semitransparent green
lines. Distribution of 9 half-periods, $\lambda_{j}$, and 9 decay lengths,
$\xi_{j}$.

We now turn to the discussion of the effect of exchange splitting energy on
the supercurrent in MJJs. So far we have used $M/M_{0}=0.5$, which yields the
experimentally observed $V_{\mathrm{ex}}=0.3$ eV. It is interesting to
investigate how a ferromagnet with a different $V_{\mathrm{ex}}$ (but
otherwise the same band structure as Ni) would affect the $w$-dependence of
$J_{1}$. In Fig. 10 we show the results for parametric variations in
$M/M_{0}$. One can see in Fig. 10 that the half-period, $\lambda_{J}$, and the
decay length, $\xi_{J}$, strongly depend on $M$. Here black points correspond
to $M/M_{0}=0.5$, and the self-consistent calculations with no rescaling
correspond to $M/M_{0}=1$. In order to understand how $\lambda_{J}$ and
$\xi_{J}$ depend on $M$, we perform the fitting procedure Eq. (14) for
different magnetic moments and plot, in Fig. 10, the supercurrent density as a
function of the rescaled thickness, $(w+\delta_{J})/2{\tilde{\lambda}}_{J}$
with ${\tilde{\lambda}}_{J}=(M_{0}/M)\,11.6$ Å. Remarkably $J_{1}$, as a
function of the rescaled thickness, collapses to the same universal curve. The
inset demonstrates that the fitted half-period, $\lambda_{J}$, is proportional
to $1/M$ showing that the oscillation period scales linearly with the inverse
of $V_{\mathrm{ex}}$ in this parameter range.222Deviations from linear regime
become significant for $M/M_{0}\gtrsim 1.7$. For the clarity of the data, we
do not show the results for $M/M_{0}>1$ in Fig. 10.

Figure 10:  $J_{1}$ for different rescaling of the magnetic moments $M/M_{0}$
as a function of 10 the thickness, $w$, and 10 the rescaled thickness,
$(w+\delta_{J})/2{\tilde{\lambda}}_{J}$, where $\delta_{J}$ is the fitting
parameter in Eq. (14) and ${\tilde{\lambda}}_{J}=(M_{0}/M)\,11.6$ Å. $J_{1}$
values are shown by circles; corresponding fittings $J_{1}^{\mathrm{fit}}(w)$
are shown by lines of the same color. (Inset) Half-period, $\lambda_{J}$,
fitted using Eq. (14) versus ${\tilde{\lambda}}_{J}$.

We now discuss the difference in crystal orientation in Nb/Ni/Nb junctions. We
have performed calculations for Nb/Ni/Nb junctions built from stacking the Ni
atomic planes in the (110) orientation instead of the (111) orientation. The
results are given in Appendix C. Qualitatively, the same physics hold for both
stacks built from (111) and (110) Ni planes. However, our calculations show
that the actual value for the period of oscillation and for the current decay
depend crucially on the details of the electronic structure of the junctions,
such as the relative crystal orientation.

Finally, we considered effect of spin-orbit coupling in SFS junctions. Spin-
orbit coupling leads to mixing of the minority and majority channels and may
change current-phase relationship. The interplay between Zeeman splitting and
spin-orbit coupling have been discussed in Ref. 45; the regime of interest is
Zeeman-field-dominated regime considered there. Indeed, we find that SOC in Ni
is much smaller than the exchange splitting because of low atomic number of
Ni. As we show in Appendix D, the SOC in Ni-based MJJ considered here does not
change qualitative picture described above but rather leads to small
quantitative changes to the Josephson current.

##  5 Conclusion

In this paper we identified two generic mechanisms for the decay of the
supercurrent with junction thickness: (i) exchange-splitting induced gap
opening for minority or majority carriers and (ii) dephasing between different
modes due to the significant quasiparticle velocity dispersion with the
transverse momentum. It was previously believed that disorder in the
ferromagnet is mainly responsible for the supercurrent decay in SFS junctions.
In the present work we have shown that band structure effects also contribute
to the critical current suppression and thus provide an upper bound for the
supercurrent in ideal (i.e. disorder-free) structure.

We found that the Nb/Ni/Nb junction is a suitable system for comparison with
the simulations because of the long mean free path in Ni relative to the
junction thickness and the quasi-ballistic nature of quasiparticle propagation
in the ferromagnet. We have found good agreement with published experimental
data for the half-period of the critical current oscillations:
$\lambda_{J}\approx 23$ Å [see Eq. (14) and text after it] versus $\approx 26$
Å in experiment, Ref. 22. We have also found that the critical current decays
exponentially with the ferromagnet thickness $w$. This is to be contrasted
with previously assumed algebraic decay based on results for the clean SFS
junctions using simple parabolic-like band structure. We believe that in
measured Nb/Ni/Nb junctions with $w\lesssim 50$ Å the mechanism (i) is likely
to be responsible for the supercurrent decay. This finding is crucial for
material and geometry optimization of MJJs and superconducting magnetic spin
valves.

Understanding the interplay of band structure effects and disorder in MJJs is
an interesting open problem. We believe that interfacial disorder due to, for
example, surface roughness will mix different $\mathbf{k}_{\parallel}$ modes
and will lead to a larger spread of half-periods. This, in turn, will further
enhance the dephasing mechanism (ii) of the supercurrent 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 in the
ferromagnet which is a significant departure from the quasi-ballistic junction
limit considered here. We think that bulk disorder would induce even more
dephasing between different modes because phase offsets in this case will
depend on different random trajectories of minority and majority carriers. We,
therefore, believe that bulk disorder will lead to even stronger decay of the
supercurrent with junction thickness, $w$.

###### Acknowledgements.

HN is much indebted to Dimitar Pashov for stimulating discussions about
developing the Questaal package. The authors express their gratitude to Mason
Thomas for the organizational help and discussions at the early stages of the
project. The authors acknowledge stimulating discussions with Norman Birge,
Anna Herr, Tom Ambrose, Nick Rizzo, and Don Miller. This work is based on
support by the U.S. Department of Energy, Office of Science through the
Quantum Science Center (QSC), a National Quantum Information Science Research
Center. HN and MvS acknowledge financial support from Microsoft Station Q via
a sponsor agreement between KCL and Microsoft Quantum. In the late stages of
this work MvS was supported by the U.S. Department of Energy, Office of
Science, Basic Energy Sciences under Award # FWP ERW7906.

##  Appendix A Mean free path estimate

In this section we provide an estimate for the mean free path
$l_{\mathrm{\scriptscriptstyle MFP}}$ in the Ni ferromagnet. Our approach is
similar to Ref. 44. We use the Kubo formula for conductivity

$\sigma_{xx}=e^{2}\sum_{n,\sigma}\tau_{n\sigma}\,\langle
v^{2}_{n\sigma}\rangle\,{\rho_{n}(E_{\mathrm{\scriptscriptstyle F}})},$

where index $n$ labels the Ni bands, $\sigma$ is the spin projection, and
$\tau_{n\sigma}$ is the scattering time in each band assumed to be momentum
independent. Mean square velocity, $\langle v^{2}_{n\sigma}\rangle$, and
density of states at the Fermi level, $\rho_{n}(E_{\mathrm{\scriptscriptstyle
F}})$, are obtained from our ab initio model calculations and are provided in
Table 1.

We use the available experimental data [74, 75] for thick Ni samples $w\gtrsim
50$ Å: the low temperature linear resistance
$1/(\sigma_{xx}^{\uparrow}+\sigma_{xx}^{\downarrow})\approx
33\,\mathrm{n}\Omega{\cdot}\mathrm{m}$, and bulk spin scattering asymmetry,
$\beta_{\mathrm{\scriptscriptstyle
F}}=(\sigma_{xx}^{\uparrow}-\sigma_{xx}^{\downarrow})/(\sigma_{xx}^{\uparrow}+\sigma_{xx}^{\downarrow})=0.14$.
Only a single $6^{\uparrow}$ band of the majority spin is occupied. Thus, one
can readily estimate mean free path for majority carriers
$l_{\mathrm{\scriptscriptstyle MFP}}^{\uparrow}=\langle
v_{\uparrow}\rangle\tau_{\uparrow}\approx 60$ Å. For the minority channel we
assume that the conductance at large thicknesses is dominated by the most
mobile band, i.e. $6^{\downarrow}$ [76]. This gives an estimate for the mean
free path $l_{\mathrm{\scriptscriptstyle MFP}}^{\downarrow}\approx 61$ Å,
similar to its exchange-split partner $6^{\uparrow}$. These estimates are
consistent with the available ARPES [77] and computational data [44]. As
mentioned in the main text, the mean free path for $6^{\uparrow}$ and
$6^{\downarrow}$ carriers exceeds the typical junction thicknesses measured
experimentally.

##  Appendix B Nb(110)/Ni(111)/Nb(110) junctions

Here we describe how coincident site lattices are constructed for the
Nb(110)/Ni(111) interface. Supercells of both Ni and Nb are separately
constructed in the following way, to make a coincident site lattice.

The primitive unit cells (fcc in the Ni case, with a 0 K lattice constant
3.515 Å, and bcc in the Nb case with a 0 K lattice constant 3.295 Å) are
rotated. For Ni, the rotation is compactly described in terms of three Euler
angles: rotation about $z$ by $\pi/4$, rotation about $x^{\prime}$ by
$\arccos(\sqrt{1/3})$, rotation about $z^{\prime\prime}$ by
$\arccos(4/\sqrt{19})$. Axes $xyz$ are shown in Fig. 1;
$x^{\prime}y^{\prime}z^{\prime}$ and
$x^{\prime\prime}y^{\prime\prime}z^{\prime\prime}$ are axes after first and
second rotation, respectively. From the first two rotations the $[1\bar{1}1]$
axis becomes the new $z^{\prime\prime}$ axis. The last rotation is needed to
make it approximately coincident with Nb. The Nb is rotated about $(1,-1,0)$
by $-\pi/2$, about $z^{\prime}$ by $\pi/4$, and about $z^{\prime}$ again by
$\arccos(5/\sqrt{43})$. After rotations both Ni $(1\bar{1}1)$ and Nb
$(\bar{1}\bar{1}0)$ planes are normal to $z$, which is the propagation
direction.

Next, superlattices must be constructed. They are generated by scaling the
primitive lattice vectors by the following integer multiples, for Ni and Nb
respectively:

$\left[{\begin{array}[]{rrr}3&0&2\\\ -1&0&4\\\
1&1&1\end{array}}\right]\\!,\qquad\left[{\begin{array}[]{rrr}-1&-4&2\\\
2&-2&0\\\ 1&1&2\end{array}}\right]\\!.$

Figure 11:  11 Top view of the Nb(110)/Ni(110) interface supercell. Ni (Nb)
atoms are shown in light blue (grey). A repetition ($3{\times}3$) of the
surface supercell is shown. The corresponding Nb(110) and Ni(110) cells
contain 2 atoms each. 11 Transverse view of the Nb(110)/Ni(110)/Nb(110)
junction with 5 layers of Ni. 11 $J_{1}$ component for the
Nb(110)/Ni(110)/Nb(110) and the Nb(110)/Ni(111)/Nb(110) junctions versus the
thickness $w$. Both currents have a decaying oscillatory behavior, with a half
period of oscillation of $\lambda_{J}\approx 30$ Å for the Ni(110) case and of
$\lambda_{J}\approx 23$ Å for the Ni(111) case.

The Ni supercell contains 14 atoms, all in a single plane perpendicular to
$z$; the Nb supercell consists of 2 planes with 10 atoms per plane [Fig. 5].
The lattice vectors transverse to $z$ are nearly coincident, but not
identically so. To render them coincident, we opt to shear the Ni by the
following linear transformation,

$\left[{\begin{array}[]{cc}0.997&0\\\ 0.050&1.028\end{array}}\right]\\!.$

This shear is a measure of the remaining mismatch of the undistorted Ni and Nb
lattices. It is close enough to unity to have a minor effect on the Ni band
structure.

Along this axis, Nb planes form an ABAB$\dots$ stacking pattern; the Ni are
stacked ABCABC$\dots$ Ni/Nb interfaces are formed by stacking varying numbers
of Ni planes on a Nb substrate, and relaxing the supercells to minimize the
total energy. This was done for a few small superlattices with 3, 4, and 5 Ni
planes, restricting the relaxation to the two Ni and two Nb frontier planes.
In this way we can build up structures of arbitrarily many Ni planes,
sandwiching unrelaxed planes between frontier planes. Thus the entire set of
structures has three families of interfaces: those with integer numbers $3n$,
$3n+1$ or $3n+2$ of Ni planes.

Figure 3 shows that normal conductance, $G$, experiences $\sim 1.3\%$
deviations from its mean $\langle G\rangle$ following the pattern of these
three families. Overall, both normal conductance and supercurrent (see Fig. 7)
are reasonably smooth functions of the thickness $w$, which suggests that the
discreteness of the lattice and the details of lattice relaxation plays a
minor role.

##  Appendix C Nb(110)/Ni(110)/Nb(110) junctions

In this section, we present results for the Nb/Ni/Nb junctions built from
stacking the Ni atomic planes with the (110) orientation instead of the (111)
orientation. Although the (110) plane orientation is not the most stable
configuration for the Nb/Ni interface and most experimental studies focus on
(111) orientation, it is interesting to ascertain whether there is an effect
of crystallographic orientation on the supercurrent. Indeed, the supercurrent
decay and half-period of oscillations with $w$ do depend on crystal
orientation.

We study the Nb(110)/Ni(110)/Nb(110) trilayer shown in Fig. 11. The inter-
plane distance between the Nb and Ni planes at the interfaces is taken to be
the average of the Nb and Ni inter-plane distances. For the sake of
simplicity, we do not perform atomic relaxations for 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 supercell. This has a
slight effect on the Nb band structure.

Calculations of the trilayers were performed self-consistently, and for
transport we reconstruct the potential by rescaling the magnetic moments by
0.5, as in the Ni(111) case. The $0$-$\pi$ and $\pi$-$0$ transitions occur
around $w\approx 11$ Å and 41 Å [shown in red in Fig. 11], yielding a half-
period for the oscillations of $\lambda_{J}\approx 30$ Å. This period is
significantly larger than $\lambda_{J}\approx 23$ Å for the
Nb(110)/Ni(111)/Nb(110) case [shown in black in Fig. 11]. The difference of
$\approx 7$ Å between the two half periods represents roughly 5 inter-plane
distances for Ni(110) [3 inter-plane distances for Ni(111)], and is not
induced by the unrelaxed atomic structures at the Nb(110)/Ni(110) interfaces.
It is noteworthy that the calculated period for the Ni(111) trilayers better
coincides with available experimental data [22] than for the Ni(110) case.

##  Appendix D Effect of spin-orbit interaction

Figure 12:  GW Fermi surfaces of bulk Ni in the $\mathbf{k}_{\parallel}$ plane
corresponding to Ni(111). 12 No SOC [same as Fig. 4]. Shown are majority band
$6^{\downarrow}$ (red) and minority bands $6^{\uparrow}$ (blue),
$3^{\downarrow}$ (dark blue), $4^{\downarrow}$ (cyan), and $5^{\downarrow}$
(green). 12 With SOC. SOC opens gaps around the small regions close to the
band crossings (circled).

In this section, we discuss the effect of spin-orbit coupling (SOC) on
Josephson current in quasi-ballistic SFS junctions. The effect of spin-orbit
interaction is two-fold: (i) SOC changes the band structure and, therefore,
may modify the dependence of the quasiparticle velocities on transverse
momenta and (ii) SOC couples spin and orbital motion and leads to spin
precession along the junction. The combination of spin precession and
scattering on non-magnetic impurities may introduce random spin-flip processes
which suppress the phase difference between minority and majority
quasiparticles. We analyze both mechanisms below and show that they represent
weak perturbations to our main results and do not change qualitative
predictions for quasi-ballistic Ni-based MJJs.

Detailed theoretical analysis of the bulk Ni band structure with SOC and
different orientations of the magnetic moment has been discussed in Ref. 69.
In Fig. 12, we compare Fermi surface of Ni with and without SOC using the same
high-level theory [46] of Fig. 4. The Questaal code can treat SOC effects, see
e.g. tutorial ‘Spin and spin orbit coupling.’ The SOC implementation is
conventional for the energy-dependent local basis set and band structure
calculations (see Sec. 2.8.2 in Ref. 47), and has also been developed for the
Green’s functions (see Supplementary Materials in Ref. 78 and Sec. 2.16 in
Ref. 47).

One may notice that the spin-orbit coupling only weakly modifies the Fermi
surface shown in Fig. 12. Indeed, apart from several pockets (close to the
degeneracy points) across the Brillouin zone, the SOC only weakly perturbs the
band structure. In fcc crystals, the $d$-orbitals are split into degenerate
$t_{\mathrm{2g}}$ and $e_{\mathrm{g}}$ levels, at the $\Gamma$ point,
according to the symmetry. The SOC splits further the degenerate
$t_{\mathrm{2g}}$ levels. From these level shifts (taken at the $\Gamma$
$\mathbf{k}_{\parallel}$-point), we estimate local SOC coupling to be of the
order of $E_{\mathrm{\scriptscriptstyle SO}}\sim 10\,$meV which is much
smaller than the exchange splitting in Ni. In Ref. 69, the corresponding SOC
strength was found to be 68meV which further corroborates our conclusions.

One can now evaluate effect of SOC coupling on quasiparticle propagation
through the junction. Given that bulk band structure of Ni is centro-
symmetric, Dresselhaus spin orbit coupling is forbidden by symmetry. The
Rashba SOC may appear due to the inversion symmetry-breaking along the
direction of the junction. However, the rapid screening of the interface
potential in the ferromagnet would limit Rashba SOC to a few atomic layers in
the junction. Finally, we believe that effect of local spin-flip processes due
to scattering on impurities is weak in the regime of interest. Indeed, here we
are considering quasi-ballistic magnetic Josephson junctions (i.e. the
thickness of the junction smaller than the mean free path
$l_{\mathrm{\scriptscriptstyle MFP}}$ in Ni). Thus, we expect impurity-induced
spin-relaxation rate to be suppressed as well.

##  Appendix E Normal scattering matrix calculation

The Questaal code is based on the LMTO technique which uses a set of electron
wave-function $\phi_{\mathrm{R}L}$ and its energy-derivative
$\dot{\phi}_{\mathrm{R}L}$ as a basis set. The wave-functions
$\phi_{\mathrm{R}L}$ are solutions of the spherical Schrödinger equation in a
sphere around a given atom at position $R$ with angular quantum numbers
$L=l,m$. The LMTO technique is an all-electron approach, with does not rely on
the use of pseudo-potential. The basis set contains core electrons as well as
‘valence’ (non-core) electrons.

The electron ground state of the systems is obtained with the DFT-LDA-ASA
framework [47]. The (energy-dependent) basis set consists of partial waves of
s, p, d character on each of the atomic sites. Core levels are integrated
separately from the valence partial waves to obtain the all-electron charge
density. However, they are not included in the secular matrix.

The calculations are converged when variations in the electron density and
total energy between the last iterations, are below $3\times 10^{-5}$
($10^{-5}$) respectively. To obtain the self-consistent charge density,
calculations were performed with a $\mathbf{k}_{\parallel}$ mesh of ($8\times
4$). As is well known, for purposes of determining the density, a finer
$\mathbf{k}_{\parallel}$ mesh is not needed since the output is a
$\mathbf{k}_{\parallel}$-independent potential for a subsequent transport
calculation.

Atomic relaxations are performed on relatively small cells with a full-
potential method [47], with periodic boundary conditions. Once the ground-
state is reached, the atomic positions of the frontier atoms at the Nb/Ni and
Ni/Nb interfaces are allowed to relax to minimize the total energy.
Convergence is achieved when all the forces on the relaxed atoms are below
$\sim 25$ mRy/bohr. Atomic relaxations are constrained to the first and second
frontier layers; these shifts are then added to the ideal geometry of the
(larger) stacked cells in the subsequent ASA calculations of transport. More
details of the calculations can be found in the Questaal tutorial ‘Nb/Ni
superlattice.’

A denser $\mathbf{k}_{\parallel}$ mesh is needed as the transmission
(reflection) probability is a quantity very sensitive to the
$\mathbf{k}_{\parallel}$-point sampling. As shown in Fig. 3 of the main text,
we used a mesh of ($24\times 24$) and a mesh of ($90\times 92$) points, giving
rise to 290 and 4142 irreducible $\mathbf{k}_{\parallel}$ points respectively.
This establishes that our transport calculations are well converged in
$\mathbf{k}_{\parallel}$.

Note that the principal layer technique (used for the transport) does not rely
on periodic boundary conditions in the (transport) direction perpendicular to
the layers (in our case the $z$-axis). Hence there is not $k_{z}$ point
sampling, the wave number $k_{z}(E)$ should be understood as a continuous
complex function of $E$.

Further details for the constructions of the family of the $3n$, $3n+1$ or
$3n+2$ of Ni planes of the Nb(110)/Ni(111)/Nb(110) junctions can be found in
the Questaal tutorial ‘Nb(110)/Ni/Nb(110) metallic trilayers.’

## References

  * Josephson [1962] B. Josephson, Possible new effects in superconductive tunnelling, Phys. Lett. 1, 251 (1962).
  * Warburton [2011] P. A. Warburton, The josephson effect: 50 years of science and technology, Phys. Educ. 46, 669 (2011).
  * Beenakker [1992] C. W. J. Beenakker, Three “universal” mesoscopic Josephson effects, in  _Transport phenomena in mesoscopic systems_, Vol. 109, edited by H. Fukuyama and T. Ando (Springer, Berlin, Heidelberg, 1992) pp. 235–253.
  * Golubov  _et al._ [2004] A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, The current-phase relation in Josephson junctions, Rev. Mod. Phys. 76, 411 (2004).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Microscopic theory of superconductivity, Phys. Rev. 106, 162 (1957).
  * Dayton  _et al._ [2018] 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, Experimental demonstration of a Josephson magnetic memory cell with a programmable $\pi$-junction, IEEE Magn. Lett. 9, 1 (2018).
  * Buzdin [2005] A. I. Buzdin, Proximity effects in superconductor-ferromagnet heterostructures, Rev. Mod. Phys. 77, 935 (2005).
  * Bergeret  _et al._ [2005] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Odd triplet superconductivity and related phenomena in superconductor-ferromagnet structures, Rev. Mod. Phys. 77, 1321 (2005).
  * Blamire and Robinson [2014] M. G. Blamire and J. W. A. Robinson, The interface between superconductivity and magnetism: understanding and device prospects, J. Condens. Matter Phys. 26, 453201 (2014).
  * Eschrig [2015] M. Eschrig, Spin-polarized supercurrents for spintronics: a review of current progress, Rep. Prog. Phys 78, 104501 (2015).
  * Fulde and Ferrell [1964] P. Fulde and R. A. Ferrell, Superconductivity in a strong spin-exchange field, Phys. Rev. 135, A550 (1964).
  * Larkin and Ovchinnikov [1965] A. I. Larkin and Y. N. Ovchinnikov, Inhomogeneous state of superconductors, Sov. Phys. JETP 20, 762 (1965).
  * Buzdin  _et al._ [1982] A. L. Buzdin, L. N. Bulaevskil, and S. V. Panyukov, Critical-current oscillations as a function of the exchange field and thickness of the ferromagnetic metal (F) in an S-F-S Josephson junction, JETP Lett. 35, 178 (1982).
  * Ryazanov  _et al._ [2001a] V. V. Ryazanov, V. A. Oboznov, A. Y. Rusanov, A. V. Veretennikov, A. A. Golubov, and J. Aarts, Coupling of two superconductors through a ferromagnet: Evidence for a $\pi$ junction, Phys. Rev. Lett. 86, 2427 (2001a).
  * Bell  _et al._ [2004] C. Bell, G. Burnell, C. W. Leung, E. J. Tarte, D.-J. Kang, and M. G. Blamire, Controllable Josephson current through a pseudospin-valve structure, Appl. Phys. Lett 84, 1153 (2004).
  * Gingrich  _et al._ [2016] 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-$\pi$ Josephson junctions containing a ferromagnetic spin valve, Nat. Phys. 12, 564 (2016).
  * Ryazanov  _et al._ [2001b] V. V. Ryazanov, V. A. Oboznov, A. V. Veretennikov, and A. Y. Rusanov, Intrinsically frustrated superconducting array of superconductor-ferromagnet-superconductor $\pi$ junctions, Phys. Rev. B 65, 020501 (2001b).
  * Kontos  _et al._ [2002] T. Kontos, M. Aprili, J. Lesueur, F. Genêt, B. Stephanidis, and R. Boursier, Josephson junction through a thin ferromagnetic layer: Negative coupling, Phys. Rev. Lett. 89, 137007 (2002).
  * Sellier  _et al._ [2003] H. Sellier, C. Baraduc, F. m. c. Lefloch, and R. Calemczuk, Temperature-induced crossover between $0$ and $\pi$ states in S/F/S junctions, Phys. Rev. B 68, 054531 (2003).
  * Robinson  _et al._ [2006] J. W. A. Robinson, S. Piano, G. Burnell, C. Bell, and M. G. Blamire, Critical current oscillations in strong ferromagnetic $\pi$ junctions, Phys. Rev. Lett. 97, 177003 (2006).
  * Khaire  _et al._ [2009] T. S. Khaire, W. P. Pratt, and N. O. Birge, Critical current behavior in Josephson junctions with the weak ferromagnet PdNi, Phys. Rev. B 79, 094523 (2009).
  * Baek  _et al._ [2017] B. Baek, M. L. Schneider, M. R. Pufall, and W. H. Rippard, Phase offsets in the critical-current oscillations of Josephson junctions based on Ni and Ni-(Ni81Fe19)xNby barriers, Phys. Rev. Applied 7, 064013 (2017).
  * Baek  _et al._ [2018] B. Baek, M. L. Schneider, M. R. Pufall, and W. H. Rippard, Anomalous supercurrent modulation in Josephson junctions with Ni-based barriers, IEEE Trans. Appl. Supercond 28, 1 (2018).
  * Aguilar  _et al._ [2020] V. Aguilar, D. Korucu, J. A. Glick, R. Loloee, W. P. Pratt, and N. O. Birge, Spin-polarized triplet supercurrent in Josephson junctions with perpendicular ferromagnetic layers, Phys. Rev. B 102, 024518 (2020).
  * Mishra  _et al._ [2021] S. S. Mishra, R. Loloee, and N. O. Birge, Supercurrent transmission through Ni/Ru/Ni synthetic antiferromagnets, Appl. Phys. Lett. 119, 172603 (2021).
  * McMillan [1968] W. L. McMillan, Theory of superconductor—normal-metal interfaces, Phys. Rev. 175, 559 (1968).
  * Wolfram [1968] T. Wolfram, Tomasch oscillations in the density of states of superconducting films, Phys. Rev. 170, 481 (1968).
  * Kulik [1970] I. O. Kulik, Macroscopic quantization and the proximity effect in S-N-S junctions, J. Exp. Theor. Phys. 30, 944 (1970).
  * Demers and Griffin [1971] J. Demers and A. Griffin, Scattering and tunneling of electronic excitations in the intermediate state of superconductors, Can. J. Phys. 49, 285 (1971).
  * Griffin and Demers [1971] A. Griffin and J. Demers, Tunneling in the normal-metal-insulator-superconductor geometry using the Bogoliubov equations of motion, Phys. Rev. B 4, 2202 (1971).
  * Entin-Wohlman [1977] O. Entin-Wohlman, Effect of a barrier at the superconducting-normal metal interface, J. Low Temp. Phys 27, 777 (1977).
  * Blonder  _et al._ [1982] G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion, Phys. Rev. B 25, 4515 (1982).
  * Furusaki and Tsukada [1991] A. Furusaki and M. Tsukada, Dc Josephson effect and Andreev reflection, Solid State Commun 78, 299 (1991).
  * Furusaki  _et al._ [1992] A. Furusaki, H. Takayanagi, and M. Tsukada, Josephson effect of the superconducting quantum point contact, Phys. Rev. B 45, 10563 (1992).
  * Furusaki [1994] A. Furusaki, DC Josephson effect in dirty SNS junctions: Numerical study, Physica B 203, 214 (1994).
  * de Jong and Beenakker [1995] M. J. M. de Jong and C. W. J. Beenakker, Andreev reflection in ferromagnet-superconductor junctions, Phys. Rev. Lett. 74, 1657 (1995).
  * Tanaka and Kashiwaya [1997] Y. Tanaka and S. Kashiwaya, Theory of Josephson effect in superconductor-ferromagnetic-insulator-superconductor junction, Physica C 274, 357 (1997).
  * Žutić and Valls [1999] I. Žutić and O. T. Valls, Spin-polarized tunneling in ferromagnet/unconventional superconductor junctions, Phys. Rev. B 60, 6320 (1999).
  * Radović  _et al._ [2003] Z. Radović, N. Lazarides, and N. Flytzanis, Josephson effect in double-barrier superconductor-ferromagnet junctions, Phys. Rev. B 68, 014501 (2003).
  * Cayssol and Montambaux [2005] J. Cayssol and G. Montambaux, Incomplete andreev reflection in a clean superconductor-ferromagnet-superconductor junction, Phys. Rev. B 71, 012507 (2005).
  * Konschelle  _et al._ [2008] F. Konschelle, J. Cayssol, and A. I. Buzdin, Nonsinusoidal current-phase relation in strongly ferromagnetic and moderately disordered sfs junctions, Phys. Rev. B 78, 134505 (2008).
  * Tzortzakakis and Flytzanis [2019] A. F. Tzortzakakis and N. Flytzanis,  _Josephson junctions with spin-orbit and spin-flip interactions_ , Master’s thesis, University of Crete (2019).
  * Demler  _et al._ [1997] E. A. Demler, G. B. Arnold, and M. R. Beasley, Superconducting proximity effects in magnetic metals, Phys. Rev. B 55, 15174 (1997).
  * Gall [2016] D. Gall, Electron mean free path in elemental metals, J. Appl. Phys 119, 085101 (2016).
  * Cheng and Lutchyn [2012] M. Cheng and R. M. Lutchyn, Josephson current through a superconductor/semiconductor-nanowire/superconductor junction: Effects of strong spin-orbit coupling and Zeeman splitting, Phys. Rev. B 86, 134522 (2012).
  * Sponza  _et al._ [2017] 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 95, 041112 (2017).
  * Pashov  _et al._ [2020] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Questaal: A package of electronic structure methods based on the linear muffin-tin orbital technique, Comput. Phys. Commun 249, 107065 (2020).
  * Faleev  _et al._ [2005] S. V. Faleev, F. Léonard, D. A. Stewart, and M. van Schilfgaarde, Ab initio tight-binding LMTO method for nonequilibrium electron transport in nanosystems, Phys. Rev. B 71, 195422 (2005).
  * Meir and Wingreen [1992] Y. Meir and N. S. Wingreen, Landauer formula for the current through an interacting electron region, Phys. Rev. Lett. 68, 2512 (1992).
  * Andersen and Jepsen [1984] O. K. Andersen and O. Jepsen, Explicit, first-principles tight-binding theory, Phys. Rev. Lett. 53, 2571 (1984).
  * Fisher and Lee [1981] D. S. Fisher and P. A. Lee, Relation between conductivity and transmission matrix, Phys. Rev. B 23, 6851 (1981).
  * Chen  _et al._ [1989] A.-B. Chen, Y.-M. Lai-Hsu, and W. Chen, Difference-equation approach to the electronic structures of surfaces, interfaces, and superlattices, Phys. Rev. B 39, 923 (1989).
  * Fujimoto and Hirose [2003] Y. Fujimoto and K. Hirose, First-principles treatments of electron transport properties for nanoscale junctions, Phys. Rev. B 67, 195315 (2003).
  * Wimmer [2008] M. Wimmer,  _Quantum transport in nanostructures: From computational concepts to spintronics in graphene and magnetic tunnel junctions_ , Ph.D. thesis, Universität Regensburg (2008).
  * Halterman and Valls [2001] K. Halterman and O. T. Valls, Proximity effects at ferromagnet-superconductor interfaces, Phys. Rev. B 65, 014509 (2001).
  * Halterman and Valls [2002] K. Halterman and O. T. Valls, Proximity effects and characteristic lengths in ferromagnet-superconductor structures, Phys. Rev. B 66, 224516 (2002).
  * Csire  _et al._ [2018] G. Csire, A. Deák, B. Nyári, H. Ebert, J. F. Annett, and B. Újfalussy, Relativistic spin-polarized KKR theory for superconducting heterostructures: Oscillating order parameter in the Au layer of Nb/Au/Fe trilayers, Phys. Rev. B 97, 024514 (2018).
  * Halterman and Valls [2004] K. Halterman and O. T. Valls, Layered ferromagnet-superconductor structures: The $\pi$ state and proximity effects, Phys. Rev. B 69, 014517 (2004).
  * Halterman  _et al._ [2007] K. Halterman, P. H. Barsic, and O. T. Valls, Odd triplet pairing in clean superconductor/ferromagnet heterostructures, Phys. Rev. Lett. 99, 127002 (2007).
  * Halterman  _et al._ [2015] K. Halterman, O. T. Valls, and C.-T. Wu, Charge and spin currents in ferromagnetic Josephson junctions, Phys. Rev. B 92, 174516 (2015).
  * Halterman and Alidoust [2016] K. Halterman and M. Alidoust, Josephson currents and spin-transfer torques in ballistic SFSFS nanojunctions, Supercond. Sci. Technol. 29, 055007 (2016).
  * Alidoust and Halterman [2020] M. Alidoust and K. Halterman, Supergap and subgap enhanced currents in asymmetric ${\mathrm{s}}_{1}{\mathrm{fs}}_{2}$ Josephson junctions, Phys. Rev. B 102, 224504 (2020).
  * Yagovtsev  _et al._ [2021] V. O. Yagovtsev, N. G. Pugach, and M. Eschrig, The inverse proximity effect in strong ferromagnet–superconductor structures, Supercond. Sci. Technol. 34, 025003 (2021).
  * Beenakker [1991] C. W. J. Beenakker, Universal limit of critical-current fluctuations in mesoscopic josephson junctions, Phys. Rev. Lett. 67, 3836 (1991).
  * van Heck  _et al._ [2014] B. van Heck, S. Mi, and A. R. Akhmerov, Single fermion manipulation via superconducting phase differences in multiterminal Josephson junctions, Phys. Rev. B 90, 155450 (2014).
  * Aguayo  _et al._ [2004] A. Aguayo, I. I. Mazin, and D. J. Singh, Why Ni3Al is an itinerant ferromagnet but Ni3Ga is not, Phys. Rev. Lett. 92, 147201 (2004).
  * Moriya [1985] T. Moriya,  _Spin fluctuations in itinerant electron magnetism_ (Springer-Verlag, Berlin, 1985).
  * Himpsel  _et al._ [1979] F. J. Himpsel, J. A. Knapp, and D. E. Eastman, Experimental energy-band dispersions and exchange splitting for ni, Phys. Rev. B 19, 2919 (1979).
  * Bünemann  _et al._ [2008] J. Bünemann, F. Gebhard, T. Ohm, S. Weiser, and W. Weber, Spin-orbit coupling in ferromagnetic nickel, Phys. Rev. Lett. 101, 236404 (2008).
  * Karlsson and Aryasetiawan [2000] K. Karlsson and F. Aryasetiawan, A many-body approach to spin-wave excitations in itinerant magnetic systems, J. Phys. Condens. Matter 12, 7617 (2000).
  * Stoutimore  _et al._ [2018] M. J. A. Stoutimore, A. N. Rossolenko, V. V. Bolginov, V. A. Oboznov, A. Y. Rusanov, D. S. Baranov, N. Pugach, S. M. Frolov, V. V. Ryazanov, and D. J. Van Harlingen, Second-harmonic current-phase relation in Josephson junctions with ferromagnetic barriers, Phys. Rev. Lett. 121, 177702 (2018).
  * Dorokhov [1984] O. N. Dorokhov, On the coexistence of localized and extended electronic states in the metallic phase, Solid State Commun 51, 381 (1984).
  * Mello  _et al._ [1988] P. A. Mello, P. Pereyra, and N. Kumar, Macroscopic approach to multichannel disordered conductors, Ann Phys 181, 290 (1988).
  * Moreau  _et al._ [2007] C. E. Moreau, I. C. Moraru, N. O. Birge, and W. P. Pratt, Measurement of spin diffusion length in sputtered Ni films using a special exchange-biased spin valve geometry, Appl. Phys. Lett 90, 012101 (2007).
  * Bass [2011] J. Bass, Giant magnetoresistance: Experiment, in  _Handbook of spin transport and magnetism_, edited by I. Zutic and E. Y. Tsymbal (CRC Press, 2011) pp. 69–94.
  * Tsymbal  _et al._ [2011] E. Y. Tsymbal, D. Pettifor, and S. Maekawa, Giant magnetoresistance: Theory, in  _Handbook of spin transport and magnetism_, edited by I. Zutic and E. Y. Tsymbal (CRC Press, 2011) pp. 95–114.
  * Petrovykh  _et al._ [1998] D. Y. Petrovykh, K. N. Altmann, H. Höchst, M. Laubscher, S. Maat, G. J. Mankey, and F. J. Himpsel, Spin-dependent band structure, Fermi surface, and carrier lifetime of permalloy, Appl. Phys. Lett 73, 3459 (1998).
  * Belashchenko  _et al._ [2015] K. D. Belashchenko, L. Ke, M. Däne, L. X. Benedict, T. N. Lamichhane, V. Taufour, A. Jesche, S. L. Bud’ko, P. C. Canfield, and V. P. Antropov, Origin of the spin reorientation transitions in (Fe1-xCox)2B alloys, Appl. Phys. Lett. 106, 062408 (2015).