%\documentclass[prl,twocolumn,amsmath,amssymb,floatfix,showpacs, galley]{revtex4-1}
%\documentclass[twocolumn, preprint]{revtex4}
% Include figure files
% bold mat
\documentclass[prl,twocolumn,floatfix,tighten,bibtex,letterpaper]{revtex4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{graphicx}
\usepackage{bm}
\usepackage{amssymb}
\usepackage{color}
\usepackage{epsfig}
\usepackage{subfigure}
%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{Version=5.50.0.2890}
%TCIDATA{}
%TCIDATA{BibliographyScheme=BibTeX}
%TCIDATA{LastRevised=Monday, July 11, 2011 00:43:05}
%TCIDATA{}
\marginparwidth=5cm
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\addtolength{\topmargin}{1.0cm}
\def\der#1#2{\frac{\partial #1}{\partial #2}}
\def \tr{{\mbox{tr~}}}
\def \ket#1{|#1\rangle}
\def \ra{{\rightarrow}}
\def \ua{{\uparrow}}
\def \da{{\downarrow}}
\def \be{\begin{equation}}
\def \ee{\end{equation}}
\def \ba{\begin{array}}
\def \ea{\end{array}}
\def \bea{\begin{eqnarray}}
\def \eea{\end{eqnarray}}
\def \nn{\nonumber}
\def \l{\left}
\def \r{\right}
\def \rr{\right}
\def \half{{1\over 2}}
\def \etal{{\it {et al}}}
\def \H{{\cal{H}}}
\def \cM{{\cal{M}}}
\def \cN{{\cal{N}}}
\def \cQ{{\cal Q}}
\def \cI{{\cal I}}
\def \cV{{\cal V}}
\def \cG{{\cal G}}
\def \bS{{\bf S}}
\def \bL{{\bf L}}
\def \bG{{\bf G}}
\def \bQ{{\bf Q}}
\def \bR{{\bf R}}
\def \br{{\bf r}}
\def \bu{{\bf u}}
\def \bq{{\bf q}}
\def \bk{{\bf k}}
\def \tJ{{\tilde{J}}}
\def \W{{\Omega}}
\def \e{{\epsilon}}
\def \lam{{\lambda}}
\def \a{{\alpha}}
\def \t{{\theta}}
\def \b{{\beta}}
\def \g{{\gamma}}
\def \D{{\Delta}}
\def \d{{\delta}}
\def \w{{\omega}}
\def \s{s}
\def \p{p}
\def \q{q}
\def \f{{\varphi}}
\def \x{{\chi}}
\def \h{{\eta}}
\def \G{{\Gamma}}
\def \z{{\zeta}}
\def \hatt{{\hat{\t}}}
\def \hn{{\bar{n}}}
\def \vk{{\bf{k}}}
\def \vq{{\bf{q}}}
\def \on{\overline{n}}
\def \intt{\int\limits}
\def \gk{{\g_{\vk}}}
\def \nd{{^{\vphantom{\dagger}}}}
\def \yd{^\dagger}
\def \av#1{{\langle#1\rangle}}
\def \csbar{\overline{\cos}}
\def\O {\Omega}
\def \summ{\sum\limits}
\def\intt{\int\limits}
\def\nbar{\overline{n}}
\def\sx{\sigma^x}
\def\sy{\sigma^y}
\def\sz{\sigma^z}
\def\tx{\tau^x}
\def\ty{\tau^y}
\def\tz{\tau^z}
\def\bra#1{\left\langle #1\right|}
\def\ket#1{\left| #1\right\rangle}
\begin{document}
\title{Unconventional Josephson signatures of Majorana bound states\\
Supplementary Material}
\maketitle
\section{Perturbative Calculation}
Let us pursue here a detailed calculation of the Josephson coupling across
the Majorana junction described in Fig. 1. We will first find
the wave functions of the Majorana states localized on each domain wall,
ignoring the existence of the other interface. We will denote these states
as $\left| L\right\rangle$ and $\left| R\right\rangle$. Next, we follow the
usual procedure for finding tight-binding states and Hamiltonians. We first
calculate the overlap matrix, $M_{\alpha\beta}=\left\langle
\alpha\right|\beta\rangle$ with $\alpha,\,\beta=L,\,R$, and the Hamiltonian
matrix within this subspace, $h_{\alpha\beta}=\left\langle \alpha\right|{%
\mathcal{H}}\left| \beta\right\rangle$. It is easy to see that the
approximate hybridization Hamiltonian is then given by
\begin{equation}
H_{maj}=M^{-1/2} h M^{-1/2}. \label{tbm}
\end{equation}
\textbf{Single Majorana solution at $x=L$.} We first solve for the
zero-energy eigenstates of the Hamiltonian (6) with parameters:
\begin{equation}
\begin{array}{c}
\Delta (x)=\Delta _{r}\Theta (x-L)+\Delta _{m}(1-\Theta (x-L)) \\
\mu (x)=\mu _{r}\Theta (x-L)+\mu _{m}(1-\Theta (x-L))%
\end{array}%
\end{equation}%
The solution has the same form on the two sides of the domain wall, but with
different parameters. We denote the side of the domain with the index $s$
being $s=r,\,m$ for right or middle. By squaring the Hamiltonian and looking
for momentum values yielding zero energy states, we find two imaginary
momenta on each side, which correspond to the spatial decay constants $%
\lambda _{s\pm }^{-1}$ given by Eq. (7). The wave function
associated with each side of the domain is:
\begin{equation}
\left\vert r\right\rangle =\Psi _{s}(x)=\mathit{R}_{s+}\psi _{s+}e^{-\frac{%
|x-L|}{\lambda _{s+}}}+\mathit{R}_{s-}\psi _{s-}e^{-\frac{|x-L|}{\lambda
_{s-}}}. \label{Psix}
\end{equation}%
with $\mathit{R}_{s\pm }$ being four complex numbers determining the
amplitude of the wave functions corresponding to the two decay lengths, and
with $\psi _{s\pm }$ being four, four-dimensional vectors, which when $\phi
_{s}=0$ are given by:
\begin{equation}
\psi _{m\pm }=\left(
\begin{array}{c}
1 \\
e^{-i\zeta _{m}} \\
\pm i \\
\mp ie^{-i\zeta _{m}}%
\end{array}%
\right) ,\,\,\psi _{r\pm }=\left(
\begin{array}{c}
1 \\
e^{\pm i\zeta _{r}} \\
i \\
-ie^{\pm i\zeta _{r}}%
\end{array}%
\right)
\end{equation}%
with $\exp (i\zeta _{s})=\frac{\mu _{s}+i\sqrt{B^{2}-\mu _{s}^{2}}}{B}$, for
$s=\ell ,\,r,\,m$. These solutions are the building blocks for each Majorana
state. In order to obtain what the wave function becomes when $\phi _{s}$
(the phases of the superconducting electrodes) deviate from zero, we can
apply the rotations: $\hat{U}_{\phi }=\exp \left( i\frac{\phi }{2}\tau
_{z}\right) $ such that:
\begin{equation}
\psi _{s\pm }^{(\phi _{s})}=\hat{U}_{\phi _{s}}\psi _{s\pm }.
\end{equation}%
Obtaining the Majorana solution follows from matching the boundary condition
of the solutions, and from them finding the coefficients $\mathit{R}_{s\pm }$%
.
To avoid the complicated expression that could arise in the most general
case of Majorana coupling, we concentrate on the case where $%
\Delta_r=\Delta_{\ell}>\Delta_m$ and $\mu_r=\mu_{\ell}=\mu$ and $\mu_m=0$.
This choice does not constitute a substantial loss of generality, and is
useful for grasping the results of our calculations. A straightforward but
rather tedious calculation leads to the following solution for the
amplitudes of the decaying waves of the right Majorana state under the above
assumptions:
\begin{equation}
\begin{array}{c}
\left(%
\begin{array}{c}
\mathit{R}_{m+} \\
\mathit{R}_{m-}%
\end{array}
\right)=\frac{2\sin\zeta_r}{1+ie^{-i\zeta_r}}\left(%
\begin{array}{c}
i\sin\left(\frac{\phi_r-\phi_m}{2}\right) \\
\cos\left(\frac{\phi_r-\phi_m}{2}\right)%
\end{array}
\right), \\
\left(
\begin{array}{c}
\mathit{R}_{r+} \\
\mathit{R}_{r-}%
\end{array}
\right)=\left(
\begin{array}{c}
-\frac{1+ie^{i\zeta_r}}{1+ie^{-i\zeta_r}} \\
1%
\end{array}
\right).%
\end{array}%
\end{equation}
By symmetry, we can infer the structure of the left Majorana, which is localized about $x=0$:
\begin{equation}
\left\vert \mathit{L}\right\rangle =\Psi _{s}^{(L)}(x)=\mathit{L}_{s+}\psi
_{s+}e^{{-\frac{|x|}{\lambda _{s+}}}}+\mathit{L}_{s-}\psi _{s-}e^{{-\frac{%
|x|}{\lambda _{s-}}}}. \label{PsixL}
\end{equation}%
The amplitudes $\mathit{L}_{s\pm }$ also depend on the phases on the left
and middle segment of the wire in a similar way:
\begin{equation}
\begin{array}{c}
\left(
\begin{array}{c}
\mathit{L}_{m+} \\
\mathit{L}_{m-}%
\end{array}%
\right) =\frac{2\sin \zeta _{\ell }}{1-ie^{i\zeta _{\ell }}}\left(
\begin{array}{c}
-i\sin \left( \frac{\phi _{\ell }-\phi _{m}}{2}\right) \\
\cos \left( \frac{\phi _{\ell }-\phi _{m}}{2}\right)%
\end{array}%
\right) , \\
\left(
\begin{array}{c}
L_{\ell +} \\
L_{\ell -}%
\end{array}%
\right) =\left(
\begin{array}{c}
-\frac{1-ie^{-i\zeta _{\ell }}}{1-ie^{i\zeta _{\ell }}} \\
1%
\end{array}%
\right)%
\end{array}%
\end{equation}
From the above results, and under the symmetric choice of parameters, we can
compute the overlap matrix, $M_{\alpha\beta}=\left\langle
\alpha\right|\beta\rangle$. Neglecting exponentially suppressed corrections,
we obtain the following form:
\begin{equation}
\begin{array}{c}
M_{\alpha\beta}=v\delta_{\alpha\beta}\left[\frac{B+\Delta_m\cos(\phi_{%
\alpha}-\phi_m)}{2(B^2-\Delta_m^2)}\right. \\
\left. +\frac{\Delta_r(B+\Delta_r)+\mu^2}{2\Delta_r(\Delta_r^2-B^2+\mu^2)}%
\right], \label{m}%
\end{array}%
\end{equation}
with $v$ being the spin-orbit velocity.
The coupling between the Majoranas could be calculated perturbatively by
considering the two domain walls juxtaposed. For instance, while the left
Majorana is an exact zero-energy eigenstate of the Hamiltonian
\[
{\mathcal{H}}_{\mathit{L}}={\mathcal{H}}_{\ell}\Theta(-x)+{\mathcal{H}}%
_m\Theta(x),
\]
the existence of the right segment of the wire perturbs this wave function,
with the perturbation potential being
\[
V_r=({\mathcal{H}}_r-{\mathcal{H}}_m)\theta(x-L).
\]
Similarly, we can write ${\mathcal{H}}={\mathcal{H}}_{\mathit{R}}+V_{\ell}$
with $V_{\ell}=({\mathcal{H}}_{\ell}-{\mathcal{H}}_m)\theta(-x)$. This
perturbation induces a hybridization matrix between the left Majorana and
the right Majorana:
\begin{equation}
h=\left(%
\begin{array}{cc}
\left\langle \mathit{L}\right|V_r\left| \mathit{L}\right\rangle=0 &
\left\langle \mathit{L}\right|V_r\left| \mathit{R}\right\rangle \\
\left\langle \mathit{R}\right|V_r\left| \mathit{L} \right\rangle &
\left\langle \mathit{R}\right|V_{\ell}\left| \mathit{R}\right\rangle=0%
\end{array}%
\right).
\end{equation}
In our case,
\[
\begin{array}{c}
V_r=\left[\left(\Delta_r\cos\phi_r-\Delta_m\cos\phi_m\right)\tau^x-\right.
\\
\left. \left(\Delta_r\sin\phi_r-\Delta_m\sin\phi_m\right)\tau^y-\mu_r\tau_z
\right]\Theta(x-L).%
\end{array}%
\]
The perturbation matrix we obtain is:
\begin{equation}
\begin{array}{c}
h=i v e^{i\nu\epsilon_{\alpha\beta}}\epsilon_{\alpha\beta}\left[%
e^{-L/\lambda^m_+}\sin\frac{\phi_r-\phi_m}{2}\sin\frac{\phi_{\ell}-\phi_m}{2}%
\right. \\
\left. +e^{-L/\lambda^m_-}\cos\frac{\phi_r-\phi_m}{2}\cos\frac{%
\phi_{\ell}-\phi_m}{2}\right].%
\end{array}%
\end{equation}
with $\nu$ an unimportant phase.
We arrive at the final answer for the Josephson coupling using Eq. (\ref{tbm}%
). The result indeed coincides with Eq. (5):
\begin{equation}
\begin{array}{c}
{\mathcal{H}}_{JM}=(2 f^{\dagger}f-1)\left[J_+ e^{-L/\lambda_+}\sin\frac{%
\phi_r-\phi_m}{2}\sin\frac{\phi_{\ell}-\phi_m}{2}\right. \\
\left. +J_-e^{-L/\lambda_-}\cos\frac{\phi_r-\phi_m}{2}\cos\frac{\phi_l-\phi_m%
}{2}\right] \\
=(2 f^{\dagger}f-1)\left(J_M\cos\frac{\phi_r-\phi_{\ell}}{2}+J_Z\cos\left(%
\frac{\phi_r+\phi_{\ell}}{2}-\phi_m\right)\right)%
\end{array}%
\end{equation}
with the constants $J_{\pm}$ being:
\begin{equation}
J_+=J_-\approx\frac{v}{\overline{M}_{rr}}.
\end{equation}
where $\overline{M}_{rr}=v\left[\frac{B}{2(B^2-\Delta_m^2)}+\frac{%
\Delta_r(B+\Delta_r)+\mu^2}{2\Delta_r(\Delta_r^2-B^2+\mu^2)}\right].$ is the
average of the overlap matrix [Eq. (\ref{m})] diagonal elements, dropping
the cosine term. The cosine term in the overlap will produce additional
harmonics of the Majorana-Josephson term but will not qualitatively change
the answer we obtained. The $J_{\pm}$ terms give rise to to the previously
explored Majorana-Josephson term, Eq. (1) and to the new zipper term,
Eq. (2).
\section{Numerical Calculation}
We now detail the procedure of our numerical calculation. In the Nambu
spinor basis $\Psi ^{T}=(\psi _{\uparrow },\psi _{\downarrow },\psi
_{\downarrow }^{\dagger },-\psi _{\uparrow }^{\dagger })$, the Bogoliubov-de
Gennes Hamiltonian for this system is
\begin{equation}
{\mathcal{H}}=v\hat{p}\sigma ^{z}\tau ^{z}-\mu \tau ^{z}+\Delta \left( \cos
\phi \tau ^{x}-\sin \phi \tau ^{y}\right) +B\sigma ^{x},
\end{equation}%
with $v$ the edge-state velocity, $\hat{p}$ the momentum, $B$ the Zeeman
energy, and $\sigma ^{a}$ and $\tau ^{a}$ Pauli matrices acting in the spin
and particle-hole sectors, respectively. We allow the chemical potential $%
\mu $, pairing amplitude $\Delta $, and superconducting phase $\phi $, to
vary spatially. In region $s$ (with $s=l,m,r$), the parameters $\left( \mu
,\Delta ,\phi \right) =\left( \mu _{s},\Delta _{s},\phi _{s}\right) $ are
constant. Without loss of generality, we assume $\phi _{m}=0$ to be a
reference of superconducting phase.
The Josephson effects in the TST junction has both bound states and
continuum contributions. In the following, we first present the procedure to
compute the exact interaction energy $E$ between two Majoranas, and then
provide the formalism to calculate the energy contribution from the
continuum.
\subsection{Bound state energy}
For TST configuration, there are two Majoranas at interfaces between
topological and trivial regions. The finite separation leads to a finite
interaction energy $E=E_{int}$ between these two Majoranas, with
spatial-dependent wave function $\Psi =\Psi \left( x\right) $ satisfying the
equation%
\begin{equation}
{\mathcal{H}}\Psi =E\Psi .
\end{equation}%
We will solve the interaction energy $E=E_{int}$ by matching the boundary
condition of the wave function.
First, we replace the momentum operator $\hat{p}$ with $-i\frac{\partial }{%
\partial x}$, and obtain the linear differential equation associated with
energy $E$%
\begin{equation}
\frac{\partial }{\partial x}\Psi \left( x\right) =\mathbf{G}_{E}\Psi \left(
x\right) ,
\end{equation}%
with $4\times 4$ matrix%
\begin{equation}
\mathbf{G}_{E}=i\frac{\mu }{v}\sigma ^{z}+\frac{\Delta }{v}\sigma ^{z}\left(
\cos \phi \tau ^{y}+\sin \phi \tau ^{x}\right) -\frac{B}{v}\sigma ^{y}\tau
^{z}+i\frac{E}{v}\sigma ^{z}\tau ^{z}.
\end{equation}%
In region $s$ (with $s=l,m,r$), the parameters $\left( \mu ,\Delta ,\phi
\right) =\left( \mu _{s},\Delta _{s},\phi _{s}\right) $ are constant, and
the matrix $G_{E}^{\left( s\right) }$ has eigensystem%
\begin{equation}
\mathbf{G}_{E}^{\left( s\right) }\vec{u}_{j}^{\left( s\right) }=\kappa
_{j}^{\left( s\right) }\vec{u}_{j}^{\left( s\right) }
\end{equation}%
with eigenvalues $\kappa _{j}^{\left( s\right) }$ and eigenvectors $\vec{u}%
_{j}^{\left( s\right) }$ for $j=1,\cdots ,4$ and $s=l,m,r$.
Then, we expand the four-component wave function $\Psi \left( x\right) $ in
terms of eigenvectors $\vec{u}_{j}^{\left( s\right) }$. We are interested in
the localized state with $E0$) and the two divergent modes ($\mathrm{Re}\kappa
_{3,4}^{\left( l\right) }<0$). Similarly, in the right region, there are two
localized modes ($\mathrm{Re}\kappa _{1,2}^{\left( r\right) }<0$) and the
two divergent modes ($\mathrm{Re}\kappa _{3,4}^{\left( r\right) }>0$). The
wave function with two localized Majoranas consists of localized modes%
\begin{equation}
\Psi \left( x\right) =\left\{
\begin{tabular}{ll}
$\sum\limits_{j=1,2}c_{j}^{\left( l\right) }e^{\kappa _{j}^{\left( l\right)
}x}\vec{u}_{j}^{\left( l\right) }$ & for $x\leq 0$ \\
$\sum\limits_{j=1,2}c_{j}^{\left( r\right) }e^{\kappa _{j}^{\left( r\right)
}\left( x-L\right) }\vec{u}_{j}^{\left( r\right) }$ & for $x\geq L$%
\end{tabular}%
\right. .
\end{equation}%
In order to match the coefficients associated with left and right regions,
we integrate the wavefunction over the middle region and obtain the condition%
\begin{equation}
\sum\limits_{j=1,2}c_{j}^{\left( r\right) }\vec{u}_{j}^{\left( r\right)
}=\Psi \left( L\right) =e^{\mathbf{G}_{E}^{\left( m\right) }L}\Psi \left(
0\right) =\sum\limits_{j=1,2}c_{j}^{\left( l\right) }e^{\mathbf{G}%
_{E}^{\left( m\right) }L}\vec{u}_{j}^{\left( l\right) },
\end{equation}%
which can be written as%
\begin{equation}
\mathbf{M}_{E}\left(
\begin{array}{c}
c_{1}^{\left( l\right) } \\
c_{2}^{\left( l\right) } \\
c_{1}^{\left( r\right) } \\
c_{2}^{\left( r\right) }%
\end{array}%
\right) =\left(
\begin{array}{c}
0 \\
0 \\
0 \\
0%
\end{array}%
\right)
\end{equation}%
with $4\times 4$ matrix%
\begin{equation}
\mathbf{M}_{E}=\left[
\begin{array}{cccc}
\left( e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}_{1}^{\left( l\right)
}\right) & \left( e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}_{2}^{\left(
l\right) }\right) & \left( -\vec{u}_{1}^{\left( r\right) }\right) & \left( -%
\vec{u}_{2}^{\left( r\right) }\right)%
\end{array}%
\right] .
\end{equation}
The necessary condition for non-zero solution is%
\begin{equation}
\det \mathbf{M}_{E}=0,
\end{equation}%
which can be used to numerically determine the interaction energy $E_{int}$.
As illustrated in Fig.~\ref{fig:Det}, the function $\det \mathbf{M}_{E}$
vanishes at $E=\pm E_{int}$. (There is a technical subtlety associated with
the fact that $\mathbf{G}_{E}^{\left( l,r\right) }$ is \emph{not} a
Hermitian matrix. For some fixed values of $E$, the eigenvalues of $\mathbf{G%
}_{E}^{\left( l,r\right) }$ have multiplicity larger than one, and the
eigenvector $\vec{u}_{j}^{\left( l,r\right) }$ might be a zero vector, which
may also lead to spurious solutions with vanishing $\det \mathbf{M}_{E}$.
This issue can be resolved by using a polynomial discriminant to identify
and remove these spurious solutions.)
\begin{figure}[tbp]
\centering
\includegraphics[width=8.7cm]{fig_Det.eps}
\caption[fig:Det]{{}The function $\det M_{E}$ vanishes when $E=\pm E_{int}$,
which can be used to numerically find the interaction energy between the
Majoranas.}
\label{fig:Det}
\end{figure}
\subsection{Continuum contribution}
We now consider the energy contribution from the continuum. The continuum
states can be characterized by the scattering matrix $\mathbf{S}_{E}=\mathbf{%
S}_{E}\left( \phi _{l},\phi _{r}\right) $, which can be computed by matching
the boundary conditions for all incoming and outgoing modes. Once we know
the scattering matrix, we can use the Fumi's sum rule to compute the
continuum contribution to the system energy \cite{Mahan00,Akkermans91}%
\begin{equation}
W\left( \phi _{l},\phi _{r}\right) =\int_{E_{gap}}^{\infty }\frac{dE}{2\pi i}%
\ln \left[ \det \left[ \mathbf{S}_{E}\left( \phi _{l},\phi _{r}\right) %
\right] \right] .
\end{equation}
The continuum contribution consists of many Fourier components%
\begin{equation}
W\left( \phi _{l},\phi _{r}\right) =\frac{1}{2}\sum_{n_{l},n_{r}=-\infty
}^{\infty }W_{n_{l},n_{r}}\cos (n_{l}\phi _{\ell }+n_{r}\phi _{r})
\end{equation}%
with $W_{n_{l},n_{r}}=W_{-n_{l},-n_{r}}$. Then conventional Josephson terms
are $J_{L/R}=W_{1,0}$ and $W_{0,1}$, and the even harmonics of the zipper
terms are $J_{Z,2n}=W_{n,n}$. In the following, we provide the formalism to
compute the scattering matrix $S_{E}\left( \phi _{l},\phi _{r}\right) .$
For energy $E>E_{gap}^{\left( l,r\right) }$, there are the propagating modes
($\mathrm{Re}\kappa _{j}^{\left( l,r\right) }=0$), with momentum $%
p_{j}^{\left( l,r\right) }=\mathrm{Im}\kappa _{j}^{\left( l,r\right) }$.
Suppose there are four incoming modes $\left( \vec{u}_{1}^{\left( l\right) },%
\vec{u}_{2}^{\left( l\right) },\vec{u}_{1}^{\left( r\right) },\vec{u}%
_{2}^{\left( r\right) }\right) $ with $p_{1,2}^{\left( l\right) }>0$ and $%
p_{1,2}^{\left( r\right) }<0$, and four outgoing modes $\left( \vec{u}%
_{3}^{\left( l\right) },\vec{u}_{4}^{\left( l\right) },\vec{u}_{3}^{\left(
r\right) },\vec{u}_{4}^{\left( r\right) }\right) $ with $p_{3,4}^{\left(
l\right) }<0$ and $p_{3,4}^{\left( r\right) }>0$. The wave function can be
written as a linear combination of all these modes%
\begin{equation}
\Psi \left( x\right) =\left\{
\begin{tabular}{ll}
$\sum\limits_{j=1,\cdots ,4}c_{j}^{\left( l\right) }e^{\kappa _{j}^{\left(
l\right) }x}\vec{u}_{j}^{\left( l\right) }$ & for $x\leq 0$ \\
$\sum\limits_{j=1,\cdots ,4}c_{j}^{\left( r\right) }e^{\kappa _{j}^{\left(
r\right) }\left( x-L\right) }\vec{u}_{j}^{\left( r\right) }$ & for $x\geq L$%
\end{tabular}%
\right. .
\end{equation}%
In order to match the coefficients associated with left and right regions,
we integrate the wavefunction over the middle region and obtain the
condition $\Psi \left( L\right) =e^{\mathbf{G}_{E}^{\left( m\right) }L}\Psi
\left( 0\right) $.
The relation between the amplitudes of incoming and outgoing modes is
\begin{equation}
\mathbf{M}_{E,in}\left(
\begin{array}{c}
c_{1}^{\left( l\right) } \\
c_{2}^{\left( l\right) } \\
c_{1}^{\left( r\right) } \\
c_{2}^{\left( r\right) }%
\end{array}%
\right) =\mathbf{M}_{E,out}\left(
\begin{array}{c}
c_{3}^{\left( l\right) } \\
c_{4}^{\left( l\right) } \\
c_{3}^{\left( r\right) } \\
c_{4}^{\left( r\right) }%
\end{array}%
\right)
\end{equation}%
with $4\times 4$ matrices%
\begin{equation}
\mathbf{M}_{E,in}=\left[
\begin{array}{cccc}
\left( e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}_{1}^{\left( l\right)
}\right) & \left( e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}_{2}^{\left(
l\right) }\right) & \left( -\vec{u}_{1}^{\left( r\right) }\right) & \left(
-\vec{u}_{2}^{\left( r\right) }\right)
\end{array}%
\right] ,
\end{equation}%
\begin{equation}
\mathbf{M}_{E,out}=\left[
\begin{array}{cccc}
\left( -e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}_{3}^{\left( l\right)
}\right) & \left( -e^{\mathbf{G}_{E}^{\left( m\right) }L}\vec{u}%
_{4}^{\left( l\right) }\right) & \left( \vec{u}_{3}^{\left( r\right)
}\right) & \left( \vec{u}_{4}^{\left( r\right) }\right)
\end{array}%
\right] .
\end{equation}%
The scattering relation is%
\begin{equation}
\mathbf{S}_{E}\left(
\begin{array}{c}
\left( p_{1}^{\left( l\right) }\right) ^{1/2}c_{1}^{\left( l\right) } \\
\left( p_{2}^{\left( l\right) }\right) ^{1/2}c_{2}^{\left( l\right) } \\
\left( p_{1}^{\left( r\right) }\right) ^{1/2}c_{1}^{\left( r\right) } \\
\left( p_{2}^{\left( r\right) }\right) ^{1/2}c_{2}^{\left( r\right) }%
\end{array}%
\right) =\left(
\begin{array}{c}
\left( -p_{3}^{\left( l\right) }\right) ^{1/2}c_{3}^{\left( l\right) } \\
\left( -p_{4}^{\left( l\right) }\right) ^{1/2}c_{4}^{\left( l\right) } \\
\left( -p_{3}^{\left( r\right) }\right) ^{1/2}c_{3}^{\left( r\right) } \\
\left( -p_{4}^{\left( r\right) }\right) ^{1/2}c_{4}^{\left( r\right) }%
\end{array}%
\right)
\end{equation}%
with scattering matrix%
\begin{equation}
\fbox{$\mathbf{S}_{E}=\mathbf{P}_{out}^{1/2}\cdot \mathbf{M}%
_{E,out}^{-1}\cdot \mathbf{M}_{E,in}\cdot \mathbf{P}_{in}^{-1/2},$}
\end{equation}%
where $\mathbf{P}_{in}=\mathrm{Diag}\left[ p_{1}^{\left( l\right)
},p_{2}^{\left( l\right) },p_{1}^{\left( r\right) },p_{2}^{\left( r\right) }%
\right] $ and $\mathbf{P}_{out}=-\mathrm{Diag}\left[ p_{3}^{\left( l\right)
},p_{4}^{\left( l\right) },p_{3}^{\left( r\right) },p_{4}^{\left( r\right) }%
\right] $. The requirement of conservation of current is%
\begin{equation}
\sum_{j}p_{j}^{\left( l\right) }\left\vert c_{j}^{\left( l\right)
}\right\vert ^{2}=\sum_{j}p_{j}^{\left( r\right) }\left\vert c_{j}^{\left(
r\right) }\right\vert ^{2},
\end{equation}%
which ensures the unitarity of the scattering matrix%
\begin{equation}
\mathbf{S}_{E}^{\dag }\mathbf{S}_{E}=I.
\end{equation}%
Hence, $\det \left[ \mathbf{S}_{E}\right] =e^{i2\delta _{E}}$ and $\frac{1}{%
2\pi i}\ln \left[ \det \left[ \mathbf{S}_{E}\left( \phi _{l},\phi
_{r}\right) \right] \right] =\frac{1}{\pi }\delta _{E}\left( \phi _{l},\phi
_{r}\right) $. Numerically, we just need to compute the quantity $\delta
_{E}\left( \phi _{l},\phi _{r}\right) $ and the integral%
\begin{equation}
W\left( \phi _{l},\phi _{r}\right) =\int_{E_{gap}}^{\infty }\frac{dE}{\pi }%
\delta _{E}\left( \phi _{l},\phi _{r}\right) .
\end{equation}%
The continuum contribution $W\left( \phi _{l},\phi _{r}\right) $ has $2\pi $
periodicity in both $\phi _{l}$ and $\phi _{r}$, with Fourier decomposition
of $W\left( \phi _{l},\phi _{r}\right) =\frac{1}{2}\sum_{n_{l},n_{r}=-\infty
}^{\infty }W_{n_{l},n_{r}}\cos (n_{l}\phi _{\ell }+n_{r}\phi _{r})$, with
Fourier coefficients of $W_{n_{l},n_{r}}$. The relevant Fourier components
are $J_{L/R}=W_{1,0}=W_{0,1}$, and $J_{Z,2n}=W_{n,n}$ for $n=1,2,\cdots $.
There is one subtle issue in the computation of the scattering matrix. There
are four propagating modes for $E>\left\vert B^{\left( l,r\right) }+\sqrt{%
\left( \mu ^{\left( l,r\right) }\right) ^{2}+\left( \Delta ^{\left(
l,r\right) }\right) ^{2}}\right\vert $, but there are two propagating modes
and two localized modes for $\left\vert B^{\left( l,r\right) }+\sqrt{\left(
\mu ^{\left( l,r\right) }\right) ^{2}+\left( \Delta ^{\left( l,r\right)
}\right) ^{2}}\right\vert >E>E_{gap}^{\left( l,r\right) }$ $=\left\vert
B^{\left( l,r\right) }-\sqrt{\left( \mu ^{\left( l,r\right) }\right)
^{2}+\left( \Delta ^{\left( l,r\right) }\right) ^{2}}\right\vert $. In the
latter case, we need to compute the effective scattering matrix that are
projected to the subspace spanned by the propagating modes.
\bibliographystyle{apsrev}
\bibliography{majorana}
\end{document}