# Supergap and subgap enhanced currents in asymmetric $\rm S_{1}FS_{2}$
Josephson junctions

Mohammad Alidoust  Department of Physics, Norwegian University of Science and
Technology, N-7491 Trondheim, Norway  Klaus Halterman  Michelson Lab, Physics
Division, Naval Air Warfare Center, China Lake, California 93555, USA

###### Abstract

We have theoretically studied the supercurrent profiles in three-dimensional
normal metal and ferromagnetic Josephson configurations, where the magnitude
of the superconducting gaps in the superconducting leads are unequal, i.e.,
$\Delta_{1}\neq\Delta_{2}$, creating asymmetric $\rm S_{1}NS_{2}$ and $\rm
S_{1}FS_{2}$ systems. Our results reveal that by increasing the ratio of the
superconducting gaps $\Delta_{2}/\Delta_{1}$, the critical supercurrent in a
ballistic $\rm S_{1}NS_{2}$ system can be enhanced by more than $100\%$, and
reaches a saturation point, or decays away, depending on the junction
thickness, magnetization strength, and chemical potential. The total critical
current in a diffusive $\rm S_{1}NS_{2}$ system was found to be enhanced by
more than $50\%$ parabolically and reaches saturation by increasing one of the
superconducting gaps. In a uniform ferromagnetic junction, the supercurrent
undergoes reversal by increasing $\Delta_{2}/\Delta_{1}$ $>1$. Through
decomposing the total supercurrent into its supergap and subgap components,
our results illustrate their crucial relative contributions to the Josephson
current flow. It was found that the competition of subgap and supergap
currents in a $\rm S_{1}FS_{2}$ junction results in the emergence of second
harmonics in the current-phase relation. In contrast to a diffusive asymmetric
Josephson configuration, the behavior of the supercurrent in a ballistic
system with $\Delta_{2}/\Delta_{1}$ $=1$ can be properly described by the
subgap current component only, in a wide range of parameter sets, including
Fermi level mismatch, magnetization strength, and junction thickness.
Interestingly, when $\Delta_{2}/\Delta_{1}$ $>1$, our results have found
multiple parameter sets where the total supercurrent is driven by the supergap
component. Therefore, our comprehensive study highlights the importance of
subgap and supergap supercurrent components in both the ballistic and
diffusive regimes. We focus on experimentally accessible material and
geometric parameters that can lead to advancements in cryogenic devices based
on Josephson junction architectures that utilize supergap currents, which are
less sensitive to temperature compared to the subgap current.

##  I introduction

When two superconductors with different macroscopic phases are weakly coupled
by proximity effects, a finite dissipationless current can flow, demonstrating
the Josephson effect B.D.Josephson . The current flow is carried through the
coherent tunneling of Cooper pairs from one superconductor (S) to the other.
The coherent nature of Cooper pairs allows for supercurrent flow through
finite-thickness normal (N) metal and ferromagnetic (F) materials (SNS and SFS
junctions). For these types of junctions, the physical quantities of interest
can have complicated variations across the structure over a wide range of
length scales due to proximity induced inhomogeneous superconductivity.

The widely accepted microscopic theory of conventional superconductors is the
mean field BCS theory, where two electrons with opposite momenta and spins
create a single boson through lattice vibrations. This theory was later
reformulated by introducing particle-hole space, which is the well-known
Bogoliubov-de Gennes (BdG) approach bdg . To microscopically study systems
containing a superconducting segment, one employs the associated BdG
Hamiltonian with pair potential $\Delta(x)$ to account for the spatially
varying superconducting correlations. The competition between superconducting
order and other phases in proximity coupled junctions can induce striking
phenomena that has attracted considerable attention over the decades Z.Radovic
; T.Karabassov ; S.Acharjee ; Z.Shomali2011 ; Hikino ; Setiawan ; jap_al ;
s1fs2_ov ; Ryazanov1 ; Ryazanov2 ; Fominov1 ; Takahashi ; H.Chakraborti ;
T.Karabassov2019 ; S.V.Bakurskiy ; Z.Shomali ; M.V.Avdeev ; E.Moen2018 ;
K.Ohnishi ; K.Kulikov ; layered ; kontos2002 ; shell2006 ; zutic ; fabian ;
Mazanik ; H.Meng ; khold . Differing approaches and a wide range of
approximations have been incorporated to study various normal and
ferromagnetic superconducting hybrids that have achieved success to describe
experimental observationsD.Culcer ; G.Tkachov ; A.V.Galaktionov ; A.G.Golubov
; Iovan1 ; Iovan2 ; L.R.Tagirov ; M.Alidoust2020 ; K.Halterman2015 ;
C.-T.Wu2018 ; E.Koshina ; C.W.J.Beenakker ; A.A.Golubov ; halterman2002 . For
instance, a recent study of superconducting (half-)metallic spin-valves has
shown good agreement between theoretical predictions and experimental
observationszep ; half ; bernard1 ; half2 .

Nevertheless, except in simple situations, it is highly challenging to obtain
analytical solutions to the BdG Hamiltonian. One particular example is
asymmetric Josephson junctions, where the pair potentials in the S regions are
different, i.e., $\Delta_{1}\neq\Delta_{2}$. There are mainly two approaches
for studying current flow in $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$
configurations. (i) The wave-function approach where one diagonalizes the BdG
Hamiltonian to obtain the wave functions and energies, which after application
of the appropriate boundary conditions, permits calculation of the subgap
bound states. To further simplify the resultant expressions, the vast majority
of works utilize the so-called Andreev approximation. (ii) The other approach
is Gorkov’s Green function technique abrikosov . Here, also one needs to
incorporate multiple simplifying assumptions for obtaining simple and solvable
equations. The best-known approximation in this approach is the quasiclassical
approximation where the Fermi energy is considered the largest energy in the
system, leading to the Eilenberger equation eilenberger . One main advantage
of this approach is that is can conveniently accommodate nonmagnetic
impurities via a white-noise scattering potential. In the presence of disorder
and nonmagnetic impurities one can integrate the Eilenberger equation over the
random quasiparticle scattering angle to arrive at the Usadel equation usadel
. This approach has been recently generalized to a spin-orbit-coupled electron
gas to study several phenomena including: the spatial distribution of spin
currents alidoust1 ; alidoust2 , the surface state of three-dimensional
topological insulatorszu1 ; zu3 , Weyl semimetalsAlidoustWS2 , and black
phosphorusAlidoustBP2 .

The former approach (i) has been used to simulate ballistic systems, where
multiple interference effects from the propagating quasiparticles strongly
influences the transport behavior of the system. It was demonstrated that this
approach, if followed analytically, can be problematic for asymmetric $\rm
S_{1}NS_{2}$ junctions even within the quasiclassical regimebagwell . The
problem becomes increasingly difficult in asymmetric $\rm S_{1}FS_{2}$
structures due to the inclusion of band spin splitting. One main issue is to
properly obtain the contribution of supergap channels to the total
supercurrent. These modes become particularly important in asymmetric
junctions bagwell ; A.V.Galaktionov ; zu3  due to the imbalance of
superconducting gaps that open up the continuum domain to states that can
carry considerable amounts of supercurrent.

In this paper, we aim to study the behavior of the supercurrent in asymmetric
three dimensional $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ Josephson junctions
in both the ballistic and diffusive regimes. Due to the asymmetry in the pair
potential “well”, three relevant energy scales play a role in the net
supercurrent response: (i) subgap energies ($\varepsilon\leq\Delta_{1}$),
which comprise the resonant Andreev bound states, (ii) supergap energies
($\Delta_{1}<\varepsilon\leq\Delta_{2}$), and (iii) energies in the continuum,
where $\varepsilon>\Delta_{2}$. We demonstrate that our microscopic numerical
approaches in the ballistic regime can adequately describe the supercurrent
flow in all three energy regimes and provides an accessible framework that
recovers previous results in various asymptotic limits for a simpler one-
dimensional quasiclassical $\rm S_{1}NS_{2}$ system bagwell . Our numerical
approaches allow for exploring realms beyond those studied in the vast
majority of the literature without imposing any limitations to Fermi level
mismatch and magnetization strength (supporting weak magnetization to a half-
metallic phase). Our results reveal that when $\Delta_{2}/\Delta_{1}$ $=1$,
the subgap component of the supercurrent that describes the resonant bound
states can properly account for the total supercurrent in a ballistic
Josephson junction, regardless of Fermi level mismatch, junction thickness,
and magnetization strength. When $\Delta_{2}/\Delta_{1}$ $>1$, we find the the
critical current can be strongly enhanced in highly asymmetric ballistic $\rm
S_{1}NS_{2}$ junctions. By means of the current-phase relations, we find that
the supergap and subgap supercurrent components can propagate in opposite
directions, and within certain regimes, the subgap supercurrent vanishes, so
that the total supercurrent arises from supergap states.

In the full proximity limit of the diffusive regime, we find that the critical
supercurrent can be enhanced by more than $50\%$ when increasing the
superconducting gap ratio to $\Delta_{2}/\Delta_{1}$ $\sim 10$, in an
asymmetric $\rm S_{1}NS_{2}$ junction. The critical supercurrent also shows an
oscillatory behavior in $\rm S_{1}FS_{2}$ junctions when increasing the
exchange field intensity, and additionally, the supercurrent undergoes
reversals as a function of $\Delta_{2}/\Delta_{1}$ for certain magnetization
strengths. By calculating the total current-phase-relation with its subgap and
supergap current components, our results reveal that the emergence of a $\sin
2\Delta\varphi$ harmonic close to a current reversal point is the consequence
of the intricate competition between the subgap and supergap currents flowing
in opposite directions.

The paper is organized as follows. In Sec. II, we have summarized the main
equations which establish the theoretical framework employed throughout the
calculations. In Secs. II.1 and II.2, we present detailed formulations of the
ballistic and diffusive regimes, respectively. In Sec. III, the main results
and findings are presented. In Secs. III.1 and III.2, we discuss the results
for the ballistic and diffusive regimes, respectively. Lastly, in Sec. IV, we
give concluding remarks.

##  II theory and model

Figure 1: (Color online). Schematic of the asymmetric Josephson junction. The
left and right superconductors have different superconducting gaps,
$\Delta_{1,2}$, and macroscopic phases, $\varphi_{l,r}$, respectively. The
superconductor-nonsuperconductor interfaces are located at $x=\pm d/2$. The
two superconductors are connected either by a normal metal (N) or ferromagnet
(F), making $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ Josephson configurations.
The subgap and supergap currents are marked by $J_{1}$ and $J_{2,3}$,
respectively, depending on the quasiparticle energy $\varepsilon$. The
superconducting phase difference is defined by
$\Delta\varphi=|\varphi_{l}-\varphi_{r}|$.

In the ballistic regime, we directly solve the microscopic Bogoliubov-de
Gennes (BdG) equations. bdg  In Appendix A we also outline a complementary
numerical method that can be employed to contrast and compare results. By
employing two distinct numerical approaches halterman2002 ; K.Halterman2015 ;
C.-T.Wu2018 ; K.Halterman_ss2016 ; AlidoustBP1 ; M.Alidoust2020 . it allows
for a comprehensive and accurate investigation into general hybrid Josephson
junctions without being limited to a narrow range of ferromagnetic exchange
fields and Fermi level differences. Moreover, both approaches produce similar
results, as expected. In other words, with these methods in the ballistic
regime, one is able to span a wide parameter space from weak magnetizations to
half-metallic systems, half ; half2 ; zep ; C.-T.Wu2018  over a wide range of
Fermi level mismatches between the S electrodes and junction insert. Of
course, limiting cases such as the quasiclassical regime can be studied as
well.AlidoustWS1 ; AlidoustBP1 ; M.Alidoust2020 ; AlidoustBP2 ; AlidoustWS2
For systems containing impurities and disorder, we make use of the Usadel
equation usadel  in the full proximity limit of the quasiclassical regime to
study the diffusive motion of quasiparticles in asymmetric $\rm S_{1}NS_{2}$
and $\rm S_{1}FS_{2}$ systems.

