# The role of Rashba spin-orbit induced spin textures in the anomalous
Josephson effect.

R. D. Monaghan<EMAIL_ADDRESS>School of Physical Sciences, The
University of Adelaide  School of Chemical and Advanced Engineering, The
University of Adelaide  G. C. Tettamanzi  School of Physical Sciences, The
University of Adelaide  School of Chemical and Advanced Engineering, The
University of Adelaide

###### Abstract

This work reports the theoretical investigation into the mechanism
underpinning the anomalous Josephson effect. The prototypical system we study
is a ballistic two-dimensional junction containing a two-dimensional Rashba
spin-orbit interaction. In this paper we demonstrate how this two-dimensional
Rashba interaction mixes the spins of adjacent transverse subbands which leads
to significant spin-asymmetry within the junction. Under an external magnetic
field, applied perpendicular to both the axis of transport and the normal
vector of the junction, the sinusoidal Josephson current can then experience
an anomalous phase shift. The role of this spin mixing in the limit of a
single sub-band is initially explored by deriving an analytical expression for
the resulting anomalous phase shift. The analysis is then extended to systems
with multiple occupied sub-bands; in this later section, starting from a
microscopic model, we derive an analytic formula for the resulting anomalous
phase shift indicating it is linear in both magnetic field and spin-orbit
strength. We then verify and validate all findings by comparing them with
numerical results evaluated by a tight-binding model.

##  I Introduction

The Josephson effect is one of the most profound results concerning the phase
of quantum objects; for a phase difference of $\varphi$ between two coupled
superconductors, separated by a region of non-superconducting material, a
sinusoidal electrical current will flow [1]. A long-standing result in the
field of Josephson junctions was that applying an appropriately large external
magnetic field to the junction could induce a $\pi$ phase shift in the
sinusoidal current, where the minimum of the free energy of the system now
occurs when $\varphi=\pi$ [2]. A natural extension of this concept is to
construct a junction where the Josephson current acquires an anomalous phase
factor $\varphi_{0}$, with $0\leq\varphi_{0}\leq 2\pi$, such that the total
current as a function of the phase difference $\varphi$ across the
superconductors is given by

$I=I_{c}\sin(\varphi+\varphi_{0})\;.$ (1)

To generate such an arbitrary phase shift it was realised that, regardless of
the physical mechanism, both time-reversal and spatial-inversion symmetry of
the system must be broken [3]. One way to do so is by constructing Josephson
junctions where the non-superconducting section of the junction is made from
material possessing a strong Rashba spin-orbit signature, and by also applying
an external magnetic field perpendicular to both the axis of transport and the
normal of the substrate. Careful experimental work over the recent years has
indeed detected this anomalous phase in systems possessing both Rashba spin-
orbit coupling and an external magnetic field [4, 5, 6, 7]. This effect is now
known as the anomalous Josephson effect. However, no clear and unified
explanation has been proposed regarding the microscopic nature of this effect
[8].

In quasi one-dimensional devices, such as the ones studied in this work, the
physical confinement can give rise to multiple transverse subbands. A well
known result is that under a two-dimensional Rashba spin-orbit interaction,
these transverse subbands can be coupled and the spins of the states mixed [9,
10]. The resulting mixing of the spin states – which we denote as the spin-
texture of the system – has been theoretically investigated previously within
the context of non-superconducting transport [11, 12, 13, 14, 15]. This mixing
has also recently been observed in experimental spectroscopy measurements on
two-dimensional Josephson junction [16]. The purpose of this work is to
explicitly investigate the effects of the spin-texture on the resulting
anomalous phase – deriving expressions for the anomalous phase as a function
of the spin-texture and, finally, to provide a more complete microscopic model
that can justify the appearance of this effect in different systems. Although
previous work has considered the role of subband mixing within the framework
of the anomalous Josephson effect (see Refs [17, 14, 15, 18]), these papers
either only consider the Fermi velocity asymmetry induced by the subband
mixing, or as in the case of Ref [14], discuss how the mixing of the spins can
act as a spin-polariser to generate the anomalous phase.

In this work, we initially study single-channel systems with two-dimensional
Rashba interactions such that we can generate analytic expressions to clearly
demonstrate the physics. The results are then extended to multi-channel
systems where it is possible, after some approximations, to make experimental
predictions; to verify these predictions, the anomalous phase shift is then
explicitly computed numerically through the use of the non-equilibrium Green
function formalism.

##  II Model for the anomalous Josephson effect

Within the literature on the anomalous Josephson effect there are currently
two leading microscopic mechanisms which can be used to explain the observed
experimental results. The first is spin-related interference effects induced
by the combination of Zeeman and Rashba spin-orbit coupling which were first
introduced by Krive et al. in Ref. [17]. However, another mechanism proposed
is orbital-related interference effects induced by the vector potential and
disorder within the system as in Refs [19, 20]. Indeed, based purely on
symmetry arguments, Rasmussen et al. showed that either spin-interference, or
orbital-interference, could induce a non-zero supercurrent at zero phase
difference [3]. In this work, our focus will be purely on the role of spin-
channel interference and, as a result, both disorder and vector potential will
be ignored within our model. Strictly speaking, this restricts our analysis to
systems where the total magnetic flux piercing the non-superconducting region
is much smaller than a flux quantum $\Phi_{0}$ such that any phase gained by
the vector potential is negligible. Similarly, to ignore disorder we only
consider systems which are purely in the ballistic limit.

###  II.1 Single particle Hamiltonian

Figure 1: Diagram of the prototypical system studied throughout this paper. It
consists of a two-dimensional material which, due to the confinement $W$ in
the $y$-axis is quasi one-dimensional along the $x$-axis. The blue material is
some semiconductor which has a Rashba spin-orbit interaction present, and is
assumed to extend to infinity in the $x$-axis. The red regions are some Type-I
superconductor which has been deposited onto the semiconducting layer such
that the left and right regions of the device contain some proximitised
superconductivity. The uncovered central region of the Josephson junction has
length $L$.

We consider the prototypical Josephson junction with a ballistic non-
superconducting central region sandwiched between two superconducting regions
depicted in Fig. 1. The blue region is some non-superconducting semi-
conducting material which possesses a two-dimensional Rashba spin-orbit
interaction. The red regions in Fig. 1 illustrate some Type-I superconducting
material which has been deposited over the semiconductor to induce semi-
infinite superconducting leads by the proximity effect [21]. Conceptually, in
our model the system consists of three regions: the ‘left’ and ‘right’
proximitised regions are those covered by a superconducting layer and are
assumed to extend to infinity, whilst the non-proximitised uncovered region is
known as the ‘central’ region. The entire device is assumed to be quasi two-
dimensional with a hard-wall boundary conditions confining the system in the
$y$-axis. We also assume that the superconducting coherence length $\xi$ is
larger than the length $L$ of the central region such that our system is in
the short junction limit [22].

We assume that an external magnetic field is applied along the $y$ axis such
that it is perpendicular to both the normal of the substrate and the direction
of transport – this is depicted in Fig. 1. We also assume that the magnetic
flux is entirely screened from the left and right regions due to the deposited
superconducting layer. As a result, the single particle Hamiltonian in the
central non-proximitised region can be written as

$H_{C}=H_{0}+H_{\alpha}+H_{Z}\;,$ (2)

where the free single-particle Hamiltonian is given by

$H_{0}=-\frac{\hbar^{2}\nabla^{2}}{2m^{*}}-\mu\;,$ (3)

and $m^{*}$ denotes the effective mass for the semiconductor being modelled,
whilst $\mu$ is the Fermi level. The two-dimensional Rashba spin-orbit term is
written as [9]

$H_{\alpha}=i\alpha\sigma_{y}\partial_{x}-i\alpha\sigma_{x}\partial_{y}\;,$
(4)

where $\alpha$ is the Rashba spin-orbit strength, and we have used the usual
Pauli matrices to span spin-space. The Zeeman term is written as

$H_{Z}=E_{Z}\sigma_{y}\;,$ (5)

where $E_{Z}=g\mu_{B}B/2$ for some material dependent $g$ factor. We have
ignored any on-site scattering potentials, placing ourselves purely in the
ballistic limit. Furthermore, as our system is two-dimensional the height $H$
of the device in the $z$-axis is small such that the flux penetrating is
negligible, $BLH\ll h/e$. As a result, it is reasonable to ignore the vector
potential [23].

For the left and right regions of the device, the non-superconducting
component of the Hamiltonian is given by

$H_{L/R}=H_{0}+H_{\alpha}\;,$ (6)

where we have omitted the Zeeman term as we assume that, due to the deposited
superconducting layer, the magnetic field has been entirely expelled by the
Meissner effect. To account for the superconducting proximity effect induced
by the deposited superconductors, we introduce a superconducting order
parameter $\Delta$ which acts microscopically to couple the electrons and
holes of the system [24]. As a result, we enlarge our basis such that the
order parameter couples the non-superconducting Hamiltonian with its time
reversed pair. Considering this, the Hamiltonian in the proximitised semi-
conducting leads takes the form

$H_{L/R}\to\begin{pmatrix}H_{0}+H_{\alpha}&i\sigma_{y}\Delta_{L/R}\\\
-i\sigma_{y}\Delta_{L/R}&-(H_{0}+H_{\alpha})^{*}\end{pmatrix}\;.$ (7)

###  II.2 Spin-textures in wavevector space

In a purely one-dimensional system, the Rashba spin-orbit term of Eq. 4 can be
simplified to the form $H_{\alpha}^{1\text{D}}=i\alpha\sigma_{y}\partial_{x}$
[25]. Within this one-dimensional system, the full Hamiltonian commutes with
the $y$ component of the electron’s spin
$\left[H_{0}+H_{\alpha}+H_{Z},\sigma_{y}\right]=0$. If we now include a
transverse dimension, such that the system has a non-zero length $L$ and width
$W$, the confinement in the $y$-axis will produce subbands which can be
labelled by an index characterising the discrete states within an infinite
potential well. Retaining the purely one-dimensional Rashba spin-orbit term,
the Hamiltonian is entirely separable in $x$ and $y$ such that both the spin
and subband indices are good quantum numbers. The resulting dispersion
relation for this system with zero external magnetic field is illustrated in
Fig. 2.a. This has been computed numerically with a tight-binding calculation
described in Ref. [26]. The colours depict the value of the corresponding
numerical eigenstate projected onto the $y$-axis, $\left<\sigma_{y}\right>$,
which we denote as the spin-projection [10]. As this is essentially the
dispersion relation of an electron waveguide with only transport in the
$x$-axis, we refer to the wavevector along the $x$-axis as simply $k$, rather
than $k_{x}$.

Figure 2:  Numerical dispersion relation computed used the tight-binding
techniques of Ref. [26]. The wavevector $k$ denotes the wavevector along the
$x$-axis. The colours of the curves indicate the corresponding numerical
eigenstate projected onto the $y$-axis, $\left<\sigma_{y}\right>$. No magnetic
field is applied to this system. (a) The dispersion relation for a quasi one-
dimensional system with a purely one-dimensional Rashba spin-orbit coupling.
(b) The dispersion relation for a quasi one-dimensional system with a
physically correct two-dimensional Rashba spin-orbit coupling. (c) The linear
states labelled for the energies around the Fermi level.

