# The Dirac paradox in 1+1 dimensions
and its realization with spin-orbit coupled nanowires

Leonid Gogin  Dipartimento di Scienza Applicata e Tecnologia del Politecnico
di Torino, I-10129 Torino, Italy  Lorenzo Rossi<EMAIL_ADDRESS>Dipartimento di Scienza Applicata e Tecnologia del Politecnico di Torino,
I-10129 Torino, Italy  Fausto Rossi  Dipartimento di Scienza Applicata e
Tecnologia del Politecnico di Torino, I-10129 Torino, Italy  Fabrizio Dolcini
Dipartimento di Scienza Applicata e Tecnologia del Politecnico di Torino,
I-10129 Torino, Italy

###### Abstract

At the interface between two massless Dirac models with opposite helicity a
paradoxical situation arises: A transversally impinging electron can seemingly
neither be transmitted nor reflected, due to the locking between spin and
momentum. Here we investigate this paradox in one spatial dimension where,
differently from higher dimensional realizations, electrons cannot leak along
the interface. We show that models involving only massless Dirac modes lead to
either no solutions or to trivial solutions to the paradox, depending on how
the helicity change across the interface is modeled. However, non trivial
scattering solutions to the paradox are shown to exist when additional massive
Dirac modes are taken into account. Although these modes carry no current for
energies within their gap, their interface coupling with the massless modes
can induce a finite and tunable transmission. Finally, we show that such
massless+massive Dirac model can be realized in suitably gated spin-orbit
coupled nanowires exposed to an external Zeeman field, where the transmission
coefficient can be controlled electrically.

##  I Introduction

Conventional semiconductor heterostructures are typically described, within
the envelope function and effective mass approximations, by a Schrödinger
Hamiltonian with a space dependent effective electron mass varying along the
growth direction and accounting for the different effective masses of the
component materials.

In the last two decades, however, it has been realized that in various
materials such as graphene, topological insulators and Weyl semimetals, the
dynamics of the conduction electrons is well captured, in physically relevant
regimes, by a (D+1)-dimensional massless Dirac electron modelgeim-novoselov ;
castroneto-review ; dassarma_review_2011 ; kane-hasan-review ; zhang-review ;
hasan_natmat_2016 ; yang_2016 , where D denotes the spatial dimension, and
“$+1$” the time dimension. These discoveries have thus spurred the interest in
the investigation of Dirac heterojunctions. Each massless Dirac cone is
characterized by a given helicity of the electron eigenstates, i.e. a sign
encoding the locking between the propagation direction and the orientation of
a “spin-like” degree of freedom, which can be a sublattice pseudospin, like in
graphene, or the actual angular momentum in topological insulators. In
particular, when a junction is formed between two Dirac materials with
opposite helicity, a paradoxical situation emerges, as sketched in Fig.1. A
right-moving electron (blue line on the left-hand side) impinging
transversally towards a spin-inactive interface can neither be transmitted nor
be reflected, due to spin conservation. The Dirac paradox has been discussed
in heterojunctions between two 3D Topological insulators, whose surface states
are governed by a 2D massless Dirac Hamiltonian. In such a case the surface
electrons turn out to “escape” the problem by leaking along the interface
surfacemurakami_2011 ; debeule_2013 ; asano_2013 . However, in a truly 1D
realization of a Dirac model such way out to an extra dimension does not exist
and the Dirac paradox becomes even more interesting. The challenging question
is whether a solution in 1D exists and, if so, whether it can be realized in
some physical system.

Figure 1:  Schematic representation of the Dirac paradox emerging at a
heterojunction between two massless Dirac models with opposite helicity. Blue
and red lines describe spin-$\uparrow$ and spin-$\downarrow$ states,
respectively. On the left-hand side of the junction, right-moving electrons
are characterized by spin-$\uparrow$ and left-moving electrons by
spin-$\downarrow$, while the opposite occurs on the right-hand side. A
spin-$\uparrow$ electron impinging from the left onto the interface can
seemingly neither be transmitted nor reflected, due to spin conservation.

In this paper we investigate the Dirac paradox in 1D and address these
problems. First, we show that, if the helicity change across the interface is
accounted for by an inhomogeneous velocity profile, the paradox has no
solution, in the sense that the continuity equation forbids the existence of
scattering states and only allows for eigenstates that involve electron
injection from both sides of the junction, which carry a vanishing current.
If, however, the interface directly introduces spin-flip processes, the
solution of the paradox is of course trivial and transmission is possible. We
then investigate whether non-trivial solutions with proper scattering states
and finite transmission exist without spin-flip processes at the interface. To
this purpose, we propose an extended model involving both massless and massive
Dirac modes. We show that, despite carrying no current for energies within
their gap, the massive modes play a crucial role in inhomogeneous problems
like the Dirac paradox. Indeed scattering states describing a transmission
from a spin-$\uparrow$ incoming massless mode to a spin-$\downarrow$ outgoing
massless mode do exist, due to the interface coupling between massless and
massive modes. Moreover, the resulting transmission coefficient is finite and
tuneable.

We then discuss the possible realization of such extended model. While
massless Dirac helical states have been proven to exist at the edges of
quantum spin Hall systemsbernevig_science_2006 ; zhang_2008 ; molenkamp-zhang
; knez_2011 ; knez_2015 , this implementation is not optimal for the Dirac
paradox in 1D. Indeed, since these states flow at the boundaries of a 2D
quantum well, an heterojunction between two such wells with opposite edge
helicity would exhibit a linear interface, whereto electrons could leak, like
in the case of heterojunctions between two 3D Topological insulators mentioned
above. However, a truly 1D implementation of helical states has been realized
with spin-orbit coupled nanowires (NWs) exposed to a magnetic
fielddepicciotto_2010 ; kouwenhoven_nanolett_2013 ; kouwenhoven_natcom_2017 ;
kouwenhoven_2012 ; liu_2012 ; heiblum_2012 ; xu_2012 ; defranceschi_2014 ;
marcus_2016 ; marcus_science_2016 ; kouwenhoven_2018 , in the regime of spin-
orbit energy much larger than the Zeeman energystreda ; vonoppen_2010 ;
dassarma_2010 ; loss_PRB_2011 ; lutchyn_2012 ; loss_PRB_2017 .

So far, this remarkable discovery has been mostly exploited in the search for
Majorana quasi-particlesvonoppen_2010 ; dassarma_2010 ; kouwenhoven_2012 ;
liu_2012 ; heiblum_2012 ; xu_2012 ; defranceschi_2014 ; marcus_2016 ;
marcus_science_2016 ; kouwenhoven_2018 ; kitaev ; alicea_review ; fujimoto ;
aguado_review . However, further interesting research areas on NWs are
fostered by the recent advances in gating techniquessasaki_2013 ; micolich ;
sasaki_2017 ; das_2019 ; guo_2021 ; sasaki_2021 , which nowadays enable one to
control the Rashba spin-orbit coupling (RSOC), both in magnitudegao_2012 ;
slomski_NJP_2013 ; wimmer_2015 ; nygaard_2016 ; sherman_2016 ;
tokatly_PRB_2017 ; loss_2018 ; goldoni_2018 ; tsai_2018 ; gao-review ;
lau_2021  and signkaindl_2005 ; meier-dil_2013 ; wang-fu_2016 ; nitta-
frustaglia . Because in a NW the helicity of the massless modes is determined
by the sign of the RSOC, a NW with two differently gated regions can represent
a truly 1D implementation of the Dirac paradox configuration. Notably, in such
an inhomogeneous setup, the massless helical modes are not sufficient to
describe the low energy physics, which turns out to be well captured by the
massless+massive Dirac model proposed here, instead. The resulting
transmission coefficient is finite and can be tuned electrically over a wide
range of values.

Our paper is organized as follows. In Sec.II we analyze the Dirac paradox with
two different models involving only massless modes. Then, in Sec.III we
introduce a model with both massless and massive Dirac modes and show how this
can yield a finite transmission coefficient depending on three parameters.
Furthermore, in Sec.IV we show that this model can be implemented in a
suitably designed setup involving spin-orbit coupled NWs. Finally, in Sec.V we
discuss our results and draw our conclusions.

##  II Massless Dirac heterojunctions

Let us thus consider a junction connecting two 1+1 dimensional massless Dirac
models

$\hat{\mathcal{H}}_{L/R}=v_{L/R}\int\hat{\Psi}^{\dagger}(x)\sigma_{z}p_{x}\,\hat{\Psi}(x)\,dx$
(1)

where $\hat{\mathcal{H}}_{L}$and $\hat{\mathcal{H}}_{R}$ denote the
Hamiltonian in the bulks on the left and on the right of the interface,
respectively, $\hat{\Psi}=(\hat{\Psi}_{\uparrow},\hat{\Psi}_{\downarrow})^{T}$
is the $2\times 1$ electron spinor field operator, $p_{x}=-i\hbar\partial_{x}$
is the momentum operator, and $\sigma_{z}$ is the Pauli matrix in spin space.
Finally $v_{L/R}$ denotes the Fermi velocity. When $v_{R}$ and $v_{L}$ have
opposite signs, the helicity changes across the interface and the Dirac
paradox emerges. The answer to the paradox, if any, heavily depends on how the
crossover from $\hat{\mathcal{H}}_{L}$ to $\hat{\mathcal{H}}_{R}$ occurs, as
we shall discuss here below considering different models.

###  II.1 Model 1: velocity sign change

The most straightforward way to implement the crossover from
$\hat{\mathcal{H}}_{L}$ to $\hat{\mathcal{H}}_{R}$ is to assume that the
entire system is characterized by an inhomogeneous velocity $v(x)$, varying
from $v_{L}$ to $v_{R}$ over a certain crossover length $\lambda$. Since the
momentum operator $p_{x}$ does not commute with an inhomogeneous velocity
profile $v(x)$, a quite natural approach is to replace their product
$p_{x}\,v$ by a half of their anticommutator, obtaining the following
Hamiltonian

$\hat{\mathcal{H}}=\int\hat{\Psi}^{\dagger}(x)\sigma_{z}\frac{\left\\{v(x)\,,p_{x}\right\\}}{2}\,\hat{\Psi}(x)\,dx\quad.$
(2)

The current operator associated to Eq.(2) is

$\hat{J}(x)={\rm e}\,v(x)\hat{\Psi}^{\dagger}(x)\sigma_{z}\hat{\Psi}(x)\quad,$
(3)

with ${\rm e}$ denoting the electron charge, whereas the Heisenberg Equation
dictated by Hamiltonian (2) reads

$\partial_{t}\hat{\Psi}=-\sigma_{z}\left(v(x)\partial_{x}\hat{\Psi}+\frac{\partial_{x}v}{2}\hat{\Psi}\right)\quad.$
(4)

Looking for stationary solutions
$\hat{\Psi}(x,t)=\hat{\Psi}_{E}(x)e^{-iEt/\hbar}$ and multiplying Eq.(4) by
$\sigma_{z}$ on the left, the equation reduces to

$\partial_{x}\hat{\Psi}=v^{-1}(x)\left(-\frac{\partial_{x}v}{2}\sigma_{0}+i\frac{E}{\hbar}\sigma_{z}\right)\hat{\Psi}$
(5)

whose formal solution is

$\displaystyle\hat{\Psi}_{E}(x)$ $\displaystyle=$
$\displaystyle\displaystyle\exp\left[-\frac{1}{2}\int_{x_{R}}^{x}\frac{\partial_{x}v}{v(x^{\prime})}dx^{\prime}\right]\times$
(6)
$\displaystyle\times\exp\left[iE\sigma_{z}\int_{x_{R}}^{x}\frac{dx^{\prime}}{\hbar
v(x^{\prime})}\right]\,\hat{\Psi}_{E}(x_{R})$

where $x_{R}$ is some arbitrary reference space point. One can now exploit
$\partial_{x}\ln|v(x)|=\partial_{x}v/v(x)$, denote $k_{E}(x)=E/\hbar v(x)$ and
write

$\hat{\Psi}_{E}(x_{R})=\frac{{u}}{\sqrt{2\pi\hbar|v(x_{R})|}}\,\hat{a}_{E}\quad,$
(7)

where ${u}$ is a position-independent $2\times 1$ spinor and $\hat{a}_{E}$ the
related energy-$E$ mode operator fulfilling
$\\{a_{E}\,,a^{\dagger}_{E^{\prime}}\\}=\delta(E-E^{\prime})$. At each energy
$E$ there are thus two independent solutions, corresponding to two mutually
orthogonal choices for the spinor $u$. Then, Eq.(6) takes the form

$\displaystyle\hat{\Psi}_{E}(x)$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2\pi\hbar|v(x)|}}\,e^{i\sigma_{z}\int_{x_{R}}^{x}k_{E}(x^{\prime})dx^{\prime}}\,u\,\hat{a}_{E}$
(8)

which straightforwardly implies that at any space point $x$, including
possible discontinuity points of $v(x)$, the following boundary condition
holds

$\sqrt{|v(x^{+})|}\hat{\Psi}_{E}(x^{+})=\sqrt{|v(x^{-})|}\hat{\Psi}_{E}(x^{-})$
(9)

where $x^{\pm}=x\pm\varepsilon$ with $\varepsilon\rightarrow 0$.

If $v(x)$ varies in magnitude from $v_{L}$ to $v_{R}$ while preserving a
positive sign, $v(x)=|v(x)|$, the Hamiltonian (2) can equivalently be
rewritten as