###  II.1 Ballistic regime

The effective Hamiltonian that describes our asymmetric Josephson junction is
given by:

$\displaystyle{\cal H}_{\rm eff}$ $\displaystyle=$ $\displaystyle\int
d^{3}r\left\\{\sum_{s}\psi_{s}^{\dagger}\left(\bm{r}\right){\cal
H}_{0}\psi_{s}\left(\bm{r}\right)\right.$ (1) $\displaystyle+$
$\displaystyle\left.\frac{1}{2}\left[\sum_{s\>s^{\prime}}\left(i\sigma_{y}\right)_{ss^{\prime}}\Delta\left(\bm{r}\right)\psi_{s}^{\dagger}\left(\bm{r}\right)\psi_{s^{\prime}}^{\dagger}\left(\bm{r}\right)+H.c.\right]\right.$
$\displaystyle-$
$\displaystyle\left.\sum_{s\>s^{\prime}}\psi_{s}^{\dagger}\left(\bm{r}\right)\left(\bm{h}\cdot\bm{\sigma}\right)_{ss^{\prime}}\psi_{s^{\prime}}\left(\bm{r}\right)\right\\},$

where $s$,$s^{\prime}$ are spin indices, $\bm{\sigma}$ are Pauli matrices, and
the exchange energy $\bm{h}$ describes the ferromagnet exchange interaction.
The kinetic part of the single particle Hamiltonian is defined as,

${\cal H}_{0}(x)=-\frac{1}{2m}\frac{\partial^{2}}{\partial
x^{2}}+\varepsilon_{\perp}-E_{F}(x),$ (2)

in which $\varepsilon_{\perp}=\frac{1}{2m}(k_{y}^{2}+k_{z}^{2})$ is the
quasiparticle energy in the $yz$ plane of Fig. 1. To accommodate the
possibility of differing bandwidths in the two junction materials, we take the
Fermi level $E_{F}(x)$ to equal $E_{FM}$ in the ferromagnet region and
$E_{FS}$ in the superconductor region. Following standard proceduresbdg , we
then utilize the generalized Bogoliubov transformation,
$\psi_{s}=\sum_{n}\left(u_{ns}\gamma_{n}+\eta_{s}v_{ns}^{\ast}\gamma_{n}^{\dagger}\right)$,
where $\eta_{s}\equiv 1(-1)$ for spin-down (up), to cast Eq. (1) in terms of
the spin-generalized BdG equationsbdg :

$\displaystyle\begin{pmatrix}{\cal H}_{0}-h_{z}&-h_{x}+ih_{y}&0&\Delta\\\
-h_{x}-ih_{y}&{\cal H}_{0}+h_{z}&\Delta&0\\\ 0&\Delta^{*}&-({\cal
H}_{0}-h_{z})&-h_{x}-ih_{y}\\\ \Delta^{*}&0&-h_{x}+ih_{y}&-({\cal
H}_{0}+h_{z})\\\ \end{pmatrix}\begin{pmatrix}u_{n\uparrow}(x)\\\
u_{n\downarrow}(x)\\\ v_{n\uparrow}(x)\\\
v_{n\downarrow}(x)\end{pmatrix}=\epsilon_{n}\begin{pmatrix}u_{n\uparrow}(x)\\\
u_{n\downarrow}(x)\\\ v_{n\uparrow}(x)\\\ v_{n\downarrow}(x)\end{pmatrix},$
(3)

where $u_{ns}$ and $v_{ns}$ are the quasiparticle and quasihole amplitudes,
respectively (for $s=\uparrow,\downarrow$). The generalized BdG technique has
been shown to provide a numerically stable framework for solving inhomogeneous
superconductivity problems K.Halterman_ss2016 ; K.Halterman2015 ; half ; half2
. For the layered Josephson junctionslayered  considered in this work, we
assume each F and S layer is infinite in the $yz$ plane and the finite layer
thicknesses extend along the $x$ axis (see Fig. 1). As a result, the BdG
equations are translationally invariant in the $yz$ plane, and become quasi-
one-dimensional in $x$.

To solve the BdG equation, we first expand K.Halterman_ss2016 ;
K.Halterman2015  the quasiparticle amplitudes in a Fourier series using a
complete set of $N$ basis functions:

$\displaystyle\psi_{n}(x)=\sqrt{\frac{2}{d}}\sum_{q=0}^{N}\sin({k_{q}x})\hat{\psi}_{q}(k_{q}),$
(4)

where
$\psi_{n}(x)=(u_{n\uparrow}(x),u_{n\downarrow}(x),v_{n\uparrow}(x),v_{n\downarrow}(x))$,
and
$\hat{\psi}_{q}=(\hat{u}_{q\uparrow},\hat{u}_{q\downarrow},\hat{v}_{q\uparrow},\hat{v}_{q\downarrow})$.
The wavevector $k_{q}=q\pi/d$ is discretized by the thickness of junction $d$.
Next, we transform the real-space BdG equations by first inserting Eq. (4)
into Eq. (3) and using orthogonality of the basis set to give:

$\displaystyle\hat{H}_{0}(q,q^{\prime})$
$\displaystyle=\frac{2}{d}\int_{0}^{d}dx\Biggl{\\{}\left(\frac{k_{q}^{2}}{2m}+\varepsilon_{\perp}-E_{F}(x)\right)$
$\displaystyle\times\sin(k_{q}x)\sin(k_{q^{\prime}}x)\Biggr{\\}},$ (5a)
$\displaystyle\hat{\Delta}({q,q^{\prime}})$
$\displaystyle=\frac{2}{d}\int_{0}^{d}dx\Delta(x)\sin(k_{q}x)\sin(k_{q^{\prime}}x),$
(5b) and $\displaystyle\hat{h}_{i}({q,q^{\prime}})$
$\displaystyle=\frac{2}{d}\int_{0}^{d}dx\,h_{i}(x)\sin(k_{q}x)\sin(k_{q^{\prime}}x).$
(5c)

Here $i=x,y,z$ and we have defined
$\hat{u}_{\sigma}=(\hat{u}_{1\sigma},\hat{u}_{2\sigma},\ldots,\hat{u}_{N\sigma})$,
$\hat{v}_{\sigma}=(\hat{v}_{1\sigma},\hat{v}_{2\sigma},\ldots,\hat{v}_{N\sigma})$.
Additional details on this solution process can be found elsewhere.khold

To compute the dc Josephson current, we numerically diagonalize the Fourier
transformed BdG equations [with matrix elements in Eqs. (5)-(5c)] to get the
eigenenergies $\epsilon_{n}$ and quasiparticle coefficients $u_{i\sigma}$,
$v_{i\sigma}$ ($i=1,...,N$). The real-space amplitudes are then obtained via
the series expansion in Eq. (4). Since we wish to determine the current-phase
relation for asymmetric Josephson junctions, the input for the pair potential
is taken to be the bulk gap, $\Delta_{1}\exp(i\varphi_{l})$, in ${S_{1}}$ and
$\Delta_{2}\exp(i\varphi_{r})$ in $S_{2}$. With this form for $\Delta(x)$, and
making use of the obtained wave functions and eigenenergies, we calculate the
charge current with the expression,

$\displaystyle J_{x}=\frac{2e}{m}\sum_{ns}{\rm Im}\left[u_{ns}\frac{\partial
u^{*}_{ns}}{\partial x}f_{n}+v_{ns}\frac{\partial v^{*}_{ns}}{\partial
x}\left(1-f_{n}\right)\right],$ (6)

where $f_{n}$ is the Fermi function. The supercurrent satisfies the
conservation law

$\displaystyle\frac{\partial J_{x}(x)}{\partial x}\\!=\\!2e\,{\rm
Im}\left\\{\Delta({x})\sum_{n}\left[u_{n\uparrow}^{*}v_{n\downarrow}+u_{n\downarrow}^{*}v_{n\uparrow}\right]\tanh\left(\frac{\epsilon_{n}}{2T}\right)\right\\}.$
(7)

Thus, within the junction region where $\Delta(x)$ vanishes, the current
density is uniform. We refer the reader to Refs. K.Halterman_ss2016, ;
K.Halterman2015,  for further details on the methods used here for calculating
the supercurrent. Also, an alternative approach to study ballistic asymmetric
$\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ systems is described in Appendix A.

###  II.2 Diffusive regime

In a system containing nonmagnetic impurities and disorder, the motion of
quasiparticles can be described by a diffusion equation because of the
scattering sources. In the quasiclassical regime, where the Fermi energy is
the largest energy scale in the system, the diffusion equation is given by the
so-called Usadel equation usadel ,

$\displaystyle
D\hat{{\bm{\nabla}}}\\{\hat{G}(\varepsilon,\mathbf{r})\hat{{\bm{\nabla}}}\hat{G}(\varepsilon,\mathbf{r})\\}+i[\varepsilon\hat{\rho}_{z},\hat{G}(\varepsilon,\mathbf{r})]=0,$
(8a)
$\displaystyle\hat{G}(\varepsilon,\mathbf{r})=\left(\begin{array}[]{cc}G^{A}&G^{K}\\\
0&G^{R}\end{array}\right),$ (8d)

where $D$ is the diffusion constant,
$\hat{{\bm{\nabla}}}\equiv(\partial_{x},\partial_{y},\partial_{z})$,
$\mathbf{r}=(x,y,z)$, and the quasiparticle energy $\varepsilon$ is measured
relative to the Fermi level. The total Green’s function
$\hat{G}(\varepsilon,\mathbf{r})$ is comprised of the advanced
$G^{A}(\varepsilon,\mathbf{r})$, retarded $G^{R}(\varepsilon,\mathbf{r})$, and
Keldysh $G^{K}(\varepsilon,\mathbf{r})$, propagators. The $4\times 4$ Pauli
matrices in Nambu space are given by $\hat{\rho_{z}}=\tau_{z}\sigma_{0}$ and
$\hat{\rho_{0}}=\tau_{0}\sigma_{0}$, in which $\tau_{i}$ and $\sigma_{i}$ are
$2\times 2$ Pauli matrices in particle-hole and spin spaces, respectively.
Throughout the following calculations, an equilibrium state is considered so
that the retarded and Keldysh components of the total Green’s function can be
obtained from the advanced component by symmetry considerations. For example,
$G^{A}(\varepsilon,\mathbf{r})=-\\{\hat{\rho}_{z}G^{R}(\varepsilon,\mathbf{r})\hat{\rho}_{z}\\}^{\dagger}$
and
$G^{K}(\varepsilon,\mathbf{r})=\\{G^{R}(\varepsilon,\mathbf{r})-G^{A}(\varepsilon,\mathbf{r})\\}\tanh(\varepsilon
k_{B}T/2)$, where $k_{B}$ is the Boltzmann constant, and the system
temperature is denoted by $T$. To simulate the asymmetric Josephson junction
shown in Fig. 1, we assume that the superconducting leads are tunnel coupled
to the ferromagnetic region and can be described by boundary_c :

$\displaystyle\zeta\hat{G}\mathbf{n}\cdot\hat{{\bm{\nabla}}}\hat{G}=[\hat{G},\hat{G}_{\text{S}}],~{}~{}G^{R}_{\text{S}}=\left(\begin{array}[]{cc}{\cal
C}&{\cal S}e^{+i\varphi}\\\ {\cal S}e^{-i\varphi}&-{\cal
C}\end{array}\right).~{}~{}$ (11)

Here $\zeta$ is the ratio of the barrier resistance to the resistivity of the
normal layer, and $\mathbf{n}$ is the unit vector normal to the interfaces.
The retarded component of the total Green’s function inside a superconducting
lead can be expressed by $\hat{G}_{\text{S}}$ so that ${\cal
C}\equiv\cosh\theta\sigma_{0}$, ${\cal S}\equiv i\sinh\theta\sigma_{y}$, and
$\theta=\text{atanh}(\Delta/\varepsilon)$. Note that in the case of asymmetric
Josephson junction, one should replace $\Delta$ and $\varphi$ by
$\Delta_{1,2}$ and $\varphi_{l,r}$, respectively, according to Fig. 1.