Importantly, a physically two-dimensional system with a one-dimensional Rashba
spin-orbit interaction is not realistic; to accurately model this system we
need to write the full two-dimensional Rashba spin-orbit interaction as

$H^{2\text{D}}_{\alpha}=i\alpha\sigma_{y}\partial_{x}-i\alpha\sigma_{x}\partial_{y}\;.$
(8)

The immediate result of the extra term containing a $\sigma_{x}$ matrix, is
that $\sigma_{y}$ no longer commutes with the Hamiltonian. This well-known
result leads to the mixing of states of opposite spin and adjacent subband
index [9]. The dispersion relation for this system is sketched in Fig. 2.b. As
spin is no longer a good quantum number, where there was previously a crossing
between odd and even numbered subbands of opposite spin there is now an
avoided crossing where the spin-projection smoothly flips.

To clarify the eventual role of this spin-texture, it is useful to linearise
the dispersion relation around the Fermi wavevectors for each state – the
Fermi wavevector for the $j$th state being denoted by $k_{j}$. As the
transport in the quasi one-dimensional structure is entirely along the
$x$-axis, we can decompose the electronic eigenstates $\ket{\psi_{j\pm}}$ into
a linear combination of left-moving, $\ket{\psi^{l}_{j\pm}}$, and right-
moving, $\ket{\psi^{r}_{j\pm}}$, components [23]

$\ket{\psi_{j\pm}}=e^{-ik_{j}x}\ket{\psi^{l}_{j\pm}}+e^{ik_{j}x}\ket{\psi^{r}_{j\pm}}\;.$
(9)

Although the subband index $j$, and spin index $\pm$ are no longer good
quantum numbers, to clarify the notation we can still denote each state by the
indices it would have in the no Rashba limit $\alpha\to 0$. We sketch the
linearised dispersion relation in Fig. 2.c where we have also labelled each
state by their three indices: subband, spin, and direction. The linearised
dispersion relation for each state is entirely characterised by three
parameters – namely, the Fermi wavevector of each state $k_{j\pm}^{l/r}$
defined such that

$E(k_{j\pm}^{l/r})-\mu=0\;,$ (10)

the Fermi velocity, which is proportional to the slope of the linearised
spectrums

$v_{j\pm}^{l/r}\equiv\frac{1}{\hbar}\left.\frac{\partial E}{\partial
k}\right|_{k=k_{j;\pm}^{l/r}}\;,$ (11)

and the spin-projection which is given by the overlap

$s_{j\pm}^{l/r}\equiv\expectationvalue{\sigma_{y}}{\psi_{j\pm}^{l/r}}\;.$ (12)

Importantly, in the absence of an external magnetic field, the system still
possesses time-reversal symmetry such that [13]

$\displaystyle E(k)$ $\displaystyle=E(-k)\;,$ (13)
$\displaystyle\ket{\psi^{l}_{j\pm}}$
$\displaystyle=\mathcal{T}\ket{\psi^{r}_{j\mp}}\;,$

where $\mathcal{T}$ is the time reversal operator. As a result of this, a pair
of useful identities which are valid in the absence of a magnetic field are
that

$v^{l}_{\pm}=-v^{r}_{\mp}\qquad s^{l}_{\pm}=-s^{r}_{\mp}\;,$ (14)

which can be observed visually in Fig. 2.c.

To now briefly study how the dispersion relation changes under an external
magnetic field applied along the spatial $y$-axis, we restrict ourselves to
systems where the magnetic field is small relative to the subband spacing.
This is not entirely restrictive, as for say a generic InAs device with a
width of $200$nm, this corresponds to magnetic fields less than $400$mT, which
is certainly within the range of magnetic fields commonly applied to these
systems [5, 7, 4, 6]. For these small magnetic fields, it is reasonable to
assume that the spin-texture and Fermi velocities are unchanged as they only
vary appreciably for energies on the scale of the subband spacing [13].
Furthermore, as the Zeeman term only contains a $\sigma_{y}$ matrix, then the
Zeeman term expressed in the basis of the linearised states
$\ket{\psi^{l/r}_{j\pm}}$ is given simply by $H_{Z}=s^{l/r}_{j\pm}E_{Z}$,
where we have substituted in the definition of spin-projection defined in Eq.
12. As a result, we can express the linearised energy spectrums as

$\varepsilon_{j\pm}^{l/r}(k)=\hbar v_{j\pm}^{l/r}k-s_{j\pm}^{l/r}E_{Z}\;.$
(15)

###  II.3 Anomalous phase induced by Andreev bound states

As introduced in Eq. 7, the order parameter acts to couple the electrons and
holes within the superconducting system. As a result of this, an electron in
the non-proximitised region of the Josephson junction which is incident on the
superconducting interface will be Andreev reflected as its time reversed
partner [27, 28]. For a Josephson junction, where the non-proximitised region
is sandwiched between two superconductors, this leads to the formation of
bound states consisting of counter-propagating electrons and holes – an object
known as an Andreev bound state [29, 30]. This process is illustrated in Fig.
3.a. It is important to remember that in this work we focus solely on the
short junction limit where the superconducting coherence length $\xi$ is much
greater than the length $L$ of the junction; in this limit the supercurrent is
entirely carried by the Andreev bound states rather than states outside of the
superconducting gap [31].

Figure 3: (a) The diagram of an Andreev reflection, where a right-moving
electron is Andreev reflected as a hole below the Fermi level. With no terms
acting on spin, the linearised dispersion relation is spin degenerate. (b)
Under a magnetic field, and with no spin-orbit coupling, the dispersion
relation is Zeeman split such that two Andreev bound states appear. Although
the resulting Andreev bound states are phase shifted, due to the opposite
spins the two phases cancel each other out so that the total system gains no
anomalous phase – illustrated by both Andreev bound states having the same
path length. (c) Under both a two-dimensional spin-orbit term and a magnetic
field, the asymmetry in the spin-projections will induce an anomalous phase –
illustrated by the different path lengths. The colours of the states
correspond to the magnitude of the spin-projection; for reference see the
colourbar in Fig. 2.

The Andreev bound states which carry the Josephson current are entirely formed
by the states within the non-proximitised central region. Hence, having
discussed the microscopic spin-texture of the states within the non-
proximitised central region in Sec. II.2, we can now discuss the ramifications
of the spin-texture on the resulting Andreev bound states. To begin, previous
literature has shown that, generically, to observe the anomalous Josephson
effect you need to introduce some asymmetry between different Andreev bound
states of the system [17, 32]. For example, if we consider the case of no
Rashba spin-orbit coupling, but a non-zero external magnetic field –
illustrated in Fig. 3.b – although we see the spin-bands shift, no asymmetry
is generated. As a result, we expect no anomalous effect. However, if we now
include the two-dimensional Rashba spin-orbit coupling, in conjunction with
the magnetic field along the $y$-axis – illustrated in Fig. 3.c – the
resulting Fermi velocity asymmetry and spin-texture, mixed with the external
magnetic field, induces significant asymmetry in the resulting pair of Andreev
bound states. Although previous theoretical work has studied this asymmetry,
they ignored the role of the spin-texture focusing instead only on the Fermi
velocities [17, 18]. The main focus of this work is to provide both a
quantitative and qualitative link between the spin-texture asymmetry induced
by the two-dimensional Rashba spin-orbit coupling and the anomalous Josephson
effect.

To understand the consequences of this asymmetry in our system, we now look to
evaluate the energy of the different Andreev bound states. A simplification of
working in the short junction limit is that each linearised state
$\psi^{l/r}_{j\pm}$ will give rise to only a single Andreev bound state [28].
Indeed, an Andreev bound state with energy $\varepsilon$ can only be formed by
matching the wavefunctions of both the linearised state $\psi^{l/r}_{j\pm}$
with energy $\varepsilon$ above the Fermi level, to both its time-reversed
hole partner, and to the exponentially decaying wavefunctions within the
superconducting leads [28]. The details are left for Appendix. A, however,
doing so leads to the following transcendental equation which the energy
$\varepsilon$ of the resulting Andreev bound state must satisfy

$e^{-i\varphi}e^{-2i\text{arccos}(\varepsilon/|\Delta|)}e^{2iL\varepsilon/\hbar
v_{j\pm}^{l/r}}e^{2iLs_{j\pm}^{l/r}E_{Z}/\hbar v_{j\pm}^{l/r}}=1\;.$ (16)

Taking the log of both sides we find a quantisation condition for the phase
acquired by the counter-propagating electrons and holes forming the Andreev
bound state

$\varphi+2\text{arccos}\left(\frac{\varepsilon}{|\Delta|}\right)-\frac{2L\varepsilon}{\hbar
v_{j\pm}^{l/r}}-\frac{2Ls_{j\pm}^{l/r}E_{Z}}{\hbar v_{j\pm}^{l/r}}=2\pi n\;,$
(17)

for some integer $n$. The term $\varphi$ is the superconducting phase
difference picked up during the Andreev reflection process, the
$2\text{arccos}(\varepsilon/|\Delta|)$ term is the phase picked up from the
evanescent states in the superconductor, whilst $2L\varepsilon/\hbar
v_{j\pm}^{l/r}$ is the phase picked up by traversing the junction of length
$L$ [33]. The interesting term is ${2Ls_{j\pm}^{l/r}E_{Z}}/{\hbar
v_{j\pm}^{l/r}}$ which is the phase picked up by the spin moving through the
magnetic field.

As we are working in the short junction limit such that $\xi_{j\pm}^{l/r}\gg
L$, where $\xi_{j\pm}^{l/r}\equiv\hbar v_{j\pm}^{l/r}/\Delta$ is the coherence
length of that subband state, then the phase picked up by traversing the
junction is negligible

$\frac{L\varepsilon}{\hbar
v_{j\pm}^{l/r}}=\frac{L}{\xi_{j\pm}^{l/r}}\frac{\varepsilon}{|\Delta|}\approx
0\;.$ (18)

As a result, solving for the allowed energies $\varepsilon$ in Eq. 17, we find
that the energy phase relation for the Andreev bound state produced via the
$\psi^{l/r}_{j\pm}$ subband state is given by

$\displaystyle\varepsilon^{r}_{j\pm}(\varphi)$
$\displaystyle=-|\Delta|\cos(\frac{\varphi}{2}-E_{Z}L\frac{s_{j\pm}^{r}}{\hbar
v_{j\pm}^{r}})\;,$ (19) $\displaystyle\varepsilon^{l}_{j\pm}(\varphi)$
$\displaystyle=|\Delta|\cos(\frac{\varphi}{2}-E_{Z}L\frac{s_{j\pm}^{l}}{\hbar
v_{j\pm}^{l}})\;.$

From these bound state energies, we can now directly compute the Josephson
current; in the short junction limit the current is entirely given by the
contribution through each discrete energy level [31]