$\hat{\mathcal{H}}=\int\hat{\Psi}^{\dagger}(x)\sigma_{z}\sqrt{v(x)}\left[p_{x}\sqrt{v(x)}\right]\,\hat{\Psi}(x)\,dx$
(10)

In this case, one finds that the transmission coefficient is always 1,
regardless of the specific values of $v_{L},v_{R}>0$, as discussed in
Ref.[peres, ].

If, however, $v(x)$ vanishes at some point $x_{0}$, like in the Dirac paradox,
the problem becomes more subtle. Indeed in such case the energy dependent
phase factor appearing in the solution Eq.(8) is well defined only if $v(x)$
vanishes as $|v(x)|=O(|x-x_{0}|^{\alpha})$ with $0<\alpha<1$. Moreover the
solution diverges as $\sim 1/\sqrt{|v(x)|}$ for $x\rightarrow x_{0}$. Yet, in
view of the condition (9), the current in Eq.(3) is finite. Denoting by
$\hat{J}_{E}^{\pm}\doteq\hat{J}_{E}(x_{0}^{\pm})$ the current operator for a
stationary solution at energy $E$ at the two sides of the point $x_{0}$ of
vanishing velocity, one straightforwardly finds from Eq.(9) that
$\hat{J}_{E}^{+}=-\hat{J}_{E}^{-}$. For stationary solutions, however, the
continuity equation requires the expectation value of the current to be
continuous and independent of the position. The only possibility is that no
current flows through the system, $\langle\hat{J}_{E}(x)\rangle\equiv
0\,\,\forall x$, implying that the spinor $u$ appearing in Eq.(8) must be
chosen to have vanishing spin along $z$, i.e. $u^{\dagger}\sigma_{z}u=0$. Up
to an overall dimensional coefficient, two independent choices are
$u_{+}=(1\,,e^{i\phi})^{T}/\sqrt{2}$ and $u_{-}=(e^{-i\phi},-1)^{T}/\sqrt{2}$,
where $\phi$ is an arbitrary phase.

As an illustrative example, consider for instance the spatially odd profile
$v(x)=-v_{F}\,\mbox{sgn}(x)\tanh^{\alpha}(|x|/\lambda)$, which describes a
velocity sign change from $v_{L}=+v_{F}$ to $v_{R}=-v_{F}$ across the
interface located at $x_{0}=0$, occurring over a lengthscale $\lambda$ and
with an exponent $0<\alpha<1$. It is straightforward to prove that the
solution (8) is spatially even. Explicitly, choosing e.g. $x_{R}=-4\lambda$ as
a reference point and taking the phase $\phi=0$ in the above spinors $u$, one
finds $k_{E}(x)\simeq-E\,\mbox{sgn}(x)/\hbar v_{F}$ for $|x|\gg\lambda$. The
two physically correct solutions of the Heisenberg Eq.(4) then read

$\hat{\Psi}_{\pm}(x,t)=\int\frac{dE}{{\sqrt{2\pi\hbar
v_{F}}}}\,\psi_{E\pm}(x)\,e^{-iEt/\hbar}\,\hat{a}_{E\pm}\quad,$ (11)

where the wavefunctions $\psi_{E\pm}$ for $|x|\gg\lambda$ take the form

$\psi_{E\pm}(x)=\frac{e^{-iE|x|/\hbar
v_{F}}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
0\end{array}\right)\pm\frac{e^{+iE|x|/\hbar
v_{F}}}{\sqrt{2}}\left(\begin{array}[]{c}0\\\ 1\end{array}\right)\quad.$ (12)

These solutions fulfill the continuity equation by carrying a vanishing
current [see Eq.(3)]. Note that the spatially even wavefunctions in Eq.(12)
involve incoming waves from both sides and cannot be scattering state
solutions. Moreover, any attempt to construct scattering states by their
linear combinations would fail and would also violate the continuity equation.

In summary, the answer to the Dirac paradox provided by model 1 is that, when
$v_{L}$ and $v_{R}$ have opposite signs, it is impossible to construct
scattering state solutions that respect the continuity equation. The
transmission coefficient cannot be properly defined. Physically correct
solutions must necessarily involve incoming waves from both sides and carry no
current, regardless of the specific magnitudes of $|v_{L}|$ and $|v_{R}|$. We
conclude this section by noticing that the model 1 only involves the
$\sigma_{z}$-component of spin [see Eq.(2)]. The space-dependent $v(x)$
changes magnitude and sign of such component. In this respect, the model is
purely scalar.

###  II.2 Model 2: spin-active interface

The second model to approach the Dirac paradox is described by the Hamiltonian

$\hat{\mathcal{H}}=v_{F}\int\hat{\Psi}^{\dagger}(x)\left(e^{-i\theta(x)\sigma_{x}/2}p_{x}\sigma_{z}\,e^{+i\theta(x)\sigma_{x}/2}\right)\hat{\Psi}(x)\,dx$
(13)

where the helicity changes sign through a counter-clockwise rotation of the
$\sigma_{z}$ spin around the $x$-axis by a space-dependent angle $\theta(x)$
varying from $\theta_{L}=0$ to $\theta_{R}=\pi$, over a certain crossover
length. Thus, differently from the purely scalar model (2), the model (13)
exploits the full ${\rm SU}(2)$ spin structure and Hamiltonian terms at two
different points do not commute in general. The current operator related to
the Hamiltonian (13) is

$\displaystyle\hat{J}(x)$ $\displaystyle=$ $\displaystyle{\rm
e}v_{F}\hat{\Psi}^{\dagger}(x)\left(e^{-i\theta(x)\sigma_{x}/2}\sigma_{z}\,e^{+i\theta(x)\sigma_{x}/2}\right)\hat{\Psi}(x)=$
(14) $\displaystyle=$ $\displaystyle{\rm
e}v_{F}\hat{\Psi}^{\dagger}(x)\left[\sigma_{z}\cos\theta(x)-\sigma_{y}\sin\theta(x)\right]\hat{\Psi}(x)\quad.$

Integrating the Heisenberg equation for the field operator

$\partial_{x}\left(e^{i\theta(x)\sigma_{x}/2}\hat{\Psi}(x)\right)=-\frac{\sigma_{z}}{v_{F}}e^{i\theta(x)\sigma_{x}/2}\partial_{t}\hat{\Psi}(x)$
(15)

around any point $x$, including possible discontinuity points of $\theta(x)$,
the following boundary condition is found

$e^{i\theta(x^{+})\sigma_{x}/2}\hat{\Psi}(x^{+})=e^{i\theta(x^{-})\sigma_{x}/2}\hat{\Psi}(x^{-})\quad,$
(16)

which in turn straightforwardly implies the continuity of the current operator
(14). In particular, for a step-like model $\theta(x<0)=0$ and
$\theta(x>0)=\pi$ of an interface located at $x_{0}=0$, Eq.(16) reduces to

$\begin{array}[]{lcl}\hat{\Psi}_{\uparrow}(0^{+})&=&i\hat{\Psi}_{\downarrow}(0^{-})\\\
\hat{\Psi}_{\downarrow}(0^{+})&=&i\hat{\Psi}_{\uparrow}(0^{-})\end{array}$
(17)

and describes a spin-flip process occurring at the interface. Differently from
model 1, the Heisenberg equation (15) does admit scattering state solutions
$\hat{\Psi}_{E,\pm}(x,t)=\exp[-iEt/\hbar]\psi_{E\pm}(x)\hat{a}_{E,\pm}/\sqrt{2\pi\hbar
v_{F}}$, where $\hat{a}_{E,\pm}$ are the energy-$E$ mode operators for
scattering from the left(+) and from right(-), respectively, and
$\psi_{E\pm}(x)$ are the related scattering wavefunctions. For instance, the
scattering state from left is given by

$\begin{array}[]{lcl}\psi_{E+}(x<0)&=&\left(\begin{array}[]{c}1\\\
0\end{array}\right)e^{iEx/\hbar v_{F}}+r\left(\begin{array}[]{c}0\\\
1\end{array}\right)e^{-iEx/\hbar v_{F}}\\\
\psi_{E+}(x>0)&=&t\left(\begin{array}[]{c}0\\\ 1\end{array}\right)e^{iEx/\hbar
v_{F}}\end{array}$ (18)

and, when inserted in Eq.(17), straightforwardly implies $r=0$ and $t=i$,
leading to a perfect transmission $T=|t|=1$.

In summary, model 2 trivially solves the Dirac paradox by simply introducing
spin-flip processes at the interface.

##  III Dirac Heterojunctions with massless and massive modes

So far, we have considered heterojunctions that purely involve Dirac massless
modes and we have obtained two opposite answers to the Dirac paradox,
depending on how the helicity change across the interface is modelled. Model
1, based on an inhomogeneous scalar velocity profile, implies that physical
solutions necessarily involve injections from both sides of the junction and
predicts no current flowing through the system, whereas model 2 “tackles” the
paradox by introducing an interface with spin-flip processes that guarantee a
perfect transmission. In this section we propose a model that, without
introducing spin-flip processes at the interface, leads to non-vanishing
transmission.

Figure 2:  The band spectrum of the massless+massive Dirac model with opposite
helicity across an interface. While blue and red curves denote the spectrum of
the massless modes, as in Fig.1, the green curves describe the spectrum of the
massive modes, characterized by a gap $2\Delta$.

Suppose that, along with the massless propagating Dirac fermions illustrated
in Fig.1, the system is also characterized by massive Dirac fermions, as
sketched by the green curves of Fig.2. Specifically, the model we consider is

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$ $\displaystyle
v_{F}\int\hat{\Psi}^{\dagger}(x)U^{\dagger}(x)\tau_{z}\sigma_{z}p_{x}\left(U(x)\,\hat{\Psi}(x)\right)\,dx+$
(19)
$\displaystyle-\frac{\Delta}{2}\int\hat{\Psi}^{\dagger}(x)(\tau_{0}-\tau_{z})\sigma_{x}\,\hat{\Psi}(x)\,dx$

where
$\hat{\Psi}=(\hat{\xi}_{\uparrow},\hat{\xi}_{\downarrow},\hat{\eta}_{\uparrow},\hat{\eta}_{\downarrow})^{T}$,
with $\hat{\xi}_{\uparrow},\hat{\xi}_{\downarrow}$ and
$\hat{\eta}_{\uparrow},\hat{\eta}_{\downarrow}$ denoting the massless and
massive fields, respectively. Here $\sigma_{0}$ and
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ denote the $2\times 2$
identity matrix and Pauli matrices acting on the spin space, whereas
$\tau_{0}$ and $\bm{\tau}=(\tau_{x},\tau_{y},\tau_{z})$ the corresponding
quantities acting on the massless-massive degree of freedom, which we shall
label as pseudospin. In the first term of Eq.(19) the $4\times 4$ matrix
$U(x)$ interpolates from $U_{L}$ on the left of the interface to its value
$U_{R}$ on the right, where $U_{L/R}$ are required to fulfill the following
properties

$\displaystyle U^{\dagger}_{L}\tau_{z}\sigma_{z}U_{L}$ $\displaystyle=$
$\displaystyle+\tau_{z}\sigma_{z}$ (20) $\displaystyle
U^{\dagger}_{R}\tau_{z}\sigma_{z}U_{R}$ $\displaystyle=$
$\displaystyle-\tau_{z}\sigma_{z}\quad,$ (21)

so that the $\hat{\xi}_{\uparrow},\hat{\xi}_{\downarrow}$ modes have helicity
$+1$ on the left of the interface and $-1$ on the right, just like in the
Dirac paradox configuration of Fig.1, whereas the opposite occurs for the
$\hat{\eta}_{\uparrow},\hat{\eta}_{\downarrow}$ modes. The simplest example of
a $U(x)$-matrix fulfilling the conditions (20)-(21) is
$U(x)=\exp[i\theta(x)\tau_{x}\sigma_{0}/2]$, where $\theta(x)$ is a space-
dependent angle describing a rotation in pseudospin space around $\tau_{x}$
from $\theta_{L}=0$ to $\theta_{R}=\pi$ and causing the helicity flip, just
like the spin-active model (13) introduces a rotation in spin space. As we
shall see below, there exists in fact a much broader set of possible choices
for $U(x)$ that turn out to describe interesting and realistic cases. The
second term in Eq.(19) describes the mass term for $\hat{\eta}_{\uparrow}$ and
$\hat{\eta}_{\downarrow}$, and we shall be interested in the energy window
$|E|<\Delta$ inside their gap, where these massive modes carry no current.

In view of Eqs.(20)-(21), the current operator related to the Hamiltonian (19)

$\displaystyle\hat{J}(x)$ $\displaystyle=$ $\displaystyle{\rm
e}v_{F}\hat{\Psi}^{\dagger}(x)U^{\dagger}(x)\tau_{z}\sigma_{z}\,U(x)\hat{\Psi}(x)$
(22)

takes opposite expressions $\hat{J}_{L/R}=\pm{\rm
e}v_{F}\hat{\Psi}^{\dagger}(x)\tau_{z}\sigma_{z}\,\hat{\Psi}(x)$ at the two
sides of the interface. However, the boundary condition

$U(x^{+})\hat{\Psi}(x^{+})=U(x^{-})\hat{\Psi}(x^{-})\quad,$ (23)