The Usadel equation together with the boundary conditions create a set of
coupled complex boundary differential equations. To solve them accurately, we
make use of a so-called Riccati parametrization to help in establishing a
stable numerical algorithm. Two auxiliary unknown 2$\times$2 matrices, i.e.,
$\gamma$ and $\tilde{\gamma}$, are defined for parameterizing the retarded
Green’s function:

$\displaystyle
G^{R}(\varepsilon,\mathbf{r})=\left(\begin{array}[]{cc}(1-\gamma\tilde{\gamma})\Gamma&2\gamma\tilde{\Gamma}\\\
2\tilde{\gamma}\Gamma&(\tilde{\gamma}\gamma-1)\tilde{\Gamma}\end{array}\right),$
(14)

in which $\Gamma=(1+\gamma\tilde{\gamma})^{-1}$ and
$\tilde{\Gamma}=(1+\tilde{\gamma}\gamma)^{-1}$. Implementing the Riccati
parameterized Green’s function, the Usadel equation, Eq. (8), in the
nonsuperconducting region of Fig. 1 reads

$\displaystyle\sum_{k}\Big{\\{}\partial_{k}^{2}\gamma-2(\partial_{k}\gamma)\tilde{\gamma}\Gamma\partial_{k}\gamma\Big{\\}}=-2i\frac{\varepsilon}{D}\gamma,$
(15a)
$\displaystyle\sum_{k}\Big{\\{}\partial_{k}^{2}\tilde{\gamma}-2(\partial_{k}\tilde{\gamma})\gamma\tilde{\Gamma}\partial_{k}\tilde{\gamma}\Big{\\}}=-2i\frac{\varepsilon}{D}\tilde{\gamma}.$
(15b)

Here we have defined $k\equiv x,y,z$ for the spatial coordinates. To account
for ferromagnetism with an arbitrary exchange field, i.e.,
$\mathbf{h}=(h_{x},h_{y},h_{z})$, one simply needs to add
$(\mathbf{h}\cdot\mathbf{\sigma})\gamma-\gamma(\mathbf{h}\cdot\mathbf{\sigma}^{\ast})$,
and
$\tilde{\gamma}(\mathbf{h}\cdot\mathbf{\sigma})-(\mathbf{h}\cdot\mathbf{\sigma}^{\ast})\tilde{\gamma}$
into the Usadel equation [Eq. (15a) and (15b), respectively]. Also, the
boundary conditions in Eq. (11) for this parametrization scheme at $x=\mp d/2$
are:

$\displaystyle\partial_{x}\gamma=\pm(2\frac{{\cal C}_{1,2}}{{\cal
S}_{1,2}}+\gamma
e^{-i\varphi_{l,r}}-\frac{e^{+i\varphi_{l,r}}}{\gamma})\frac{{\cal
S}_{1,2}\gamma}{\zeta},~{}~{}~{}$ (16a)
$\displaystyle\partial_{x}\tilde{\gamma}=\pm(2\frac{{\cal C}_{1,2}}{{\cal
S}_{1,2}}+\tilde{\gamma}e^{+i\varphi_{l,r}}-\frac{e^{-i\varphi_{l,r}}}{\tilde{\gamma}})\frac{{\cal
S}_{1,2}\tilde{\gamma}}{\zeta}.$ (16b)

Finally, the charge current density in the equilibrium state is given by

${\bf J}(\mathbf{r})=\int
d\varepsilon\text{Tr}\left\\{\rho_{z}\left[\hat{G}(\varepsilon,\mathbf{r})\hat{{\bm{\nabla}}}\hat{G}(\varepsilon,\mathbf{r})\right]^{K}\right\\},$
(17)

where ‘Tr’ represents the trace operator. To obtain the total charge current
flowing across the junction shown in Fig. 1, one performs an spatial
integration over the charge flow component perpendicular to the junction
interfaces, namely, $J(x)=\int dy\int dz\;{\bf J}_{x}(\mathbf{r})$. Due to the
charge conservation law, $J(x)$ is a constant within the nonsuperconducting
region of Fig. 1.

##  III results and discussions

In the diffusive regime, all lengths are normalized by the superconducting
coherence length in the left superconductor, $\xi_{S}=\sqrt{\hbar
D/|\Delta_{1}|}$, and energies are scaled by the superconducting gap of the
left superconducting electrode at zero temperature, $|\Delta_{1}|$. For the
ballistic regime, unless otherwise indicated, all lengths are measured in
units of $k^{-1}_{F}$, where $k_{F}$ is the Fermi wavevector in the S regions.
We also have the dimensionless zero-temperature coherence length
$k_{F}\xi_{0}=(2/\pi)(E_{FS}/\Delta_{1})$, and fix $k_{F}\xi_{0}=100$. As in
the diffusive regime, energies are normalized by $\Delta_{1}$. The measure of
mismatch between the Fermi levels is given by the ratio
$\Lambda=E_{FM}/E_{FS}$.halterman2002 ; halterman2004 ; halterman2005
Throughout this paper, we assume $\Delta_{2}>\Delta_{1}$, and we consider a
uniform magnetization oriented along the $z$ direction so that
$\mathbf{h}=h_{z}=h$. Dimensionless units are implied with $\hbar=k_{B}=1$.

To gain a detailed view of the supercurrent profile in asymmetric Josephson
junctions, we divide the supercurrent into its three constituent parts.
According to Fig. 1, these parts consist of (i) the subgap supercurrent
$J_{1}$, with energies less than $\Delta_{1}$, (ii) the supergap supercurrent
$J_{2}$, with energies $\Delta_{1}<\varepsilon<\Delta_{2}$, and (iii) the
supergap supercurrent $J_{3}$, for scattering states with energies larger than
$\Delta_{2}$. When determining the supercurrent via Eq. (6) or Eq. (17), we
divide the energy integrals into three parts:
$J_{1}=\int_{0}^{\Delta_{1}}\square~{}d\varepsilon$,
$J_{2}=\int_{\Delta_{1}}^{\Delta_{2}}\square~{}d\varepsilon$, and
$J_{3}=\int_{\Delta_{2}}^{\infty}\square~{}d\varepsilon$. For all approaches,
the critical supercurrent is calculated in the usual way by finding the
maximum of the supercurrent within a phase difference interval of
$\Delta\varphi\in[0,2\pi]$, namely, $J_{c}=\max[J(\Delta\varphi)]$.

###  III.1 Ballistic $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ Josephson
junctions

To study the supercurrent profile in ballistic Josephson configurations with
asymmetric superconducting gaps, we employ the quantum particle-in-a-box
formalism described above K.Halterman_ss2016 . We emphasize that the
alternative approach, described in Appendix A, produces similar results and
has been used to study various systems, including Rashba-Dresselhaus spin-
orbit coupled, type-II Weyl semimetals, and black phosphorusAlidoustWS1 ;
AlidoustBP1 ; M.Alidoust2020 ; AlidoustBP2 ; AlidoustWS2 . Our comprehensive
numerical treatment was also found to agree with certain asymptotic limits,
including previous resultsbagwell  that considered quasiclassical one-
dimensional $\rm S_{1}NS_{2}$ junctions within the short and long junction
limits. These limitations were mainly imposed so that analytical solutions
could be found. Our fully microscopic approach however does not suffer from
these limitations and allows for investigations into asymmetric Josephson
configurations with more complicated band structures K.Halterman_ss2016 ;
K.Halterman2015 ; M.Alidoust2020 ; AlidoustBP1 ; AlidoustWS1 ; AlidoustBP2 .
In what follows, our results cover a broad range of geometrical and material
parameters, including junction thickness, Fermi energies, superconducting gap
ratio, and magnetization strength. With regards to normalization schemes, the
supercurrent density $J$ is normalized by $J_{0}\equiv en_{e}v_{F}$, where
$n_{e}$ is the bulk electron density and $v_{F}$ is the Fermi velocity. For
clarity, plots involving the supercurrent are also scaled by $10^{-2}$.

To begin, we plot the critical supercurrent as a function of
$\Delta_{2}/\Delta_{1}$ in Fig. 2. The strength of the magnetization in the
central junction region [see Fig. 1] varies as $h/\Delta_{1}=0,1,5,10,15$,
corresponding to both $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ systems. Each
panel in Figs. 2(a)-2(c) examines a different junction thickness, with
$k_{F}d_{F}=5,50,200$, respectively. Considering the nonmagnetic cases first
($h=0$), it is seen that the correlation between the two superconducting leads
(and thus the supercurrent), decays by increasing the junction thickness. This
effect becomes more pronounced when now considering uniformly magnetized
ferromagnets, as the pair-breaking exchange field in the magnet tends to
induce damped oscillations in the Cooper pair wave function with a
characteristic decay that goes as $1/h$. This causes $J_{c}$ to become
vanishingly small for $k_{F}d_{F}=200$ and $h/\Delta_{1}>5$. Note that the
results in Fig. 2 have no Fermi level mismatch ($\Lambda=1$), which when
present can amplify the supercurrent significantly, as will be seen below. One
pronounced feature seen in Fig. 2(a) for the short junction limit (i.e.,
$k_{F}d_{F}=5$), is the enhancement of the critical supercurrent by more than
$100\%$ when increasing the gap ratio to $\Delta_{2}/\Delta_{1}$ $\approx 25$,
for $h/\Delta_{1}<5$. This enhancement is diminished as the junction thickness
increases whereas the maximum enhancement of the critical current occurs at
significantly lower ratios of $\Delta_{2}/\Delta_{1}$ $\approx 1$. Also, as
$h$ increases, the critical current maximum gets shifted to larger
$\Delta_{2}/\Delta_{1}$ ratios and becomes relatively insensitive to changes
in gap asymmetry beyond $\Delta_{2}/\Delta_{1}\approx 10$.

Figure 2: (Color online). The critical supercurrent $J_{c}$ as a function of
superconducting gap ratio $\Delta_{2}/\Delta_{1}$ for three values of the
normalized junction thickness: $k_{F}d_{F}=5,50,200$, and five values of
magnetization strength $h/\Delta_{1}=0,1,5,10,15$. There is no Fermi level
mismatch ($\Lambda=1$).

We next investigate the effects of Fermi level mismatch, characterized by the
ratio of the Fermi levels in the two regions: $\Lambda=E_{FM}/E_{FS}$. For F/S
junctions, it was previously found gap ; halterman2002 ; halterman2004 ;
halterman2005  that the characteristic damped oscillations of the singlet pair
correlations within the ferromagnet become drastically modified for
$\Lambda<1$. When there is Fermi level mismatch between the junction layers,
the energy gap of the system gap  tends to close, as revealed in signatures of
the density of states gap ; halterman2004 ; halterman2005 . In Fig. 3, we plot
total critical supercurrent and its subgap and supergap components against the
Fermi energy ratio $\Lambda$ for a $\rm S_{1}FS_{2}$ junction with
$h/\Delta_{1}=15$. Two cases of gap anisotropy are shown in Fig. 3: (b)
$\Delta_{2}/\Delta_{1}$ $=5$ and (c) $\Delta_{2}/\Delta_{1}$ $=10$. The
symmetric case ($\Delta_{2}/\Delta_{1}$ $=1$) is also shown in panel (a) for
comparison purposes. The subgap critical current ($J_{c1}$) and the supergap
($J_{c2},J_{c3}$) critical currents are defined as
$J_{ci}\equiv|J_{i}(\Delta\varphi_{c})|$, where $\Delta\varphi_{c}$ is the
phase angle that leads to largest magnitude of the total supercurrent, i.e.,
$|J(\Delta\varphi_{c})|=J_{c}$. The results reveal that the total critical
supercurrent oscillates as a function of $\Lambda$. The overall profile of
$J_{c}$ remains approximately the same for all gap ratios
$\Delta_{2}/\Delta_{1}$, with an overall amplification of the magnitude of the
supercurrent when gap anisotropy is present. For the symmetric case
$\Delta_{2}/\Delta_{1}$ $=1$ in Fig. 3(a), the $J_{c2}$ supergap current
component vanishes and $J_{c3}$ is negligibly small. Thus, nearly the entire
contribution to the critical supercurrent comes from the $J_{c1}$ subgap
component. This is a general feature that arises in short Josephson junctions
with symmetric gaps, and it is clear that the Andreev bound states are the
dominate mechanism for supercurrent flow. Increasing the superconducting gap
ratio to $\Delta_{2}/\Delta_{1}$ $=5$, opens up the transport channel for
states between $\Delta_{1}$ and $\Delta_{2}$, and thus the $J_{c2}$ supergap
current component can contribute considerably to the total supercurrent.
Increasing the asymmetry further to $\Delta_{2}/\Delta_{1}$ $=10$, Fig. 3(c)
shows that the $J_{c2}$ contribution becomes even greater. In general, as
$\Lambda\rightarrow 0$, the number of available states for supercurrent flow
declines to zero. This also is true for the other extreme case of mismatch
with $J_{c}=0$ as $\Lambda\rightarrow\infty$. For the continuum of states with
energies exceeding $\Delta_{2}$, the quasiparticles are no longer confined to
the gap regions and lose phase-coherence, resulting in $J_{c3}$ being small
relative the other supercurrent components.