$I=-\frac{2e}{\hbar}\sum_{j,l/r,\pm}\left[\tanh(\frac{\varepsilon_{j\pm}^{l/r}}{2k_{B}T})\partialderivative{\varepsilon_{j\pm}^{l/r}}{\varphi}\right]\;,$
(20)

where $k_{B}T$ is the thermal energy of the system. Utilising the time-
reversal symmetry discussed in Eq. 13 then we can show that

$\varepsilon_{j\pm}^{r}=-\varepsilon_{j\mp}^{l}\;,$ (21)

such that Eq. 20 can be rewritten to include only right-moving states as

$\displaystyle I=\frac{e|\Delta|^{2}}{\hbar
k_{B}T}\sum_{j}\left[\cos(\frac{\varphi}{2}-E_{Z}L\frac{s_{j+}^{r}}{\hbar
v_{j+}^{r}})\sin(\frac{\varphi}{2}-E_{Z}L\frac{s_{j+}^{r}}{\hbar
v_{j+}^{r}})\right]$ (22) $\displaystyle+\frac{e|\Delta|^{2}}{\hbar
k_{B}T}\sum_{j}\left[\cos(\frac{\varphi}{2}-E_{Z}L\frac{s_{j-}^{r}}{\hbar
v_{j-}^{r}})\sin(\frac{\varphi}{2}-E_{Z}L\frac{s_{j-}^{r}}{\hbar
v_{j-}^{r}})\right]$

where we have assumed that we are working near the critical temperature
($k_{B}T\sim|\Delta|$) to simplify the analytical expression. This is not an
essential approximation; however, it allows us to recover a sinusoidal
Josephson current even in the ballistic regime [22]. Using a double angle
trigonometric identity this is simplified further to

$\displaystyle I=\frac{e|\Delta|^{2}}{2\hbar
k_{B}T}\sum_{j}\left[\sin(\varphi-2E_{Z}L\frac{s_{j+}^{r}}{\hbar
v_{j+}^{r}})\right.$ (23)
$\displaystyle+\left.\sin(\varphi-2E_{Z}L\frac{s_{j-}^{r}}{\hbar
v_{j-}^{r}})\right]\;.$

Writing the current as a sum of sinusoids, each with its own phase offset, as
in Eq. 23, allows us to explicitly evaluate the anomalous phase $\varphi_{0}$
by phasor addition. The resulting current can be expressed as

$I=I_{C}\sin(\varphi+\varphi_{0})\;,$ (24)

where $I_{C}$ is the critical current and the anomalous phase is given by

$\varphi_{0}=\atan\left[\frac{\sum_{j}\left(\sin\frac{2E_{Z}Ls_{j+}^{r}}{\hbar
v_{j+}^{r}}+\sin\frac{2E_{Z}Ls_{j-}^{r}}{\hbar
v_{j-}^{r}}\right)}{\sum_{j}\left(\cos\frac{2E_{Z}Ls_{j+}^{r}}{\hbar
v_{j+}^{r}}+\cos\frac{2E_{Z}Ls_{j-}^{r}}{\hbar v_{j-}^{r}}\right)}\right]\;.$
(25)

We can identify the phase factor ${2E_{Z}Ls_{j\pm}^{r}}/{\hbar v_{j\pm}^{r}}$
in this expression from the quantisation condition of Eq. 17 as the phase
gained by right-moving electrons with spin-index $\pm$ as they pass through
the magnetic field. As we can observe from this expression, the presence of an
anomalous phase is solely due to asymmetry in the Fermi velocities and the
spin-projections – namely, there will be an anomalous phase when
$v^{r}_{j+}\neq v^{r}_{j-}$ or when $s^{r}_{j+}\neq-s^{r}_{j-}$.

##  III Anomalous phase in single-channel systems

In this section the anomalous phase shift is evaluated for systems with only a
single occupied subband. The benefit of this system is due to only a pair of
spin states being occupied, the resulting physics is simplified and can
generally be carried out analytically. We consider a Josephson junction with
length $L$ and width $W$ as depicted in Fig. 1. We utilise a two-dimensional
Rashba spin-orbit term as given by Eq. 8.

###  III.1 Effective 1D Hamiltonian

We begin with the same Hamiltonian as outlined in Sec II.1, however, we assume
the Fermi level $\mu$ is such that only a single transverse mode, and hence
two spin channels, are occupied. As this system is only quasi one-dimensional,
rather than truly one-dimensional, it is difficult to extract analytical
results. As a result, we look to integrate out the transverse dimension $y$ to
construct an effective one-dimensional Hamiltonian $H^{\text{eff}}$ which can
be analytically diagonalised. To this end, and following the work of Ref.
[13], the Hamiltonian is first partitioned into a free term $H_{0}$, and a
perturbing term $V$

$H_{C}=H_{C}^{0}+V\;,$ (26)

where the free term is given by

$H_{C}^{0}=-\frac{\hbar^{2}\nabla^{2}}{2m^{*}}-\mu+i\alpha\sigma_{y}\partial_{x}+E_{Z}\sigma_{y}\;,$
(27)

and the perturbation is given by

$V=-i\alpha\sigma_{x}\partial_{y}\;.$ (28)

As $H_{C}^{0}$ is both separable in $x$, $y$, and commutes with $\sigma_{y}$,
then we can write the energy eigenstates as

$\Phi_{j\pm}(x,y)=e^{ik_{x}x}\phi_{j}(y)\chi_{\pm}\;,$ (29)

where

$\phi_{j\pm}(y)=\sqrt{\frac{2}{W}}\sin\left(\frac{j\pi y}{W}\right)\;,$ (30)

are orthogonal eigenstates of $-\hbar^{2}\partial_{y}^{2}/(2m^{*})$ that
satisfy the hard-wall boundary conditions the confinement along the $y$-axis
imposes. The eigenenergies of these transverse states are given by

$E_{j}=\frac{\hbar^{2}\pi^{2}j^{2}}{2mW^{2}}\;.$ (31)

To integrate out the $y$-dimension, a projection operator $P_{0}$ that
projects onto the basis of $\phi_{j}$ is utilised. It is defined such that it
takes a generic function $f(x,y)$ and expresses it in the basis of $\phi_{j}$:

$P_{0}f(x,y)\equiv\sum_{i}\phi_{i}(y)\left[\int\phi_{i}^{*}(y)f(x,y)\;\differential{y}\right]\;.$
(32)

With the use of this projection operator, to first order in the perturbation
the effective Hamiltonian is given by [34]

$H^{\text{eff}}\equiv P_{0}HP_{0}\;.$ (33)

Writing the resulting eigenvalue equation as

$H^{\text{eff}}\Psi_{\pm}(x,y)=\varepsilon\Psi_{\pm}(x,y)\;,$ (34)

we can left multiply by $\phi_{i}^{*}$ and integrate over $y$ to obtain the
matrix equation

$\left<H^{\text{eff}}\right>_{ij}\psi_{j\pm}(x)=\varepsilon\psi_{i\pm}(x)\;,$
(35)

where

$\left<H^{\text{eff}}\right>_{ij}\equiv\int\phi_{i}^{*}H^{\text{eff}}\phi_{j}\differential
y\;,$ (36)

and

$\psi_{i\pm}(x)\equiv\int\phi_{i}^{*}\Psi_{\pm}(x,y)\differential y\;.$ (37)

As only the lowest transverse subband is occupied, we can truncate the
expansion of Eq. 35 to only include the first two transverse subbands. In
doing so, the only term which mixes the transverse subbands in Eq. 36 is given
by

$\displaystyle\left<H^{\text{eff}}\right>_{12}=-\left<H^{\text{eff}}\right>_{21}$
$\displaystyle=-i\alpha\int\phi_{1}^{*}\sigma_{x}\partial_{y}\phi_{2}\differential
y$ (38) $\displaystyle=i\sigma_{x}\frac{8\alpha}{3W}$ $\displaystyle\equiv
i\sigma_{x}\eta$

where $\eta\equiv 8\alpha/(3W)$ is the geometry dependent measure of the
coupling between different transverse subbands of opposite spin. The diagonal
terms are also given by

$\displaystyle\left<H^{\text{eff}}\right>_{11}$
$\displaystyle=\frac{\hbar^{2}k^{2}}{2m^{*}}+E_{1}-\mu-k\alpha\sigma_{y}+E_{Z}\sigma_{y}\;,$
(39) $\displaystyle\left<H^{\text{eff}}\right>_{22}$
$\displaystyle=\frac{\hbar^{2}k^{2}}{2m^{*}}+E_{2}-\mu-k\alpha\sigma_{y}+E_{Z}\sigma_{y}\;,$

where again we are dropping subscripts to simply denote the momentum along the
$x$-axis by $k$. Focusing on the lowest energy solution – as only a single
subband is assumed to be occupied – the energy of the ground state
$\varepsilon_{\pm}$ can be solved from Eq. 35. In doing so, two bands are
found with energy given by

$\varepsilon_{\pm}(k)=\frac{\hbar^{2}k^{2}}{2m^{*}}+\frac{5E_{1}}{2}-\mu+\sqrt{\left(\frac{3E_{1}}{2}\pm
E_{Z}\pm\alpha k\right)^{2}+\eta^{2}}$ (40)

We plot an example dispersion relation using Eq. 40 in Fig. 4. Using Eq. 11
and Eq. 12, the spin-projections and the group velocity as a function of the
wavevector are given by [13]

$\displaystyle s_{1\pm}(k)$ $\displaystyle=\frac{\frac{3E_{1}}{2}\pm\alpha
k}{\sqrt{\left(\frac{3E_{1}}{2}\pm\alpha k\right)^{2}+\eta^{2}}}\;,$ (41)
$\displaystyle v_{1\pm}(k)$ $\displaystyle=\frac{\hbar
k}{m^{*}}+\frac{\alpha}{\hbar}\frac{\frac{3E_{1}}{2}\pm\alpha
k}{\sqrt{\left(\frac{3E_{1}}{2}\pm\alpha k\right)^{2}+\eta^{2}}}\;.$

The spin-projections $s_{1\pm}(k)$ are also plotted in Fig. 4. A useful
diagnostic is that by setting $\eta=0$, such that we are ignoring the overlap
between adjacent transverse subbands, then Eq. 41 states that $s_{1\pm}=\pm 1$
as expected.

Figure 4: Analytic dispersion relation computed using Eq. 40 for a two-
dimensional material with width $W=300\,$nm, spin-orbit strength
$\alpha=3\times 10^{-11}$\text{\,}\mathrm{e}\mathrm{V}\mathrm{m}$$, and no
external magnetic field $E_{Z}=0$. The dashed grey line is when $\eta=0$ such
that there is no coupling between adjacent transverse subbands – we do not
colour these lines as to simplify the figure. The qualitative similarity of
this figure to the numerical bandstructure of Fig. 2.b indicates the quality
of the analytic dispersion relation.

###  III.2 Anomalous phase

The general expression for the anomalous Josephson phase given by Eq. 25 is
written as a sum over all of occupied subbands; in systems with only a single
occupied subband it can be simplified to