obtained from integration of the Heisenberg equation around any point $x$,
guarantees that the current is in fact continuous for any $U(x)$. In
particular, adopting again a step-like model $U(x)=U_{L}\,{\rm
H}(-x)+U_{R}\,{\rm H}(x)$ for an interface located at $x_{0}=0$, with ${\rm
H}(x)$ denoting the Heaviside function, the field $\Psi$ fulfills the
interface boundary condition

$\hat{\Psi}(0^{+})=\mathsf{M}\,\hat{\Psi}(0^{-})\quad,$ (24)

where

$\mathsf{M}=U^{-1}_{R}U_{L}$ (25)

is the transfer matrix, which must fulfill

$\displaystyle\mathsf{M}^{\dagger}\tau_{z}\sigma_{z}\mathsf{M}$
$\displaystyle=$ $\displaystyle-\tau_{z}\sigma_{z}\hskip
28.45274pt\mbox{\small({\it Requirement \\#1})}$ (26)

as a straightforward consequence of Eqs.(20)-(21). Note that Eq.(24) implies
that the field $\hat{\Psi}$ is discontinuous, as is customary for Dirac models
in the presence of a $\delta(x)$-term, which in this case originates from
$p_{x}U(x)$ term in the Hamiltonian (19).

Importantly, in order to avoid trivial solutions to the Dirac paradox like in
model (13), we require that the model (19) does not directly introduce any
spin-flip process at the interface. This leads to impose another requirement
on the transfer matrix Eq.(24), namely that $\mathsf{M}$ is diagonal in spin
space, i.e.

$\begin{array}[]{l}\mathsf{M}\,\,\,\mbox{must involve}\\\ \mbox{only
$\sigma_{0}$ and $\sigma_{z}$}\end{array}\hskip 28.45274pt\mbox{\small({\it
Requirement \\#2})}\quad.$ (27)

It can be shown (see Appendix A for details) that the most general matrix
fulfilling the requirements Eqs.(26)-(27) has the following form in the
$\tau\otimes\sigma$ basis

$\displaystyle\mathsf{M}=\left(\begin{array}[]{cccc}i\beta_{\uparrow}\,e^{i(\nu_{\uparrow}-\gamma_{\uparrow})}&0&(1-i\beta_{\uparrow})e^{i(\nu_{\uparrow}+\chi_{\uparrow})}&0\\\
0&i\beta_{\downarrow}e^{i(\nu_{\downarrow}-\gamma_{\downarrow})}&0&(1-i\beta_{\downarrow})e^{i(\nu_{\downarrow}+\chi_{\downarrow})}\\\
(1+i\beta_{\uparrow})e^{i(\nu_{\uparrow}-\chi_{\uparrow})}&0&-i\beta_{\uparrow}e^{i(\nu_{\uparrow}+\gamma_{\uparrow})}&0\\\
0&(1+i\beta_{\downarrow})e^{i(\nu_{\downarrow}-\chi_{\downarrow})}&0&-i\beta_{\downarrow}e^{i(\nu_{\downarrow}+\gamma_{\downarrow})}\end{array}\right)$
(32)

and depends on 8 parameters, namely 4 real parameters
$\chi_{\sigma},\gamma_{\sigma},\beta_{\sigma},\nu_{\sigma}$ for each spin
sector $\sigma=\uparrow,\downarrow$. The vanishing entries in Eq.(32) encode
the decoupling of the two spin sectors dictated by Eq.(27).

###  III.1 Scattering states

Let us now focus on $E=0$, i.e. on the middle of the massive energy gap, and
build up scattering state solutions on both sides of the junction, namely

$\displaystyle\hat{\Psi}(x<0)=$ (33) $\displaystyle=$
$\displaystyle\hat{a}_{L\uparrow}\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\
0\end{array}\right)e^{ik_{0}x}\,+\hat{b}_{L\downarrow}\left(\begin{array}[]{c}0\\\
1\\\ 0\\\
0\end{array}\right)e^{-ik_{0}x}+\frac{\hat{c}_{L}}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
0\\\ -i\\\ 1\end{array}\right)e^{\kappa_{0}x}$ (46)

and

$\displaystyle\hat{\Psi}(x>0)=$ (47) $\displaystyle=$
$\displaystyle\hat{a}_{R\uparrow}\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\
0\end{array}\right)e^{-ik_{0}x}+\hat{b}_{R\downarrow}\left(\begin{array}[]{c}0\\\
1\\\ 0\\\
0\end{array}\right)e^{ik_{0}x}+\frac{\hat{c}_{R}}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
0\\\ -i\\\ 1\end{array}\right)e^{-\kappa_{0}x}$ (60)

where $k_{0}=0$, $\kappa_{0}=\Delta/\hbar v_{F}$. Here $\hat{a}_{L\uparrow}$
and $\hat{a}_{R\uparrow}$ are incoming operators describing a propagating mode
impinging from the left(L) and from the right(R) of the interface,
respectively, whereas $\hat{b}_{L\downarrow}$ and $\hat{b}_{R\downarrow}$ are
outgoing operators for modes propagating to the left and to the right,
respectively. Note that in the Dirac paradox configuration (see Fig.2)
incoming states and outgoing states have opposite spin, namely spin-$\uparrow$
and spin-$\downarrow$, respectively. Furthermore in Eqs.(LABEL:PsiL) and
(LABEL:PsiR) $\hat{c}_{L}$ and $\hat{c}_{R}$ describe evanescent modes on the
left- and on the right-hand side of the interface. Importantly, because they
are massive, their spinors are characterized by both spin components.
Introducing Eqs.(LABEL:PsiL) and (LABEL:PsiR) into Eq.(24) and using Eq.(32),
one can write

$\left(\begin{array}[]{c}\hat{b}_{L\downarrow}\\\ \hat{b}_{R\downarrow}\\\
\hat{c}_{L}\\\ \hat{c}_{R}\end{array}\right)=\begin{pmatrix}&\mathsf{S}&\\\
&\tilde{\mathsf{S}}&\end{pmatrix}\left(\begin{array}[]{c}\hat{a}_{L\uparrow}\\\
\hat{a}_{R\uparrow}\end{array}\right)$ (61)

where $\mathsf{S}$ denotes the Scattering Matrix returning the outgoing
propagating modes

$\displaystyle\mathsf{S}=\frac{ie^{-i\Delta\chi}}{(1-i\beta_{\uparrow})(1+i\beta_{\downarrow})}\times$
$\displaystyle\begin{pmatrix}e^{i\Delta\nu}+e^{-i\Delta\gamma}\beta_{\uparrow}\beta_{\downarrow}&i(e^{i(\gamma_{\downarrow}-\nu_{\uparrow})}\beta_{\downarrow}-e^{i(\gamma_{\uparrow}-\nu_{\downarrow})}\beta_{\uparrow})\\\
&\\\
i(e^{i(\nu_{\uparrow}-\gamma_{\downarrow})}\beta_{\downarrow}-e^{i(\nu_{\downarrow}-\gamma_{\uparrow})}\beta_{\uparrow})&e^{-i\Delta\nu}+e^{i\Delta\gamma}\beta_{\uparrow}\beta_{\downarrow}\end{pmatrix}$

with $\Delta\chi\doteq\chi_{\uparrow}-\chi_{\downarrow}$,
$\Delta\nu\doteq\nu_{\uparrow}-\nu_{\downarrow}$ and
$\Delta\gamma\doteq\gamma_{\uparrow}-\gamma_{\downarrow}$, whereas

$\displaystyle\tilde{\mathsf{S}}=\frac{\sqrt{2}\,e^{-i\chi_{\uparrow}}}{1-i\beta_{\uparrow}}\begin{pmatrix}\beta_{\uparrow}e^{-i\gamma_{\uparrow}}&ie^{-i\nu_{\uparrow}}\\\
&\\\
ie^{+i\nu_{\uparrow}}&\beta_{\uparrow}e^{i\gamma_{\uparrow}}\end{pmatrix}$
(63)

is the matrix yielding the evanescent modes.

In Eq.(61), setting $\hat{a}_{R\uparrow}\rightarrow 0$ yields a scattering
state with injection from left, while a scattering state with injection from
right is obtained for $\hat{a}_{L\uparrow}\rightarrow 0$. Thus, differently
from model 1 in Eq.(2), the model in Eq.(19) does allow for scattering
solutions. The transmission coefficient $T_{0}=|t_{0}|^{2}$, obtained from the
off-diagonal entries of the Scattering Matrix (LABEL:Smat), reads

$T_{0}=\frac{\beta_{\uparrow}^{2}+\beta_{\downarrow}^{2}-2\beta_{\uparrow}\beta_{\downarrow}\cos\varphi}{(1+\beta_{\uparrow}^{2})(1+\beta_{\downarrow}^{2})}\quad,$
(64)

and depends on the three parameters $\beta_{\uparrow}$, $\beta_{\downarrow}$
and $\varphi=\Delta\gamma+\Delta\nu$. Moreover, at variance with model 2, here
the transfer matrix in Eqs.(24) and (32) is diagonal in spin, and the
transmission from a spin-$\uparrow$ incoming state to a spin-$\downarrow$
outgoing state occurs indirectly, i.e. thanks to the presence of the
evanescent massive modes at the interface.

To illustrate this effect, let us consider e.g. the scattering from left,
setting $\hat{a}_{R\uparrow}\rightarrow 0$ in Eq.(61). For simplicity we can
also set all phases to zero ($\gamma_{\sigma}=\chi_{\sigma}=\nu_{\sigma}=0$)
in Eqs.(32), (LABEL:Smat) and (63). We first focus on the case
$\beta_{\uparrow}=0$, where the scattering state resulting from
Eqs.(61)-(LABEL:Smat)-(63) is sketched on the left-hand side of Fig.3(a): The
blue (red) wiggy line describes the incoming spin-$\uparrow$ state (outgoing
spin-$\downarrow$ states), while an evanescent wave (green solid line) is
present only for $x>0$. Its role is elucidated on the right hand side of
Fig.3(a), which is a graphical representation of Eq.(24) where the non-
vanishing components of such a scattering state are connected across the
interface by the transfer matrix entries (black lines). When the massless
spin-$\uparrow$ state propagates towards the interface from the left, the
transfer matrix Eq.(32) connects it through the entry $\mathsf{M}_{31}=1$ to
its massive evanescent partner with the same spin located across the
interface, represented by a green dashed box, with the thick solid lines
inside it denoting its two spin components. Because such a mode is massive,
inside the gap it always exhibits both spin components [see third term in
Eq.(LABEL:PsiR)]. Thus, its spin-$\downarrow$ component is also present and is
connected through the transfer matrix entry
$\mathsf{M}_{42}=1+i\beta_{\downarrow}$ to its spin-$\downarrow$ massless
partner, which describes the reflected wave propagating to the left of the
junction. Finally, the latter is also coupled, through the entry
$\mathsf{M}_{22}=i\beta_{\downarrow}$, to the massless spin-$\downarrow$ state
outgoing to the right of the junction. Thus, despite the interface connects
only states with the same spin on the two sides, the presence of an evanescent
massive mode exhibiting both spin components leads to an effective spin-flip
transmission between massless modes.

Let us now consider the case $\beta_{\downarrow}=0$. In this case the
scattering state resulting from the solution Eqs.(61)-(LABEL:Smat)-(63)
exhibits evanescent modes on both sides of the junction, as sketched in the
left hand side of Fig.3(b). The scheme on the right hand side of Fig.3(b)
illustrates the related Eq.(24). While the entry $\mathsf{M}_{31}$ is modified
to $\mathsf{M}_{31}=1+i\beta_{\uparrow}$, a connection
$\mathsf{M}_{33}=-i\beta_{\uparrow}$ opens up across the junction between the
two spin-$\uparrow$ components of the massive modes. In turn, their
corresponding spin-$\downarrow$ components are connected through the entries
$\mathsf{M}_{24}=\mathsf{M}_{42}=1$ to the spin-$\downarrow$ massless modes
across the junction, thereby inducing again spin-flipped reflection and
transmission.

Figure 3:  For each panel, the left-hand side sketches the scattering state
wavefunction in the case of injection from the left, resulting from
Eqs.(61)-(LABEL:Smat)-(63) for $\hat{a}_{R\uparrow}\rightarrow 0$. Blue and
red wiggy lines describe spin-$\uparrow$ and spin-$\downarrow$ propagating
massless states, respectively, whereas solid green lines describe the
evanescent wave of the massive mode. The right-hand side of each panel is a
graphical representation of Eq.(24), where black lines represent the transfer
matrix entries connecting the non vanishing components of such scattering
state. (a) The case with $\beta_{\uparrow}=0$. The evanescent mode is present
only on the right side of the interface. Here $\mathsf{M}_{31}=1$,
$M_{42}=1+i\beta_{\downarrow}$ and $\mathsf{M}_{22}=i\beta_{\downarrow}$. (b)
The case with $\beta_{\downarrow}=0$. In this case the evanescent modes are
present on both sides of the junction. Here
$\mathsf{M}_{31}=1+i\beta_{\uparrow}$, $M_{33}=-i\beta_{\uparrow}$ and
$\mathsf{M}_{24}=\mathsf{M}_{42}=1$. In all cases, despite the transfer matrix
only connects states with the same spin, the presence of the evanescent modes
of the massive field enables a spin-flip transmission between the propagating
modes.

The general case, where both $\beta_{\uparrow}$ and $\beta_{\downarrow}$ are
non vanishing, is a combination of the two elementary cases and yields the
transmission coefficient (64). Note that in the limit where both
$\beta_{\uparrow}\rightarrow 0$ and $\beta_{\downarrow}\rightarrow 0$, the
transmission coefficient (64) vanishes. This can also be understood by
realizing that in such limit the transfer matrix (32) reduces to
$\mathsf{M}=\tau_{x}\sigma_{0}$, yielding the boundary conditions

$\left\\{\begin{array}[]{lcl}\hat{\xi}_{\sigma}(0^{+})&=&\hat{\eta}_{\sigma}(0^{-})\\\
\hat{\eta}_{\sigma}(0^{+})&=&\hat{\xi}_{\sigma}(0^{-})\end{array}\right.\hskip
28.45274pt\sigma=\uparrow,\downarrow$ (65)

so that e.g. a massless mode incoming from the left towards the interface is
completely transformed into its massive evanescent mode partner across the
interface (with the same spin), which carries no current.

In summary, although massive modes do not carry any current inside the gap,
their presence is important in inhomogeneous problems because they may
localize at the interfaces. In particular they are crucial in the Dirac
paradox, for they provide an indirect coupling between the two spin channels
that would be otherwise uncoupled by the interface transfer matrix. This leads
to an effective spin-flip transmission of the massless propagating modes.
Moreover, in contrast with the models 1 and 3 discussed in Sec.II, here the
transmission coefficient is tunable from 0 to 1 through the 3 knobs
$\beta_{\uparrow},\beta_{\downarrow}$ and $\varphi$. This is one of the main
results of our paper.

##  IV Spin-orbit coupled nanowires

In this section we shall describe how the model presented in Sec.III can be
realized with spin-orbit coupled NWs, under suitable circumstances. First, we
shall briefly recall how these systems, when exposed to an external magnetic
field, can host helical states described by Dirac massless fermions, as well
as gapped Dirac states. Then, focussing on energies inside the gap opened up
by the magnetic field, we shall explicitly derive the effective low-energy
model for these systems. Finally, we shall consider an inhomogeneous spin-
orbit coupling profile that, in suitable regimes, realizes the Dirac paradox
configuration involving both massless and massive modes, just like in the
model proposed above.

###  IV.1 The NW Hamiltonian and its low energy limit

We consider a single-channel semiconductor NW deposited on a substrate. For
NWs like InSb or InAs, the structural inversion asymmetry can lead to quite
strong RSOCnilsson_2009 ; deng_2012 ; kouwenhoven_PRL_2012 ; ensslin_2010 ;
gao_2012 ; joyce_2013 , which can further be tuned with appropriate gating
techniques. Furthermore we assume that a uniform magnetic field is applied
parallel to the nanowire axis, denoted by $x$, while the substrate plane will
be identified as $x$-$z$.

We shall adopt a widely used model to describe the NWstreda ; vonoppen_2010 ;
dassarma_2010 ; loss_PRB_2011 ; lutchyn_2012 ; loss_PRB_2017 , whose main
ingredients are summarized here below, while details are reported in Appendix
B for the sake of completeness. The second-quantized NW Hamiltonian consists
of three terms
$\hat{\mathcal{H}}_{NW}=\hat{\mathcal{H}}_{kin}+\hat{\mathcal{H}}_{R}+\hat{\mathcal{H}}_{Z}$
and can be written as
$\hat{\mathcal{H}}_{NW}=\int\hat{\Phi}^{\dagger}(x)\,H_{NW}(x)\,\hat{\Phi}(x)\,dx$.
Here
$\hat{\Phi}(x)=(\hat{\Phi}_{\uparrow}(x)\,,\,\hat{\Phi}_{\downarrow}(x))^{T}$
is the electron spinor field, with $\uparrow,\downarrow$ corresponding to spin
projections along $z$, and

$H_{NW}(x)=\frac{p_{x}^{2}}{2m^{*}}\sigma_{0}-\frac{\alpha}{\hbar}p_{x}\sigma_{z}\,\,-h_{x}\sigma_{x}$
(66)

contains the kinetic term characterized by an effective mass $m^{*}$, the
Rashba term with a RSOC $\alpha$, and the Zeeman term describing the coupling
$h_{x}=g\mu_{B}B_{x}/2$ with the external magnetic field
$\mathbf{B}=(B_{x},0,0)$, with $\mu_{B}$ denoting the Bohr magneton and $g$
the NW Landé factor. The model is characterized by two energy scales, namely
the spin-orbit energy

$E_{SO}=\frac{m^{*}\alpha^{2}}{2\hbar^{2}}\quad,$ (67)

and the Zeeman energy

$E_{Z}=|h_{x}|$ (68)

For definiteness, we shall henceforth assume $h_{x}>0$ and identify
$h_{x}=E_{Z}$. The spin-orbit wavevector

$k_{SO}=\frac{\sqrt{{2m^{*}E_{SO}}}}{\hbar}=\frac{|\alpha|m^{*}}{\hbar^{2}}\quad,$
(69)

and the Zeeman wavevector

$k_{Z}=\sqrt{{2m^{*}E_{Z}}}\,/\hbar\quad.$ (70)

are the wavevectors associated to such energies. Diagonalizing the model in
momentum space, one obtains a spectrum characterized by two energy bands
separated at $k=0$ by a gap $2E_{Z}$ centered around the midgap energy $E=0$.
Depending on the ratio of $E_{Z}$ to $E_{SO}$, the qualitative behavior of
these bands is different. In fact, two regimes can be identified: (a) for
$E_{Z}>2E_{SO}$ (Zeeman-dominated regime) both bands have a minimum at $k=0$
taking values $E_{\pm}^{\rm min}=\pm E_{Z}$, whereas (b) for $E_{Z}<2E_{SO}$
(Rashba-dominated regime) the upper band still has a minimum $E_{+}^{\rm
min}=+E_{Z}$ at $k=0$, while the lower band has a local maximum at $k=0$ and
acquires two degenerate minima $E_{-}^{\rm min}=-E_{SO}-E_{Z}^{2}/4E_{SO}$ at
$k=\pm k_{SO}\sqrt{1-E_{Z}^{2}/4E_{SO}^{2}}$.

In the following we shall focus on the deep Rashba-dominated regime ($E_{Z}\ll
2E_{SO}$), illustrated in Fig.4, and analyze the energy window inside the
magnetic gap ($|E|\ll E_{Z}$), highlighted by the dashed box. As is well
known, in this range the NW propagating eigenstates are helicalstreda ;
vonoppen_2010 ; dassarma_2010 ; loss_PRB_2011 ; lutchyn_2012 ; loss_PRB_2017 ;
depicciotto_2010 ; kouwenhoven_nanolett_2013 ; kouwenhoven_natcom_2017 : Their
dispersion relation is well described by a linear behavior near the Fermi
points $k\simeq\pm 2k_{SO}$, while their spin orientation, mainly dictated by
the Rashba term, is locked to the propagation direction. For $\alpha>0$,
right-moving electrons near the right Fermi point $k\simeq+2k_{SO}$ are
characterized by spin-$\uparrow$, while left-moving electrons near the left
Fermi point $k\simeq-2k_{SO}$ have spin-$\downarrow$ (see Fig.4). The opposite
occurs if $\alpha<0$. The dynamics of these low energy propagating modes,
which we shall denote by $\hat{\xi}_{\uparrow}$, $\hat{\xi}_{\downarrow}$, is
thus described by a massless Dirac Hamiltonian, whose helicity is determined
by the sign of the RSOC $\alpha$,

$s_{\alpha}=\mbox{sgn}(\alpha)\quad.$ (71)

This suggests that a junction between two NW regions with opposite values of
RSOC realizes the Dirac paradox configuration.

Figure 4:  In the deep Rashba-dominated regime ($2E_{SO}\gg E_{Z}$) and for
energies $|E|\ll E_{Z}$, a NW exposed to an external magnetic field exhibits
helical states near the Fermi points $\pm 2k_{SO}$. Their spin orientation is
locked to the propagation direction and is determined by the sign of the RSOC
$\alpha$. The spin orientation is shown for the case $\alpha>0$. The low
energy massless and massive modes are highlighted with the same colors as in
Fig.4.

However, as highlighted by the green lines in Fig.4, the NW also exhibits low-
energy gapped modes near $k=0$, which shall be denoted as
$\hat{\eta}_{\uparrow}$ and $\hat{\eta}_{\downarrow}$. Notably, these modes
turn out to behave as massive Dirac fermions with a mass term $\Delta=E_{Z}$.
In problems involving homogeneous NWs these modes are dropped because they are
not normalizable. However, as observed in Sec.III, in inhomogeneous problems
such as the Dirac paradox configuration they describe evanescent waves that,
despite carrying no current, ensure the wavefunction matching at the
interface. For these reasons, the effective low energy theory capturing the
physical properties of the Dirac paradox configuration realized with NWs is a
Dirac model involving both massless and massive modes.

To develop such effective theory describing low energy excitations $|E|\ll
E_{Z}\ll 2E_{SO}$, we assume that the ground state is the Fermi sea where all
NW states below the midgap energy ($E<0$) are occupied, and we perform an
expansion near the points $k\simeq\pm 2k_{SO}$ and $k\simeq 0$. It is possible
to show (details can be found in Appendix B), that the low energy excitations
of the NW Hamiltonian are equivalent to low energy excitations of the
massless+massive Dirac model

$\displaystyle\hat{\mathcal{H}}_{NW}=\sum_{q=-\infty}^{+\infty}(\hat{\xi}_{q\uparrow}^{\dagger}\
\hat{\xi}_{q\downarrow}^{\dagger})\begin{pmatrix}\hbar s_{\alpha}v_{so}q&&0\\\
0&&-\hbar
s_{\alpha}v_{SO}q\end{pmatrix}\left(\begin{array}[]{c}\hat{\xi}_{q\uparrow}\\\
\hat{\xi}_{q\downarrow}\end{array}\right)$ (74)
$\displaystyle\displaystyle+\sum_{q=-\infty}^{+\infty}(\hat{\eta}_{q\uparrow}^{\dagger}\
\hat{\eta}_{q\downarrow}^{\dagger})\begin{pmatrix}-\hbar
s_{\alpha}v_{so}q&&-E_{Z}\\\ -E_{Z}&&\hbar
s_{\alpha}v_{SO}q\end{pmatrix}\left(\begin{array}[]{c}\hat{\eta}_{q\uparrow}\\\
\hat{\eta}_{q\downarrow}\end{array}\right)\hskip 28.45274pt$ (77)

where $v_{SO}=\hbar k_{SO}/m^{*}$. Introducing the low-energy fields
($\sigma=\uparrow,\downarrow$)

$\hat{\xi}_{\sigma}(x)=\frac{1}{\sqrt{\Omega}}\sum_{q}\hat{\xi}_{q,\sigma}e^{iqx}\hskip
8.5359pt,\hskip
8.5359pt\hat{\eta}_{\sigma}(x)=\frac{1}{\sqrt{\Omega}}\sum_{q}\hat{\eta}_{q,\sigma}e^{iqx}$
(78)

that physically describe excitations varying over lengthscales much longer
than the spin-orbit length $l_{SO}=k_{SO}^{-1}$, the NW Hamiltonian can be
expressed as

$\displaystyle\hat{\mathcal{H}}_{NW}=$ $\displaystyle=$ $\displaystyle\int
dx\hat{\Psi}^{\dagger}(x)\left(s_{\alpha}v_{SO}\tau_{z}\sigma_{z}p_{x}-\frac{E_{Z}}{2}(\tau_{0}-\tau_{z})\sigma_{x}\right)\hat{\Psi}(x)$

where
$\hat{\Psi}(x)=(\hat{\xi}_{\uparrow},\hat{\xi}_{\downarrow},\hat{\eta}_{\uparrow},\hat{\eta}_{\downarrow})^{T}$
is a 4-component spinor field. One can now realize the connection between the
NW Hamiltonian (LABEL:HNW-lowene) and the model introduced in Sec.III in
Eq.(19). Indeed, identifying $v_{SO}\rightarrow v_{F}$ and
$E_{Z}\rightarrow\Delta$, Eq.(LABEL:HNW-lowene) describes one side of the
junction model (19), where the sign $s_{\alpha}$ of the RSOC [see Eq.(71)]
implements the condition Eq.(20) or (21) and determines which side of the
junction is described.

Finally, the original field $\hat{\Phi}$ can be expressed in terms of the
Dirac slowly varying modes $(\hat{\xi},\hat{\eta})$ and the fast oscillating
plane waves related to the midgap Dirac points, as follows

$\left(\begin{array}[]{c}\hat{\Phi}_{\uparrow}(x)\\\
\hat{\Phi}_{\downarrow}(x)\end{array}\right)=\left(\begin{array}[]{c}e^{+2is_{\alpha}k_{SO}x}\,\hat{\xi}_{\uparrow}(x)+\hat{\eta}_{\uparrow}(x)\\\
e^{-2is_{\alpha}k_{SO}x}\,\hat{\xi}_{\downarrow}(x)+\hat{\eta}_{\downarrow}(x)\end{array}\right)\quad.$
(80)

###  IV.2 The case of inhomogeneous RSOC

Because the helicity of the NW low energy massless modes is determined by the
sign of the RSOC, one can envisage a setup where two different NW portions are
characterized by values of $\alpha$ with opposite signs, as illustrated in
Fig.5. Indeed the huge advances of gating techniques enable one to realize
different gate potentials to various portions of the nanowiressasaki_2013 ;
micolich ; sasaki_2017 ; das_2019 ; guo_2021 ; sasaki_2021 , thereby locally
varying the magnitude and even the sign of the RSOCgao_2012 ; slomski_NJP_2013
; wimmer_2015 ; nygaard_2016 ; sherman_2016 ; tokatly_PRB_2017 ; loss_2018 ;
goldoni_2018 ; tsai_2018 ; gao-review ; lau_2021 ; kaindl_2005 ; meier-
dil_2013 ; wang-fu_2016 ; nitta-frustaglia . The overall system can thus be
described by a inhomogeneous spin-orbit coupling $\alpha(x)$ and the
Hamiltonian (66) is generalized tosanchez_2006 ; sanchez_2008 ; sherman_2011 ;
sherman_2013 ; sherman_2017 ; loss_EPJB_2015 ; dolcini-rossi_PRB_2018 ;
lorenzo_2020 ; rossi_EPJB

$H(x)=\frac{p_{x}^{2}}{2m^{*}}\sigma_{0}-\frac{\left\\{\alpha(x),p_{x}\right\\}}{2\hbar}\sigma_{z}\,\,-h_{x}\sigma_{x}$
(81)

where the anticommutator is necessary since $p_{x}$ does not commute with the
space-dependent RSOC. In particular, as an elementary building block, one can
consider a step-like RSOC profile $\alpha(x)=\alpha_{L}\,{\rm
H}(x_{0}-x)+\alpha_{R}\,{\rm H}(x-x_{0})$ describing an interface located at
$x_{0}$ between two regions with RSOC equal to $\alpha_{L}$ and $\alpha_{R}$.
In such a configuration one can straightforwardly derive the following
matching conditionslorenzo_2020

$\begin{cases}\hat{\Phi}_{\uparrow}(x_{0}^{-})=\hat{\Phi}_{\uparrow}(x_{0}^{+})\\\
\hat{\Phi}_{\downarrow}(x_{0}^{-})=\hat{\Phi}_{\downarrow}(x_{0}^{+})\\\
\partial_{x}\hat{\Phi}_{\uparrow}(x_{0}^{-})=\partial_{x}\hat{\Phi}_{\uparrow}(x_{0}^{+})-i\frac{m^{*}}{\hbar^{2}}(\alpha_{R}-\alpha_{L})\hat{\Phi}_{\uparrow}(x_{0})\\\
\partial_{x}\hat{\Phi}_{\downarrow}(x_{0}^{-})=\partial_{x}\hat{\Phi}_{\downarrow}(x_{0}^{+})+i\frac{m^{*}}{\hbar^{2}}(\alpha_{R}-\alpha_{L})\hat{\Phi}_{\downarrow}(x_{0})\,.\end{cases}$
(82)

This provides all the ingredients for a concrete implementation of the Dirac
paradox. In order to be more realistic, we shall consider a three-region
configuration where the RSOC varies as

$\alpha(x)=\left\\{\begin{array}[]{lcll}+\alpha>0&\mbox{for}&x<-L/2&\mbox{\small(region
1)}\\\ 0&\mbox{for}&|x|<L/2&\mbox{\small(region 2)}\\\
-\alpha<0&\mbox{for}&x>+L/2&\mbox{\small(region 3)}\end{array}\right.$ (83)

where the two outer regions 1 and 3 with opposite RSOC are both assumed in the
deep Rashba-dominated regime ($2E_{SO}\gg E_{Z}$), and are separated by the
central crossover region 2 with length $L$ and with vanishing RSOC, i.e. in
the Zeeman-dominated regime (see Fig.5).

Figure 5: By gating different portions of the NW with metallic electrodes an
inhomogeneous RSOC like (83) can be realized. Inside the Zeeman gap, induced
by an external magnetic field applied along the nanowire axis, the helical
states in the two outer regions have opposite helicity and thus realize the
Dirac paradox configuration.

Applying the general interface condition (82) to the two interfaces
$x_{1}=-L/2$ and $x_{2}=+L/2$ of the piecewise constant profile (83), one can
match the NW wavefunction in the three regions and obtain the solution for the
NW scattering problem with standard techniquesdatta-book  (see Appendix C).
Although the resulting transmission coefficient is numerically exact and
available for arbitrary values of $E$, $E_{SO}$ and $E_{Z}$, it is not quite
amenable. However, in the energy window $|E|\ll E_{Z}\ll 2E_{SO}$ where the
Dirac paradox emerges, an analytical expression can be gained from the
effective low energy model. To this purpose, one can insert the expression
(80) for the field $\hat{\Phi}$ in the outer Rashba-dominated regions into the
interface condition (82) and obtain the low energy boundary conditions at the
left interface $x_{1}=-L/2$

$\left\\{\begin{array}[]{l}e^{-ik_{SO}L}\hat{\xi}_{\uparrow}(x_{1}^{-})+\hat{\eta}_{\uparrow}(x_{1}^{-})=\hat{\Phi}_{\uparrow}(x_{1}^{+})\\\
e^{+ik_{SO}L}\hat{\xi}_{\downarrow}(x_{1}^{-})+\hat{\eta}_{\downarrow}(x_{1}^{-})=\hat{\Phi}_{\downarrow}(x_{1}^{+})\\\
+ik_{SO}\left[e^{-ik_{SO}L}\hat{\xi}_{\uparrow}(x_{1}^{-})-\hat{\eta}_{\uparrow}(x_{1}^{-})\right]=\partial_{x}\hat{\Phi}_{\uparrow}(x_{1}^{+})\\\
-ik_{SO}\left[e^{+ik_{SO}L}\hat{\xi}_{\downarrow}(x_{1}^{-})-\hat{\eta}_{\downarrow}(x_{1}^{-})\right]=\partial_{x}\hat{\Phi}_{\downarrow}(x_{1}^{+})\end{array}\right.$
(84)

and at the right interface $x_{2}=+L/2$

$\left\\{\begin{array}[]{l}\hat{\Phi}_{\uparrow}(x_{2}^{-})=e^{-ik_{SO}L}\hat{\xi}_{\uparrow}(x_{2}^{+})+\hat{\eta}_{\uparrow}(x_{2}^{+})\\\
\hat{\Phi}_{\downarrow}(x_{2}^{-})=e^{+ik_{SO}L}\hat{\xi}_{\downarrow}(x_{2}^{+})+\hat{\eta}_{\downarrow}(x_{2}^{+})\\\
\partial_{x}\hat{\Phi}_{\uparrow}(x_{2}^{-})=-ik_{SO}\left[e^{-ik_{SO}L}\hat{\xi}_{\uparrow}(x_{2}^{+})-\hat{\eta}_{\uparrow}(x_{2}^{+})\right]\\\
\partial_{x}\hat{\Phi}_{\downarrow}(x_{2}^{-})=+ik_{SO}\left[e^{+ik_{SO}L}\hat{\xi}_{\downarrow}(x_{2}^{+})-\hat{\eta}_{\downarrow}(x_{2}^{+})\right]\end{array}\right.\quad,$
(85)

where, consistently with the low energy limit, we have neglected the
derivatives $\partial_{x}\hat{\xi}$ and $\partial_{x}\hat{\eta}$ of the slowly
varying fields with respect to the term proportional to $k_{SO}$, since they
are characterized by wavevectors $|q|\ll k_{SO}$.

In the central region 2, where only the Zeeman coupling is present, the field
$\hat{\Phi}$ can be expressed as a linear combination of propagating and
evanescent waves that are eigenfunctions of $\sigma_{x}$, so that for
$|x|<L/2$

$\displaystyle\hat{\Phi}(x)$ $\displaystyle=$
$\displaystyle\frac{\hat{h}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
1\end{array}\right)e^{ik_{2,E}x}+\frac{\hat{g}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
1\end{array}\right)e^{-ik_{2,E}x}$ (90) $\displaystyle+$
$\displaystyle\frac{\hat{d}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-1\end{array}\right)e^{\kappa_{2,E}x}+\frac{\hat{f}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-1\end{array}\right)e^{-\kappa_{2,E}x}$ (95)

where $\hat{h}_{E},\hat{g}_{E},\hat{d}_{E}$ and $\hat{f}_{E}$ are mode
operators, whereas $k_{2,E}=k_{Z}\sqrt{1+E/E_{Z}}$,
$\kappa_{2,E}=k_{Z}\sqrt{1-E/E_{Z}}$ and $k_{Z}$ is given in Eq.(70).
Inserting Eq.(90) into Eqs.(84)-(85), one can obtain the link between the
fields in the outer Rashba-dominated regions

$\displaystyle\left(\begin{array}[]{c}\hat{\xi}_{\uparrow}(L/2)\\\
\hat{\xi}_{\downarrow}(L/2)\\\ \hat{\eta}_{\uparrow}(L/2)\\\
\hat{\eta}_{\downarrow}(L/2)\end{array}\right)=\mathsf{M}_{E}\left(\begin{array}[]{c}\hat{\xi}_{\uparrow}(-L/2)\\\
\hat{\xi}_{\downarrow}(-L/2)\\\ \hat{\eta}_{\uparrow}(-L/2)\\\
\hat{\eta}_{\downarrow}(-L/2)\end{array}\right)$ (104)

where the transfer matrix $\mathsf{M}_{E}$ depends on the energy $E$ and on
the size $L$ of the central region through two dimensionless parameters
$k_{Z}L$ and $k_{SO}L$. Details about the derivation of $\mathsf{M}_{E}$ can
be found in the Appendix C. As an illustrative example, here we shall focus on
the midgap value ($E=0$), which in fact well represents the entire low energy
range $|E|\ll E_{Z}$. Moreover, since in the deep Rashba-dominated regime
$k_{Z}L\ll k_{SO}L$, one can keep $k_{SO}L$ finite and consider $k_{Z}L$ as a
small parameter, performing an expansion of $\mathsf{M}_{E=0}$ in its powers.
Neglecting orders $O((k_{Z}L)^{4})$ one obtains

$\mathsf{M}_{0}\simeq\begin{pmatrix}ik_{SO}L/2&A&(1-ik_{SO}L/2)e^{ik_{SO}L}&B\\\
A^{*}&-ik_{SO}L/2&B^{*}&(1+ik_{SO}L/2)e^{-ik_{SO}L}\\\
(1+ik_{SO}L/2)e^{-ik_{SO}L}&-B&-ik_{SO}L/2&A^{*}e^{2iLk_{SO}}\\\
-B^{*}&(1-ik_{SO}L/2)e^{ik_{SO}L}&Ae^{-2ik_{SO}L}&ik_{SO}L/2\end{pmatrix}$
(105)

where

$\displaystyle A$
$\displaystyle=i\frac{-6+k_{SO}L(k_{SO}L+6i)}{12k_{SO}L}e^{2iLk_{SO}}(k_{Z}L)^{2}$
(106) $\displaystyle B$
$\displaystyle=-i\frac{(k_{SO}L)^{2}+6}{12k_{SO}L}e^{ik_{SO}L}(k_{Z}L)^{2}$
(107)

The 8 entries of the transfer matrix (105) containing $A$ and $B$ couple
spin-$\uparrow$ to spin-$\downarrow$ components. Notably, such terms are of
the order $O((k_{Z}L)^{2})$ and in the regime $k_{Z}L\ll 1$ can be neglected
with respect to the other terms, which are $O(1)$ with respect to the variable
$k_{Z}L$. Then, the transfer matrix reduces to

$\mathsf{M}_{0}\simeq\begin{pmatrix}i\frac{k_{SO}L}{2}&0&(1-i\frac{k_{SO}L}{2})e^{ik_{SO}L}&0\\\
0&-i\frac{k_{SO}L}{2}&0&(1+i\frac{k_{SO}L}{2})e^{-ik_{SO}L}\\\
(1+i\frac{k_{SO}L}{2})e^{-ik_{SO}L}&0&-i\frac{k_{SO}L}{2}&0\\\
0&(1-i\frac{k_{SO}L}{2})e^{ik_{SO}L}&0&i\frac{k_{SO}L}{2}\end{pmatrix}\quad.$
(108)

The expression (108) has precisely the form Eq.(32) of the transfer matrix of
the massless+massive Dirac model described in Sec.III, when setting
$\beta_{\uparrow}=-\beta_{\downarrow}=k_{SO}L/2$,
$\chi_{\uparrow}=-\chi_{\downarrow}=k_{SO}L$ and
$\gamma_{\uparrow}=\gamma_{\downarrow}=\nu_{\uparrow}=\nu_{\downarrow}=0$.
Thus, in the regime $k_{Z}L\ll 1$, where the central region is much shorter
than the Zeeman wavelength $l_{Z}=k_{Z}^{-1}$ characterizing the wavefunction
(90) at $E=0$, the transfer matrix is diagonal in spin and becomes independent
of the Zeeman energy $E_{Z}$. Yet, $\mathsf{M}_{0}$ couples massless to
massive modes and still depends on $k_{SO}L$. This parameter, which represents
the ratio of the crossover region $L$ to the spin-orbit length
$l_{SO}=k_{SO}^{-1}$, may be finite because of the deep Rashba-dominated
regime $k_{Z}\ll k_{SO}$.

In turn, the transmission coefficient related to Eq.(108) can be obtained from
the general formula Eq.(64),

$T_{0}=\frac{(k_{SO}L)^{2}}{\left(1+(k_{SO}L/2)^{2}\right)^{2}}\quad,$ (109)

and varies over the entire range $T_{0}\in[0,1]$ as a function of $k_{SO}L$,
as shown in Fig.6. In particular, while for small values $k_{SO}L\ll 1$ the
transmission is low, $T_{0}\sim(k_{SO}L)^{2}$, for finite values of $k_{SO}L$
we observe from Fig.6 that $T_{0}$ increases, and a perfect transmission
$T_{0}=1$ is obtained for $k_{SO}L=2$. Then, for large values of $k_{SO}L$ the
transmission decreases again as $T_{0}\sim 16/(k_{SO}L)^{2}$. The ratio of the
spin-orbit length $l_{SO}=k_{SO}^{-1}$ to the distance $L$ is thus the
parameter controlling the value of $T_{0}$.

Figure 6: The transmission coefficient (109), plotted as a function of
$k_{SO}L$, covers the entire range $T_{0}\in[0,1]$.

###  IV.3 Transmission coefficient in the case of InSb

For definiteness, we consider here an implementation with a InSb NW with
effective electron mass $m^{*}=0.015m_{e}$. Two different portions of the NW
are supposed to be gated by differently biased metals inducing opposite RSOC
values, as previously sketched in Fig.5, and are separated by a crossover
region $L=100{\rm nm}$ where the RSOC is negligible. In Fig.7(a) the solid
curves display the midgap transmission coefficient $T_{0}=T_{E=0}$ as a
function of the spin-orbit energy $E_{SO}$, for different values of the Zeeman
energy $E_{Z}$, obtained from the numerically exact solution of the model (81)
with the profile (83) (see Appendix C for technical details). Moreover the
dashed curve describes the analytical result (109) obtained from the low-
energy limit in the Rashba-dominated regime of the outer regions, i.e. the
massless+massive Dirac model. As one can see, for $E_{SO}\rightarrow 0$, the
exact transmission coefficient tends to 1, regardless of the value of $E_{Z}$,
since all three regions become equal in such a limit. However, for each Zeeman
energy value, when $E_{SO}$ is sufficiently large to enter the deep Rashba-
dominated regime ($2E_{SO}\gg E_{Z}$), all solid curves are well reproduced by
the low-energy limit Eq.(109) (dashed curve), which is independent of $E_{Z}$.
This is thus the regime where the NW gap states are helical and the setup
realizes the Dirac paradox. Despite the absence of spin-flip processes at the
interface, the transmission coefficient is non-vanishing because of the
mediation of the massive modes. In Fig.7(b) the same quantities as in panel
(a) are shown, with a zoom in the range of spin-orbit energy values up to
$E_{SO}=0.5{\rm meV}$, which is the realistic range presently reachable.
Correspondingly, the range of Zeeman energy values $E_{Z}$ ensuring a deep
Rashba dominated regime for the external gated regions is $E_{Z}<0.1{\rm
meV}$. This implies that the linear conductance $G_{0}$, straightforwardly
connected to transmission coefficient through the relation $G_{0}=({\rm
e}^{2}/h)T_{E=0}$, is tunable from low to high values with varying the spin-
orbit energy, which can electrically be done through the gate voltage.

Figure 7:  The Dirac paradox configuration realized with a InSb NW setup where
two outer gated regions are characterized by opposite RSOC and the central
region has a width $L=100\,{\rm nm}$ (see Fig.5). The midgap transmission
coefficient ($E=0$), obtained from the numerically exact solution of the model
(81) with the profile (83), is plotted as a function of the spin orbit energy,
for different values of the external magnetic field $E_{Z}$ (solid curves). At
each value of $E_{Z}$, when the Rashba-dominated regime ($2E_{SO}\gg E_{Z}$)
is reached, the various solid curves all tend to the dashed curve describing
the result Eq.(109), obtained in the low energy limit from the effective
massless+massive Dirac model. Panel (b) is a zoom of panel (a) in the regime
of spin-orbit values that are realistic with present gating techniques.

##  V Conclusions

In this paper we have analyzed the Dirac paradox, illustrated in Fig.1, where
an electron impinging towards an interface can seemingly neither be
transmitted nor reflected. In particular, we have focussed on the interesting
case of one spatial dimension. Indeed, differently from higher dimensional
realizations such as heterojunctions between two 3D Topological Insulators
where electrons can leak along the interface surface, in 1D electrons do not
have a “way out” to escape the paradox. We have first analyzed models that
purely involve massless modes. The first model Eq.(2), where the helicity
change across the interface is accounted for by a spatially inhomogeneous
velocity, leads to conclude that the paradox has no actual solution, namely it
is not possible to build up a scattering state solution. Indeed physical
solutions must necessarily involve electron injection from both sides and are
characterized by a vanishing current. In contrast, the second model Eq.(13),
where the helicity change occurs through a rotation of the spin across the
interface, provides a trivial solution to the paradox, for it directly
introduces spin-flip processes at the interface, which lead to a perfect
transmission.

Then, we have proposed a model, Eq.(19), involving both massless and massive
Dirac modes [see Fig.2] and we have shown that it leads to a non-trivial
solution of the Dirac paradox, even when the transfer matrix forbids spin-flip
processes at the interface. This is possible because of the massive modes
that, despite carrying no current for energies inside their gap, always
exhibit both spin components. Thus, a massless-massive coupling at the
interface indirectly enable an incoming massless electron impinging with
spin-$\uparrow$ to get transmitted as a massless electron with
spin-$\downarrow$ (see Fig.3). Properly defined scattering state solutions
thus exist, and the transmission coefficient depends in general on three
parameters [see Eq.(64)].

Moreover, in Sec.IV, we have discussed how such model can be implemented in
spin-orbit coupled NWs exposed to an external magnetic field, whose midgap
states are characterized by massless modes near the Fermi points $k\sim\pm
2k_{SO}$ and massive modes near $k\sim 0$ [see Fig.4]. The massless modes are
helical in the deep Rashba-dominated regime ($2E_{SO}\gg E_{Z}$) and their
spin orientation is determined by the sign of the RSOC. Because the latter can
be tuned by state-of-the-art gating techniques, a NW with two regions
characterized by opposite RSOC values, as shown in Fig.5, is a suitable
candidate to realize the Dirac paradox configuration in one spatial dimension.
We have shown that the low energy limit of such inhomogeneous NW model is
precisely a particular case of the proposed model (19). The resulting
transmission coefficient Eq.(109) varies over the full range $T_{0}\in[0,1]$
(see Fig.6) as a function of the parameter $k_{SO}L$, where $L$ is the
distance between the two differently gated regions and $k_{SO}$ is the spin-
orbit wavevector that is directly controlled by the RSOC [see Eq.(69)].
Focussing on the specific case of an inhomogeneous InSb NW, we have determined
from model (81) the exact transmission coefficient, which in general depends
both on the spin-orbit and the Zeeman energies [solid curves of Fig.7].
However, whenever the Rashba-dominated regime is reached, the setup realizes
the Dirac paradox configuration. Then, the transmission coefficient is well
captured by the low energy limit result (109) (dashed curve of Fig.7) obtained
from the proposed massless+massive Dirac model and only depends on the spin-
orbit energy $E_{SO}$. Because $E_{SO}$ can be controlled via the gate bias
coupled to the NW, the transmission coefficient and the linear conductance are
electrically tuneable. These results thus represent a conceptual advance in
the understanding of Dirac heterojunctions and pave the way to fruitful
applications of the helical states realized in spin-orbit coupled NWs.

###### Acknowledgements.

Discussions with Björn Trauzettel, Christoph Fleckenstein and Niccolò Traverso
Ziani in the early stage of this work are acknowledged.

##  Appendix A Derivation of the transfer matrix Eq.(32)

In this Appendix we provide details about the derivation of the transfer
matrix (32), i.e. the most general $4\times 4$ matrix fulfilling the
requirements (26) and (27). We first observe that the former requirement (26)
straightforwardly stems from Eq.(21) and the property
$(U_{R}^{\dagger})^{-1}=(U_{R}^{-1})^{\dagger}$, which imply that $U_{R}^{-1}$
fulfills Eq.(21) as well. When taking into account Eq.(20) and the definition
(25), the condition Eq.(26) follows. Second, we observe that the requirement
(27) can equivalently be formulated by requiring that $\mathsf{M}$ must only
involve the combinations $\sigma_{\uparrow}=(\sigma_{0}+\sigma_{z})/2$ and
$\sigma_{\downarrow}=(\sigma_{0}-\sigma_{z})/2$, i.e. $\mathsf{M}$ must have
the form

$\mathsf{M}=\mathsf{M}_{\uparrow}\sigma_{\uparrow}+\mathsf{M}_{\downarrow}\sigma_{\downarrow}\quad,$
(110)

where $\mathsf{M}_{\uparrow,\downarrow}$ are $2\times 2$ matrices acting on
the massless-massless pseudospin space and fulfilling

$\mathsf{M}_{\sigma}^{\dagger}\tau_{z}\mathsf{M}_{\sigma}=-\tau_{z}\hskip
28.45274pt\sigma=\uparrow,\downarrow$ (111)

as a consequence of Eq.(26) and of the properties
$\sigma^{2}_{\uparrow,\downarrow}=\sigma_{\uparrow,\downarrow}$ and
$\sigma_{\uparrow}\sigma_{\downarrow}=[\sigma_{\uparrow},\sigma_{z}]=[\sigma_{\downarrow},\sigma_{z}]=0$.
For each spin sector $\sigma=\uparrow,\downarrow$, the requirement Eq.(111)
imposed on a generic $2\times 2$ complex matrix

$\mathsf{M}_{\sigma}=\begin{pmatrix}a_{\sigma}&b_{\sigma}\\\
c_{\sigma}&d_{\sigma}\end{pmatrix}$ (112)

implies that $|c_{\sigma}|^{2}-|a_{\sigma}|^{2}=1$,
$|b_{\sigma}|^{2}-|d_{\sigma}|^{2}=1$ and
$a^{*}_{\sigma}b_{\sigma}=c^{*}_{\sigma}d_{\sigma}$. These conditions
straightforwardly imply the following expression

$\mathsf{M}_{\sigma}=e^{i\nu_{\sigma}}\left(\begin{array}[]{cc}i\beta_{\sigma}e^{-i\gamma_{\sigma}}&(1-i\beta_{\sigma})\,e^{i\chi_{\sigma}}\\\
(1+i\beta_{\sigma})\,e^{-i\chi_{\sigma}}&-i\beta_{\sigma}e^{i\gamma_{\sigma}}\end{array}\right)\quad,$
(113)

which also fulfills the properties
$\mathsf{M}^{-1}_{\sigma}(\beta_{\sigma},\chi_{\sigma},\nu_{\sigma},\gamma_{\sigma})=\mathsf{M}_{\sigma}(\beta_{\sigma},\chi_{\sigma},-\nu_{\sigma},-\gamma_{\sigma})$
and ${\rm det}(\mathsf{M}_{\sigma})=-\exp[2i\nu_{\sigma}]$. Inserting the two
independent matrices $\mathsf{M}_{\uparrow}$ and $\mathsf{M}_{\downarrow}$
given in Eq.(113) into Eq.(110), the transfer matrix $\mathsf{M}$ in the
$\tau\otimes\sigma$ basis takes the form given in Eq.(32).

Finally, an explicit expression can be given for $U_{L}$ and $U_{R}$ as well.
The requirement (20) can always be fulfilled by choosing for $U_{L}$ the form

$U_{L}=\tau_{0}\sigma_{0}\quad.$ (114)

Then, the expression for $U_{R}=\mathsf{M}^{-1}$ following from Eq.(25) can
straightforwardly be obtained from Eq.(32) by exploiting the property
$\mathsf{M}^{-1}(\bm{\beta},\bm{\chi},\bm{\nu},\bm{\gamma})=\mathsf{M}(\bm{\beta},\bm{\chi},-\bm{\nu},-\bm{\gamma})$,
where each bold symbols denotes the pair of related parameters, e.g.
$\bm{\beta}=(\beta_{\uparrow},\beta_{\downarrow})$.

##  Appendix B details about the NW Hamiltonian and its low energy limit

For the sake of completeness, we provide here some details about the NW
Hamiltonian described in Sec.IV.1. Denoting by $\Omega$ the total NW length
and re-expressing the field in terms of its Fourier modes
$\hat{C}_{k}=(\hat{c}_{k\uparrow},\ \hat{c}_{k\downarrow})^{T}$

$\left(\begin{array}[]{c}\hat{\Phi}_{\uparrow}(x)\\\
\hat{\Phi}_{\downarrow}(x)\end{array}\right)=\frac{1}{\sqrt{\Omega}}\sum_{k}e^{ikx}\left(\begin{array}[]{c}\hat{c}_{k\uparrow}\\\
\hat{c}_{k\downarrow}\end{array}\right)\quad,$ (115)

the NW Hamiltonian $\hat{\mathcal{H}}_{NW}$ is compactly rewritten in terms of
a $2\times 2$ matrix $H_{NW}(k)$, i.e.
$\hat{\mathcal{H}}_{NW}=\sum_{k}\hat{C}_{k}^{\dagger}H_{NW}(k)\hat{C}_{k}$. In
turn, this also highlights the energy scales involved in the problem. In
particular, the first two terms acquire the form

$\displaystyle\hat{\mathcal{H}}_{kin}+\hat{\mathcal{H}}_{R}=$ $\displaystyle=$
$\displaystyle\sum_{k}\hat{C}_{k}^{\dagger}\left(\frac{\hbar^{2}}{2m^{*}}(k\sigma_{0}-s_{\alpha}k_{SO}\sigma_{z})^{2}-E_{SO}\sigma_{0}\right)\hat{C}_{k}$

and describe two parabolic spin bands that are lowered by the spin-orbit
energy (67) and horizontally shifted by the spin-orbit wavevector (69) with
the sign (71) of the RSOC determining whether the shift is positive or
negative in $k$-axis. Assuming $h_{x}>0$ for definiteness, the Zeeman term is
rewritten as

$\displaystyle
H_{Z}=-E_{Z}\sum_{k}\hat{C}_{k}^{\dagger}\sigma_{x}\hat{C}_{k}\quad,$ (117)

where $E_{Z}$ is the Zeeman energy given in Eq.(68). Summing up
Eqs.(LABEL:eq:H_rk) and (117) the diagonalization of the resulting $H_{NW}(k)$
is straightforward. Denoting $\varepsilon^{0}_{k}=\hbar^{2}k^{2}/2m^{*}$, the
spectrum consists of two energy bands

$E_{\pm}(k)=\,\varepsilon^{0}_{k}\pm\sqrt{E_{Z}^{2}+\alpha^{2}k^{2}}\quad,$
(118)

separated at $k=0$ by a gap $2E_{Z}$ centered around the midgap energy $E=0$.
The eigenfunctions related to the spectrum (118) are
$\psi_{k,\pm}(x)=w_{k,\pm}\exp[ikx]/\sqrt{\Omega}$. They describe plane waves
with spinors

$\displaystyle w_{k,-}=\left(\begin{array}[]{c}\cos\frac{\theta_{k}}{2}\\\ \\\
\sin\frac{\theta_{k}}{2}\end{array}\right)\hskip
14.22636ptw_{k,+}=\left(\begin{array}[]{c}-\sin\frac{\theta_{k}}{2}\\\ \\\
\cos\frac{\theta_{k}}{2}\end{array}\right)\,,\hskip 14.22636pt$ (125)

whose spin orientation
$\mathbf{n}(k)\equiv\left(\sin\theta_{k}\,,0\,,\,\cos\theta_{k}\right)$ lies
on the $xz$-plane and depends on the wavevector $k$, forming with the $z$-axis
an angle $\theta_{k}\in[0,\pi]$ defined through

$\left\\{\begin{array}[]{lcl}\cos\theta_{k}&=&\displaystyle\frac{\alpha
k}{\sqrt{E_{Z}^{2}+\alpha^{2}k^{2}}}\\\
\sin\theta_{k}&=&\displaystyle\frac{{E_{Z}}}{\sqrt{E_{Z}^{2}+\alpha^{2}k^{2}}}\end{array}\right.\quad.$
(126)

Furthermore, for energies $|E|<E_{Z}$, the model also exhibits evanescent wave
solutions $\tilde{\psi}_{\kappa,\pm}(x)=\tilde{w}_{\kappa,\pm}\exp[\kappa
x]/\sqrt{\Omega}$. They describe plane waves with spinors

$\displaystyle\tilde{w}_{\kappa,\pm}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\mp\exp\left[\pm
i\mbox{arctan}\left(\alpha\kappa/\sqrt{E_{Z}^{2}-(\alpha\kappa)^{2}}\right)\right]\\\
\\\ 1\end{array}\right)\hskip 14.22636pt$ (130)

with energy
$E_{\pm}=-\varepsilon^{0}_{\kappa}\pm\sqrt{E_{Z}^{2}-(\alpha\kappa)^{2}}$.
While these solutions are not normalizable in a homogeneous NW, they must be
taken into account in the inhomogeneous RSOC problem.

Let us now focus on the regime ($|E|\ll E_{Z}\ll 2E_{SO}$) and derive an
effective low energy NW Hamiltonian.

Expansion near $k=\pm 2k_{SO}$. In the deep Rashba-dominated regime ($E_{Z}\ll
2E_{SO}$), one finds that, up to $O\big{(}(E_{Z}/2E_{SO})^{2}\big{)}$,

$\displaystyle E_{-}(k)\approx 0\,\Leftrightarrow\,k\approx\pm 2k_{SO}$ (131)
$\displaystyle\cos\theta_{k=\pm 2k_{SO}}\approx\pm s_{\alpha}\quad,$ (132)

so that the spinors (125) of the lower band propagating modes near $k\sim\pm
2k_{SO}$ reduce to eigenstates of $\sigma_{z}$, $(1,0)^{T}$ or $(0,1)^{T}$,
depending on the sign $s_{\alpha}$ of the RSOC [see Eq.(71)]. To extract the
low energy Hamiltonian governing their dynamics, let us consider, for
instance, $\alpha>0$ like in Fig.4, and focus e.g. on the vicinity of the
right Fermi point $+2k_{SO}$. Setting $k=2k_{SO}+q$ and performing an
expansion of Eqs.(LABEL:eq:H_rk) and (117) for $|q|\ll k_{SO}$, one obtains

$\displaystyle\left.\hat{\mathcal{H}}_{NW}\right|_{k\simeq+2k_{SO}}\simeq$
(135) $\displaystyle\simeq$ $\displaystyle\sum_{|q|\ll
k_{SO}}\hat{C}^{\dagger}_{2k_{SO}+q}\left(\begin{array}[]{cc}\hbar
v_{SO}q&-E_{Z}\\\ -E_{Z}&8E_{SO}\end{array}\right)\hat{C}_{2k_{SO}+q}$
$\displaystyle\simeq$ $\displaystyle\sum_{|q|\ll k_{SO}}\hbar
v_{SO}q\,\hat{c}^{\dagger}_{2k_{SO}+q,\uparrow}\hat{c}_{2k_{SO}+q,\uparrow}\quad,$
(136)

where $v_{SO}=\hbar k_{SO}/m^{*}$. The last line of Eq.(136) follows from the
fact that, while the spin-$\uparrow$ band is characterized by a low-energy
$\hbar v_{F}q$, the spin-$\downarrow$ band has a large energy $8E_{SO}$ much
above the magnetic gap. The weak Zeeman energy $E_{Z}\ll 2E_{SO}$ cannot
couple them, so that in the low energy sector $|E|\ll E_{Z}$ only the
spin-$\uparrow$ states matter. One can proceed in a similar manner near the
$-2k_{SO}$ Fermi point, obtaining that only the spin-$\downarrow$ states
matter, proving that the states are helical. Repeating the same calculation
for $\alpha<0$ one obtains the opposite helicity. From Eq.(136) the massless
propagating low energy excitations ($|q|\ll k_{SO}$) are thus described by the
set of operators

$\left\\{\begin{array}[]{lcl}\hat{\xi}_{q\uparrow}&\doteq&\hat{c}_{2s_{\alpha}k_{SO}+q,\uparrow}\\\
&&\\\
\hat{\xi}_{q\downarrow}&\doteq&\hat{c}_{-2s_{\alpha}k_{SO}+q,\downarrow}\end{array}\right.$
(137)

where $s_{\alpha}$ is given by Eq.(71).

Expansion near $k=0$. In the low energy range $|E|\ll E_{Z}$ there are also
gapped (i.e. massive) modes, related to the upper and lower bands for $k\sim
0$ (see Fig.4). Performing an expansion of Eqs.(LABEL:eq:H_rk) and (117) in
$q=k$ (with $|q|\ll k_{SO}$) and introducing the new set of operators

$\left\\{\begin{array}[]{lcl}\hat{\eta}_{q\uparrow}&=&\hat{c}_{q\uparrow}\\\
&&\\\
\hat{\eta}_{q\downarrow}&=&\hat{c}_{q\downarrow}\end{array}\right.\quad,$
(138)

one obtains the low energy expression

$\displaystyle\left.\hat{\mathcal{H}}_{NW}\right|_{k\simeq 0}$
$\displaystyle\simeq-\sum_{|q|\ll k_{SO}}\hbar
s_{\alpha}v_{SO}\,q\,(\hat{\eta}_{q\uparrow}^{\dagger}\
\hat{\eta}_{q\downarrow}^{\dagger})\sigma_{z}\left(\begin{array}[]{c}\hat{\eta}_{q\uparrow}\\\
\hat{\eta}_{q\downarrow}\end{array}\right)$ (141) $\displaystyle-$
$\displaystyle E_{Z}\\!\\!\sum_{|q|\ll
k_{SO}}\\!\\!(\hat{\eta}_{q\uparrow}^{\dagger}\
\hat{\eta}_{q\downarrow}^{\dagger})\,\sigma_{x}\left(\begin{array}[]{c}\hat{\eta}_{q\uparrow}\\\
\hat{\eta}_{q\downarrow}\end{array}\right)\quad.$ (144)

Summing up Eqs.(136) and (144) one obtains a low-energy NW Hamiltonian.
Moreover, one can observe that such model shares the same low-energy physics
as the Dirac model given in Eq.(77), obtained by removing the constraints on
wave vector $q$, which can therefore be regarded to as the effective low
energy model for the NW.

##  Appendix C The scattering problem for the inhomogeneous NW with the
profile (83)

The solution of the scattering problem for the model (81) with the piecewise
constant profile (83) can be obtained from the expression of the electron
field operator in the three regions. For an energy $E$ within the magnetic gap
($|E|<E_{Z}$) one has

$\displaystyle\hat{\Phi}_{E}(x)=\left\\{\begin{array}[]{lcl}\displaystyle\hat{a}_{LE}w_{k_{E},-}e^{ik_{E}x}+\hat{b}_{LE}w_{-k_{E},-}e^{-ik_{E}x}+\hat{c}_{LE}\tilde{w}_{\kappa_{E},s_{E}}e^{\kappa_{E}x}&&x<-L/2\\\
&&\\\ \displaystyle\frac{\hat{h}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
1\end{array}\right)e^{ik_{2,E}x}+\frac{\hat{g}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
1\end{array}\right)e^{-ik_{2,E}x}+\frac{\hat{d}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-1\end{array}\right)e^{\kappa_{2,E}x}+\frac{\hat{f}_{E}}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-1\end{array}\right)e^{-\kappa_{2,E}x}&&|x|<L/2\\\ &&\\\
\displaystyle\hat{a}_{RE}w_{-k_{E},-}e^{-ik_{E}x}+\hat{b}_{RE}w_{k_{E},-}e^{ik_{E}x}+\hat{c}_{RE}\tilde{w}_{-\kappa_{E},s_{E}}e^{-\kappa_{E}x}&&x>+L/2\\\
\end{array}\right.$ (158)

where

$\displaystyle k_{E}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2m^{*}}}{\hbar}\sqrt{E+2E_{SO}+\sqrt{4EE_{SO}+4E_{SO}^{2}+E_{Z}^{2}}}\hskip
28.45274pt$ (159) $\displaystyle\kappa_{E}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2m^{*}}}{\hbar}\sqrt{-E-2E_{SO}+\sqrt{4EE_{SO}+4E_{SO}^{2}+E_{Z}^{2}}}\hskip
28.45274pt$ (160) $\displaystyle s_{E}$ $\displaystyle=$
$\displaystyle\mbox{sgn}(E+E_{Z}^{2}/4E_{SO})$ (161) $\displaystyle k_{2,E}$
$\displaystyle=$ $\displaystyle k_{Z}\sqrt{1+E/E_{Z}}\hskip
28.45274pt\kappa_{2,E}=k_{Z}\sqrt{1-E/E_{Z}}\quad,\hskip 28.45274pt$ (162)

while the spinors $w_{\pm k_{E},-}$ and $\tilde{w}_{\pm\kappa_{E},s_{E}}$ are
given in Eqs.(125) and (130), respectively.

Imposing the boundary conditions (82) to the field (158), one expresses the
outgoing operators $\hat{b}_{L/R\,E}$ in terms of the operators $\hat{a}_{LE}$
and $\hat{a}_{RE}$ describing the modes incoming from the left and from the
right region, respectively. The transmission amplitudes $t_{E}$ and
$t^{\prime}_{E}$ are then obtained through the relations
$\hat{b}_{RE}=t_{E}\hat{a}_{LE}$ and
$\hat{b}_{LE}=t^{\prime}_{E}\hat{a}_{RE}$. The resulting transmission
coefficient $T_{E}=|t_{E}|^{2}=|t^{\prime}_{E}|^{2}$ is numerically exact and
is plotted in the solid curves of Fig.7 as a function of the spin orbit
energy, at the midgap energy $E=0$ and for different values of the external
magnetic field $E_{Z}$.

However, as observed in Sec.IV, an analytical expression for the transmission
coefficient can be obtained in the low energy limit (dashed curve in Fig.7),
where the inhomogeneous NW physics is well captured by the effective
massless+massive Dirac theory. Such an expression directly follows from the
transfer matrix (104) connecting the massless and massive fields of the outer
Rashba-dominated regions, which can be obtained as follows. Inserting Eq.(90)
into the low energy boundary conditions Eqs.(84)-(85), the latter can be re-
expressed in a matrix form as

$\displaystyle\mathsf{P}\left(\begin{array}[]{c}\hat{\xi}_{\uparrow}(-L/2)\\\
\hat{\xi}_{\downarrow}(-L/2)\\\ \hat{\eta}_{\uparrow}(-L/2)\\\
\eta_{\downarrow}(-L/2)\end{array}\right)$
$\displaystyle=\mathsf{V}(-L/2)\left(\begin{array}[]{c}\hat{h}\\\ \hat{g}\\\
\hat{d}\\\ \hat{f}\end{array}\right)$ (171)
$\displaystyle\mathsf{V}(L/2)\left(\begin{array}[]{c}\hat{h}\\\ \hat{g}\\\
\hat{d}\\\ \hat{f}\end{array}\right)$
$\displaystyle=\mathsf{Q}\left(\begin{array}[]{c}\hat{\xi}_{\uparrow}(L/2)\\\
\hat{\xi}_{\downarrow}(L/2)\\\ \hat{\eta}_{\uparrow}(L/2)\\\
\hat{\eta}_{\downarrow}(L/2)\end{array}\right)\quad,$ (180)

where the energy dependence of the operators has been dropped to make the
notation lighter. Here

$\mathsf{P}=\begin{pmatrix}e^{-ik_{SO}L}&0&1&0\\\ 0&e^{ik_{SO}L}&0&1\\\
k_{SO}e^{-ik_{SO}L}&0&-k_{SO}&0\\\
0&-k_{SO}e^{ik_{SO}L}&0&k_{SO}\end{pmatrix}\quad,$ (181)
$\mathsf{Q}=\begin{pmatrix}e^{-ik_{SO}L}&0&1&0\\\ 0&e^{ik_{SO}L}&0&1\\\
-k_{SO}e^{-ik_{SO}L}&0&k_{SO}&0\\\
0&k_{SO}e^{ik_{SO}L}&0&-k_{SO}\end{pmatrix}$ (182)

and

$\displaystyle\mathsf{V}(x)=\frac{1}{\sqrt{2}}\times$
$\displaystyle\begin{pmatrix}e^{ik_{2,E}x}&e^{-ik_{2,E}x}&e^{\kappa_{2,E}x}&e^{-\kappa_{2,E}x}\\\
e^{ik_{2,E}x}&e^{-ik_{2,E}x}&-e^{\kappa_{2,E}x}&-e^{-\kappa_{2,E}x}\\\
ik_{2,E}e^{ik_{2,E}x}&-ik_{2,E}e^{-ik_{2,E}x}&\kappa_{2,E}e^{\kappa_{2,E}x}&-\kappa_{2,E}e^{-\kappa_{2,E}x}\\\
ik_{2,E}e^{ik_{2,E}x}&-ik_{2,E}e^{-ik_{2,E}x}&-\kappa_{2,E}e^{\kappa_{2,E}x}&\kappa_{2,E}e^{-\kappa_{2,E}x}\\\
\end{pmatrix}$

The transfer matrix $\mathsf{M}_{E}$ appearing in Eq.(104) can thus
straightforwardly be obtained as
$\mathsf{M}_{E}=\mathsf{Q}^{-1}\mathsf{V}(L/2)\mathsf{V}^{-1}(-L/2)\mathsf{P}$.
In particular, setting the energy to the midgap value $E=0$ and expanding in
the parameter $k_{Z}L$ one obtains Eq.(105), up to $O((k_{Z}L)^{4})$ terms.

## References

## References

  * (1) A. K. Geim, and K.S. Novoselov, Nat. Mater. 6, 183 (2007).
  * (2) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * (3) S. Das Sarma, S. Adam, E. H. Hwang, E. Rossi, Rev. Mod. Phys. 83, 407 (2011).
  * (4) M. Z. Hasan, and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (5) X.-L. Qi, and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (6) S. Jia, S.-Y. Xu, and M. Z. Hasan, Nat. Mater. 15, 1140 (2016).
  * (7) S.A. Yang, SPIN 6, 1640003 (2016).
  * (8) R. Takahashi, and S. Murakami, Phys. Rev. Lett. 107, 166805 (2011)
  * (9) Ch. De Beule, and B. Partoens, Phys. Rev. B 87, 115113 (2013)
  * (10) T. Habe, and Y. Asano, Phys. Rev. B 88, 155442 (2013)
  * (11) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * (12) C. Liu, T. L. Hughes, X.-L. Qi, K. Wang, and S.-C. Zhang, Phys. Rev. Lett. 100, 236601 (2008).
  * (13) M. König, H. Buhmann, L. W. Molenkamp, T. Hughes, C.-X. Liu, X.-L. Qi, and S.-C. Zhang, J. Phys. Soc. Jpn 77, 031007 (2008).
  * (14) I. Knez, R.R. Du , and G. Sullivan, Phys. Rev. Lett. 107, 136603 (2011).
  * (15) L. Du, I. Knez, G. Sullivan, and R.R. Du, Phys. Rev. Lett. 114, 096802 (2015).
  * (16) C. H. L. Quay, T. L. Hughes, J. A. Sulpizio, K. W. Pfeiffer, K. W. Baldwin, K. W. West, D. Goldhaber-Gordon, and R. De Picciotto, Nat. Phys. 6, 336 (2010).
  * (17) I. van Weperen, S. R. Plissard, E.P.A.M. Bakkers, S.M. Frolov, and Leo P. Kouwenhoven, Nano Lett. 13, 387 (2013).
  * (18) J. Kammhuber, M.C. Cassidy, F. Pei, M.P. Nowak, A. Vuik, Ö. Gül, D. Car, S.R. Plissard, E.P.A.M. Bakkers, M. Wimmer, and L.P. Kouwenhoven, Nat. Comm. 8, 478 (2017).
  * (19) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, L. P. Kouwenhoven, Science 336, 1003 (2012).
  * (20) L.P. Rokhinson, X. Liu, and J. K. Furdyna, Nat. Phys. 8, 795 (2012).
  * (21) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nat. Phys. 8, 887 (2012).
  * (22) M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Nano Lett. 12, 6414 (2012).
  * (23) E.J.H. Lee, X. Jiang, M. Houzet, R. Aguado, C. M. Lieber, and S. De Franceschi, Nat. Nanotech. 9, 79 (2014).
  * (24) S. M. Albrecht, A. P. Higginbotham, M. Madsen, F. Kuemmeth, T. S. Jespersen, J. Nygård, P. Krogstrup, and C. M. Marcus, Nature 531, 206 (2016).
  * (25) M.T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus, Science 354, 1557 (2016).
  * (26) Ö. Gül , H. Zhang, J.D.S. Bommer, M.W. A. de Moor, D. Car, S.R. Plissard, E.P.A.M. Bakkers, A. Geresdi, K. Watanabe, T. Taniguchi, and L.P. Kouwenhoven, Nat. Nanotech. 13, 192 (2018).
  * (27) P. Strěda, and P. Šeba, Phys. Rev. Lett. 90, 256601 (2003).
  * (28) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * (29) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * (30) C. Kloeffel, M. Trif, and D. Loss, Phys. Rev. B 84, 195314 (2011).
  * (31) M. Cheng, and R. M. Lutchyn, Phys. Rev. B 86, 134522 (2012).
  * (32) P. Szumniak, D. Chevallier, D. Loss, and J. Klinovaja, Phys. Rev. B 96, 041401(R) (2017).
  * (33) A. Y. Kitaev, Phys. Usp. 44, 131 (2001).
  * (34) J. Alicea, Rep. Progr. Phys. 75, 076501 (2012).
  * (35) M. Sato, and S. Fujimoto J. Phys. Soc. Jpn 85, 072001 (2016).
  * (36) R. Aguado, La Rivista del Nuovo Cimento, 40, 523 (2017).
  * (37) S. Sasaki, K. Tateno, G. Zhang, H. Suominen, Y. Harada, S. Saito, A. Fujiwara, T. Sogawa, and K. Muraki, Appl. Phys. Lett. 103, 213502 (2013).
  * (38) A. M. Burke, D. J. Carrad, J. G. Gluschke, K. Storm, S. Fahlvik Svensson, H. Linke, L. Samuelson, and A. P. Micolich, Nano Lett. 15, 2836 (2015).
  * (39) K. Takase, Y. Ashikawa, G. Zhang, K. Tateno, and S. Sasaki, Sci. Rep. 7, 930 (2017).
  * (40) S.R. Das, in Nanoelectronics: Devices, Circuits and Systems, edited by B.K. Kaushik (Elsevier, Amsterdam, 2019), Chap.11, p. 355.
  * (41) Y. Yin, Z. Zhang, H. Zhong, C. Shao, X. Wan, C. Zhang, J. Robertson, and Y. Guo, ACS Appl. Mat. and Interfaces 13, 3387 (2021).
  * (42) K. Takase, K. Tateno, and S. Sasaki, Appl. Phys. Lett. 119, 013102 (2021).
  * (43) D. Liang and X. P.A. Gao, Nano Lett. 12, 3263 (2012).
  * (44) B. Slomski, G. Landolt, S. Muff, F. Meier, J. Osterwalder, and J.H. Dil, New J. Phys. 15, 125031 (2013).
  * (45) I. van Weperen, B. Tarasinski, D. Eeltink, V. S. Pribiag, S. R. Plissard, E.P.A.M. Bakkers, L.P. Kouwenhoven, and M. Wimmer, Phys. Rev. B 91, 201413(R) (2015).
  * (46) Z. Scherübl, G. Fülöp, M. H. Madsen, J. Nygård, and S. Csonka, Phys. Rev. B 94, 035444 (2016).
  * (47) J. R. Bindel, M. Pezzotta, J. Ulrich, M. Liebmann, E. Y. Sherman, and M. Morgenstern, Nat. Phys. 12, 920 (2016).
  * (48) J. Borge and I. V. Tokatly, Phys. Rev. B 96, 115445 (2017).
  * (49) Ch. Kloeffel, M. J. Rančić, and D. Loss, Phys. Rev. B 97, 235422 (2018).
  * (50) P. Wójcik, A. Bertoni, and G. Goldoni, Phys. Rev. B 97, 165401 (2018).
  * (51) H.Tsai, S. Karube, K. Kondou, N. Yamaguchi, Y. Otani, Sci. Rep. 8, 5564 (2018).
  * (52) K. Premasiri and X. P. A. Gao, J. Phys. Cond. Matt. 31, 193001 (2019).
  * (53) D. Shcherbakov, P. Stepanov, S. Memaran, Y. Wang, Y. Xin, J. Yang, K. Wei, R. Baumbach, W. Zheng, K. Watanabe, T. Taniguchi, M. Bockrath, D. Smirnov, T. Siegrist2, W. Windl, L. Balicas, C.N. Lau, Sci. Adv. 7, eabe2892 (2021).
  * (54) O. Krupin, G. Bihlmayer, K. Starke, S. Gorovikov, J. E. Prieto, K. Döbrich, S. Blügel, and G. Kaindl, Phys. Rev. B 71, 201403(R) (2005).
  * (55) B. Slomski, G. Landolt, S. Muff, F. Meier, J. Osterwalder, and J.H. Dil, New J. Phys. 15, 125031 (2013).
  * (56) W. Wang, X.M. Li, J.Y. Fu, J. Magn. Magn. Mat. 411, 84 (2016).
  * (57) F. Nagasawa, A.A. Reynoso, J.P. Baltanás, D. Frustaglia, H. Saarikoski, and J. Nitta, Phys. Rev. B 98, 245301 (2018).
  * (58) N.M.R. Peres, J. Phys.: Condens. Matter 21, 095501 (2009).
  * (59) H. A. Nilsson, Ph. Caroff, C. Thelander, M. Larsson, J. B. Wagner, L.-E. Wernersson, L. Samuelson, and H. Q. Xu, Nano Lett. 9, 3151 (2009).
  * (60) M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Nano Lett. 12, 6414 (2012).
  * (61) S. Nadj-Perge, V. S. Pribiag, J.W.G. van den Berg, K. Zuo, S. R. Plissard, E. P. A. M. Bakkers, S. M. Frolov, and L. P. Kouwenhoven, Phys. Rev. Lett. 108, 166801 (2012).
  * (62) P. Roulleau, T. Choi, S. Riedi, T. Heinzel, I. Shorubalko, T. Ihn, and K. Ensslin, Phys. Rev. B 81, 155449 (2010).
  * (63) H. J. Joyce, C. J. Docherty, Q. Gao, H. H. Tan, C. Jagadish, J. Lloyd-Hughes, L. M. Herz and M. B. Johnston, Nanotechnology 24, 214006 (2013).
  * (64) D. Sánchez, L. Serra, Phys. Rev. B 74 153313 (2006).
  * (65) D. Sánchez, L. Serra, and M.-S. Choi, Phys. Rev. B 77, 035315 (2008).
  * (66) M. M. Glazov, and E. Ya. Sherman, Phys. Rev. Lett. 107, 156602 (2011).
  * (67) A. F. Sadreev, and E. Y. Sherman, Phys. Rev. B 88, 115302 (2013).
  * (68) M. Modugno, E. Ya. Sherman, and V. V. Konotop, Phys. Rev. A 95, 063620 (2017).
  * (69) J. Klinovaja, and D. Loss, Eur. Phys. J. B 88, 62 (2015).
  * (70) F. Dolcini and F. Rossi, Phys. Rev. B 98, 045436 (2018).
  * (71) L. Rossi, F. Dolcini, and F. Rossi, Phys. Rev. B 101, 195421 (2020).
  * (72) L. Rossi, F. Dolcini, and F. Rossi, Eur. Phys. J. Plus 135, 597 (2020).
  * (73) S. Datta, Electronic Transport in Mesoscopic Systems, (Cambridge University Press, Cambridge, 1995)