Figure 3: (Color online). Total critical supercurrent with its supergap
$J_{c2,c3}$ and subgap $J_{c1}$ components as a function of the Fermi level
mismatch ratio $\Lambda$. The labels show the different parameter values used
for the calculations. The junction thickness is set to $k_{F}d_{F}=5$ and the
ferromagnetic exchange field is $h/\Delta_{1}=15$.

In Fig. 4, we illustrate how $\Lambda$ affects the behavior of the subgap
($J_{c1}$) and the supergap ($J_{c2},J_{c3}$) critical current components as
functions of $\Delta_{2}/\Delta_{1}$. Figures 4(a)-4(e) correspond to the
nonmagnetic $\rm S_{1}NS_{2}$ case ($h=0$), while the right set 4(f) and 4(g)
exhibits the critical current behavior for $\rm S_{1}FS_{2}$ junctions with
$h/\Delta_{1}=15$. Considering first the nonmagnetic case, the left set of
panels shows that for a given gap ratio $\Delta_{2}/\Delta_{1}$, changing
$\Lambda$ modulates the critical current, consistent with the findings shown
in Fig. 3. When $\Lambda$ $=3$, the subgap supercurrent is strongly suppressed
and the supergap current dominates the net behavior of the critical current
throughout the wide range of $\Delta_{2}/\Delta_{1}$ considered. Indeed,
Andreev bound states with energies less than $\Delta_{1}$ play little, if any
role in the establishment of a supercurrent. For $\Lambda$ $=2,4$, the subgap
components all have the same trends, including the subgap supercurrent
$J_{c1}$, which exceeds the supergap currents for $\Delta_{2}/\Delta_{1}$
$\lesssim 5$, after which it decays considerably. The supergap current
$J_{c2}$ on the other hand, rapidly increases in the region
$\Delta_{2}/\Delta_{1}$ $\lesssim 25$, surpassing $J_{c1}$ at
$\Delta_{2}/\Delta_{1}$ $\approx 10$, before eventually leveling out close to
the total critical current curve. Thus, for $\Delta_{2}/\Delta_{1}$ $\gtrsim
25$, only energies that fall within $\Delta_{1}<\epsilon<\Delta_{2}$ are
needed when calculating the critical supercurrent response, while the Andreev
bound states with $\epsilon<\Delta_{1}$ and scattering states with
$\epsilon>\Delta_{2}$ can be neglected. It is evident that there is an
intricate and nontrivial relationship between the Fermi level mismatch and the
critical current. For situations where there is no mismatch $\Lambda=1$, Fig.
4(d) shows that the influence of $J_{c1}$ in this regime is significant, and
cannot be neglected for most gap ratios. In particular, for
$\Delta_{2}/\Delta_{1}$ $\leq 50$, the $J_{c1}$ component exceeds all other
components before slowly decaying at higher values of $\Delta_{2}/\Delta_{1}$.
In Fig. 4(e), where $\Lambda=0.5$, the crossover point occurs at the much
smaller $\Delta_{2}/\Delta_{1}$ $\approx 6$, indicating that both $J_{c1}$ and
$J_{c2}$ must be accounted for, even for moderate $\Delta_{2}/\Delta_{1}$
ratios. These results indicate that when characterizing the supercurrent
decomposition, the degree of Fermi level mismatch and gap mismatch play an
important role in which quasiparticle energies contribute to the supercurrent
response. Also, a notable feature in Figs. 4(a) and 4(b) is that $J_{c2}$
sometimes exceeds the total supercurrent. This arises mainly due to the subgap
$J_{1}$ and supergap $J_{2}$ currents flowing in opposite directions. This
important point shall be discussed further below.

Figure 4: (Color online). Critical supercurrent $J_{c}$ with its subgap
$J_{c1}$, and supergap $J_{c2},J_{c3}$ components shown as functions of the
superconducting gap ratio $\Delta_{2}/\Delta_{1}$. The junction thickness is
set at $k_{F}d_{F}=5$ for various values of
$\Lambda=E_{\text{FM}}/E_{\text{F}}=0.5,1,2,3,4$. The magnetization is set to
zero in (a)-(e) while $h/\Delta_{1}=15$ in panels (f)-(j).

Next, when the central layer possesses a uniform magnetization, the adjacent
panels in Figs. 4(f)-4(j) reveal a clear modification to all three components
of the critical current. For each of the four cases of Fermi level mismatch,
the component $J_{c1}$, has an extremely slow decay for gap ratios within
$\Delta_{2}/\Delta_{1}\gtrsim 5$. Thus within this regime, and for
quasiparticle energies with $\epsilon>5\Delta_{1}$, the critical current is
insensitive to the relative gap ratios characterizing the superconducting
leads. When $\Delta_{2}/\Delta_{1}=1$, corresponding to the commonly used
scenario of no asymmetry in the gaps, the supergap component $J_{c2}$
vanishes. In Figs. 4(f) and 4(h), when $\Lambda=4$ and $\Lambda=2$,
respectively, the subgap component $J_{c1}$ dominates all other critical
current components for relatively small gap ratios
$\Delta_{2}/\Delta_{1}\lesssim 2$. Both of the two components $J_{c1}$ and
$J_{c2}$ contribute equally to the critical current for
$\Delta_{2}/\Delta_{1}\approx 5$, with $J_{c2}$ the largest contributor for
larger gap asymmetry. The picture changes considerably when $\Lambda=3$, as
Fig. 4(g) illustrates that the conventional Andreev bound states characterized
by $J_{c1}$ play a minor role in the net supercurrent behavior when
$\Delta_{2}/\Delta_{1}\gtrsim 2$. When the Fermi levels are the same in each
segment of the junction [Fig. 4(i)], the subgap component becomes more
influential with $J_{c1}$ and $J_{c2}$ crossing at
$\Delta_{2}/\Delta_{1}\approx 4.5$. Finally, in Fig. 4(j), we consider
$E_{FM}<E_{F}$ and $\Lambda=0.5$. In this regime, both the subgap and supergap
components contribute nearly equally for most gap ratios
$\Delta_{2}/\Delta_{1}\gtrsim 15$, while less asymmetry again has $J_{c1}$ the
larger of the components. It should be noted that the enhancement of each of
the current components against $\Delta_{2}/\Delta_{1}$ occurs within a smaller
interval compared to the $\rm S_{1}NS_{2}$ junction.

Figure 5: (Color online). Critical supercurrent with its supergap and subgap
components as a function of uniform magnetization strength $h/\Delta_{1}$. In
the left column panels (a)-(e) the superconducting gap ratio is set to
$\Delta_{2}/\Delta_{1}$ $=1$ while in the right column panels (f)-(j)
$\Delta_{2}/\Delta_{1}$ $=5$.

In conventional $\rm S_{1}FS_{2}$ junctions, the exchange field induces damped
oscillations in the pair potential and the Josephson current exhibits
oscillations as a function of the ferromagnet exchange field.buzzed  To see
how this picture may change for asymmetric junctions, we present in Fig. 5 the
behavior of the critical supercurrent and its components versus the normalized
exchange field $h/\Delta_{1}$. The study covers the full range from
nonmagnetic, $h/\Delta_{1}=0$, to half-metallic, $h/\Delta_{1}\approx 160$,
where only one spin band is available. Each panel corresponds to one of the
five different Fermi level mismatch ratios $\Lambda$ that are considered (as
labeled). To further clarify the importance of supergap and subgap
supercurrents when there is gap asymmetry, we have considered
$\Delta_{2}/\Delta_{1}$ $=1$ in Figs. 5(a)-5(e) and set
$\Delta_{2}/\Delta_{1}$ $=5$ in Figs. 5(f)-5(j). The results in Figs.
5(a)-5(e) show that for symmetric junctions ($\Delta_{2}/\Delta_{1}$ $=1$),
the subgap supercurrent $J_{c1}$ is the dominant current component over the
entire range of magnetization strengths and Fermi level mismatches. The
supergap contribution $J_{c2}$ to the total critical current of course
vanishes entirely, while the $J_{c3}$ component arising from states in the
continuum is negligible. Hence, when the gap asymmetry vanishes
($\Delta_{2}/\Delta_{1}$ $=1$), the behavior of the total supercurrent can be
well captured through $J_{c1}$ only. It then suffices to take only the subgap
current components to accurately account for the general features of the
supercurrent, including its reversal in certain regions where it displays
cusps for a given $\Lambda$ and $h/\Delta_{1}$ (see, e.g.,
$h/\Delta_{1}\approx 30,110$ in Fig. 5(d)). Examining Figs. 5(a)-5(e), it is
clear that the critical current is nonmonotonic with a modulation that
strongly depends on $\Lambda$. Interestingly, tuning the Fermi level mismatch
to $\Lambda=3$ creates a situation where increasing the exchange field beyond
the first cusp at $h/\Delta_{1}\approx 90$ results in a dramatic rise in the
supercurrent response. Indeed, Fig. 5(b) demonstrates that compared to a
nonmagnetic junction ($h=0$), using a half-metallic insert
($h/\Delta_{1}\approx 160$) causes $J_{c}$ to nearly double.

If gap asymmetry is now introduced into the system
($\Delta_{2}/\Delta_{1}~{}=5$), Figs. 5(f)-5(j) show the emergence of the
$J_{c2}$ component, which can at times make sizable contributions to the total
supercurrent. It is evident that the cusps in $J_{c}$ where the current
reverses, can be ascribed to the cumulative effects of $J_{c1}$ and $J_{c2}$,
regardless of the Fermi level ratio $\Lambda$. Therefore, when
$\Delta_{2}/\Delta_{1}$ $\neq 1$, the subgap current component alone is unable
to provide an accurate and complete picture of the total supercurrent. Further
details on the origins of both the subgap and supergap supercurrent components
in terms of their discrete energy spectra and energy-resolved supercurrents is
given in Appendix B.

Figure 6: (Color online). Critical supercurrent with its subgap and supergap
components as a function of the normalized junction thickness $k_{F}d_{F}$ for
an asymmetric $\rm S_{1}FS_{2}$ junction. The gap asymmetry for both panels
(a) and (b) corresponds to $\Delta_{2}/\Delta_{1}=5$ and the normalized
exchange field is set to $h/\Delta_{1}=15$. In (a) the Fermi level mismatch is
set to $\Lambda=2$, while in (b) the relative Fermi levels are the same
($\Lambda=1$).

Figure 7: (Color online). Critical supercurrent with its subgap and supergap
components as a function of junction thickness for differing chemical
potential ratios $\Lambda$ $=0.5,1,2,3,4$. The superconducting gap ratio and
exchange field strength are set to $\Delta_{2}/\Delta_{1}$ $=5$ and
$h/\Delta_{1}=15$, respectively.