$\varphi_{0}=E_{Z}L\left(\frac{s^{r}_{1+}}{\hbar
v^{r}_{1+}}+\frac{s^{r}_{1-}}{\hbar v^{r}_{1-}}\right)\;.$ (42)

If the asymmetry of the spin-texture is now ignored such that we naïvely
simply set $s^{r}_{1+}=1$ and $s^{r}_{1-}=-1$, then the anomalous phase is
solely a function of the Fermi velocity asymmetry

$\varphi_{0}\underset{\begin{subarray}{c}s^{r}_{1+}=+1\\\
s^{r}_{1-}=-1\end{subarray}}{\longrightarrow}E_{Z}L\left(\frac{1}{\hbar
v^{r}_{1+}}-\frac{1}{\hbar v^{r}_{1-}}\right)\;.$ (43)

This equation is not new, having appeared in the literature previously to
model the anomalous phase in quasi one-dimensional Josephson junctions [17,
18]. It is a good validation tool to note that ignoring the spin-texture our
equation resolves to the standard result in the literature.

To highlight the explicit role of the spin-texture, we will now analytically
solve for the anomalous phase using both Eq. 42, which contains the spin-
texture asymmetry, and the conventional Eq. 43, which ignores the spin-
texture. To this end, we consider a Josephson junction with a length of
$200\,\text{nm}$, a width of $20\,\text{nm}$, and a magnetic field is applied
along the $y$-axis with a Zeeman strength of $100\mu\text{eV}$. To compute the
anomalous phase, the spin-projection and group velocities of the states at the
Fermi level must be determined. To determine the Fermi wavevectors $k_{1\pm}$,
which can then be substituted into Eq. 41 to determine the spin-projection and
Fermi velocities, the roots of the effective dispersion relation derived in
Eq. 40 must be found. Although this is equivalent to solving for the roots of
a quartic, and hence has an analytic form, it is not illuminating to write
down the resulting expressions as they are particularly cumbersome.
Regardless, substituting in these roots into Eq. 41 provides an analytic form
for the Fermi velocity and spin-projection at the Fermi level. Substituting
these into Eq. 42 and Eq. 43 will then provide analytic forms for the
anomalous phase. Sweeping both the strength of the Rashba spin-orbit strength
$\alpha$ and the Fermi level $\mu$, the results are plotted in Fig. 5.a and
Fig. 5.b respectively.

Figure 5: The anomalous phase in the Josephson junction as a function of the
spin-orbit coupling strength. As the effective one-dimensional model for the
dispersion derived in Sec.III.1 is only correct to first order in the subband
coupling, the spin-orbit strength is kept such that the energy within the
spin-orbit energy scale does not exceed the subband energy scale
$m^{*}\alpha^{2}/\hbar^{2}<E_{1}$. (a) The case of two-dimensional Rashba
spin-orbit coupling with subband mixing. (b) The case ignoring the spin-
texture, and considering only Fermi velocity asymmetry. (c) The asymmetry
parameter computed using Eq. 44. The Fermi level $\mu$ has been fixed to the
value $(E_{1}+E_{2})/2$ for each curve.

By comparing Fig. 5.a and Fig. 5.b, it is clear that the case with the spin-
texture results in a significantly larger anomalous phase. The equation for
anomalous phase in Eq. 42 indicates that the driver of the anomalous phase is
the asymmetry between the two spin-bands of the system. To this end, we define
an asymmetry parameter for some quantity $x_{\pm}$ as [17]

$\lambda(x_{\pm})\equiv\absolutevalue{\frac{x_{+}-x_{-}}{x_{+}+x_{-}}}\;.$
(44)

Figure 5.c compares the asymmetry of $s^{r}_{1\pm}$ with $v^{r}_{1\pm}$ and
their ratio $s^{r}_{1\pm}/v^{r}_{1\pm}$. Fig. 5.c indicates that the asymmetry
within the device is almost entirely due to the spin-texture rather than the
Fermi velocity asymmetry. This is not unsurprising as although increasing the
spin-orbit coupling causes the Fermi velocities to deviate slightly, the spin-
projections can experience a full sign change. These results indicate that it
is the spin-texture which is the largest driver of the anomalous Josephson
effect.

##  IV Anomalous phase in multi-channel systems

In Sec. III, having only a single subband occupied allowed for a number of
analytic simplifications. However, in systems containing multiple occupied
subbands, and with a two-dimensional Rashba spin-orbit term mixing states, it
is not possible to evaluate the dispersion relation analytically. As the
generic expression for the anomalous phase in Eq. 25 requires the spin-
projection and group velocity of every state at the Fermi level, then it is
also not possible to derive analytic expressions for the resulting anomalous
phase. Fortunately, there are some approximations which can be made when
multiple transverse subbands are occupied. For example, for a large enough
Fermi level relative to the subband spacing, the linearised dispersions will
appear almost parallel. As a result, we can replace the Fermi velocity of each
state by simply the velocity of the lowest subband $v_{1}$; this is a
reasonable upper bound as in the absence of spin-orbit interactions this would
be the supremum of the set of Fermi velocities. By making this replacement,
the phase gained by the Andreev bound state moving through the magnetic field,
defined in Eq. 17, is now given by

$\frac{2Ls_{j\pm}^{r}E_{Z}}{\hbar
v_{j\pm}^{r}}\to\frac{2Ls_{j\pm}^{r}E_{Z}}{\hbar
v_{1}}\equiv\theta_{B}s_{j,\pm}^{r}\;,$ (45)

where we have pulled out all constant coefficients into a magnetic field
induced phase factor $\theta_{B}$ [18]. This simplification explicitly ignores
Fermi velocity asymmetry, allowing us to focus entirely on the spin-texture
asymmetry; even for the case of a single occupied band, we demonstrated in
Fig. 5.c that this is a reasonable simplification to make. This can then be
substituted into the expression for the total phase shift given by Eq. 25

$\varphi_{0}=\atan\left[\frac{\sum_{j}\left(\sin\theta_{B}s_{j+}^{r}+\sin\theta_{B}s_{j-}^{r}\right)}{\sum_{j}\left(\cos\theta_{B}s_{j+}^{r}+\cos\theta_{B}s_{j-}^{r}\right)}\right]\;.$
(46)

Another approximation we can make is that in multi-channel systems, the Fermi
wavelength of the electron within the lowest subband is significantly smaller
than the confined width of the system $\lambda_{1}\ll W$. Combining this with
the fact that we are considering systems where the Zeeman energy is smaller
than the subband spacing, $E_{Z}<\Delta E\sim\hbar^{2}/(2mW^{2})$, then we can
bound the magnetic phase $\theta_{B}$ by

$\displaystyle\theta_{B}\equiv\frac{2LE_{Z}}{\hbar v_{1}^{r}}<\frac{2L\Delta
E}{\hbar v_{1}^{r}}\sim\frac{L}{W}\frac{\lambda_{1}}{W}\;,$ (47)

where we have substituted in the Fermi wavelength of the lowest subband
electron $\lambda_{j}\equiv h/mv_{j}$. By making the multi-channel
approximation that $\lambda_{1}\ll W$, and implicitly assuming that the length
of the device is not significantly greater than its width, we find that
$\theta_{B}\ll 1$. Utilising this result, the anomalous phase in Eq. 46 can be
Taylor expanded to first order in the magnetic field strength to write

$\displaystyle\varphi_{0}$
$\displaystyle\sim\atan\left[\frac{\theta_{B}}{2N}\sum_{j}^{N}\left(s_{j+}^{r}+s_{j-}^{r}\right)\right]\;,$
(48)
$\displaystyle\sim\frac{\theta_{B}}{2N}\sum_{j}^{N}\left(s_{j+}^{r}+s_{j-}^{r}\right)\;,$
$\displaystyle\equiv\left<s_{y}\right>\theta_{B}\;,$

where we have defined the average spin-projection over all $N$ occupied
subbands as

$\left<s_{y}\right>\equiv\frac{1}{2N}\sum_{j}^{N}\left(s_{j+}^{r}+s_{j-}^{r}\right)\;.$
(49)

Eq. 48 is one of the main results of this work – namely, an expression for the
anomalous phase which is linear in the external magnetic field, and directly
proportional to the spin-texture asymmetry. Importantly, the only physical
mechanism used to derive this result was the inclusion of the spin-texture
induced by the two-dimensional Rashba spin-orbit coupling.

###  IV.1 Average spin-projection in multi-channel systems

To understand the microscopic origin that gives rise to a non-zero average
spin-projection, $\left<s_{y}\right>\neq 0$, it is worth considering some
generic features of the transverse subbands of the non-superconducting region
of the Josephson junction. Indeed, we can view the two-dimensional Rashba
spin-orbit coupling as one term, $\alpha\sigma_{y}p_{x}$, which shifts the
wavevectors of the different spin-states by some amount $\pm k_{\alpha}$,
where [35]

$k_{\alpha}=\frac{m\alpha}{\hbar^{2}}\;,$ (50)

and another term, $\alpha\sigma_{x}p_{y}$, which couples the transverse
subbands and mixes the spins. To illustrate this, the dispersion relation for
a system with $N=7$, no magnetic field, and varying spin-orbit strengths is
plotted in Fig. 6. As in Fig. 2, this was computed using the tight-binding
calculation described in Ref. [26], however in this case we are considering a
wider range of spin-orbit values. The output eigenstates of this numerical
calculation can then be used to extract the spin-projection of that state by
solving Eq. 12 numerically.

Figure 6: Dispersion relation for a system with seven occupied subbands –
totalling fourteen different spin states – with a two-dimensional Rashba spin-
orbit coupling term and no external magnetic field. In each plot the spin-
orbit coupling strength is set such that: (a) $k_{\alpha}=0.1k_{1}$ (b)
$k_{\alpha}=k_{1}$ (c) $k_{\alpha}=2k_{1}$. Note that unlike Fig. 2, here we
are considering a wide range of spin-orbit strengths

From analysing Fig. 6, we observe that even for small spin-orbit strengths
there is significant mixing between the spin-projections of right-moving
states. However, when $k_{\alpha}\sim k_{1}$ the coupling between the
different spin states is great enough that the spin-projection of each state
is locked to the sign of its wavevector i.e. in this limit essentially every
right-moving state is spin-up. This is strongly reminiscent of the band-
structure of one-dimensional nanowires in the topological insulator limit
where instead of a two dimensional spin-orbit interaction, a magnetic field
parallel to the nanowire provides the coupling between spin states [36].
Interestingly, when the spin-orbit strength is increased such that
$k_{\alpha}\gg k_{1}$, the different spin bands are so far shifted from each
other that they effectively decouple. This results in two almost separate
parabolas which appear to have no subband coupling. As a result, the right-
moving states are either well-defined spin-up or spin-down states.

To investigate this point further, by considering the same system as in Fig.
6, we can extract the spin-projections at the Fermi level to compute the
average spin-projection. We plot the corresponding average spin-projection in
Fig. 7.a. The net result is a function that initially increases somewhat
linearly before reaching a maximum at the point $k_{\alpha}=k_{1}$, then
decreasing to effectively zero at $k_{\alpha}=2k_{1}$. The sharp jumps
observed in the average are due to discrete changes in the number $N$ of
right-moving states within the system which strongly perturbs the mean.