Figure 8: (Color online). Supergap ($J_{2,3}$) and subgap ($J_{1}$) current
components as a function of the superconducting phase difference
$\Delta\varphi$. (a)-(d): The superrcurent profile for a $\rm S_{1}NS_{2}$
junction. (e)-(h): The supercurrent profile for a $\rm S_{1}FS_{2}$ junction
with $h/\Delta_{1}=15$. The junction thickness is set fixed to $k_{F}d_{F}=5$,
and the relative Fermi levels are equal ($\Lambda=1$).

The study of the critical current dependence on the ferromagnet thickness in
$\rm S_{1}FS_{2}$ Josephson junctions has been extensively investigated both
theoreticallybuzzed  and experimentallykontos2002 ; shell2006 . Proximity
effects arising from the coupling of the outer superconducting banks and
ferromagnet leads to oscillations of the pair amplitude in the ferromagnet.
Under certain conditions, these oscillations can cause the ground state of the
system to transition to a state with $\Delta\varphi=\pi$. During these
transitions, the supercurrent reverses direction and reveals itself as cusps
near the minima of the critical current. Unfortunately, the study of transport
in clean ferromagnetic Josephson junctions with Fermi level mismatch and gap
asymmetry from a purely microscopic perspective is lacking. Therefore, to
address these deficiencies, we present in Fig. 6, the critical current and its
associated components as a function of normalized junction thickness
$k_{F}d_{F}$. The junction possesses a gap asymmetry of
$\Delta_{2}/\Delta_{1}=5$, and two different values of the Fermi level
mismatch parameter $\Lambda=1,2$ are considered. Since the damped oscillations
in the ferromagnet are governed by the spin-split Fermi wavevectors there,
having $E_{F}$ vary across different segments of the Josephson junction can
modify the oscillatory period of the pair amplitude. The spin splitting in the
ferromagnet introduces the length scale $\xi_{F}$ set by the difference in the
spin up and spin down Fermi wavevectors,
$\xi_{F}\propto(k_{F\uparrow}-k_{F\downarrow})^{-1}$. Accordingly, the cusps
are found to repeat in intervals of $\pi\xi_{F}\approx\pi
E_{F}\sqrt{\Lambda}/h$, which for $h/\Delta_{1}=15$ and $\Lambda=1,2$,
corresponds to $\pi\xi_{F}\approx 32.9,46.9$, respectively. This is seen when
comparing Figs. 6(a) and 6(b). Note that the damped oscillations in the
critical current with thickness have also been observed in the critical
temperature for ballistic spin valves.zep

As seen in Fig. 6, for a given $\Lambda$, the supergap component $J_{c2}$
contributes the most for small ferromagnet thicknesses. To explore this
further, Fig. 7 displays the critical current and its components over a
narrower range of thicknesses for the $\rm S_{1}FS_{2}$ configuration. To
reveal how the oscillations and magnitude of the supercurrent changes with
variations in the Fermi level ratios, a broad range of mismatch parameters
$\Lambda$ is considered. As seen in Figs. 7(a)-7(e), the $J_{c3}$ component is
negligible, as scattering states again contribute little to the supercurrent
response. The total supercurrent oscillates as a function of $k_{F}d_{F}$, but
in contrast to Fig. 6, these oscillations are over the much smaller Fermi
length scale. These small-scale oscillations are neglected in quasiclassical
treatments where atomic scale features are eliminated. Moreover, increasing
$\Lambda$ is shown to decrease the period of oscillations in Figs. 3(a)- 3(e),
which is consistent with the corresponding increases in the ratios of the
Fermi wavevectors in the ferromagnet and superconductor regions
$k_{FM}/k_{F}$. For Fermi level ratios corresponding to $\Lambda=4$ and
$\Lambda=3$, Figs. 3(a) and 3(b) illustrate that while $J_{c2}$ has gradual
variations as the thickness changes, the subgap component $J_{c1}$ has
pronounced oscillations that periodically vanish (or nearly so) for thin
ferromagnets. Therefore certain thicknesses can be chosen that result in
maximal supercurrent flow that is almost entirely comprised of supergap
$J_{c2}$ states. As the difference between Fermi levels lessens, Figs. 3(c)
and 3(d) show that the oscillatory behavior dampens out for thicker
ferromagnets. Finally in Fig. 3(e), we find that the critical current mainly
declines rapidly as the normalized thickness increases and then levels off to
greatly diminished values.

We now present the current-phase relations for the total supercurrent and its
components at specific parameter values corresponding to points of interest
found in earlier critical supercurrent plots. In Figs. 8(a)-8(d), we have set
the exchange field $h$ to zero and consider a $\rm S_{1}NS_{2}$ configuration
with varying degrees of gap asymmetry. As seen in Fig. 8(a), when there is no
gap asymmetry ($\Delta_{2}/\Delta_{1}=1$), the $J_{3}$ component is slightly
discernible, as states in the continuum contribute little to the supercurrent
in short junctions. Also, the $J_{2}$ supergap component vanishes, as expected
for the symmetric case. Thus the current in this case is carried nearly
entirely by bound states. Examining Figs. 8(b)-8(d), it is evident that by
increasing the asymmetry through the superconducting gap ratio,
$\Delta_{2}/\Delta_{1}$, a supergap current $J_{2}$ emerges, and begins to
play a more impactful role in the total supercurrent response. In contrast,
the subgap current $J_{1}$ accordingly becomes weaker, as its contribution to
at least 50% of the total supercurrent is limited to relatively narrow phase
differences, $130^{\circ}\lesssim\Delta\varphi\lesssim 160^{\circ}$. Despite
the significant changes to $J_{1}$ and $J_{2}$, the overall supercurrent,
$J_{\rm tot}$, increases, but retains its overall profile. Increasing
$\Delta_{2}/\Delta_{1}$ also is seen to enhance the supercurrent overall,
while shifting the supercurrent peaks to smaller $\Delta\varphi$. The
enhancement of the supercurrent follows in part from the broadening of the
discrete energy states that occurs for larger ratios $\Delta_{2}/\Delta_{1}$
[see Appendix B]. The microscopic numerical results presented in Figs.
8(a)-8(d) are also consistent with one-dimensional quasiclassical
modelsbagwell . In Appendix B, Fig. 12 reveals the interplay between the bound
and scattering states related to Fig. 8, and are discussed in terms of the
energy-resolved and phase-resolved current density.

Figure 9: (Color online). Total current-phase relation with its supergap and
subgap components for several key exchange-field values around a current
reversal point. The superconducting gap ratio is set to
$\Delta_{2}/\Delta_{1}$ $=5$, the junction thickness is $k_{F}d_{F}=5$, and
there is no mismatch in Fermi levels ( $\Lambda$ $=1$).

Next, upon incorporating a uniform magnetization of $h/\Delta_{1}=15$, Figs.
8(e)-8(h) illustrate how this leads to a drastic modification to the profiles
of the current-phase relations. For $\rm S_{1}FS_{2}$ junctions, we find a
slight drop to the overall net supercurrent response for the whole range of
gap asymmetries considered. Increasing the gap asymmetry tends to enhance the
overall supercurrent, but the general current-phase profile remains relatively
unchanged, with minimal change in the peak locations. One noticeable
difference however occurs for $\Delta_{2}/\Delta_{1}=2,5$, where the
supercurrent undergoes a current reversal at $\varphi\approx 162^{\circ}$. We
find that similar to the symmetric $\rm S_{1}NS_{2}$ junction, Fig. 8(e) shows
a small contribution from $J_{3}$ arising from scattering states
($\epsilon>\Delta_{1}$) when there is gap symmetry. For extreme gap asymmetry
$\Delta_{2}/\Delta_{1}=100$, shown in in Fig. 8(h), the supergap component
$J_{2}$ is the main contributor to the overall current, and we see that
although $J_{1}$ has broadened compared to the nonmagnetic case, it never
exceeds approximately 40% of the total critical current.

Figure 10: (Color online). Critical supercurrent in a diffusive $\rm
S_{1}FS_{2}$ Josephson junction as a function of $\Delta_{2}/\Delta_{1}$. The
junction thickness is fixed at $d=0.8\xi_{S}$ and various values of the
normalized exchange field are considered: $h/\Delta_{1}=0$, $0.4$, $0.8$,
$1.2$, $1.6$, $2.6$, $2.8$, $3.0$, $3.5$, $4.0$, $4.5$. Figure 11: (Color
online). The three components of supercurrent at four different values of the
gap anisotropy ratio: $\Delta_{2}/\Delta_{1}$=$1.0$, $1.5$, $2.0$, and $2.6$.
In (a)-(d), the exchange field $h$ is set to zero, and in (e)-(h), we set to
$h=2.6\Delta_{1}$.

To elaborate on the damped oscillations in the pairing correlations that
induce $0-\pi$ transitions responsible for supercurrent reversals and cusps in
the critical current profile of asymmetric junctions, we have studied the
exchange field dependence of the current-phase relations in Fig. 9. The total
supercurrent $J_{\rm tot}$ is shown alongside its constitutive components
$J_{1,2,3}$. As a representative parameter set, we have chosen the first
crossover state in Fig. 5(i), which occurs around $h/\Delta_{1}\approx 32$.
Therefore, in Figs. 9(a)-9(e) we set $\Delta_{2}/\Delta_{1}$ $=5$, $\Lambda$
$=1$, and $k_{F}d_{F}=5$ and vary the normalized magnetization according to
$h/\Delta_{1}=26,28,31,34,36$, respectively. As seen, the overall profile and
magnitude of the subgap supercurrent component $J_{1}$ is relatively
unaffected as the exchange field increases. The supergap current $J_{2}$
however experiences a transition in which it goes from entirely positive for
$h/\Delta_{1}=26$ to negative for most phase differences $\Delta\varphi$ when
$h/\Delta_{1}=36$. This results in an overall suppression of the total
supercurrent, illustrating the influential control of the $J_{2}$ component on
the total supercurrent. Therefore, adjacent to a $0-\pi$ crossover, the
supergap $J_{2}$ and subgap $J_{1}$ current components propagate in opposite
directions, creating a competing situation. Below, we shall see that a similar
occurrence takes place near the $0-\pi$ transition point in the diffusive
regime. It should be noted that the supercurrent in the ballistic regime
involves the superposition of quasiparticle trajectories undergoing normal and
Andreev reflections at the two superconductor interfaces. Due to the
microscopic method used, length scales as small as the Fermi wavelength are
included, permitting the capture of Friedel-like oscillations,sipr  which
emerge as highly oscillatory signatures in the supercurrent response for small
normalized layer thicknesses $k_{F}d_{F}$. The resonant modes are also highly
sensitive to the relevant geometrical and material parameters such as the
junction length, Fermi level, magnetization strength, and superconducting gap
ratio, all of which intricately combine to give the observed small
oscillations in e.g., Figs. 7 and 8. Nonetheless, the final conclusions and
central findings are clearly independent of these subfeatures.

###  III.2 Diffusive $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ Josephson
junctions

We now consider the supercurrent response in asymmetric diffusive Josephson
junctions. To properly capture quasiparticles with energies deemed relevant to
the net supercurrent, we have set an energy cutoff of
$\varepsilon_{max}=25\Delta_{1}$ when performing the integration over
quasiparticle energies. We have also used representative values for the
junction thickness and interface opacity, with $d=0.8\xi_{S}$ and $\zeta=4$,
respectively. The maximum supercurrent flow as a function of
$\Delta_{2}/\Delta_{1}$ is shown in Fig. 10. The maximum of critical
supercurrent occurs when the magnetization in the junction vanishes. As seen,
the critical supercurrent enhances by $50\%$ when
$\Delta_{2}/\Delta_{1}\approx 10$ and $h=0$. Upon increasing the strength of
the uniform exchange field $h$, the critical supercurrent becomes suppressed.
For weak ferromagnets with exchange energies corresponding to $h\sim
2.6\Delta_{1}$, the supercurrent undergoes a sign reversal as the gap
asymmetry $\Delta_{2}/\Delta_{1}$ is varied. By further increasing
$h/\Delta_{1}$ to $2.8,3.0,4.0,4.5$, the overall critical supercurrent is
enhanced and returns to its previous monotonic growth as a function of
$\Delta_{2}/\Delta_{1}$. In contrast to the other normalized exchange field
strengths, the case with $h/\Delta_{1}=2.6$ exhibits a clear supercurrent
suppression for $\Delta_{2}/\Delta_{1}\lesssim 7$. Increasing the exchange
field to $h/\Delta_{1}=2.8,3.0$, causes the supercurrent to have a short-lived
enhancement, and a subsequent cusp at $\Delta_{2}/\Delta_{1}\approx 3$, before
slowly declining as the gap ratio increases. To create a stable numerical
scheme, we have introduced a small imaginary part $\delta=0.001$ to the energy
of the quasiparticles, $\varepsilon$. This imaginary part can act as a source
of inelastic scattering, which if increased, can wash out the dominant and
important parts of the curves, such as the overall supercurrent response and
associated $0$-$\pi$ transitions. For a relatively small imaginary part, as is
considered throughout our numerical study, the algorithm can at times
introduce insignificant artifacts that show up as small oscillations in the
supercurrent.