Figure 7: (a) Average spin-projection $\left<s_{y}\right>$ for the system
defined in Fig. 6 as a function of the Rashba spin-orbit coupling strength.
(b) The average spin-projection of three systems where the number of occupied
subbands differs – every other parameter is kept fixed. The spin-orbit
wavevector has been restricted to $0\leq k_{\alpha}\leq 0.2k_{1}$ such that we
are best observing only the ‘linear’ region.

These qualitative features can be understood from the dispersion relation for
the states depicted in Fig. 6. When $k_{\alpha}\sim k_{1}$ we expect to
observe a maximum average spin-projection as all right-moving states have
positive spin-index. Similarly, when the spin-orbit strength is increased
further, the spin bands effectively decouple such that the spin-projections of
the occupied spin-up and spin-down states cancel – hence the average should
decreases towards zero as observed.

A particular feature of Fig. 7.a which we want to focus on is that for $0\leq
k_{\alpha}\ll k_{1}$, the average spin-projection appears to be a linear
function of the spin-orbit coupling strength

$\left<s_{y}\right>\sim\frac{k_{\alpha}}{k_{1}}\;.$ (51)

To understand this, we first note that what perturbs the average spin-
projection is not the mixing between occupied states – as mixing between
occupied states does not change the sum of the occupied spin-projections due
to pairwise cancellation – but rather unoccupied higher subbands which mix
with an occupied state. This is most clear in the case of a single occupied
subband studied in Sec. III.1, where it was the coupling with the unoccupied
$j=2$ subband with the occupied $j=1$ subband which drove the effect. We leave
the full description for Appendix. B, however, the reason for the apparent
linear increase is that the highest occupied spin-down state has its spin
fully flipped whenever the spin-orbit strength wavevector $k_{\alpha}$ is an
integer multiple of $k_{1}/N$. Effectively, as a function of the spin-orbit
strength, the total sum of the spin-projections increases by two at regular
intervals of $k_{1}/N$. In the limit as $N\to\infty$, this leads to the
apparent linear increase in the average spin-projection. To demonstrate this,
in Fig. 7.b we plot the average spin-projection for systems with different
numbers of occupied states $N$ and observe that the average spin-projection
progressively becomes more linear with $N\to\infty$.

###  IV.2 Anomalous phase

The expression for the anomalous phase in multi-channel systems given by Eq.
48 states that the anomalous phase is equal to the average spin-projection
multiplied by the phase factor induced by the magnetic field $\theta_{B}$.
From the results of Sec. IV.1, this allows us to make some qualitative
predictions on the resulting anomalous phase. Most interestingly, for Rashba
spin-orbit strengths such that $k_{\alpha}<k_{1}$, the resulting average spin-
projection can be written approximately as $\left<s_{y}\right>\approx
k_{\alpha}/k_{1}$. As a result, within this regime the total anomalous phase
is given by the formula

$\varphi_{0}=\frac{k_{\alpha}\theta_{B}}{k_{1}}\;.$ (52)

By substituting in the relation between the Fermi wavevector and velocity
$k_{1}\sim mv_{1}/\hbar$, and the expression for $\theta_{B}$ given in Eq. 45,
then the anomalous phase is given by

$\varphi_{0}=\frac{2\alpha E_{Z}L}{\hbar^{2}v_{1}^{2}}\;.$ (53)

Up to a factor of two, Eq. 53 is the same expression as derived in Ref. [37]
within the context of ballistic junctions in the long junction limit
$L\gg\xi$. In that work, the model was not microscopic but instead
phenomenological hence it is interesting to note that the results are in such
good agreement.

More generically, as the spin-orbit strength is increased further such that
$\alpha\sim\hbar v_{1}$ then the average spin-projection reaches a maximum
value close to unity such that

$\varphi_{0}\sim\theta_{B}\;,$ (54)

where small changes in the spin-orbit strength have little to no effect on the
resulting anomalous phase. Increasing the spin-orbit strength further such
that $\alpha\gg\hbar v_{1}$, the average spin-projection should appear to
decrease linearly as the different spin bands decouple. These heuristic
expectations on the anomalous phase can be conveyed by the piecewise function

$\varphi_{0}\sim\begin{cases}\theta_{B}{\alpha}/({\hbar
v_{1}})\;,&\alpha\ll\hbar v_{1}\\\ \theta_{B}\;,&\alpha\sim\hbar v_{1}\\\
\theta_{B}\left[2-{\alpha}/({\hbar v_{1}})\right]\;,&\alpha\gg\hbar v_{1}\\\
0\;,&\alpha\to\infty\end{cases}\;$ (55)

##  V Numerical results

Having developed a mechanism for the anomalous Josephson effect specifying the
role of the spin texture, we now look to compare the predictions with
numerical calculations which solve for the anomalous phase directly. To do so
the Hamiltonian of Eq. 2 is cast into a tight binding model and the anomalous
phase is solved for using the non-equilibrium Green’s function formalism
(NEGF).

For the remainder of this section the electron effective mass $m^{*}$ is set
to $0.023m$, and the temperature set such that $2k_{B}T=|\Delta|$ and hence,
as usual, we are working close to the critical temperature; both are assumed
constant throughout the system. The magnitude of the superconducting order
parameter is set to $1\mu$eV – although this is reasonably small, it is simply
to enforce the short junction limit for every numerical result and does not
qualitatively change the results. Similarly, a temperature close to the
critical temperature was used simply to retain a sinusoidal Josephson current
which best illustrates the anomalous phase shift.

###  V.1 Numerical methods

We consider a discrete system with lattice spacing $a$ and a length of $Na$
and a width of $Ma$ for integers $N$ and $M$. Within a tight-binding model,
the Hamiltonian for the non-superconducting central region of the device,
given by Eq. 2, is written as

$\displaystyle H_{C}=$ $\displaystyle\sum_{(i,j)\in
C}\psi_{i,j}^{\dagger}\left[4t+E_{Z}\sigma_{y}-\mu\right]\psi_{i,j}$ (56)
$\displaystyle-\sum_{\langle
i,j\rangle}\left[\psi_{i,j}^{\dagger}\left(t+\frac{i\alpha}{2a}\sigma_{y}\right)\psi_{i+1,j}\right.$
$\displaystyle\qquad+\left.\psi_{i,j}^{\dagger}\left(t-\frac{i\alpha}{2a}\sigma_{x}\right)\psi_{i,j+1}+\text{h.c.}\right]\;,$

where the operator
$\psi_{i,j}\equiv\begin{pmatrix}\psi_{i,j;\uparrow}&\psi_{i,j;\downarrow}\end{pmatrix}^{T}$
annihilates an electron at the lattice site $(i,j)$ [38, 39, 40]. The hopping
parameter $t\equiv\hbar/(2m^{*}a^{2})$ factors in both the effective mass
$m^{*}$ and the lattice spacing $a$ of the simulation. Written in this
quadratic form, we can extract the corresponding first quantised Hamiltonians
for the central region as an $(2NM\times 2NM)$ block tri-diagonal matrix

$h_{C}=\begin{bmatrix}\varepsilon_{1}&t_{1}&&&\\\
t_{1}^{\dagger}&\varepsilon_{2}&t_{2}&&\\\ &\ddots&\ddots&\ddots&&\\\
&&t_{N-1}^{\dagger}&\varepsilon_{N-1}&t_{N}\\\
&&&t_{N}^{\dagger}&\varepsilon_{N}\end{bmatrix}\;,$ (57)

where

$\displaystyle\varepsilon_{j}$
$\displaystyle=\text{TriDiag}\left[-t+\frac{i\alpha}{2a}\sigma_{y},\,4t-\mu+E_{Z}\sigma_{y},\,-t-\frac{i\alpha}{2a}\sigma_{y}\right]\;,$
(58) $\displaystyle t_{j}$
$\displaystyle=\text{Diag}\left[-t+\frac{i\alpha}{2a}\sigma_{x}\right]\;,$

are themselves $(2M\times 2M)$ banded matrices [41].

For the superconducting sections, as indicated in Eq. 7, the degrees of
freedom are doubled to include the holes within the system

$\Psi=\begin{pmatrix}\psi\\\ \psi^{*}\end{pmatrix}\;,$ (59)

where $\psi$ is a vector containing the annihilation operator at every lattice
site. The superconducting order parameter $\Delta$ now acts to couple these
two subspaces [42]

$\displaystyle H_{L/R}$ $\displaystyle=\Psi^{\dagger}h_{L/R}\Psi$ (60)
$\displaystyle=\begin{pmatrix}\psi^{\dagger}&\psi^{T}\end{pmatrix}\begin{pmatrix}h_{0}+h_{\alpha}&i\sigma_{y}\Delta_{L/R}\\\
-i\sigma_{y}\Delta_{L/R}&-(h_{0}+h_{\alpha})^{*}\end{pmatrix}\begin{pmatrix}\psi\\\
\psi^{*}\end{pmatrix}$

Using the discretised first-quantised Hamiltonians the Green’s functions for
the central region at energy $\varepsilon$ can be determined by solving the
corresponding Dyson’s equation [24]

$\mathcal{G}(\varepsilon)=[(\varepsilon+i\kappa)\mathbb{I}-h_{C}-\Sigma_{L}-\Sigma_{R}]^{-1}\;,$
(61)

where $\Sigma_{L/R}$ are the self-energy terms due to the coupling from the
leads, and $\kappa$ is an infinitesimal energy perturbation. By assuming that
the leads are spatially homogenous along the $x$-axis, then the self-energies
can be solved using an efficient recursive process [43]. Furthermore, assuming
both leads are at thermal equilibrium at temperature $T$, a necessary quantity
known as the lesser self-energy, $\Sigma^{<}$, can be computed as [38]

$\displaystyle\Sigma^{<}$
$\displaystyle=f(\varepsilon)(\Sigma_{L}^{\dagger}-\Sigma_{L})+(\Sigma_{R}^{\dagger}-\Sigma_{R})f(\varepsilon)\;,$
(62) $\displaystyle\equiv\Sigma^{<}_{L}+\Sigma^{<}_{R}$

where $f(\varepsilon)$ is the usual Fermi function.

From these quantities, the current can be computed by [38]

$I=e\int\left(\mathcal{G}\Sigma^{<}_{L}+\mathcal{G}\Sigma^{<}\mathcal{G}^{\dagger}\Sigma_{L}^{\dagger}\right)\differential\omega\;.$
(63)

Another important quantity is the density of states within the device which is
computed by [24]

$D=i\Tr\left[\mathcal{G}-\mathcal{G}^{\dagger}\right]\;,$ (64)

where the trace is over both spatial and spin indices.

###  V.2 Single-mode system

In Sec. III analytic formulas for the anomalous phase shift in single band
systems were derived. We now want to explicitly evaluate the anomalous phase
numerically for these single-band systems to compare with the analytic
formulas. To do so, the Fermi level $\mu$ of the numerical simulations will be
tuned such that it lies in-between the first and second transverse subband
energy. To keep the electron density constant as the Rashba spin-orbit
interaction is increased, the Fermi level is modified as [18]

$\mu\to\mu-\frac{\hbar^{2}k_{\alpha}^{2}}{2m^{*}}\;.$ (65)

The lattice spacing was set to $0.1\,\mu$m, and the width and length to
$0.8\,\mu$m and $12\,\mu$m respectively. It is important to reiterate that
this is not a one-dimensional simulation, but rather a two-dimensional
simulation with a tuned Fermi level such that only a single band is occupied.

To confirm that only a single subband is occupied, the normalised density of
states $D$ as a function of the phase difference $\varphi$ across the
Josephson junction is plotted in Fig. 8. In each plot the spin-orbit coupling
strength $\alpha$ is varied; to break the spin degeneracy, a constant magnetic
field with strength $E_{Z}=15\,\mu$eV is also applied.

Figure 8: Normalised density of states computed using Eq. 64 for a Josephson
junction with width $800\;$nm and length $12\mu$m as a function of both energy
and phase difference $\varphi$ across the superconducting leads. The Fermi
level $\mu$ has been tuned such that only a single subband is occupied. As
stated in the text, a small magnetic field with Zeeman strength of
$E_{Z}=15\,\mu$eV is applied to induce a phase shift between the Andreev bound
states. As the Rashba spin-orbit coupling strength $\alpha$ is increased the
asymmetry between the two spin-bands becomes evident.

We can compare the peaks of the density of states with the equations for the
Andreev bound state energies given by Eq. 19. For the case of no spin-orbit
coupling ($\alpha=0$) such that the states have trivial spin-projection, then
the magnetic field identically shifts each Andreev bound state in the opposite
direction by an amount given by $\theta_{B}$ – this is essentially the process
depicted in Fig. 3.b. However, increasing the strength of the Rashba spin-
orbit interaction leads to a mixing of the spin-projections. This causes the
amount each band is shifted to become asymmetric – a process illustrated in
Fig. 3.c. This is taken to the extreme when $\alpha=1\times 10^{-11}$eVm where
both Andreev bound states are visually coincident. This indicates that at this
spin-orbit coupling strength, the formerly spin-down state is so strongly
coupled with the higher order spin-up subband that its spin has effectively
flipped. Due to this spin-projection flipping, both occupied states behave
identically under the external magnetic field. Note that this change in
direction of the phase shift cannot be explained solely by Fermi velocity
asymmetry, as the sign of the Fermi velocity of the right-moving state is
always positive.

Figure 9: (a) The Josephson current corresponding to Fig. 8 evaluated using
Eq. 63. (b) The extracted anomalous phase for the numerical NEGF calculation
overlaid over the results from the analytic band structure calculation given
by solving Eq. 42. (c) The anomalous phase as a function of the external
magnetic field. (d) The anomalous phase as a function of the length $L$ of the
system. As the length is quantised in units of the lattice spacing $a$, we
only plot a discrete set of data points; unlike in (c) where a much higher
resolution could be simulated. Note that the ratio between the magnetic field
$B$ and the length $L$ of the device was held fixed which explains why the
plots (c) and (d) appear identical. In each plot the Rashba spin-orbit
coupling strength has units of $10^{-11}$eVm – this information was removed to
safe space within the figure.

The Josephson currents corresponding to each trace of Fig. 8 are plotted in
Fig. 9.a by solving Eq. 63 numerically. Fig. 9.a clearly demonstrates that as
the spin-orbit coupling strength is increased, the sinusoidal Josephson
current experiences a phase shift. Although not the purpose of this work, we
also observe that the critical current – the maximum supercurrent – increases
slightly with the spin-orbit coupling strength, as predicted in Ref. [17].
This anomalous phase shift can then be extracted and plotted as a function of
the Rashba spin-orbit strength, which is shown in Fig. 9.b. Overlaid on top is
the analytic result for the anomalous phase shift given by Eq. 42. For small
spin-orbit strengths the agreement between the two results are excellent,
however, for larger spin-orbit strengths a slight discrepancy is visible; this
is likely due to the fact that the analytic result was only computed
perturbatively up to first order in the spin-orbit strength.

Another computational test for the analytic equation given by Eq. 42 is by
numerically evaluating the anomalous phase as a function of both the Zeeman
strength $E_{Z}$, and the length $L$ of the device. The results for this are
plotted in Fig. 9.c and Fig 9.d respectively. The plots here are strikingly
linear, which is unsurprising given that the numerical simulations are
strictly within the short-junction limit and that the Zeeman energy is
significantly smaller than the subband spacing – the two requirements for the
linear behaviour discussed in Sec. II.3.

###  V.3 Multi-mode system

In this section, we want to compare with the results of Sec. IV concerning
systems with multiple occupied subbands. To do so numerically, the lattice
spacing $a$ is set to $0.25\mu$m and the width and length to $10\,\mu$m and
length $20\,\mu$m respectively. Although this constitutes a rather large
device, the small subband spacing this induces makes the numerical
calculations easier. It should be noted that although this small subband
spacing imposes a strict energy-scale for the other parameters such that we
can observe the anomalous physics – on the order of $\mu$eV – it does not
qualitatively change the results.

We tune the Fermi level such that there are twelve occupied states ($N=12$).
As in the single-band case, to illustrate that there are multiple occupied
states, we can plot the density of states for the Josephson junction as a
function of the phase difference $\varphi$ across the two superconductors. As
before, to break the spin degeneracy of the system a magnetic field with
strength $E_{Z}=15\mu$eV is applied. This result is illustrated in Fig. 10.
Note that rather than referring to the spin-orbit strength $\alpha$ directly,
it is simpler to refer to the spin-orbit wavevector $k_{\alpha}$ defined in
Eq. 9.

Figure 10: Normalised density of states computed using Eq. 64 for a Josephson
junction with width $10\mu$m and length $20\mu$m as a function of both energy
and phase difference $\varphi$ across the superconducting leads. The Fermi
level is tuned such that there are twelve occupied subbands. As stated in the
text, a small magnetic field of Zeeman strength $E_{Z}=15\;\mu$eV is applied
to induce a relative phase shift between the Andreev bound states.

Making the multi-channel approximations discussed in Sec. IV.1, the energy of
each Andreev bound states derived in Eq. 19 can be written as

$\displaystyle\varepsilon^{r}_{j\pm}(\varphi)$
$\displaystyle=-|\Delta|\cos(\frac{\varphi}{2}-{s_{j\pm}^{r}}\theta_{B})\;,$
(66) $\displaystyle\varepsilon^{l}_{j\pm}(\varphi)$
$\displaystyle=-\varepsilon^{r}_{j\mp}(\varphi)\;,$

such that the amount each Andreev bound state is shifted by is entirely
determined by the spin-projection of its constituent linearised electron
state, $s_{j\pm}^{r}$. Considering this, the striking observation of Fig. 10
is that when the spin-orbit strength is increased such that $k_{\alpha}\sim
k_{1}$, there appears to be only a single state phase shifted to the right,
the rest are phase shifted towards the left. Qualitatively, this agrees with
Fig. 6.b where we found that at this spin-orbit coupling strength, almost
every occupied right-moving state has a positive spin-projection. As a result,
almost every Andreev bound state will be shifted in the same direction under
the external magnetic field. Similarly, when the spin-orbit strength is
increased further such that $k_{\alpha}=1.5k_{1}$, the normalised density of
states appears similar to the system with no spin-orbit coupling. This is
again predicted by Fig. 6.c where we found the spin-bands decouple such that
both spin-up and spin-down states are again occupied. Due to the occupation of
both states with positive and negative spin-projection, there will be Andreev
bound states present which are shifted in either direction as observed.

Figure 11: (a) Josephson current computed using Eq. 63 for the system outlined
in Sec. V.3. (b) The corresponding anomalous phase extracted from the
Josephson current. The blue markers indicate the non-equilibrium Green
function (NEGF) calculation, whilst the overlaid red markers indicate the
dispersion relation calculation for the average spin-projection outlined in
Sec. IV.1. (c) The anomalous phase as a function of the spin-orbit coupling
strength computed from the NEGF calculation for three system with different
numbers of occupied subbands. To change the number of occupied subbands, every
parameter remained the same except the Fermi level which was manually
adjusted.

Solving for the current given by Eq. 63, we plot the Josephson current
corresponding to each figure within Fig. 10 in Fig. 11.a. In Fig. 11.a we
observe that the sinusoidal Josephson relation experiences a phase shift as
the spin-orbit strength is increased, however, at larger values the phase
shift appears to decreases. To study this behaviour, we extract the phase
shift from the Josephson current and plot them in Fig. 11.b. Overlaid on top
of the NEGF calculations, we have also plotted the average spin-projection
$\left<s_{y}\right>$ for the same physical system numerically computed using
the dispersion relation calculation of Sec. IV.1. It is important to note that
although both curves are computational, the NEGF results are solving the
dynamics of the system such that we can explicitly extract the anomalous
phase, whilst the dispersion relation calculation is simply computing the
average spin-projection of the occupied right-moving states. The fact that
they are similar is strong evidence that the spin-texture is the predominant
effect in the anomalous Josephson effect for ballistic systems.

An important feature of the predicted anomalous phase was that for spin orbit
strengths where $k_{\alpha}<k_{1}$, the anomalous phase is predicted to be
linear in the spin-orbit strength. To check this, in Fig. 11.c we plot the
anomalous phase extracted from the NEGF calculations within this regime for
systems with different number of occupied states. Overall, we observe that the
behaviour is highly linear, with the slope appearing to approach unity as the
number of subbands increases; effectively identical to the prediction made in
Fig. 7.b.

Considering the other predictions made in Sec. IV.2, for spin-orbit strengths
where $k_{\alpha}\sim k_{1}$ we observe in Fig. 11.b a rounding of the
anomalous phase such that it reaches a plateau. In this regime the anomalous
phase is effectively independent of the value of the spin-orbit strength
$\varphi\sim\theta_{B}$ as predicted. Similarly, for larger spin-orbit
strengths we observe that the numerical anomalous phase decays linearly before
decaying to zero as predicted. It is important to note that these qualitative
results were highly robust to variations in the system such as varying the
number $N$ of occupied states or changing lattice parameters such as the
lattice spacing $a$ or effective mass $m^{*}$. However, the exact nature of
the anomalous phase, most particularly the discrete jumps observed as the
spin-orbit coupling strength was increased – due to the number of occupied
states within the system changing – did differ slightly between computational
samples.

##  VI Conclusions and discussion

In this work, we have studied the anomalous Josephson effect in quasi one-
dimensional ballistic structures with both two-dimensional Rashba spin-orbit
coupling and an external magnetic field. We constructed a microscopic model
for the resulting dynamics by focusing solely on the role of the spin-texture
which forms due to the coupling between adjacent transverse subbands of
opposite spin. This microscopic model found great agreement when compared to
computational non-equilibrium Green function calculations. We note our results
are only applicable for a subset of experimental devices; we are working
purely in the ballistic regime, are studying only the short junction limit,
and have ignored the vector potential. The purpose of this was to focus solely
on the role of spin-texture asymmetry, which represents an important area of
study within the field of proximitised nanowire devices.