To gain a comprehensive picture of the supercurrent in a diffusive asymmetric
Josephson junction, we have plotted the components $J_{1,2,3}(\Delta\varphi)$,
and total supercurrent $J(\Delta\varphi)$ in Fig. 11. Various levels of gap
asymmetry are shown corresponding to $\Delta_{2}/\Delta_{1}=1.0,1.5,2.0,2.6$.
The nonmagnetic $\rm S_{1}NS_{2}$ case is shown in Figs. 11(a)-11(d), whereas
Figs. 11(e)-11(h) correspond to a $\rm S_{1}FS_{2}$ junction. As seen for both
cases, when $\Delta_{2}=\Delta_{1}$, the $J_{2}(\Delta\varphi)$ component
vanishes, as expected. Also, there is an overall reduction of the supercurrent
magnitudes in the magnetic case. Other than their magnitudes, the current-
phase relations of the $J_{2}(\Delta\varphi)$ components are similar for $\rm
S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ junctions with asymmetric superconducting
gaps (albeit with different signs). The subgap supercurrent
$J_{1}(\Delta\varphi)$ behaves slightly differently for $\rm S_{1}NS_{2}$ and
$\rm S_{1}FS_{2}$ junctions. In the latter case, when $\Delta_{2}=\Delta_{1}$,
the supercurrent deviates strongly from the usual sinusoidal relation. In
contrast to the $\rm S_{1}NS_{2}$ case, the component $J_{3}(\Delta\varphi)$
shows a $\sim\sin 2\Delta\varphi$ relation around $\Delta_{2}\approx
2\Delta_{1}$, and changes sign for larger values. The competition between
these components results in the total supercurrents shown in the far right
panels. As is apparent, unlike the strikingly different responses of
$J_{1,2,3}(\Delta\varphi)$ for $\Delta_{2}>\Delta_{1}$, the total supercurrent
changes uniformly except when transitioning from $\Delta_{2}=\Delta_{1}$ to
$\Delta_{2}>\Delta_{1}$ for the $\rm S_{1}FS_{2}$ junction with
$h=2.6\Delta_{1}$. This variation results in the form of the current-phase
relation changing from $\sim\sin 2\Delta\varphi$ to $\sim\sin\Delta\varphi$.
Note that the total supercurrent response in Fig. 11(h) has the form $\sim\sin
2\Delta\varphi$, which appears due to the competition between the subgap
supercurrent states comprising $J_{1}(\Delta\varphi)$ and the scattering
states which embody $J_{3}(\Delta\varphi)\sim\sin\Delta\varphi$. The
competition originates from the opposite propagation directions of the
$J_{1}(\Delta\varphi)$ and $J_{3}(\Delta\varphi)$ current components. These
findings thus complement the ballistic results that found many instances where
supergap modes must be accounted for appropriately to obtain accurate and
reliable results. Further insight into the supergap and subgap responses are
presented in Appendix C, where the energy-resolved and phase-resolved
supercurrent density is analyzed.

##  IV Conclusions

By employing complementary numerical approaches in the ballistic and diffusive
regimes, we have performed a comprehensive study of supercurrent flow through
asymmetric $\rm S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ Josephson junctions where
the superconducting gap in the $\rm S_{1,2}$ regions are unequal, i.e.,
$\Delta_{2}\neq\Delta_{1}$. In the ballistic regime, we have directly solved
the Bogoliubov de-Gennes Hamiltonian that allows for exploring a parameter
space with a wide range of energy and length scales, whereas when impurities
and disorder are present, we make use of the full proximity limit of the
quasiclassical regime. Our results found that for asymmetric junctions with
$\Delta_{2}/\Delta_{1}$ $\approx 25,10$, the critical supercurrent can be
enhanced by more than $100\%$ and $50\%$ in the ballistic and diffusive $\rm
S_{1}NS_{2}$ cases, respectively. Our results in the ballistic cases reveal
that when $\Delta_{2}/\Delta_{1}$ $=1$, the subgap current is the main
contributor to total supercurrent. Introducing an imbalance to the
superconducting gap ratio $\Delta_{2}/\Delta_{1}$ $>1$, the supergap currents
were discovered to play key roles and for certain parameter values were the
main contributors to the total supercurrent. Through our investigations of
asymmetric junctions (with $\Delta_{2}/\Delta_{1}$ $>1$), the current phase
relations with their supergap and subgap current components were explored
around $0-\pi$ current crossover points. We found that the emergence of second
harmonics in the current-phase relations of $\rm S_{1}FS_{2}$ junctions is a
direct consequence of the competition between subgap and supergap current
components with opposite flow directions. It was shown in an earlier
workbagwell  that supergap currents are relatively insensitive to temperature
compared to the subgap component, as the former originates from coherent
evanescent modes in the continuum, whereas the latter is carried through
resonant bound states. Therefore, the findings of this paper should serve to
stimulate experiments that pave the way for designing new superconducting
devices that utilize robust supergap currents. The asymmetric $\rm
S_{1}NS_{2}$ and $\rm S_{1}FS_{2}$ structures studied here can apply to
Josephson configurations where the amplitude of the superconducting gaps might
fluctuate independently when the system is subject to a strong external
magnetic field or high temperatures near the critical temperature.
Furthermore, the enhancement of the critical supercurrent due to
$\Delta_{1}\neq\Delta_{2}$ can be beneficial in magnetic Josephson junctions
that suffer from weakened currents in the presence of ferromagnetism.

###### Acknowledgements.

M.A. is supported by Iran’s National Elites Foundation (INEF). K.H. is
supported in part by ONR and a grant of HPC resources from the DOD HPCMP.

Figure 12: (Color online). Total supercurrent density mappings as a function
of phase difference $\Delta\varphi$ and energy $\varepsilon$ for a $\rm
S_{1}NS_{2}$ Josephson junction with varying levels of gap asymmetry: (a)
$\Delta_{2}/\Delta_{1}=1$, (b) $\Delta_{2}/\Delta_{1}=2$, (c)
$\Delta_{2}/\Delta_{1}=5$, and (d) $\Delta_{2}/\Delta_{1}=100$. The Fermi
levels are equal throughout the system ($\Lambda=1$).

Figure 13:  (Color online). The energy-resolved supercurrent density
$J_{\varepsilon}$ for an $\rm S_{1}FS_{2}$ junction with
$\Delta_{2}/\Delta_{1}=5$, and Fermi level mismatch $\Lambda=3$. Three
different values of the normalized exchange energy $h/\Delta_{1}$ are
considered: (a) $h/\Delta_{1}=0$, (b) $h/\Delta_{1}=45$, and (c)
$h/\Delta_{1}=90$. The macroscopic phase difference chosen in each case
corresponds to $\Delta\varphi_{c}$, where $\Delta\varphi_{c}$ is the phase
angle that leads to the largest magnitude of the total supercurrent. In panels
(a) and (b) $\Delta\varphi_{c}=106^{\circ}$, and for (c) we have
$\Delta\varphi_{c}=125^{\circ}$.

Figure 14: (Color online). Total supercurrent density mappings as a function
of phase difference $\Delta\varphi$ and energy $\varepsilon$ for an asymmetric
Josephson junction with $\Delta_{2}/\Delta_{1}=5$. The Fermi level mismatch
parameter is set to $\Lambda=3$. Panels (a)-(h) depict differing normalized
exchange fields $h/\Delta_{1}$ corresponding to $0,45,70,90,110,130,140$, and
$160$, respectively.

Figure 15:  (Color online). The discrete energy spectrum for an asymmetric
Josephson junction with $\Delta_{2}/\Delta_{1}=5$, and Fermi level mismatch
corresponding to $\Lambda=3$. The normalized exchange fields considered
correlate with Fig. 14 above.

##  Appendix A Alternative wave-function approach: ballistic regime

The method outlined in Sec. II.1 provides an effective way to solve Josephson
junction systems with limited approximations in the ballistic regime. There
also exists wave-function approaches that can provide exact solutions, one of
which we outline belowAlidoustWS1 ; AlidoustBP1 ; M.Alidoust2020 ; AlidoustBP2
; AlidoustWS2 .

To simulate the low-energy physics of heterostructures in the presence of a
magnetization with arbitrary direction, $\mathbf{h}=(h_{x},h_{y},h_{z})$, one
employs an effective single-particle Hamiltonian:

$H=\frac{1}{2}\int
d\mathbf{p}~{}\hat{\psi}^{{\dagger}}(\mathbf{p})\left[\frac{\mathbf{p}^{2}}{2m}+{\bm{\sigma}}\cdot\mathbf{h}\right]\hat{\psi}(\mathbf{p}),$
(18)

where $\mathbf{p}=(p_{x},p_{y},p_{z})$ is the momentum and $m$ is the
effective mass of a charged particle. The associated field operator in spin
space is given by
$\hat{\psi}=(\psi_{\uparrow},\psi_{\downarrow})^{\mathrm{T}}$ and
${\bm{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector comprised of
Pauli matrices. The spin-singlet superconductivity in the BCS scenario can be
described by the following electron-electron amplitudes:

$\Delta\langle\psi_{\uparrow}^{\dagger}\psi_{\downarrow}^{\dagger}\rangle+\text{H.c.}~{}.$
(19)

Accounting for the electron-electron amplitudes in the BdG formalism, the low-
energy Hamiltonian in spin-Nambu space reads:

${\cal
H}(\mathbf{p})=\left(\begin{array}[]{cc}H(\mathbf{p})-\mu\hat{1}&\hat{\Delta}\\\
\hat{\Delta}^{\dagger}&-H^{\dagger}(-\mathbf{p})+\mu\hat{1}\end{array}\right),$
(20)

in which $\hat{\Delta}$ is the superconducting gap
$\hat{\Delta}\equiv(\Delta,-\Delta)e^{i\varphi_{l,r}}$, and $\mu$ is the
chemical potential. The field operators in the rotated spin-Nambu space are
given by
$\hat{\psi}=(\psi_{\uparrow},\psi_{\downarrow},\psi_{\downarrow}^{{\dagger}},-\psi_{\uparrow}^{{\dagger}})^{\mathrm{T}}$.
The continuity equation for charged carriers is written:

$\begin{split}\frac{\partial\rho_{\text{c}}}{\partial
t}=\lim\limits_{\mathbf{r}\rightarrow\mathbf{r}^{\prime}}\sum\limits_{\sigma\tau\sigma^{\prime}\tau^{\prime}}\frac{1}{i}\Big{[}\psi^{\dagger}_{\sigma\tau}(\mathbf{r}^{\prime}){\cal
H}_{\sigma\tau\sigma^{\prime}\tau^{\prime}}(\mathbf{r})\psi_{\sigma^{\prime}\tau^{\prime}}(\mathbf{r})\\\
-\psi^{\dagger}_{\sigma\tau}(\mathbf{r}^{\prime}){\cal
H}_{\sigma\tau\sigma^{\prime}\tau^{\prime}}^{\dagger}(\mathbf{r}^{\prime})\psi_{\sigma^{\prime}\tau^{\prime}}(\mathbf{r})\Big{]},\end{split}$
(21)

where ${\cal H}_{\sigma\tau\sigma^{\prime}\tau^{\prime}}$ is the component
form of Eq. (20) and $\sigma,\tau$ indices label the spin and particle-hole
degrees of freedom, respectively. In a situation where no sink or source of
charge is present, the time variation of charge density vanishes:
$\partial_{t}\rho_{\text{c}}\equiv~{}0$. Accounting for the current
conservation law, the charge current density can be expressed by:

$\mathbf{J}=\int
d\mathbf{r}\Big{\\{}\hat{\psi}^{\dagger}(\mathbf{r})\overrightarrow{{\cal
H}}(\mathbf{r})\hat{\psi}(\mathbf{r})-\hat{\psi}^{\dagger}(\mathbf{r})\overleftarrow{{\cal
H}}(\mathbf{r})\hat{\psi}(\mathbf{r})\Big{\\}},$ (22)

where the real-space Hamiltonian ${\cal H}(\mathbf{r})$ is obtained by
substituting
$\mathbf{p}\equiv~{}{-im^{-1}(\partial_{x},\partial_{y},\partial_{z})}$ in Eq.
(20). The arrows indicate the specific wave functions that the Hamiltonian
operates on. To obtain the supercurrent, one computes the current density
perpendicular to the interfaces, in our geometry shown in Fig. 1, $J_{x}$, and
integrates over the junction cross section in the $y$ direction:
$J(\Delta\varphi)=J_{0}\int_{-W/2}^{+W/2}dyJ_{x}(x,y,\Delta\varphi)$. Here
$J_{0}=2e|\Delta|/\hbar$, and $e$ is the electron charge. Upon diagonalizing
${\cal H}(\textbf{p})$ in Eq. (20), one obtains the electronic wave functions
$\hat{\psi}_{l,m,r}(\textbf{p})$ within the left ($l$), middle ($m$), and
right ($r$) regions. Next, the wave functions are matched at the left
$\hat{\psi}_{l}$=$\hat{\psi}_{m}|_{x=0}$ and the right boundaries
$\hat{\psi}_{m}$=$\hat{\psi}_{r}|_{x=d}$. The continuity conditions must also
be satisfied: $(\partial_{\textbf{p}}{\cal
H}_{l})_{\textbf{r}}\hat{\psi}_{l}$=$(\partial_{\textbf{p}}{\cal
H}_{m})_{\textbf{r}}\hat{\psi}_{m}|_{x=0}$, $(\partial_{\textbf{p}}{\cal
H}_{m})_{\textbf{r}}\hat{\psi}_{m}$=$(\partial_{\textbf{p}}{\cal
H}_{r})_{\textbf{r}}\hat{\psi}_{r}|_{x=d}$. The index r indicates a switch to
real space after taking the derivatives in momentum space. It is important to
note that we apply no simplifying assumptions and approximations to the wave
functions in the numerical calculations. This however results in highly
complicated and lengthy expressions for the wave functions and supercurrent.
We therefore are only able to evaluate them numerically.

##  Appendix B energy dispersion and energy-resolved supercurrent density:
ballistic regime

When calculating the supercurrent via Eq. (6), all positive energy states
within an energy cutoff $\epsilon_{c}$ are summed over, with $\epsilon_{c}$
chosen sufficiently large so that including any additional states has no
distinguishable effects on the results. It is important to note that when
summing the quantum states for the supercurrent, the summation implicitly
includes an integration over the continuum of states with transverse energy
$\epsilon_{\perp}$ [see Eq. (2)]. To isolate the supercurrent contribution at
supergap and subgap energies, it is beneficial to extract the supercurrent
response as a function of the quasiparticle energy $\epsilon$. This procedure
involves calculating the supercurrent for each quasiparticle trajectory with
associated energy $\epsilon_{\perp}$. All trajectories are then summed over to
arrive at the supercurrent for a given energy.

This procedure gives the results shown in Fig. 12, where the supercurrent is
mapped out as a function of energy and phase difference. For concreteness, we
take the parameters used in Figs. 8(a)-8(d), where a broad range of gap
asymmetries were considered. Note the emergence of the $J_{2}$ supercurrent in
the current phase relations can be accounted for in Figs. 12(b)-12(d) where
the current carrying states get shifted upwards into the supergap region and
broaden with increased $\Delta_{2}/\Delta_{1}$. This also leads to an
amplification of the $J_{2}$ component, and in turn the total supercurrent.

Figure 16: (Color online). Color map of the supercurrent density as a function
of the normalized quasiparticle energy $\varepsilon/\Delta_{1}$, and
superconducting phase difference $\Delta\varphi$ in the diffusive regime. The
junction parameter values are set in each panel to those of Fig. 11, except
now we have (a) $\rm S_{1}NS_{2}$: $h=0$, $\Delta_{2}/\Delta_{1}=1$, (b) $\rm
S_{1}NS_{2}$: $h=0$, $\Delta_{2}/\Delta_{1}=2.6$, (c) $\rm S_{1}FS_{2}$:
$h=2.6\Delta_{1}$, $\Delta_{2}/\Delta_{1}=1$, and (d) $\rm S_{1}FS_{2}$:
$h=2.6\Delta_{1}$, $\Delta_{2}/\Delta_{1}=2.6$.

The energy-resolved supercurrent for an asymmetric Josephson junction at fixed
phase is shown in Fig. 13. We take a relative ratio of
$\Delta_{2}/\Delta_{1}=5$, and mismatch in Fermi levels corresponding to
$\Lambda=3$. Figure 13(a) corresponds to a nonmagnetic junction while Figs.
13(b) and 13(c) describe a ferromagnetic junction with differing exchange
field strengths, as shown. For consistency, these system parameters correlate
with Fig. 5(g). To correlate with the critical current components, from Fig.
5(g), it is seen that for $h/\Delta_{1}=0$ the supergap component $J_{2}$
dominates, while for $h/\Delta_{1}=90$, the subgap $J_{1}$ component does. The
intermediate case of $h/\Delta_{1}=45$ has the subgap $J_{1}$ and supergap
$J_{2}$ components contributing nearly equally to the critical current. This
behavior is accounted for in the energy dependence of the supercurrent
presented in Figs. 13(a)-13(c), where we take
$\Delta\varphi=\Delta\varphi_{c}$ in each case to ensure that it gives the
critical current shown in Fig. 5(g). The shift in current carrying states is
evident as the exchange energy increases, until eventually reversing direction
for $h/\Delta_{1}=90$ where the subgap bound states dominate.

To give a more comprehensive view of the energy-resolved supercurrent we
present in Fig. 14, the normalized supercurrent as a function of the
normalized energy $\epsilon/\Delta_{1}$ and phase difference $\Delta\varphi$.
Eight different exchange fields are considered, and the system parameters
again coincide with the critical current study in Fig. 5(g). As seen in Figs.
14(a)-14(h), the supercurrent profiles exhibit varying amounts of mixtures
between the subgap and supergap states. As the exchange field increases, the
supercurrent flow evolves to contain states that have oppositely directed flow
in both the subgap ($\epsilon\leq\Delta_{1}$) and supergap
($\Delta_{1}\leq\epsilon\leq\Delta_{2}$) regions. Remarkably, increasing the
exchange field can increase the supercurrent, including the half-metallic
limit [Fig. 14(h)], where the supercurrent can exceed the supercurrent found
in Fig. 14(a) for the nonmagnetic case. This can have important consequences
for devices that utilize the spin degree of freedom in Josephson junction
systems. Note that the various admixtures of subgap and supergap supercurrents
exhibited here can be directly correlated with the critical current signatures
found in Fig. 5(g).

To delve further into the phase dependence of the supercurrent flow, we next
present in Fig. 15 the quasiparticle energy spectra $\epsilon_{n}$ for each of
the cases shown in Fig. 14. We consider the $\epsilon_{\perp}=0$ mode, noting
that other transverse modes have similar behavior, and the cumulative effect
serves to only broaden the overall supercurrent profile. From the diagrams, it
is clear that the current vanishes at $\Delta\varphi=0^{\circ}$ and
$\Delta\varphi=180^{\circ}$, where
$\partial\epsilon_{n}/\partial(\Delta\varphi)=0$. The cusps in the energy
dispersion are consistent with Fig. 14, where the current at certain energies
becomes reversed. We also find that as the exchange field $h$ increases, the
additional branches of the energy dispersion which emerge increase in
separation before coalescing at high exchange fields. Thus, although the
supercurrent is a cumulation of quasiparticle amplitudes and energies, the
energy spectrum alone gives valuable insight into the transport properties of
asymmetric Josephson junctions.

##  Appendix C energy-resolved supercurrent density: diffusive regime

When impurity scattering dominates, we turn to the diffusive regime. In Fig.
16, we show the supercurrent response as a function of macroscopic phase
difference $\Delta\varphi$ and normalized energy $\epsilon/\Delta_{1}$. In
Figs. 16(a) and 16(b), a nonmagnetic junction $h=0$ is considered. In Fig.
16(a), the junction is symmetric $\Delta_{2}/\Delta_{1}=1$, while in Fig.
16(b), an asymmetric junction with $\Delta_{2}/\Delta_{1}=2.6$ is shown. The
bottom set of panels corresponds to a weak ferromagnet junction with
$h/\Delta_{1}=2.6$. For the nonmagnetic case, Figs. 16(a) and 16(b) show how
the supercurrent with energies above the gap $\Delta_{1}$ get shifted by an
amount corresponding to the gap asymmetry $\Delta_{2}/\Delta_{1}=2.6$. In the
bottom set of panels, the exchange field $h=2.6\Delta_{1}$ introduces an
additional energy scale that redistributes the supercurrent response. In the
symmetric case Fig. 16(c) demonstrates how the presence of magnetism induces a
supercurrent reversal for energies centered around the gap $\Delta_{1}$. When
the junction becomes asymmetric, the mutual effects of the exchange field and
gap asymmetry lead to an enhancement of the supercurrent density at larger
energies around $\Delta_{2}$ [Fig. 16(d)].

## References

  * (1) B. D. Josephson, Possible new effect in superconducting tunneling, Phys. Lett. 1, 251 (1962).
  * (2) P.G. de Gennes, Superconductivity of Metals and Alloys, (Addison-Wesley, reading, Massachussetts, 1989).
  * (3) V. V. Ryazanov, V. A. Oboznov, A. V. Veretennikov, and A. Yu. Rusanov, Intrinsically frustrated superconducting array of superconductor-ferromagnet-superconductor $\pi$ junctions, Phys. Rev. B 65, 020501(R) (2001).
  * (4) V. V. Ryazanov, V. A. Oboznov, A. Yu. 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 (2001).
  * (5) Z. Radovic, N. Lazarides, and N. Flytzanis, Josephson effect in double-barrier superconductor-ferromagnet junctions, Rev. B 68, 014501 (2003).
  * (6) T. Kontos, M. Aprili, J. Lesueur, F. Genet, B. Stephanidis, and R. Boursier, Josephson junction through a thin ferromagnetic layer: negative coupling, Phys. Rev. Lett. 89, 137007 (2002).
  * (7) Ya. V. Fominov, A. F. Volkov, and K. B. Efetov, Josephson effect due to the long-range odd-frequency triplet superconductivity in junctions with Neel domain walls, Phys. Rev. B 75, 104509 (2007).
  * (8) S. Takahashi, S. Hikino, M. Mori, J. Martinek, and S. Maekawa, Supercurrent pumping in Josephson junctions with a half-metallic ferromagnet, Phys. Rev. Lett. 99, 057003 (2007).
  * (9) K. Halterman, O. T. Valls, Layered ferromagnet-superconductor structures:The $0-\pi$ state and proximity effects, Phys. Rev. B 69, 014517 (2004).
  * (10) K. Halterman, O. T. Valls, Proximity effects at ferromagnet-superconductor interfaces, Phys. Rev. B 65, 014509 (2001).
  * (11) V. Shelukhin, A. Tsukernik, M. Karpovski, Y. Blum, K. B. Efetov, A. F. Volkov, T. Champel, M. Eschrig, T. Löfwander, G. Schön, and A. Palevski, Observation of periodic $\pi$-phase shifts in ferromagnet-superconductor multilayers, Phys. Rev. B73, 174506 (2006).
  * (12) T. Karabassov, A. V. Guravova, A. Yu. Kuzin, E. A. Kazakova, S. Kawabata, B. G. Lvov and A. S. Vasenko, Anomalous current-voltage characteristics of SFIFS Josephson junctions with weak ferromagnetic interlayers, Beilstein J.Nanotechnol. 11, 252 (2020).
  * (13) S. Acharjee, U. D. Goswami, Spin transport and Spin Tunnelling Magneto-Resistance (STMR) of F—NCSC—F spin valve, J. Mag. Mag. Mat. 495, 165844 (2020).
  * (14) Z. Shomali, M. Zareyan, and W. Belzig, Spin supercurrent in Josephson contacts with noncollinear ferromagnets, New J. Phys. 13, 083033 (2011).
  * (15) S. Hikino, Magnetization reversal by tuning Rashba spin–orbit interaction and Josephson phase in a ferromagnetic Josephson junction, J. Phys. Soc. Jpn. 87, 074707 (2018).
  * (16) F. Setiawan, A. Stern, and E. Berg, Topological superconductivity in planar Josephson junctions: Narrowing down to the nanowire limit, Phys. Rev. B 99, 220506(R) (2019).
  * (17) M. Alidoust, and K. Halterman, Proximity Induced Vortices and Long-Range Triplet Supercurrents in Ferromagnetic Josephson Junctions and Spin Valves, J. Appl. Phys. 117, 123906 (2015).
  * (18) E. Moen and O. T. Valls, Quasiparticle conductance in spin valve Josephson structures, Phys. Rev. B 101, 184522 (2020).
  * (19) H. Chakraborti, S. Deb, R. Schott, V. Thakur, A. Chatterjee, S. Yadav, R. K. Saroj, A. Wieck, S. M. Shivaprasad, K. Das Gupta and S. Dhar, Coherent transmission of superconducting carriers through a $\sim 2\mu m$ polar semiconductor, Supercond. Sci. Technol. 31 085007 (2018).
  * (20) T. Karabassov, V. S. Stolyarov, A. A. Golubov, V. M. Silkin, V. M. Bayazitov, B. G. Lvov, and A. S. Vasenko, Competitive 0 and $\pi$ states in S/F/S trilayers: Multimode approach, Phys. Rev. B 100, 104502 (2019).
  * (21) S. V. Bakurskiy, A . A. Neilo, N. V. Klenov, I. I. Soloviev and M. Yu. Kupriyanov, Dynamic properties of asymmetric double Josephson junction stack with quasiparticle imbalance, Nanotechnology 30, 324004 (2019).
  * (22) Z. Shomali, R. Asgari, Spin transfer torque and exchange coupling in Josephson junctions with ferromagnetic superconductor reservoirs, J. Phys.: Condens. Matter 32, 035806 (2020).
  * (23) M. V. Avdeev and Yu. N. Proshin, Long-range spin-singlet proximity effect for a Josephson system with a single-crystal ferromagnet due to its band-structure features, Phys. Rev. B 97, 100502(R) (2018).
  * (24) E. Moen and O. T. Valls, Spin-split conductance and subgap peak in ferromagnet/superconductor spin valve heterostructures, Phys. Rev. B 98, 104512 (2018).
  * (25) K. Ohnishi, S. Komori, G. Yang, K.-R. Jeon, L. A. B. Olde Olthof, X. Montiel, M. G. Blamire, and J. W. A. Robinson, Spin-transport in superconductors, Appl. Phys. Lett. 116, 130501 (2020).
  * (26) K. Kulikov, D. Sinha, Yu. M. Shukrinov, and K. Sengupta, Josephson junctions of Weyl and multi-Weyl semimetals, Phys. Rev. B 101, 075110 (2020).
  * (27) T. Vezin, C. Shen, J. E. Han, and I. Zutic, Enhanced spin-triplet pairing in magnetic junctions with s-wave superconductors, Phys. Rev. B 101, 014515 (2020).
  * (28) A. Costa and J. Fabian, Anomalous Josephson Hall effect charge and transverse spin currents in superconductor/ferromagnetic-insulator/superconductor junctionsPhys. Rev. B 101, 104508 (2020).
  * (29) A.A. Mazanik, I.R. Rahmonov, A.E. Botha, and Yu.M. Shukrinov, Analytical Criteria for Magnetization Reversal in a $\varphi_{0}$ Josephson Junction, Phys. Rev. Applied 14, 014003 (2020).
  * (30) H. Meng, Y. Ren, J. E. Villegas, and A. I. Buzdin, Josephson current through a ferromagnetic bilayer: Beyond the quasiclassical approximation, Phys. Rev. B 100, 224514 (2019).
  * (31) C. W. J. Beenakker, Universal Limit of Critical-Current Fluctuations in Mesoscopic Josephson Junctions, Phys. Rev. Lett. 67, 3836 (1991).
  * (32) E. Koshina and V. Krivoruchko, Spin polarization and $\pi$ phase state of the Josephson contact: Critical current of mesoscopic SFIFS and SFIS junctions, Phys. Rev. B 63, 224515 (2001).
  * (33) A. A. Golubov, M. Yu. Kupriyanov, and Ya. V. Fominov, Critical current in SFIFS junctions, Pis’ma Zh. Eksp. Teor. Fiz. 75, 223 (2002) [JETP Lett. 75, 190 (2002)].
  * (34) K. Halterman and O.T. Valls, Proximity effects and characteristic lengths in ferromagnet-superconductor structures, Phys. Rev. B 66, 224516 (2002).
  * (35) M. Alidoust, Critical supercurrent and $\varphi_{0}$ state for probing a persistent spin helix, Phys. Rev. B 101, 155123 (2020).
  * (36) K. Halterman, O. T. Valls, and C.-T. Wu, Charge and spin currents in ferromagnetic Josephson junctions, Phys. Rev. B 92, 174516 (2015).
  * (37) C.-T. Wu and K. Halterman, Spin transport in half-metallic ferromagnet-superconductor junctions, Phys. Rev. B 98, 054518 (2018).
  * (38) G. Tkachov, Topological Insulators, Taylor and Francis, Abingdon, 2016.
  * (39) A. V. Galaktionov, Supersymmetric Hamiltonian solutions simulated by Andreev bound states, Phys. Rev. B 101, 134501 (2020).
  * (40) A. G. Golubov, M. Yu. Kupriyanov, and E. llichev, The current-phase relation in Josephson junctions, Rev. Mod. Phys. 76, 411 (2004).
  * (41) A. Iovan and V. M. Krasnov, Signatures of the spin-triplet current in a Josephson spin valve: A micromagnetic analysis, Phys. Rev. B 96, 014511 (2017).
  * (42) O. M. Kapran, A. Iovan, T. Golod, V. M. Krasnov, Observation of the dominant spin-triplet supercurrent in Josephson spin valves with strong Ni ferromagnets, Phys. Rev. Research 2, 013167 (2020).
  * (43) L. R. Tagirov, M. Yu. Kupriyanov, V. N. Kushnir, and A. Sidorenko, Superconducting Triplet Proximity and Josephson Spin Valves, Functional Nanostructures and Metamaterials for Superconducting Spintronics. NanoScience and Technology. Springer, Cham (2018).
  * (44) D. Culcer, A. C. Keser, Y. Li, G. Tkachov, Transport in two-dimensional topological materials: recent developments in experiment and theory, 2D Materials 7, 022007 (2020).
  * (45) M. Alidoust, K. Halterman, and O.T. Valls, Zero-Energy Peak and Triplet Correlations in Nanoscale SFF Spin-Valves, Phys. Rev. B 92, 014508 (2015).
  * (46) K. Halterman and M. Alidoust, Half-Metallic Superconducting Triplet Spin Valve, Phys. Rev. B 94, 064503 (2016).
  * (47) A. Srivastava, L. A. B. Olde Olthof, A. Di Bernardo, S. Komori, M. Amado, C. Palomares-Garcia, M. Alidoust, K. Halterman, M. G. Blamire, J. W. A. Robinsonn, Magnetization-control and transfer of spin-polarized Cooper pairs into a half-metal manganite, Phys. Rev. Applied 8, 044008 (2017).
  * (48) M. Alidoust and K. Halterman, Half-Metallic Superconducting Triplet Spin Multivalve, Phys. Rev. B 97, 064517 (2018).
  * (49) A. A. Abrikosov, L. P. Gorkov, I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics, (Edited by A. Silverman, Dover Publications, 1963).
  * (50) G. Eilenberger, Transformation of Gorkov’s equation for type-II superconductors into transport-like equations, Z. Phys. 214, 195 (1968).
  * (51) K. D. Usadel, Generalized Diffusion Equation for Superconducting Alloys, Phys. Rev. Lett. 25, 507 (1970).
  * (52) M. Alidoust and K. Halterman, Spontaneous edge accumulation of spin currents in finite-size two-dimensional diffusive spin–orbit coupled SFS heterostructures, New J. Phys. 17, 033001 (2015).
  * (53) M. Alidoust and K. Halterman, Long-range spin-triplet correlations and edge spin currents in diffusive spin–orbit coupled SNS hybrids with a single spin-active interface, J. Phys: Cond. Matt. 27, 235301 (2015).
  * (54) A. Zyuzin, M. Alidoust, D. Loss, Josephson junction through a disordered topological insulator with helical magnetization, Phys. Rev. B 93, 214502 (2016).
  * (55) M. Alidoust and H. Hamzehpour, Spontaneous supercurrent and $\varphi_{0}$ phase shift parallel to magnetized topological insulator interfaces, Phys. Rev. B 96, 165422 (2017).
  * (56) M. Alidoust, Self-biased current, magnetic interference response, and superconducting vortices in tilted Weyl semimetals with disorder, Phys. Rev. B 98, 245418 (2018).
  * (57) M. Alidoust, M. Willatzen, A.-P. Jauho, Fraunhofer response and supercurrent spin switching in black phosphorus with strain and disorder, Phys. Rev. B 98, 184505 (2018).
  * (58) O. Šipr and B.L. Györffy, J. Phys.: Condens. Matter 7, 5239 (1995).
  * (59) L. F. Chang, P. F. Bagwell, Ballistic Josephson-current flow through an asymmetric superconductor/normal-metal/superconductor junction, Phys. Rev. B 49, 15853 (1994).
  * (60) K. Halterman, M. Alidoust, Josephson currents and spin-transfer torques in ballistic SFSFS nanojunctions, Supercond. Sci. Technol. 29, 055007 (2016).
  * (61) M. Alidoust, M. Willatzen, A.-P. Jauho, Strain-engineered Majorana zero energy modes and $\varphi_{0}$ Josephson state in black phosphorus, Phys. Rev. B 98, 085414 (2018).
  * (62) M. Alidoust and K. Halterman, Evolution of Pair Correlation Symmetries and Supercurrent Reversal in Tilted Weyl Semimetals, Phys. Rev. B 101, 035120 (2020).
  * (63) A.I. Buzdin, L.N. Bulaevskii, 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, Pis’ma Zh. Eksp. Teor. Fiz. 35, 147 (1982) [JETP Lett. 35, 178 (1982)].
  * (64) A. V. Zaitsev, Sov. Phys. JETP 59, 1015 (1984); M. Y. Kuprianov and V. F. Lukichev, Sov. Phys. JETP 67, 1163 (1988).
  * (65) K. Halterman and O.T. Valls, Energy gap of ferromagnet-superconductor bilayers, Physica C 397 151-15 (2003).
  * (66) K. Halterman and O.T. Valls, Local density of states and order parameter configurations in layered ferromagnet-superconductor structures, Physica C 420, 111 (2005).
  * (67) K. Halterman and O.T. Valls, Stability of $\pi$-junction configurations in ferromagnet-superconductor heterostructures Phys. Rev. B 70, 104516 (2004).