For the case of systems with only a single occupied subband we analytically
derived a closed-form expression – albeit requiring the analytic solution of a
quartic polynomial – for the anomalous phase produced. This formula was very
similar to previous results in the literature, however, it explicitly included
a term which factors in the role of spin mixing. Indeed, we found that the
asymmetry in this spin-mixing is the largest driver of the anomalous phase,
rather than asymmetry in the Fermi velocities as has been the main subject of
discussion in the literature.

For the case of systems with multiple occupied subbands, we explicitly ignored
the Fermi velocity asymmetry to solely focus on the role of the spin-texture
within the system. Based on this microscopic model, we can list four
theoretical predictions:

  1. 1.

The anomalous phase should be linear in the external magnetic field for Zeeman
energies weaker than the subband splitting.

  2. 2.

For weaker spin-orbit strengths, $\alpha<\hbar v_{F}$, the anomalous phase
should be linear in the spin-orbit strength.

  3. 3.

For an appropriately large spin-orbit strength ($\alpha\sim\hbar v_{F}$), and
under a constant magnetic field, the anomalous phase should parabolically
reach a maximum value.

  4. 4.

In the limit of large spin-orbit strengths, the anomalous phase is predicted
to decay to zero.

Each of these predictions were directly observed in numerical calculations
computing using the non-equilibrium Green function formalism. As a result, it
appears that for ballistic systems the spin-texture induced by the Rashba
spin-orbit coupling is the critical microscopic quantity which produces the
anomalous Josephson effect.

Interestingly, there has been experimental evidence for the first three
theoretical predictions [7, 5, 6]. Perhaps most striking is our prediction for
the anomalous phase to briefly saturate as a function of the spin-orbit
strength when $\alpha\sim\hbar v_{F}$, where similar behaviour was observed in
Ref. [7]. In that work, they noted that the origin of that observation was
unknown. On the possibility of experimentally observing the final prediction –
namely, that the anomalous phase should vanish in the large spin-orbit
coupling limit – we note that for InAs systems that it is experimentally
possible to tune the spin-orbit coupling strength within two orders of
magnitude hence this limit may be observable for certain systems [44].

Recently the anomalous Josephson effect has been experimentally linked to the
Josephson diode effect [6]. Hence in future work we look to explore explicitly
this link between the spin-texture induced by Rashba spin-orbit coupling and
critical current asymmetry.

## Acknowledgments

We acknowledge helpful discussions with both Francesco Giazotto and Elia
Strambini concerning the anomalous Josephson effect. We also acknowledge
discussions with Tyler Whittaker on the formation of spin-textures in semi-
conducting nanowires.

##  Appendix A Andreev bound state formation

Considering the Josephson junction of length $L$ illustrated in Fig. 1, the
interface between the three regions occur at $\pm L/2$. Focusing on the non-
superconducting central region, after factoring out the fast oscillations on
the order of $k_{j}$, the linearised eigenenergies of Eq. 2 have the form

$\varepsilon(k)=\hbar v_{j\pm}^{l/r}k-s_{j\pm}^{l/r}E_{Z}\;,$ (67)

We also introduce a spinor $\chi^{l/r}_{j\pm}$ defined such that

${\left(\chi^{l/r}_{j\pm}\right)}^{\dagger}\sigma_{y}{\chi^{l/r}_{j\pm}}={s^{l/r}_{j\pm}}\;,$
(68)

As a result, the travelling wave solution for a left or right moving electron
with energy $\varepsilon$ in the normal region is given by

$\ket{\psi^{l/r}_{j\pm}(x)}_{\text{electron}}=\chi^{l/r}_{j\pm}e^{i\left[(\varepsilon+s_{j\pm}^{l/r}E_{Z})/\hbar
v_{j\pm}^{l/r}\right]x}\;,$ (69)

the hole states are then given by

$\ket{\psi^{l/r}_{j\pm}(x)}_{\text{hole}}=i\sigma_{y}\chi^{l/r}_{j\pm}e^{-i\left[(\varepsilon+s_{j\pm}^{l/r}E_{Z})/\hbar
v_{j\pm}^{l/r}\right]x}\;,$ (70)

The same can be done for the proximitised superconducting regions by
evaluating the eigenstates of Eq. 7. However, as we are interested in energies
within the superconducting gap, $\varepsilon<|\Delta|$, then the only allowed
states are evanescent states which decay into the superconductors. The
resulting eigenstates in the left (L) and right (R) superconductors are given
by [28]

$\displaystyle\ket{\psi^{l/r}_{j\pm}(x)}_{L}$
$\displaystyle=\begin{pmatrix}(\varepsilon-i\Lambda)\chi^{l/r}_{j\pm}\\\
\Delta_{L}^{*}(i\sigma_{y})\chi^{l/r}_{j\pm}\end{pmatrix}e^{-\Lambda(x+L/2)/\hbar
v^{l/r}_{j\pm}}\;,$ (71) $\displaystyle\ket{\psi^{l/r}_{j\pm}(x)}_{R}$
$\displaystyle=\begin{pmatrix}(\varepsilon-i\Lambda)\chi^{l/r}_{j\pm}\\\
\Delta_{R}^{*}(i\sigma_{y})\chi^{l/r}_{j\pm}\end{pmatrix}e^{-\Lambda(x-L/2)/\hbar
v^{l/r}_{j\pm}}\;,$

where

$\Lambda\equiv\sqrt{|\Delta|^{2}-\varepsilon^{2}}\;.$ (72)

and we note that we are now working in a particle-hole basis. Within this
basis, the total wavefunction for the non-superconducting central (C) region,
consisting of counter-propagating electrons and holes, is given by the
superposition

$\displaystyle\ket{\psi^{l/r}_{j\pm}(x)}_{\text{C}}=A\begin{pmatrix}\chi^{l/r}_{j\pm}\\\
0\end{pmatrix}e^{i\left[(\varepsilon+s_{j\pm}^{l/r}E_{Z})/\hbar
v_{j\pm}^{l/r}\right]x}$ (73) $\displaystyle+B\begin{pmatrix}0\\\
i\sigma_{y}\chi^{l/r}_{j\pm}\end{pmatrix}e^{-i\left[(\varepsilon+s_{j\pm}^{l/r}E_{Z})/\hbar
v_{j\pm}^{l/r}\right]x}\;.$

By matching the wavefunctions at $x=\pm L/2$ – a necessary condition for the
formation of a bound state – we find that

$\frac{\varepsilon-i\Lambda}{\varepsilon+i\Lambda}=\frac{\Delta_{L}}{\Delta_{R}}e^{-2iL\varepsilon/\hbar
v_{j\pm}^{l/r}}e^{-2iLs_{j\pm}^{l/r}E_{Z}/\hbar v_{j\pm}^{l/r}}\;.$ (74)

Utilising some identities, this can be rewritten as

$e^{-2i\text{acos}(\varepsilon/|\Delta|)}=e^{i\varphi}e^{-2iL\varepsilon/\hbar
v_{j\pm}^{l/r}}e^{-2iLs_{j\pm}^{l/r}E_{Z}/\hbar v_{j\pm}^{l/r}}\;,$ (75)

as in the maintext.

##  Appendix B Linear average spin-projection

To understand why the average spin-projection is approximately linear in the
spin-orbit coupling strength, we first note that what perturbs the average
spin-projection is not the mixing between occupied states – as any mixing does
would not change the total sum of the occupied spin-projections – but rather
unoccupied higher subbands which mix with an occupied state. Hence the
critical parameter to determine is at what spin-orbit coupling strengths does
some unoccupied state couple with an occupied state. To determine at what
specific spin-orbit value a higher subband mixes with an occupied state, it is
useful to briefly return to the case of a one-dimensional Rashba spin-orbit
coupling where there is no overlap between adjacent subbands. In this case,
for a system with $N$ states, the Fermi-wavevectors of each subband can be
computed analytically and are given by

$k_{j\pm}^{r}=\frac{\pi}{W}\sqrt{N^{2}-(j-1)^{2}}\pm k_{\alpha}\;,$ (76)

As a result, the $j^{\text{th}}$ spin-up state overlaps with the
$(j+q)^{\text{th}}$ spin-down at the Fermi level when the spin-orbit strength
is tuned such that

$k_{\alpha}=\frac{\pi}{2W}\left[\sqrt{N^{2}-(j-1)^{2}}-\sqrt{N^{2}-(j+q-1)^{2}}\right]\;,$
(77)

Pulling out a factor or $k_{1}=N\pi/W$, and Taylor expanding the roots, we can
write the Rashba wave-vector required to overlap the states at the Fermi level
as

$\frac{k_{\alpha}}{k_{1}}\sim\frac{qj}{2N^{2}}$ (78)

Returning to the proper two-dimensional Rashba interaction – such that
overlapping bands form avoided crossings – the spin-orbit strengths defined in
Eq. 78 now define critical values where we expect to observe a spin-flip
between the $j^{\text{th}}$ state and its $q^{\text{th}}$ neighbour. To
demonstrate this, we plot the different spin-projections of each occupied
right-moving state as the spin-orbit coupling strength is increased a system
with $N=15$ occupied subbands (thirty occupied spin states) in Fig. 12.a.
Overlaid on top is the critical spin-orbit strength defined in Eq. 78.
Visually it provides a good approximation to when the spin-projections flip
sign.

Figure 12: (a) The spin-projection of each state $s_{j}$ plotted as a function
of the spin-orbit coupling strength. The black regions indicate where a spin-
projection has changed signs. Overlaid in green is the critical spin-orbit
strength – given by Eq. 78 – where we naïvely expect a sign flip to occur. (b)
The average spin-projection of the system defined in (a).

Unoccupied higher subband states always couple initially with the highest
occupied state – namely, the $j=N$ state. As a result, using Eq. 78 the
highest occupied subband will be fully mixed with the $N+q$ subband when
$k_{\alpha}\sim qk_{1}/2N$. Although this mixing likely obeys a complicated
sigmoid as in Eq. 41 for the single subband case, we can consider an ideal
system where the states begin mixing initially when $k_{\alpha}=0$, are fully
mixed when $k_{\alpha}\sim qk_{1}/2N$, and doubling this, have experienced a
full sign change in the spin-projection when $k_{\alpha}\sim qk_{1}/N$. As a
result, as the Rashba wavevector $k_{\alpha}$ is increased, every multiple of
${k_{1}}/{N}$ the total spin of the right-moving occupied states increases by
two. We can approximate these jumps by the use of the Heaviside step function
such that

$\left<s_{y}\right>\sim\frac{1}{2N}\sum_{j=1}^{N}2\Theta\left(\frac{k_{\alpha}}{k_{1}}-\frac{j}{N}\right)\;,$
(79)

In the limit as $N\to\infty$ we can replace the discrete sum with an integral
– note, we have replaced the discrete $j/N$ with the continues variable $z$

$\displaystyle\left<s_{y}\right>$
$\displaystyle\underset{N\to\infty}{\longrightarrow}\frac{1}{2}\int_{0}^{1}2\Theta\left(\frac{k_{\alpha}}{k_{1}}-z\right)\;dz\;,$
(80) $\displaystyle=\frac{k_{\alpha}}{k_{1}}\;.$

To illustrate the effectiveness of this approximation, in Fig. 12.b we plot
the average spin-projection for the $N=15$ system shown in Fig. 12.a. We
visually observe a very good agreement with the linear approximation.

## References

  * Josephson [1962] B. D. Josephson, Possible new effects in superconductive tunnelling, Physics Letters 1, 251 (1962).
  * Bulaevskii  _et al._ [1977] L. N. Bulaevskii, V. V. Kuzii, and A. A. Sobyanin, Superconducting system with weak coupling to the current in the ground state, JETP Lett. (USSR) (Engl. Transl.); (United States) 25:7 (1977).
  * Rasmussen  _et al._ [2016] A. Rasmussen, J. Danon, H. Suominen, F. Nichele, M. Kjaergaard, and K. Flensberg, Effects of spin-orbit coupling and spatial symmetries on the Josephson current in SNS junctions, Physical Review B 93, 155406 (2016).
  * Assouline  _et al._ [2019] A. Assouline, C. Feuillet-Palma, N. Bergeal, T. Zhang, A. Mottaghizadeh, A. Zimmers, E. Lhuillier, M. Eddrie, P. Atkinson, M. Aprili, and H. Aubin, Spin-Orbit induced phase-shift in Bi2Se3 Josephson junctions, Nature Communications 10, 126 (2019).
  * Strambini  _et al._ [2020] E. Strambini, A. Iorio, O. Durante, R. Citro, C. Sanz-Fernández, C. Guarcello, I. V. Tokatly, A. Braggio, M. Rocci, N. Ligato, V. Zannier, L. Sorba, F. S. Bergeret, and F. Giazotto, A Josephson phase battery, Nature Nanotechnology 15, 656 (2020).
  * Reinhardt  _et al._ [2023] S. Reinhardt, T. Ascherl, A. Costa, J. Berger, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, J. Fabian, D. Kochan, C. Strunk, and N. Paradiso, Link between supercurrent diode and anomalous Josephson effect revealed by gate-controlled interferometry (2023), arxiv:2308.01061 [cond-mat] .
  * Mayer  _et al._ [2020] W. Mayer, M. C. Dartiailh, J. Yuan, K. S. Wickramasinghe, E. Rossi, and J. Shabani, Gate controlled anomalous phase shift in Al/InAs Josephson junctions, Nature Communications 11, 212 (2020).
  * Shukrinov [2022] Y. M. Shukrinov, Anomalous Josephson effect, Physics-Uspekhi 65, 317 (2022).
  * Moroz and Barnes [1999] A. V. Moroz and C. H. W. Barnes, Effect of the spin-orbit interaction on the band structure and conductance of quasi-one-dimensional systems, Physical Review B 60, 14272 (1999).
  * Bercioux and Lucignano [2015] D. Bercioux and P. Lucignano, Quantum transport in Rashba spin–orbit materials: A review, Reports on Progress in Physics 78, 106001 (2015).
  * Governale and Zülicke [2002] M. Governale and U. Zülicke, Spin accumulation in quantum wires with strong Rashba spin-orbit coupling, Physical Review B 66, 073311 (2002).
  * Murani  _et al._ [2017] A. Murani, A. Chepelianskii, S. Guéron, and H. Bouchiat, Andreev spectrum with high spin-orbit interactions: Revealing spin splitting and topologically protected crossings, Physical Review B 96, 165415 (2017).
  * Park and Yeyati [2017] S. Park and A. L. Yeyati, Andreev spin qubits in multichannel Rashba nanowires, Physical Review B 96, 125416 (2017).
  * Reynoso  _et al._ [2008] A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Anomalous Josephson Current in Junctions with Spin Polarizing Quantum Point Contacts, Physical Review Letters 101, 107001 (2008).
  * Reynoso  _et al._ [2012] A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Spin-orbit-induced chirality of Andreev states in Josephson junctions, Physical Review B 86, 214519 (2012).
  * Tosi  _et al._ [2019] L. Tosi, C. Metzger, M. F. Goffman, C. Urbina, H. Pothier, S. Park, A. L. Yeyati, J. Nygård, and P. Krogstrup, Spin-Orbit Splitting of Andreev States Revealed by Microwave Spectroscopy, Physical Review X 9, 011010 (2019).
  * Krive  _et al._ [2004] I. V. Krive, L. Y. Gorelik, R. I. Shekhter, and M. Jonson, Chiral symmetry breaking and the Josephson current in a ballistic superconductor–quantum wire–superconductor junction, Low Temperature Physics 30, 398 (2004).
  * Yokoyama  _et al._ [2014] T. Yokoyama, M. Eto, and Y. V. Nazarov, Anomalous Josephson effect induced by spin-orbit interaction and Zeeman effect in semiconductor nanowires, Physical Review B 89, 195407 (2014).
  * Mironov  _et al._ [2015] S. V. Mironov, A. S. Mel’nikov, and A. I. Buzdin, Double Path Interference and Magnetic Oscillations in Cooper Pair Transport through a Single Nanowire, Physical Review Letters 114, 227001 (2015).
  * Zuo  _et al._ [2017] K. Zuo, V. Mourik, D. B. Szombati, B. Nijholt, D. J. van Woerkom, A. Geresdi, J. Chen, V. P. Ostroukh, A. R. Akhmerov, S. R. Plissard, D. Car, E. P. A. M. Bakkers, D. I. Pikulin, L. P. Kouwenhoven, and S. M. Frolov, Supercurrent Interference in Few-Mode Nanowire Josephson Junctions, Physical Review Letters 119, 187704 (2017).
  * De Gennes [1964] P. G. De Gennes, Boundary Effects in Superconductors, Reviews of Modern Physics 36, 225 (1964).
  * Beenakker [1992] C. W. J. Beenakker, Three “Universal” Mesoscopic Josephson Effects, in  _Transport Phenomena in Mesoscopic Systems_, Springer Series in Solid-State Sciences, edited by H. Fukuyama and T. Ando (Springer, Berlin, Heidelberg, 1992) pp. 235–253.
  * van Heck  _et al._ [2017] B. van Heck, J. I. Väyrynen, and L. I. Glazman, Zeeman and spin-orbit effects in the Andreev spectra of nanowire junctions, Physical Review B 96, 075404 (2017).
  * Rammer [2007] J. Rammer,  _Quantum Field Theory of Non-equilibrium States_ (Cambridge University Press, Cambridge, 2007).
  * Pedder  _et al._ [2016] C. J. Pedder, T. Meng, R. P. Tiwari, and T. L. Schmidt, Dynamic response functions and helical gaps in interacting Rashba nanowires with and without magnetic fields, Physical Review B 94, 245414 (2016).
  * Mireles and Kirczenow [2001] F. Mireles and G. Kirczenow, Ballistic spin-polarized transport and Rashba spin precession in semiconductor nanowires, Physical Review B 64, 024426 (2001).
  * Andreev [1964] A. F. Andreev, The thermal conductivity of the intermediate state in superconductors, Zhurnal Eksperimental’noj I Teoreticheskoj Fiziki 46 (1964).
  * Sauls [2018] J. A. Sauls, Andreev bound states and their signatures, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20180140 (2018).
  * Kulik [1969] I. O. Kulik, Macroscopic Quantization and the Proximity Effect in S-N-S Junctions, Soviet Journal of Experimental and Theoretical Physics 30, 944 (1969).
  * Bardeen and Johnson [1972] J. Bardeen and J. L. Johnson, Josephson Current Flow in Pure Superconducting-Normal-Superconducting Junctions, Physical Review B 5, 72 (1972).
  * Beenakker [1991] C. W. J. Beenakker, Universal limit of critical-current fluctuations in mesoscopic Josephson junctions, Physical Review Letters 67, 3836 (1991).
  * Nesterov  _et al._ [2016] K. N. Nesterov, M. Houzet, and J. S. Meyer, Anomalous Josephson effect in semiconducting nanowires as a signature of the topologically nontrivial phase, Physical Review B 93, 174502 (2016).
  * Bagwell [1992] P. F. Bagwell, Suppression of the Josephson current through a narrow, mesoscopic, semiconductor channel by a single impurity, Physical Review B 46, 12573 (1992).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–Wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).
  * Manchon  _et al._ [2015] A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, New perspectives for Rashba spin–orbit coupling, Nature Materials 14, 871 (2015).
  * von Oppen  _et al._ [2017] F. von Oppen, Y. Peng, and F. Pientka, Topological superconducting phases in one dimension, in  _Topological Aspects of Condensed Matter Physics: Lecture Notes of the Les Houches Summer School: Volume 103, August 2014_, edited by C. Chamon, M. O. Goerbig, R. Moessner, and L. F. Cugliandolo (Oxford University Press, 2017) p. 0.
  * Buzdin [2008] A. Buzdin, Direct Coupling Between Magnetism and Superconducting Current in the Josephson ${\varphi}_{0}$ Junction, Physical Review Letters 101, 107005 (2008).
  * Jauho  _et al._ [1994] A.-P. Jauho, N. S. Wingreen, and Y. Meir, Time-dependent transport in interacting and noninteracting resonant-tunneling systems, Physical Review B 50, 5528 (1994).
  * Martín-Rodero  _et al._ [1994] A. Martín-Rodero, F. J. García-Vidal, and A. Levy Yeyati, Microscopic theory of Josephson mesoscopic constrictions, Physical Review Letters 72, 554 (1994).
  * Liu and Chan [2010] J.-F. Liu and K. S. Chan, Anomalous Josephson current through a ferromagnetic trilayer junction, Physical Review B 82, 184533 (2010).
  * Sriram  _et al._ [2019] P. Sriram, S. S. Kalantre, K. Gharavi, J. Baugh, and B. Muralidharan, Supercurrent interference in semiconductor nanowire Josephson junctions, Physical Review B 100, 155431 (2019).
  * Zeng  _et al._ [2003] Z. Y. Zeng, B. Li, and F. Claro, Electronic transport in hybrid mesoscopic structures: A nonequilibrium Green function approach, Physical Review B 68, 115319 (2003).
  * Sancho  _et al._ [1984] M. P. L. Sancho, J. M. L. Sancho, and J. Rubio, Quick iterative scheme for the calculation of transfer matrices: Application to Mo (100), Journal of Physics F: Metal Physics 14, 1205 (1984).
  * Wickramasinghe  _et al._ [2018] K. S. Wickramasinghe, W. Mayer, J. Yuan, T. Nguyen, L. Jiao, V. Manucharyan, and J. Shabani, Transport properties of near surface InAs two-dimensional heterostructures, Applied Physics Letters 113, 262104 (2018).