# Experimental Demonstration of Quantum Tunneling in IBM Quantum Computer

Narendra N. Hegade1 E-mail<EMAIL_ADDRESS>Bikash K. Behera2 &
Prasanta K. Panigrahi2 E-mail<EMAIL_ADDRESS>author.E-mail<EMAIL_ADDRESS>
###### Abstract

Quantum computers are the promising candidates for simulation of large quantum
systems, which is a daunting task to perform in a classical computer. Here we
report the experimental realization of quantum tunneling of a single particle
through different types of potential barriers by simulating it in the IBM
quantum computer, which here acts as a universal quantum simulator. We
consider two-qubit and three-qubit systems for visualizing the tunneling
process illustrating its unique quantum nature. We clearly observe the
tunneling and oscillations of the particles in a step potential and in double-
and multi-well potentials through our experimental results. The proposed
quantum circuits and simulational techniques used here can be extended for
observing the tunneling phenomena for multi-particle systems in different
shaped potentials.

Department of Physics, National Institute of Technology Silchar, Silchar,
788010, India

Department of Physical Sciences, Indian Institute of Science Education and
Research Kolkata, Mohanpur, 741246, West Bengal, India

Keywords:Quantum simulation, Quantum tunneling, IBM Quantum Experience

##  1 Introduction

Quantum simulation is one of the problems that a quantum computer could
perform more efficiently than a classical computer as it provides a
significant improvement in computational resources [1, 2, 3, 4, 5, 6, 7, 8,
9]. It has been applied to a wide range of areas of physics like quantum many-
body theory [10, 11, 12, 13], quantum entanglement [14, 15], quantum phase
transitions [16, 17], molecular physics [18, 19, 20, 21, 22] etc. Algorithms
have been used in simulating many quantum field theoretic problems [23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34] where Hamiltonian of the system splits
into kinetic and potential energy operators which are then simulated using
Trotter’s formula [35, 36]. Experimental realizations of quantum simulation
have already been made in systems like NMR [10, 37, 38, 11, 15, 39, 12, 20],
ion-trap [40, 41, 42, 43, 44], atomic [16, 45] and photonic [46, 47] quantum
computers. The current status of this field can be found out from these recent
papers [48, 49, 50, 51, 52, 53, 54, 55, 56, 57].

Quantum tunneling acts as one of the exciting and unique fundamental phenomena
in quantum mechanics. It has been observed in superconducting Cooper pairs
[58] and in modern technologies such as narrow p-n junctions [59] and scanning
tunneling microscope [60]. Important science puzzles like lattice quantum
chromodynamics can be solved using this tunneling simulation approach [61]. A
number of digital simulation on quantum tunneling has been performed on
classical computers [62, 63] and photonic systems [64]. This type of
simulation has remained untested in a quantum computer due to the requirement
of large number of ancillary qubits and quantum gates. Recently, an algorithm
proposed by Sornborger [65] illustrates the simulation with no ancillary
qubits and a small number of quantum gates which motivates the possibility of
simulating in today’s quantum computer consisting of a few number of qubits.
Feng _et al._ [66] have demonstrated the tunneling effect for two-qubit and
three-qubit systems using NMR quantum computer. Ostrowski [67] has also
explicated this process on a rectangular potential by digital simulation.
Here, in the present work, we illustrate the simulation of quantum tunneling
using IBM’s 5-qubit quantum computer “IBM Q 5 Tenerife” and 14-qubit quantum
computer “IBM Q 14 Melbourne”. Using two-qubit and three-qubit systems, CNOT
gates, a set of Hadamard and controlled phase gates, we were able to simulate
the tunneling process of a single particle in a step potential, double-well
potential and multi-well potential. We utilize the IBM quantum experience’s
QISKit, to simulate the tunneling Hamiltonian, using which a number of
research works have been performed [68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94].

##  2 Results

###  2.1 Theoretical Protocol.

A digital quantum simulator is a controllable quantum system which can be used
to efficiently simulate the dynamics of any other quantum system with local
interactions. The digital quantum simulation [95, 96, 97, 98, 99, 100, 101,
102] consists of three main steps: initial state preparation, time evolution,
and measurement. In this work, we are particularly interested in the time
evolution of the system, so our main idea is to decompose the evolution of the
system in terms of single and two-qubit gates. The Schr$\ddot{o}$dinger’s
equation, for a single particle moving in an one-dimensional space, is
expressed as

$i\frac{\partial}{\partial t}|\psi(x,t)\rangle=\hat{H}|\psi(x,t)\rangle$ (1)

where $\hat{H}=\hat{K}+\hat{V}$, $\hat{K}$ and $\hat{V}$ are kinetic and
potential energy operators respectively. Here we set the value of $\hbar$ to 1
throughout the manuscript. The time evolution of the wave function of the
system can be given as

$\displaystyle|\psi(x,t+\Delta t)\rangle=$ $\displaystyle e^{-i\hat{H}\Delta
t}|\psi(x,t)\rangle$ $\displaystyle=$ $\displaystyle
e^{-i(\hat{K}+\hat{V})\Delta t}|\psi(x,t)\rangle$ (2)

Using first-order Suzuki-Trotter’s formula [103, 104, 105], the exponential
operator can be decomposed up to order $O(\Delta t)$ as follows.

$e^{-i(\hat{K}+\hat{V})\Delta t}\approx e^{-i\hat{K}\Delta
t}e^{-i\hat{V}\Delta t}$ (3)

We discretize the continuous coordinate space x on a lattice (with spacing
$\Delta l$) within the boundary region ($0<x<L$) with a periodic boundary
condition $\psi(x+L,t)=\psi(x,t)$. The wave function then can be mapped to a
n-qubit register as,

$|\psi(x,t)\rangle\rightarrow\sum_{k=0}^{2^{n}-1}\psi(x_{k},t)|k\rangle$ (4)

Here $|k\rangle$ represents the particle location corresponding to binary
number k, and $x_{k}=(k+\frac{1}{2})\Delta l$, $\Delta l=\frac{L}{2^{n}}$. The
mapping given in expression (4) can be a good approximation for large value n.
As an example, for n=3 qubits, the binary representation of a wavefunction on
a 1D lattice is represented in Fig. 1. Here the wave function mapped to the
3-qubit register can be written as,

$\displaystyle|\psi(x,t)\rangle$
$\displaystyle=\psi(x_{0},t)|0\rangle+\psi(x_{1},t)|1\rangle+\psi(x_{2},t)|2\rangle...+\psi(x_{7},t)|7\rangle$

where, $\sum_{k=0}^{7}|\psi(x_{k},t)|^{2}=1$, and
$x_{k}=(k+\frac{1}{2})\frac{L}{8}$. The quantum states $|0\rangle$,
$|1\rangle$,…,$|7\rangle$ are mapped to 3-qubit quantum states as
$|000\rangle$, $|001\rangle$,…,$|111\rangle$ respectively.

Figure 1: Binary representation of a wavefunction on a 1D lattice for 3-qubit
case. The 3-qubit quantum states $|000\rangle$, $|001\rangle$,…,$|111\rangle$
represent the binary form of the states $|0\rangle$, $|1\rangle$,…,$|7\rangle$
respectively.

Implementing potential energy operator: The potential energy operator plays a
key role for introducing and studying various types of potential structures.
Here we implement step-well, double-well and multi-well potentials with the
help of single qubit rotation operators with no ancillary qubits. For step
potential, we apply a single-qubit Z-rotation gate on the highest order qubit
[65].

$e^{-iV\Delta t}=e^{-iv{\sigma_{z}}^{n-1}\Delta t}=e^{-iv\sigma_{z}\Delta
t}\otimes I\otimes I...,$ (6)

where $v$ is a parameter, $\sigma_{z}$ is the Pauli $z$-matrix, and $n-j$
indicates the application of the operator ($e^{-iv\sigma_{z}\Delta t}$) on the
$j$th qubit. The double-well potential can be implemented by applying the
Z-rotation gate on the second highest order qubit.

$e^{-iV\Delta t}=e^{-iv{\sigma_{z}}^{n-2}\Delta t}=I\otimes
e^{-iv\sigma_{z}\Delta t}\otimes I...$ (7)

Similarly, by applying Z-rotation operator on the next highest order qubits we
can implement multi-well potentials. It can be observed that just one single
qubit operation can reduce the complexity in the quantum circuit by replacing
a large number of gates and ancillary qubits.

$e^{-iV\Delta t}=I\otimes I\otimes e^{-iv\sigma_{z}\Delta t}...$ (8)

Implementing kinetic energy operator: Here the kinetic energy operator
$\hat{K}$ can be expressed in terms of momentum operator as,
$\hat{K}=\frac{{\hat{p}}^{2}}{2m}$. For finding the quantum circuit for
kinetic energy operator, let’s discretize the wave function of momentum as

$\ket{\chi(p,t)}=\sum_{l=0}^{2^{n}-1}\chi(p_{l},t)\ket{l}$ (9)

where $\chi(p,t)$ is the wave function in the momentum operator
representation. The eigen values of momentum $p_{l}$ are given by

$p_{l}=\begin{cases}\frac{2\pi}{2^{n}}l&0\leqslant l\leqslant 2^{n-1}\\\
\frac{2\pi}{2^{n}}(2^{n-1}-l)&2^{n-1}<l<2^{n}\end{cases}$ (10)

In the momentum representation, the diagonal operator $\hat{P}$ is written as

$\hat{P}=\sum_{l=0}^{2^{n-1}}\frac{2\pi}{2^{n}}l\ket{l}\bra{l}+\sum_{l=2^{n-1}+1}^{2^{n}-1}\frac{2\pi}{2^{n}}(2^{n-1}-l)\ket{l}\bra{l}$
(11)

The kinetic energy operator is diagonal in momentum representation. It can be
written in the coordinate representation by using quantum Fourier transform as
follows,

$\displaystyle e^{-i(\hat{K}+\hat{V})\Delta t}\approx$
$\displaystyle(QFT)e^{-i\frac{\hat{P}^{2}}{2m}\Delta
t}(QFT^{-1})e^{-i\hat{V}\Delta t}$ $\displaystyle\approx$
$\displaystyle(QFT)D(QFT^{-1})P$ (12)

where

$\displaystyle QFT=\frac{1}{\sqrt{2^{n}}}\sum_{l,k=0}^{2^{n}-1}e^{\frac{2\pi
ilk}{2^{n}}}|l\rangle\langle k|$ $\displaystyle
D=e^{-i\frac{\hat{P}^{2}}{2m}\Delta t}$ $\displaystyle P=e^{-i\hat{V}\Delta
t}$ (13)

The equivalent quantum circuit for the Fourier transformation operator $QFT$
can be realized using a series of Hadamard and controlled-phase gates [106].
Thus after a small time interval $\Delta t$, the time evolution of the system
can be implemented by using Eqs. (2.1), (3), (4) and (2.1).

$\sum_{k=0}^{2^{n}-1}\psi(x_{k},t+\Delta
t)|k\rangle=(QFT)D(QFT^{-1})P\sum_{k=0}^{2^{n}-1}\psi(x_{k},t)|k\rangle$ (14)

The explicit construction of quantum circuit for $QFT$, D, $QFT^{-1}$ and P
are detailed in Methods section.

###  2.2 Experimental Procedures and Results.

The quantum circuit for one-time step evolution using two qubits is shown in
Fig. 11. The step potential is implemented by the quantum operation
$P=e^{-iv\sigma_{z}^{1}\Delta t}$ acting on the highest order qubit.
Similarly, the double-well potential is implemented by the quantum operation
$P=e^{-iv\sigma_{z}^{0}\Delta t}$ acting on the lowest order qubit. In our
experiment, we set the parameter $v$=0 for free particle, and $v=50$ for step
potential, double-well potential and $v=10$ for multi-well potential. We set
the time interval $\Delta t=0.1$. Mass of the particle is taken to be 0.5.

Figure 2: Free particle ($v$=0) probability distribution. Case (a): Simulated
result obtained from “IBM Q QASM simulator” for free particle evolution using
two-qubit system. Initially the particle was confined at $\ket{00}$ state, the
probability distribution of the particle slowly spreads out after six time
steps. Case (b): Experimentally obtained result from “ibmqx4” for probability
distribution of free particle. Similar observation has been obtained as case
(a). The horizontal axis represents the number of time steps and the vertical
axis represents the corresponding probabilities. Figure 3: (a) and (b) are the
real and imaginary parts of the density matrix elements for the initial state
$\ket{00}$ of the free particle obtained from the “IBM Q QASM simulator”; (c)
and (d) are the experimentally reconstructed density matrix elements obtained
from “ibmqx4” chip. Figure 4: (a) and (b) are the real and imaginary parts of
the density matrix elements for the final state (i.e after six time steps) of
the free particle obtained from the QASM simulator; (c) and (d) are obtained
from “ibmqx4” device.

In Fig. 2, the time evolution of a free particle with potential $v$ = 0 is
depicted. Initially, the particle was confined at $\ket{00}$ state. It can be
seen that after a number of time steps, the particle probability distribution
spreads over other basis states. Both the theoretical and experimental results
are illustrated in Fig. 2 (a) and (b) cases respectively. To compare the
result obtained from the quantum processor “ibmqx4” and the simulated result
from “IBM Q QASM simulator” quantum state tomography is performed. In both the
cases, the number of shots taken were 8192. The density matrix elements for
the initial state $\ket{00}$ and the final state after six time steps are
depicted in Fig. 3 and 4 respectively. The experimental fidelity for the
initial state is $96.79\%$, and for the final state after six time step is
$93.83\%$.

Figure 5: Particle probability distribution in a step-well potential ($v$=50).
Case (a,b): Simulated results obtained from “IBM Q QASM simulator”
illustrating the single particle dynamics in a step potential. Case (c,d):
Experimentally observed result from “IBM Q 14 Melbourne” showing the hopping
of the particle in a step potential of finite length. Case (e,f):
Experimentally observed result from “ibmqx4” depicting the hopping of the
particle in a step potential of finite length. Here, the barriers are located
at $\ket{10}$ and $\ket{11}$ states, From the probability distribution in all
the cases, it can be observed that, as time evolves, the particle which is
initially confined at $\ket{00}$ state hops to the $|01\rangle$ state, also
the particle has some non zero probability to be found at classically
forbidden region.

In Fig. 5, we illustrate the tunneling phenomena of a single particle in a
finite length step potential of $v$=50. The potential barriers are situated at
$\ket{10}$ and $\ket{11}$ states. At time, t=0, the particle is located at
$\ket{00}$ state, as time evolves we can observe the probability distribution
of the particle spreads to the state $\ket{01}$. Notably, it can be observed
that after six time steps, the particle has some non-zero probability to be
found at classically forbidden region $\ket{10}$ and $\ket{11}$. Which clearly
confirms the tunneling of particle through the potential barrier. We performed
this experiment on “IBM Q QASM simulator”, “ibmqx4” and “IBM Q 14 Melbourne”.
Both the theoretical and experimental results are explicated in Fig. 5. In all
the cases, the number of shots taken to perform the experiment is 8192.

Figure 6: Particle probability distribution in double-well potential ($v$=50).
Case (a) : Theoretically calculated result of quantum tunneling in double-well
potential using two-qubit system. Initially the particle was confined at
$\ket{01}$ state. After six time step the particle is found to slowly tunnel
to the $\ket{11}$ state. Case (b): Experimentally observed result of quantum
tunneling for 6 time steps. Figure 7: The density matrix elements for the
initial state $\ket{01}$ of a particle in double-well potential. (a) and (b)
are the real and imaginary parts of the density matrix obtained from QASM
simulator, (c) and (d) are from “IBM Q 14 Melbourne”. Figure 8: (a) and (b)
are the real and imaginary parts of the density matrix elements for the final
state (i.e after six time steps) of the particle in double-well potential
obtained from the QASM simulator; (c) and (d) are obtained from “IBM Q 14
Melbourne”. Figure 9: Particle probability distributions for ten time steps in
a multi-well potential ($v$ = 10). The potential wells are at four sites,
$|000\rangle$, $|010\rangle$, $|100\rangle$ and $|110\rangle$. Initially, the
particle is confined at $\ket{100}$ state, as the time evolves the particle
tunnels through the barriers situated at $\ket{001},\ket{011},\ket{101}$ and
$\ket{111}$ states. After 10 time steps, the particle is most probable to be
found at $|010\rangle$ state.

The tunneling of particle in a double-well potential ($v$=50) is shown in Fig.
6. The wells are located at $\ket{01}$ and $\ket{11}$ states and the barriers
are at $\ket{00}$ and $\ket{10}$ states. At time, t=0, the particle is
confined in the potential well situated at $\ket{01}$ state, as the time
evolves we can clearly observe the tunneling of particle from the potential
well located at $\ket{01}$ state to the well at $\ket{11}$ state. Both from
the theoretical and experimental data (Fig. 6 (a) and (b)), it is observed
that the particle oscillates between the two wells which signifies tunneling
through the potential barriers at the other two places. To compare the result
obtained from the quantum processor “IBM Q 14 Melbourne” and the ideal QASM
simulator we performed quantum state tomography. To perform the experiment,
the number of shots are chosen to be 8192. The density matrix elements for the
initial state $\ket{01}$ and the final state after six time steps are depicted
in Fig. 7 and 8 respectively. The experimental fidelity for the initial state
is $95.61\%$, and for the final state after six time step is $95.18\%$. Fig. 9
shows the tunneling of the particle in a multi-well potential implemented on
the “IBM Q QASM simulator”. The wells are situated at
$\ket{000},\ket{010},\ket{100}$ and $\ket{110}$ positions. Initially, at time
t=0, the particle is trapped inside one of the above wells situated at
$\ket{100}$ state, as time evolves the particle is tunneled through the
barriers situated at $\ket{001},\ket{011},\ket{101}$ and $\ket{111}$ states.
At different time steps, the tunneling of the particle through all the
potential barriers is observed. After 10 time steps, the particle is most
likely to be found at $|010\rangle$ state.

##  3 Discussion

To conclude, we have experimentally demonstrated here the quantum tunneling
phenomena of a single particle in a step potential, double-well potential and
multi-well potential. We have designed the equivalent quantum circuit for the
Hamiltonian of the given system in the real quantum processors “ibmqx4” and
“IBM Q 14 Melbourne”. We have shown the architecture of this processor with
important device parameters including gate errors and readout errors. We have
illustrated the tunneling process by running the quantum circuit for six time
steps in two-qubit case and ten time steps in three-qubit case respectively.
After comparing the theoretical and experimental results it is concluded that
the tunneling process has been properly carried out with IBM’s quantum
processors. The basic circuit mechanism used here for demonstrating quantum
tunneling can be extended to observe quantum tunneling of multi-particle
systems in different shaped potential wells.

##  4 Methods

Figure 10: Chip layout of 5-qubit (ibmqx4) and 14-qubit (IBM Q 14 Melbourne)
quantum processors.

We investigate quantum simulation of the tunneling process using the IBM’s
5-qubit quantum processor “ibmqx4” and 14-qubit quantum processor “IBM Q 14
Melbourne”, whose layout is depicted in Fig. 10. The connectivities between
the qubits in both the processors are shown. The devices are stored in a
temperature order of mK. Tables 1 & 2 illustrate the experimental parameters
of the devices, where frequency, coherence time ($T_{1}$), relaxation time
($T_{2}$), gate error and readout errors (RE) of each qubit are listed.

Figure 11: Quantum circuit for two qubit system. Case (A): The quantum circuit
for one time step ($\Delta t=0.1$) evolution of tunneling Hamiltonian using
two qubits is depicted. Case (B): The inverse quantum Fourier transform
$QFT^{-1}$ is built using two Hadamard gates and one controlled-rotation gate
$U_{1}(\pi/2)$. Case (C): The decomposition of diagonal operator D is
represented, which requires two single qubit $R_{z}$ rotation gates and one
controlled-rotation gate $U_{1}(\pi^{2}/10)$.

Circuit construction for two-qubit simulation. The explicit construction of
$QFT^{-1}$, D and P are illustrated in Fig. 11. $QFT^{-1}$ is prepared by two
Hadamard gates and one controlled-phase gate,
$QFT^{-1}=H_{0}C-U_{1}{(\pi/2)_{10}}H_{1}$, where $H_{i}$ denotes the
application of Hadamard gate on the $q[i]$ qubit, $C-U_{1}{(\theta)_{ij}}$
denotes controlled-phase gate of angle $\theta$, where $q[i]$ is the control
qubit and $q[j]$ is the target qubit, $U_{1}(\theta)=[[1,0],[0,e^{i\theta}]]$.
For two-qubit simulation, the diagonal elements of $\hat{P}$ are calculated to
be 0, $\pi/2$, $\pi$ and $-\pi/2$ (See Eq. (11)). Then, the kinetic evolution
operator is given by

$\displaystyle e^{-i\hat{K}\Delta t}=$ $\displaystyle
QFT^{-1}e^{-i\hat{P^{2}}\Delta t}QFT$ $\displaystyle=$ $\displaystyle
QFT^{-1}\Phi_{10}Z_{0}Z_{1}QFT$ (15)

The diagonal operator D can be expressed as a product of the following
operators,

$D=\Phi_{10}Z_{0}Z_{1}$ (16)

where $Z_{1}=e^{-i\gamma c_{0}\sigma_{z}^{0}\Delta t}$, $Z_{0}=e^{-i\gamma
c_{1}\sigma_{z}^{1}\Delta t}$ and $\phi_{10}=e^{-i\gamma
c_{2}diag(1,1,1,-1)_{10}\Delta t}$. Here, $\phi_{10}$ has been prepared with
the $C-U_{1}(\pi^{2}/10)_{10}$ gate. Here $Z_{i}$ means application of Z
operation on $q[i]$ qubit. The constant values in (16) are obtained to be
$\gamma=\frac{\pi^{2}}{8}$, $c_{0}=-1$, $c_{1}=-4$ and $c_{2}=4$ [65]. The
time step is taken to be $\Delta t=0.1$. The quantum circuit for D operator
can be given as,
$D=R_{z}(-\pi^{2}/40)_{1}R_{z}(-\pi^{2}/10)_{0}C-U_{1}(\pi^{2}/10)_{10}$,
where $R_{z}(\theta)_{i}$ is the rotation operator about Z-axis by an angle
$\theta$ on the $q[i]$th qubit (See Fig. 11). The $R_{z}(\theta)$ can be
expressed using standard “ibmqx4” gates as
$R_{z}(\theta)=HU_{3}(\theta,\phi=-\frac{\pi}{2},\lambda=\frac{\pi}{2})H$,
where
$U_{3}(\theta,\phi,\lambda)=[[cos(\theta/2),-e^{i\lambda}sin(\theta/2)],[e^{i\phi}sin(\theta/2),e^{i(\lambda+\phi)}cos(\theta/2)]]$.
The potential operator P is prepared by the application of
$HU_{3}(\theta=5,\phi=-\frac{\pi}{2},\lambda=\frac{\pi}{2})H$ operator. It has
been applied on the qubits $q[0]$ and $q[1]$ in case of step potential and
double well potential respectively. The controlled phase gates used in the
experiment, are directly applied using the codes available in QISKit (See
supplementary Information).

Figure 12: Quantum circuit for three-qubit system. Case(A): The quantum
circuit for one time step ($\Delta t=0.1$) evolution of tunneling Hamiltonian
using three qubits is depicted. Case(B): The inverse quantum Fourier transform
QFT-1 is built using three Hadamard gates and three controlled phase gates of
$\pi/2$ and $\pi/4$ angle. Case(C): The kinetic energy operator D is
decomposed into rotation and controlled-phase gates of
$-\pi^{2}/160,-\pi^{2}/40,-\pi^{2}/10,\pi^{2}/10$ and $\pi^{2}/20$ angles.

Circuit construction for three-qubit simulation. The explicit construction of
$QFT^{-1}$, D and P are illustrated in Fig. 12. $QFT^{-1}$ is prepared by
three Hadamard gates and three controlled-phase gates,
$QFT^{-1}=H_{0}C-U_{1}(\pi/2)_{10}H_{1}C-U_{1}(\pi/4)_{20}C-U_{1}(\pi/2)_{21}H_{2}$.
For three-qubit simulation, the diagonal elements of $\hat{P}_{p}$ are
calculated to be 0, $\pi/4$, $\pi/2$, $3\pi/4$, $\pi$, $-\pi/4$, $-\pi/2$ and
$-3\pi/4$ (See Eq. (11)). The decomposition of D operator is as follows,
$D=R_{z}(-\pi^{2}/160)_{2}R_{z}(-\pi^{2}/40)_{1}R_{z}(-\pi^{2}/10)_{0}C-U_{1}(\pi^{2}/10)_{21}C-U_{1}(\pi^{2}/20)_{20}C-U_{1}(-\pi^{2}/40)_{10}$.
The diagonal operator D can be expressed as a product of the following
operators, $Z_{0}=e^{-i\gamma c_{0}\sigma_{z}^{0}\Delta t}$,
$Z_{1}=e^{-i\gamma c_{1}\sigma_{z}^{1}\Delta t},Z_{2}=e^{-i\gamma
c_{2}\sigma_{z}^{2}\Delta t},\phi_{21}=e^{-i\gamma
c_{3}diag(1,1,1,-1)_{21}\Delta t},\phi_{20}=e^{-i\gamma
c_{4}diag(1,1,1,-1)_{20}\Delta t},\phi_{10}=e^{-i\gamma
c_{5}diag(1,1,1,-1)_{10}\Delta t}$. The constant values are obtained to be
$\gamma=-\frac{\pi^{2}}{32\sqrt{2}}$, $c_{0}=-1.42$, $c_{1}=-5.66$ and
$c_{2}=-22.63$, $c_{3}=-22.63$, $c_{4}=11.31$ and $c_{5}=-5.66$ [65]. It is to
be noted that the time step is taken to be $\Delta t=0.1$. The quantum circuit
for D operator can be given as,
$D=R_{z}(-\pi^{2}/40)_{1}R_{z}(-\pi^{2}/10)_{0}C-U_{1}{(\pi^{2}/10)_{10}}$,
where $R_{z}$ is the rotation operator about Z-axis (See Fig. 12). The
potential operator P prepared by
$HU_{3}(\theta,\phi=-\frac{\pi}{2},\lambda=\frac{\pi}{2})H$ operator is
applied on the qubit $q[2]$ in case of multi-well potential. Similarly, the
controlled phase gates used here are directly applied using the codes
available in QISKit (See supplementary Information).

Qubits | Freq||(GHz) |  $T^{\dagger}_{1}$ ($\mu$s) |  $T^{\ddagger}_{2}$ ($\mu$s) |  $U_{1}$ err |  $U_{2}$ err $(10^{-3})$ |  $U_{3}$ err $(10^{-3})$ | RE⟂ $(10^{-2})$
---|---|---|---|---|---|---|---
Q0 | 5.25 | 48.23 | 23.49 | 0 | 0.86 | 1.72 | 7.53
Q1 | 5.29 | 53.58 | 10.58 | 0 | 1.55 | 3.09 | 9.15
Q2 | 5.35 | 38.68 | 12.52 | 0 | 1.89 | 3.78 | 3.58
Q3 | 5.43 | 48.83 | 13.94 | 0 | 2.23 | 4.47 | 4.00
Q4 | 5.17 | 46.94 | 7.37 | 0 | 1.80 | 3.61 | 8.78

$||$ Frequency, $\dagger$ Coherence time, $\ddagger$ Relaxation time, $\perp$
Readout Error.

Multi-Qubit Gate | Gate Error | Multi-Qubit Gate | Gate Error
---|---|---|---
$C_{1}X_{0}$ | 0.03488 | $C_{3}X_{2}$ | 0.0866
$C_{2}X_{0}$ | 0.03158 | $C_{3}X_{4}$ | 0.06682
$C_{2}X_{1}$ | 0.05546 | $C_{4}X_{2}$ | 0.0842

Table 1: Experimental parameters of the device “ibmqx4” (information collected
and experimented on Date: 15.02.2019).

Qubits | Freq||(GHz) |  $T^{\dagger}_{1}$ ($\mu$s) |  $T^{\ddagger}_{2}$ ($\mu$s) |  $U_{1}$ err |  $U_{2}$ err $(10^{-3})$ |  $U_{3}$ err $(10^{-3})$ | RE⟂ $(10^{-2})$
---|---|---|---|---|---|---|---
Q0 | 5.10 | 64.45 | 22.04 | 0 | 1.62 | 3.24 | 3.47
Q1 | 5.24 | 56.76 | 81.54 | 0 | 31.89 | 63.78 | 8.83
Q2 | 5.03 | 74.67 | 111.63 | 0 | 2.92 | 5.85 | 11.82
Q3 | 4.89 | 61.06 | 53.99 | 0 | 1.78 | 3.56 | 25.23
Q4 | 5.03 | 56.47 | 40.94 | 0 | 1.90 | 3.79 | 3.41
Q5 | 5.07 | 25.51 | 57.04 | 0 | 2.08 | 4.15 | 7.22
Q6 | 4.92 | 62.04 | 54.84 | 0 | 1.59 | 3.18 | 3.62
Q7 | 4.97 | 45.74 | 84.04 | 0 | 1.26 | 2.52 | 4.38
Q8 | 4.74 | 55.26 | 77.63 | 0 | 1.54 | 3.09 | 7.25
Q9 | 4.96 | 37.90 | 56.65 | 0 | 3.01 | 6.02 | 27.92
Q10 | 4.94 | 55.71 | 69.20 | 0 | 2.22 | 4.44 | 4.21
Q11 | 5.01 | 59.16 | 76.46 | 0 | 1.43 | 2.86 | 10.68
Q12 | 4.76 | 69.46 | 109.99 | 0 | 2.90 | 5.81 | 4.18
Q13 | 4.97 | 25.46 | 39.96 | 0 | 7.74 | 15.48 | 11.45

$||$ Frequency, $\dagger$ Coherence time, $\ddagger$ Relaxation time, $\perp$
Readout Error.

Multi-Qubit Gate | Gate Error | Multi-Qubit Gate | Gate Error | Multi-Qubit Gate | Gate Error
---|---|---|---|---|---
$C_{1}X_{0}$ | 0.03624 | $C_{7}X_{8}$ | 0.02562 | $C_{5}X_{6}$ | 0.06062
$C_{1}X_{2}$ | 0.03289 | $C_{9}X_{8}$ | 0.05133 | $C_{1}2X_{2}$ | 0.06609
$C_{2}X_{3}$ | 0.03919 | $C_{9}X_{1}0$ | 0.06855 | $C_{5}X_{9}$ | 0.05182
$C_{4}X_{3}$ | 0.03584 | $C_{1}1X_{3}$ | 0.06448 | $C_{1}3X_{1}$ | 0.1436
$C_{4}X_{1}0$ | 0.0272 | $C_{1}1X_{1}0$ | 0.04068 | $C_{6}X_{8}$ | 0.03824
$C_{5}X_{4}$ | 0.05023 | $C_{1}1X_{1}2$ | 0.04357 | $C_{1}3X_{1}2$ | 0.04426

Table 2: Specifications of the parameters of each qubit in the “IBM Q 14
Melbourne” quantum computer (information collected and experimented on Date:
13.02.2019)

## References

  * [1] Feynman, R. P.  Simulating physics with computers.  _Int. J. Theor. Phys._ 21, 467–488 (1982).
  * [2] Lloyd, S.  Universal Quantum Simulators.  _Science_ 273, 5278 (1996).
  * [3] Abrams, D. S. & Lloyd, S.  Simulation of many-body Fermi systems on a universal quantum computer.  _Phys. Rev. Lett._ 79, 2586–2589 (1997).
  * [4] Zalka, C.  Simulating quantum systems on a quantum computer.  _Proc. R. Soc. Lond. A_ 454, 313–322 (1998).
  * [5] Berry, D. _et al._ Efficient quantum algorithms for simulating sparse Hamiltonians.  _Commun. Math. Phys._ 270, 359–371 (2007).
  * [6] Deutsch, D. & Jozsa, R.  Rapid solutions of problems by quantum computation.  _Proc. R. Soc. Lond. A_ 439, 553 (1992).
  * [7] Simon, D. R.  On the power of quantum computation.  _SIAM J. Comput._ 26, 1474–1483 (1997).
  * [8] Gerjuoy, E.  Shor’s factoring algorithm and modern cryptography. An illustration of the capabilities inherent in quantum computers.  _Am. J. Phys._ 73, 523–540 (2005).
  * [9] Grover, L. K.  Quantum Mechanics helps in searching for a needle in a haystack.  _Phys. Rev. Lett._ 79, 325 (1997).
  * [10] Tseng, C. H. _et al._ Quantum simulation of a three-body-interaction Hamiltonian on an NMR quantum computer.  _Phys. Rev. A_ 61, 012302 (1999).
  * [11] Negrevergne, C. _et al._ Liquid-state NMR simulations of quantum many-body problems.  _Phys. Rev. A_ 71, 032344 (2005).
  * [12] Peng, X. H., Zhang, J. F., Du, J. F. & Suter, D.  Quantum simulation of a system with competing two- and three-body interactions.  _Phys. Rev. Lett._ 103, 140501 (2009).
  * [13] Peng, G. R., Xu, G. F. & Long, G. L.  Experimental realization of nonadiabatic holonomic quantum computation.  _Phys. Rev. Lett._ 110, 190501 (2013).
  * [14] Zhang, J. F., Wei, T. C. & Laflamme, R.  Experimental quantum simulation of entanglement in many-body systems.  _Phys. Rev. Lett._ 107, 010501 (2011).
  * [15] Peng, X. H., Du, J. F. & Suter, D.  Quantum phase transition of ground-state entanglement in a Heisenberg spin chain simulated in an NMR quantum computer.  _Phys. Rev. A_ 71, 012307 (2005).
  * [16] Edwards, E. E. _et al._ Quantum simulation and phase diagram of the transverse-field Ising model with three atomic spins.  _Phys. Rev. B_ 82, 060412 (2010).
  * [17] Zhang, J. F. _et al._ Direct observation of quantum criticality in Ising spin chains.  _Phys. Rev. A_ 79, 012305 (2009).
  * [18] Aspuru-Guzik, A., Dutoi, A. D., Love, P. J. & Head-Gordon, M.  Simulated quantum computation of molecular energies.  _Science_ 309, 5741 (2005).
  * [19] Lanyon, B. P. _et al._ Towards quantum chemistry on a quantum computer.  _Nat. Chem._ 2, 106 (2010).
  * [20] Du, J. F. _et al._ NMR implementation of a molecular hydrogen quantum simulation with adiabatic state preparation.  _Phys. Rev. Lett._ 104, 030502 (2010).
  * [21] Kassal, I., Whitfield, J. D., Perdomo-Ortiz, A., Yung, M. H. & Aspuru-Guzik, A.  Simulating Chemistry using Quantum Computers.  _Annu. Rev. Phys. Chem._ 62, 185-207 (2011).
  * [22] Yung, M. H., Whitfield, J. D., Boixo, S., Tempel, D. G. & Aspuru-Guzik, A.  Introduction to quantum algorithms for physics and chemistry.  _arXiv preprint arXiv:1203.1331_ (2012).
  * [23] Büchler, H. P. _et al._ Atomic quantum simulator for lattice gauge theories and ring exchange models.  _Phys. Rev. Lett._ 95, 040402 (2005).
  * [24] Zohar, E. & Reznik, B.  Confinement and lattice quantum-electrodynamic electric flux tubes simulated with ultracold atoms.  _Phys. Rev. Lett._ 107, 275301 (2011).
  * [25] Cirac, J. I. _et al._ Cold atom simulation of interacting relativistic quantum field theories.  _Phys. Rev. Lett._ 105, 190403 (2010).
  * [26] Cirac, L. _et al._ An optical-lattice-based quantum simulator for relativistic field theories and topological insulators.  _New J. Phys._ 14, 015007 (2012).
  * [27] Kapit, E. & Mueller, E.  Optical-lattice Hamiltonians for relativistic quantum electrodynamics.  _Phys. Rev. A_ 83, 033625 (2011).
  * [28] Bermudez, A. _et al._ Wilson fermions and axion electrodynamics in optical lattices.  _Phys. Rev. Lett._ 105, 190404 (2010).
  * [29] Lepori, L. _et al._ (3+1) massive Dirac fermions with ultracold atoms in frustrated cubic optical lattices.  _Europhys. Lett._ 92, 50003 (2010).
  * [30] Maeda, K. _et al._ Simulating dense QCD matter with ultracold atomic boson-fermion mixtures.  _Phys. Rev. Lett._ 103, 085301 (2009).
  * [31] Casanova, J. _et al._ Quantum simulation of quantum field theories in trapped ions.  _Phys. Rev. Lett._ 107, 260501 (2011).
  * [32] Doucot, B. _et al._ Discrete non-Abelian gauge theories in Josephson junction arrays and quantum computation.  _Phys. Rev. B_ 69, 214501 (2004).
  * [33] Johanning, M. _et al._ Quantum simulations with cold trapped ions.  _J. Phys. B_ 42, 154009 (2009).
  * [34] Bender, J. _et al._ Digital quantum simulation of lattice gauge theories in three spatial dimensions.  _New J. Phys._ 20, 093001 (2018).
  * [35] Jordan, S. P. _et al._ Quantum algorithms for quantum field theories.  _Science_ 336, 1130–1133 (2012).
  * [36] Byrnes, T. & Yamamoto, Y.  Simulating lattice gauge theories on a quantum computer.  _Phys. Rev. A_ 73, 022328 (2006).
  * [37] Somaroo, S. _et al._ Quantum simulations on a quantum computer.  _Phys. Rev. Lett._ 82, 5381–5383 (1999).
  * [38] Khitrin, A. K. & Fung, B. M.  NMR simulation of an eight-state quantum system.  _Phys. Rev. A_ 64, 032306 (2001).
  * [39] Brown, K. R. _et al._ Limitations of quantum simulation examined by a pairing Hamiltonian using nuclear magnetic resonance.  _Phys. Rev. Lett._ 97, 050504 (2006).
  * [40] Friedenauer, A. _et al._ Simulating a quantum magnet with trapped ions.  _Nat. Phys._ 4, 757–761 (2008).
  * [41] Gerritsma, R. _et al._ Quantum simulation of the Dirac equation.  _Nature_ 463, 68–71 (2010).
  * [42] Gerritsma, R. _et al._ Quantum simulation of the Klein paradox with trapped ions.  _Phys. Rev. Lett._ 106, 060503 (2011).
  * [43] Lanyon, B. P. _et al._ Universal digital quantum simulation with trapped ions.  _Science_ 334, 57–61 (2011).
  * [44] Lanyon, B. P. _et al._ Towards quantum chemistry on a quantum computer.  _Nat. Chem._ 2, 106–111 (2010).
  * [45] Kinoshita, T. _et al._ Observation of a one-dimensional Tonks-Girardeau gas.  _Phys. Rev. B_ 82, 060412 (2010).
  * [46] Ma, X. S. _et al._ Quantum simulation of the wavefunction to probe frustrated Heisenberg spin systems.  _Nat. Phys._ 7, 399–405 (2011).
  * [47] Kassal, I. _et al._ Simulating chemistry using quantum computers.  _Annu. Rev. Phys. Chem._ 62, 185–207 (2011).
  * [48] Childs, A. M. _et al._ Toward the first quantum simulation with quantum speedup  _PNAS_ 115, 9456-9461 (2018).
  * [49] Boixo, S. _et al._ Characterizing quantum supremacy in near-term devices.  _Nat. Phys._ 14, 595-600 (2018).
  * [50] Preskill, J.  Quantum Computing in the NISQ era and beyond.  _Quantum_ 2, 79 (2018).
  * [51] Nguyen, T. L. _et al._ Towards Quantum Simulation with Circular Rydberg Atoms.  _Phys. Rev. X_ 8, 011032 (2018).
  * [52] Bermejo-Vega, J. _et al._ Architectures for Quantum Simulation Showing a Quantum Speedup.  _Phys. Rev. X_ 8, 021010 (2018).
  * [53] Kivlichan, I. D. _et al._ Quantum Simulation of Electronic Structure with Linear Depth and Connectivity.  _Phys. Rev. Lett._ 120, 110501 (2018).
  * [54] Hempel, C. _et al._ Quantum Chemistry Calculations on a Trapped-Ion Quantum Simulator.  _Phys. Rev. X_ 8, 031022 (2018).
  * [55] Noh, C. & Angelakis, D. G.  Quantum simulations and many-body physics with light.  _Rep. Prog. Phys._ 80, 016401 (2016).
  * [56] Georgescu, I. M. _et al._ Quantum simulation.  _Rev. Mod. Phys._ 86, 153 (2014).
  * [57] Cirac, J. I. & Zoller, P.  Goals and opportunities in quantum simulation.  _Nat. Phys._ 8, 264–266 (2012).
  * [58] Josephson, B. D.  Possible new effects in superconductive tunnelling.  _Phys. Lett._ 1, 251 (1962).
  * [59] Esaki, L. & Miyahara, Y.  A new device using the tunneling process in narrow p-n junctions.  _Solid-State Electron._ 1, 13 (1960).
  * [60] Binnig, G. & Rohrer, H.  Scanning tunneling microscopy.  _IBM J. Res. Dev._ 44, 279–293 (2000).
  * [61] Kogut, J. B.  The lattice gauge-theory approach to quantum chromodynamics  _IBM J. Res. Dev._ 55, 775–836 (1983).
  * [62] Smith, K. & Blaylock, G.  Simulations in quantum tunneling  _Am. J. Phys._ 85, 763 (2017).
  * [63] Hartman, T. E.  Tunneling of a wave packet  _J. Appl. Phys._ 33, 3427–3433 (1962).
  * [64] Steinberg, A. M., Kwiat, P. G. & Chiao, R. Y.  Measurement of the single-photon tunneling time  _Phys. Rev. Lett._ 71, 708–711 (1993).
  * [65] Sornborger, A. T.  Quantum simulation of tunneling in small systems.  _Sci. Rep._ 2, 597 (2012).
  * [66] Feng, G.-R., Lu, Y., Hao, L., Zhang, F.-H. & Long, G.-L.  Experimental simulation of quantum tunneling in small systems.  _Sci. Rep._ 3, 2232 (2013).
  * [67] Ostrowski, M.  Quantum simulation of the tunnel effect  _Bull. Pol. Acad. Tech. Sci._ 63, 42 (2015).
  * [68] García-Martín, D., & Sierra, G.  Five Experimental Tests on the 5-Qubit IBM Quantum Computer.  _J. App. Math. Phys._ 6, 1460 (2018).
  * [69] Behera, B. K., Banerjee, A. & Panigrahi, P. K.  Experimental realization of quantum cheque using a five-qubit quantum computer.  _Quantum Inf. Process._ 16, 312 (2016).
  * [70] Majumder, A., Mohapatra, S. & Kumar, A.  Experimental Realization of Secure Multiparty Quantum Summation Using Five-Qubit IBM Quantum Computer on Cloud.  _arXiv preprint arXiv:1707.07460_ (2017).
  * [71] Sisodia, M., Shukla, A., Thapliyal, K. & Pathak, A.  Design and experimental realization of an optimal scheme for teleportion of an n-qubit quantum state.  _Quantum Inf. Process._ 16, 292 (2017).
  * [72] Dash, A., Rout, S., Behera, B. K. & Panigrahi, P. K.  A Verification Algorithm and Its Application to Quantum Locker in IBM Quantum Computer.  _arXiv preprint arXiv:1710.05196_ (2017).
  * [73] Deffner, S.  Demonstration of entanglement assisted invariance on IBM’s quantum experience.  _Heliyon_ 3, e00444 (2016).
  * [74] Huffman, E. & Mizel, A.  Violation of noninvasive macrorealism by a superconducting qubit: Implementation of a Leggett-Garg test that addresses the clumsiness loophole.  _Phys. Rev. A_ 95, 032131 (2017).
  * [75] Vishnu, P. K., Joy, D., Behera, B. K. & Panigrahi, P. K.  Experimental demonstration of non-local controlled-unitary quantum gates using a five-qubit quantum computer.  _Quantum Inf. Process._ 17, 274 (2018).
  * [76] Alsina, D. & Latorre, J. I.  Experimental test of Mermin inequalities on a five-qubit quantum computer.  _Phys. Rev. A_ 94, 012314 (2016).
  * [77] Yalçinkaya, İ. & Gedik, Z.  Optimization and experimental realization of quantum permutation algorithm.  _Phys. Rev. A_ 96, 062339 (2017).
  * [78] Ghosh, D., Agarwal, P., Pandey, P., Behera, B. K. & Panigrahi, P. K.  Automated Error Correction in IBM Quantum Computer and Explicit Generalization.  _Quantum Inf. Process._ 17, 153 (2018).
  * [79] Kandala, A. _et al._ Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  _Nature_ 549, 242–246 (2017).
  * [80] Satyajit, S., Srinivasan, K., Behera, B. K. & Panigrahi, P. K.  Nondestructive discrimination of a new family of highly entangled states in IBM quantum computer.  _Quantum Inf. Process._ 17, 212 (2018).
  * [81] Alvarez-Rodriguez, U., Sanz, M., Lamata, L. & Solano, E.  Quantum Artificial Life in an IBM Quantum Computer.  _arXiv preprint arXiv:1711.09442_ (2017).
  * [82] Schuld, M., Fingerhuth, M. & Petruccione, F.  Implementing a distance-based classifier with a quantum interference circuit.  _Europhys. Lett._ 119, 60002 (2017).
  * [83] Sisodia, M. _et al._ Experimental realization of nondestructive discrimination of Bell states using a five-qubit quantum computer.  _Phys. Lett. A_ 381, 3860–3874 (2017).
  * [84] Tannu, S. S. & Qureshi, M. K.  A Case for Variability-Aware Policies for NISQ-Era Quantum Computers.  _arXiv preprint arXiv:1805.10224_ (2018).
  * [85] Wootton, J. R.  Benchmarking of quantum processors with random circuits.  _arXiv preprint arXiv:1806.02736_ (2018).
  * [86] Harper, R. & Flammia, S.  Fault tolerance in the IBM Q Experience.  _arXiv preprint arXiv:1806.02359_ (2018).
  * [87] Aggarwal, D., Raj, S., Behera, B. K. & Panigrahi, P. K.  Application of quantum scrambling in Rydberg atom on IBM quantum computer.  _arXiv preprint arXiv:1806.00781_ (2018).
  * [88] Srinivasan, K., Satyajit, S., Behera, B. K. & Panigrahi, P. K.  Efficient quantum algorithm for solving travelling salesman problem: An IBM quantum experience.  _arXiv preprint arXiv:1805.10928_ (2018).
  * [89] Dash, A., Sarmah, D., Behera, B. K. & Panigrahi, P. K.  Exact search algorithm to factorize large biprimes and a triprime on IBM quantum computer.  _arXiv preprint arXiv:1805.10478_ (2018).
  * [90] Roffe, J., Headley, D., Chancellor, N., Horsman, D. & Kendon, V.  Protecting quantum memories using coherent parity check codes.  _Quantum Sci. Technol._ 3, 3 (2018).
  * [91] Plesa, M.-I., & Mihai, T.  A New Quantum Encryption Scheme.  _Adv. J. Grad. Res._ 3, 1 (2018).
  * [92] Manabputra, Behera, B. K. & Panigrahi, P. K.  A Simulational Model for Witnessing Quantum Effects of Gravity Using IBM Quantum Computer.  _arXiv preprint arXiv:1806.10229_ (2018).
  * [93] Jha, R., Das, D., Dash, A., Jayaraman, S., Behera, B. K. & Panigrahi, P. K.  A Novel Quantum N-Queens Solver Algorithm and its Simulation and Application to Satellite Communication Using IBM Quantum Experience.  _arXiv preprint arXiv:1806.10221_ (2018).
  * [94] Gangopadhyay, S., Manabputra, Behera, B. K. & Panigrahi, P. K.  Generalization and Demonstration of an Entanglement Based Deutsch-Jozsa Like Algorithm Using a 5-Qubit Quantum Computer.  _Quantum Inf. Process._ 17, 160 (2018).
  * [95] Buluta, I. & Nori, F.  Quantum Simulators.  _Science_ 326, 108–111 (2009).
  * [96] Mezzacapo, A. _et al._ Digital Quantum Simulation of the Holstein Model in Trapped Ions.  _Phys.Rev.Lett._ 109, 200501 (2012).
  * [97] Heras, U. L. _et al._ Digital Quantum Simulation of Spin Systems in Superconducting Circuits.  _Phys. Rev. Lett._ 112, 200501 (2014).
  * [98] Barends, R. _et al._ Digital quantum simulation of fermionic models with a superconducting circuit.  _Nat. Commun._ 6, 7654 (2015).
  * [99] Garcia-Alvarez, L. _et al._ Digital Quantum Simulation of Minimal AdS/CFT.  _Phys. Rev. Lett._ 119, 040501 (2017).
  * [100] Lanyon, B. P. _et al._ Universal Digital Quantum Simulation with Trapped Ions.  _Science_ 332, 57–61 (2011).
  * [101] Langford, N. K. _et al._ Experimentally simulating the dynamics of quantum light and matter at deep-strong coupling.  _npj Quantum Information_ 8, 1715 (2017).
  * [102] Benenti, G. & Strini, G.  Quantum simulation of the single-particle Schrödinger equation.  _Am. J. Phys._ 76, 657 (2008).
  * [103] Yoshida, H.  Construction of higher-order symplectic integrators.  _Phys. Lett. A_ 150, 262–268 (1990).
  * [104] Sornborger, A. T. & Stewart, E. D.  Higher-order methods for simulations on quantum computers.  _Phys. Rev. A_ 60, 1956–1965 (1999).
  * [105] Hatano, N. & Suzuki, M.  Quantum Annealing and Other Optimization Methods.  _Springer, Heidelberg (2005)_.
  * [106] Coppersmith, D.  An approximate Fourier transform useful in quantum factoring.  _IBM Res. Rep._ (1994).

##  5 Data Availability statement

The data collected and the codes used in the experiment are available in the
Supplementary Information. Request for any materials should be addressed to
P.K.P. (pprasanta@iiserkol.ac.in).

## Acknowledgements

N.N.H. acknowledges the hospitality of Indian Institute of Science Education
and Research Kolkata during the project work. B.K.B. acknowledges the
financial support of IISER-K Institute fellowship. We acknowledge the support
of IBM Quantum Experience for using the quantum processors. The views
expressed are those of the authors and do not reflect the official position of
IBM or the IBM quantum experience team.

## Author contributions

Theoretical analysis and design of quantum circuits were done by N.N.H and
B.K.B. Simulation of quantum circuit, collection and analysis of data were
done by N.N.H. The project was supervised by B.K.B. P.K.P. has thoroughly
reviewed the manuscript. N.N.H. and B.K.B. have completed the project under
the guidance of P.K.P.

## Competing interests

The authors declare no competing financial as well as non-financial interests.
Readers are welcome to comment on the online version of the paper.
Correspondence and requests for materials should be addressed to P.K.P.
(pprasanta@iiserkol.ac.in).

## Supplementary Information: Experimental Demonstration of Quantum Tunneling
in IBM Quantum Computer

For simulating quantum tunneling, we used QISKit to take both simulation and
experimental results. The python code for the same is as follows.

⬇

1#Quantum simulation of tunneling in double-well potential

# Import the QISKit SDK

3import numpy as np

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

5from qiskit import Aer, execute

from math import pi

7from qiskit import IBMQ

from qiskit.tools.monitor import job_monitor

9from qiskit.tools.visualization import plot_histogram

11# Set your API Token.

# IBMQ.enable_account(’API Token’)

13

# Create a Quantum Register with 2 qubits.

15q = QuantumRegister(2)

# Create a Classical Register with 2 bits.

17c = ClassicalRegister(2)

# Create a Quantum Circuit

19qc = QuantumCircuit(q, c)

21# Circuit Implementation

# Initial state

23qc.x(q[1])

# Time step=1

25# Time evolution of potential energy part

qc.h(q[1])

27qc.u3(5, -pi/2, pi/2, q[1])

qc.h(q[1])

29# 2-qubit Inverse Quantum Fourier Transform

qc.h(q[0])

31qc.cu1(pi/2, q[1], q[0])

qc.h(q[1])

33# Time evolution of kinetic energy part

qc.h(q[0])

35qc.u3(-(pi**2)/40, -pi/2, pi/2, q[0])

qc.h(q[0])

37qc.h(q[1])

qc.u3(-(pi**2)/10, -pi/2, pi/2, q[1])

39qc.h(q[1])

qc.cu1((pi**2)/10, q[1], q[0])

41# 2-qubit Quantum Fourier Transform

qc.h(q[1])

43qc.cu1(-pi/2, q[1], q[0])

qc.h(q[0])

45# Time evolution of potential energy part

qc.h(q[1])

47qc.u3(5, -pi/2, pi/2, q[1])

qc.h(q[1])

49

# Measurement

51qc.measure(q, c)

# Choose backend

53backend = IBMQ.get_backend(’ibmq_qasm_simulator’)

# backend = IBMQ.get_backend(’ibmqx4’)

55# backend = IBMQ.get_backend(’ibmqx2’)

# backend = IBMQ.get_backend(’ibmq_16_melbourne’)

57shots = 8192 # Number of shots to run the program (experiment); maximum is
8192 shots.

max_credits = 10 # Maximum number of credits to spend on executions.

59job_exp = execute(qc, backend=backend, shots=shots, max_credits=max_credits)

job_monitor(job_exp)

61result = job_exp.result()

counts_exp = result.get_counts(qc)

63print(result.get_counts(qc))

plot_histogram(counts_exp)

doublewell_tunneling.py Figure 13: Free particle probability distribution: we
have used the qubits q[0] and q[1] from IBM Q 5 Tenerife (ibmqx4) for our
experiment, No. of shots=8192. Figure 14: Probability distribution for
particle in step potential: we have used the qubits q[0], q[1] from ibmqx4 and
q[7], q[8] from IBM Q 14 Melbourne for our experiment, No. of shots=8192.
Figure 15: Probability distribution for particle in double-well potential: we
have used the qubits q[7] and q[8] from IBM Q 14 Melbourne for our experiment,
No. of shots=8192.


# Josephson Junction through a Disordered Topological Insulator with Helical
Magnetization

Alexander Zyuzin  Mohammad Alidoust  Daniel Loss  Department of Physics,
University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland

###### Abstract

We study supercurrent and proximity vortices in a Josephson junction made of
disordered surface states of a three-dimensional topological insulator with a
proximity induced in-plane helical magnetization. In a regime where the
rotation period of helical magnetization is larger than the junction width, we
find supercurrent $0$-$\pi$ crossovers as a function of junction thickness,
magnetization strength, and parameters inherent to the helical modulation and
surface states. The supercurrent reversals are associated with proximity
induced vortices, nucleated along the junction width, where the number of
vortices and their locations can be manipulated by means of the
superconducting phase difference and the parameters mentioned above.

###### pacs:

74.50.+r, 73.20.-r, 73.63.-b,

##  I Introduction

Topological insulators (TIs) form a new state of matter which has offered
novel prospects towards topological superconducting spintronics and
topological quantum computation pankratov ; kane_rmp ; zhang_rmp . The TIs
provide a platform to observe quantum relativistic phenomena, stemming from
strong spin orbit interaction, such as the quantum spin Hall effect, spin-
momentum locking, and manipulation of Dirac fermions (see reviews kane_rmp ;
zhang_rmp ). The surface of three-dimensional TI in the presence of time
reversal symmetry hosts metallic helical states, i.e., for each momentum on
the Fermi surface of the surface states, the spin has a rigidly defined
direction, transverse to the momentum. Interestingly, the interplay of
superconductivity and magnetism at the surface of a TI may give rise to
topological superconductivity and Majorana fermions which has been receiving
strong interest both experimentally and theoretically Volovik ; exp_ti2a ;
exp_ti2b ; exp_ti2c ; exp_ti2c2 ; exp_ti4a ; exp_ti4b ; exp_ti4c .

To the best of our knowledge, the theoretical studies on superconducting TI
structures have considered fully ballistic surface states. However,
experimentally realistic systems inevitably contain impurity scattering
processes that can play a key role in the actual quantum transport through the
surface states exp_ti4a ; burkov1 ; Schwab ; exp_ti2b ; Kurter_TI_exp . A
theory for the disordered limit of ‘superconducting’ TI with the possibility
of inclusion of ‘magnetism’ is still lacking. Therefore, with the rapid
experimental progress in TI heterostructures such a theoretical framework is
becoming more essential for describing and predicting important physical
phenomena in this field. The quasiclassical formalism is capable of providing
such a theoretical technique brgt_rmp1 ; brgt_rmp2 ; brgt_rmp3 .

Figure 1:  (Color online). Schematic of a superconductor (S) - topological
insulator (TI) - superconductor junction with a helical magnetization pattern.
The magnetization vector follows a helical pattern given by
$\mathbf{h}(\mathbf{r})=h_{0}(\cos Qy,\sin Qy,0)$. The junction plane resides
in the $xy$ plane so that the S-TI interfaces lie in the $y$ direction at
$x=0,d_{F}$. The junction has a length and width of $d_{F}$ and $W_{F}$,
respectively. The superconducting electrodes are connected to the diffusive
surface states of the TI through tunneling barriers.

It is well understood that the ground state of a uniformly magnetized 3D
metallic superconductor - ferromagnet - superconductor (SFS) junction can be
reverted from $0$ to $\pi$ by varying the magnetization strength, F layer
thickness, and the system temperature, for a review see Refs. brgt_rmp1 ;
brgt_rmp2 ; brgt_rmp3 . Such a $0$-$\pi$ transition, first predicted
theoretically buz_0_pi_0a ; buz_0_pi_0b ; buz_0_pi_0c , is now well
established experimentally Ryazanov_0a ; Ryazanov_0b ; Ryazanov_0c . The
$0$-$\pi$ crossover is important for applications in quantum computations and
for the development of ultrafast switches in superconducting spintronics and
functional nanodevices brgt_rmp1 ; brgt_rmp2 ; brgt_rmp3 . On the other hand,
it is known that the free energy of the SFS Josephson junction in the presence
of the spin-orbit interaction in the F region has a minimum at phase
difference $\phi=\phi_{0},~{}(\neq 0,\pi)$, phi_jncta ; phi_jnctb ; phi_jnctc
; phi_jnctd ; phi_jncte  (the so-called $\phi_{0}$-junctions). Signatures of
the Josephson $\phi_{0}$-junction based on a nanowire quantum dot were
recently observed experimentally in Ref. phi0exp . We also note that a
$\phi_{0}$ junction may be also achieved in a different situation where the
superconducting phase is random along the junction width zyuzin .

It was shown that the critical supercurrent shows no reversal in Josephson
junctions made of ballistic surface states of a TI (S-TI-S) with proximity
induced ‘uniform’ magnetization tanaka_sfs . This follows from the fact that
the TI surface states do not respond to a static and uniform in-plane magnetic
(or exchange) field due to the momentum-spin locking of Dirac fermions. In-
plane magnetic fields via Zeeman effect simply shift the surface Dirac cones
in momentum space. Thus, the current-phase relation has a phase shift of
$\phi_{0}$, proportional to the junction thickness and the Zeeman energy.

Here we generalize a quasiclassical approach for disordered surface states of
three-dimensional TIs in the presence of $s$-wave superconductivity and in-
plane magnetization with an arbitrary pattern. We derive the Eilenberger
equation eiln  which describes the system from a fully ballistic to a weakly
disordered regime. We also derive the Usadel equation usadel , governing the
transport of quasiparticles in fully disordered TI surface states where
diffusive motion dominates the ballistic one.

We first employ the Usadel approach for a S-TI-S junction with uniform
magnetization and show that similarly to the ballistic limit the critical
supercurrent is featureless and an exchange field parallel to the S-TI
interfaces transfers into the superconducting phase difference across the
junction. We then study the supercurrent in a diffusive 2D S-TI-S with an in-
plane helical magnetization, illustrated in Fig. 1. Interestingly, we find
that multiple supercurrent reversals and proximity vortices can appear by
simply manipulating the helical magnetization parameters and varying the
junction thickness.

The rest of the article is organized as follows. In the next section (Sec. II)
we describe our model and discuss main results. In Sec. III we present the
detailed analytical derivations of the Eilenberger and Usadel equations as
well as the supercurrent through TI, and finally in Sec. IV we present
concluding remarks.

##  II  Main Results

The equation for the Green function of quasiparticle at the surface of 3D TI
reads:

$\displaystyle\left(\begin{array}[]{cc}-i\omega_{n}+\hat{H}(\mathbf{r})&0\\\
0&i\omega_{n}+\hat{\sigma}^{y}\hat{H}^{*}(\mathbf{r})\hat{\sigma}^{y}\end{array}\right)\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})$
(3)
$\displaystyle=\delta(\mathbf{r}-\mathbf{r}^{\prime})+\frac{1}{\pi\nu\tau}\check{G}(\omega_{n};\mathbf{r},\mathbf{r})\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime}),$
(4)

where the Hamiltonian
$\hat{H}(\mathbf{r})=-i\alpha({\bm{\nabla}}\times\mathbf{e}_{z})\cdot\hat{\bm{\sigma}}+\mathbf{h}(\mathbf{r})\cdot\hat{{\bm{\sigma}}}-\mu$
describes the surface states of the TI with a proximity induced in-plane
exchange field
$\mathbf{h}(\mathbf{r})=(h_{x}(\mathbf{r}),h_{y}(\mathbf{r}),0)$. Here,
$\alpha$ is the Fermi velocity, $\mu$ is the chemical potential,
$\mathbf{e}_{z}$ is an unit vector normal to the surface of TI, and
$\hat{\bm{\sigma}}$ is a vector comprised of the Pauli matrices acting on the
spin degree of freedom. We assume $\hbar=k_{B}=1$ throughout the paper. The
Green function is averaged over a nonmagnetic scattering potential
$V(\mathbf{r})$, which is assumed Gaussian $\langle
V(\mathbf{r})V(\mathbf{r}^{\prime})\rangle=\delta(\mathbf{r}-\mathbf{r}^{\prime})/\pi\nu\tau$,
where $\tau$ is the mean free time of particles in the disordered system and
$\nu=\mu/2\pi\alpha^{2}$ is the density of states per spin at the Fermi level
of the normal state of the TI. The matrix structure of the Green function in
the rotated particle-hole and spin basis has the following form:

$\displaystyle\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=\left(\begin{array}[]{cc}-\hat{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&-\hat{F}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})i\hat{\sigma}^{y}\\\
-i\hat{\sigma}^{y}\hat{F}^{{\dagger}}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&\hat{\sigma}^{y}\hat{\bar{G}}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\hat{\sigma}^{y}\end{array}\right),$

where $\omega_{n}=\pi T(2n+1)$ is the Matsubara frequency, $T$ is the
temperature, and $n\in\mathbb{Z}$. The check symbol $\check{1}$ represents
$4\times 4$ matrices in the particle-hole and spin spaces, while the hat
symbol $\hat{1}$ defines $2\times 2$ matrices. Notice that the off-diagonal
components of the Green function describe the penetration of Cooper pairs into
the surface states.

In order to solve Eq. (3) we employ the quasiclassical approximation which
results in the Eilenberger equation eiln  for the surface channels:

$\displaystyle\frac{\alpha}{2}\Big{\\{}\hat{{\bm{\eta}}},{\bm{\nabla}}\check{g}\Big{\\}}=\bigg{[}\check{g},\omega_{n}\hat{\tau}^{z}+i\mathbf{h}\cdot\hat{{\bm{\sigma}}}\hat{\tau}^{z}+i\mu\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F}+\frac{\langle\check{g}\rangle}{\tau}\bigg{]},~{}~{}~{}~{}$
(7)

where we have performed a Fourier transformation of the Green function Eq.
(II) with respect to the relative space arguments and then taken the integral
over $\xi_{p}=\alpha p-\mu$ which results in
$\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{F})=\int\frac{d\xi_{p}}{\pi
i}\check{G}(\omega_{n};\mathbf{R},\mathbf{p})$. Here,
$\mathbf{n}_{F}=\mathbf{p}_{F}/|\mathbf{p}_{F}|$ is an unit vector in the
direction of momentum $\mathbf{p}_{F}$ at Fermi level $\mu=\alpha p_{F}$,
$\hat{{\bm{\eta}}}\equiv(-\sigma^{y},\sigma^{x})$ is the vector of Pauli
matrices and $\hat{\tau}^{z}$ is a Pauli matrix acting in the particle-hole
space. The disorder potential scatter quasiparticles in random directions in
the momentum space. Therefore, one can integrate the quasiclassical Green
function over all possible directions of quasiparticle momentum
$\langle\check{g}\rangle$. To find a solution to Eq. (7), we expand the Green
function through Pauli matrices
$\check{g}=(\hat{g}^{\prime}\hat{1}+\hat{\mathbf{g}}^{\prime}\cdot\hat{{\bm{\eta}}}+\hat{g}^{\prime}_{z}\hat{\sigma}^{z})/2$,
where
$\hat{\mathbf{g}}^{\prime}=(\hat{g}^{\prime}_{x},\hat{g}^{\prime}_{y},0)$. The
spin structure of $\check{g}$ in the limit of $|\mathbf{h}|\ll\mu$ is defined
by the conduction band projector $(1+\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F})/2$.
We find that the main contributing components of $\check{g}$ commute with
$\hat{\eta}\mathbf{n}_{F}$. Hence, the $\hat{g}^{\prime}_{z}$ component is
smaller than $\hat{g}^{\prime}$ and $\hat{g}^{\prime}_{x,y}$ by a factor of
$\sim\mathrm{max}(\frac{1}{\tau},|\mathbf{h}|)/\mu\ll 1$ and can be neglected.

Figure 2:  (Color online). Normalized critical supercurrent through diffusive
TI surface states with an induced helical magnetization. Panel ($a$) shows the
critical current as a function of rotation parameter $q$ at three different
values of junction length: $d_{F}=2.0\xi_{S},2.5\xi_{S},3.5\xi_{S}$. Panel
($b$) shows the critical current as a function of junction length $d_{F}$ for
$q=0.1\pi,0.4\pi$, and $0.7\pi$. The junction ground state oscillates between
$0$ and $\pi$ superconducting phase difference by varying $q,d_{F},h_{0}$.

In the diffusive limit the Green function can be expanded up to the first two
terms of 2D harmonics, namely
$\hat{g}^{\prime}(\omega_{n};\mathbf{R},\mathbf{n}_{F})=\hat{g}_{s}(\omega_{n};\mathbf{R})+\mathbf{n}_{F}\cdot\hat{\mathbf{g}}_{a}(\omega_{n};\mathbf{R})$,
where the zero harmonic is isotropic and its amplitude is larger than the
first harmonic. We substitute this expansion into Eq. (7) and perform an
integration over momentum directions. By taking a spin trace we arrive at
$\hat{\mathbf{g}}_{a}=-2\alpha\tau\hat{g}_{s}\hat{{\bm{\nabla}}}\hat{g}_{s}$,
which results in the Usadel equation usadel :

$\displaystyle
D\hat{{\bm{\nabla}}}(\hat{g}_{s}\hat{{\bm{\nabla}}}\hat{g}_{s})=[\omega_{n}\hat{\tau}^{z},\hat{g}_{s}].$
(8)

Here, $D=\alpha^{2}\tau$ denotes the diffusion coefficient in the TI. We have
defined a covariant derivative $\hat{{\bm{\nabla}}}$ so that
$\hat{{\bm{\nabla}}}X={\bm{\nabla}}X+\frac{i}{\alpha}(h_{x}\mathbf{e}_{y}-h_{y}\mathbf{e}_{x})[\hat{\tau}^{z},X]$,
in which $\mathbf{e}_{x,y}$ are unit vectors in the $x,y$ directions. Contrary
to the zero harmonics of the Green function expansion, the first harmonic
satisfies a non-diffusive equation due to the fast spin relaxation time
proportional to the mean free time at the surface. This finding is consistent
with that of the electron spin and charge densities at the TI surface in the
absence of superconductivity burkov1 ; Schwab . We note that the Usadel
equation for an electron gas with spin orbit coupling in the presence of a
Zeeman field was derived in Refs. huzet_bergeret, . Here we however focus on
the disordered surface states of a topological insulator which are described
by the Dirac Hamiltonian. We utilize the derived Usadel equation to study the
supercurrent and proximity induced vortices’ profile in a 2D Josephson
junction with helical magnetization.

Using the definition of quasiclassical Green functions, we obtain an
expression for the charge current density in the diffusive limit:

$\mathbf{J}=-\frac{\pi
i\sigma_{N}}{4e}T\sum_{n}\mathrm{Tr}[\hat{\tau}^{z}\hat{g}_{s}\hat{{\bm{\nabla}}}\hat{g}_{s}],$
(9)

where $\sigma_{N}=2e^{2}\nu D$ is the normal-state conductivity. To calculate
the charge supercurrent flow across the junction shown in Fig. 1, one needs to
find a solution to Eq. (8) in the TI region and match the Green function at
S-TI boundaries. Here, we assume that the TI surface is sandwiched between two
spin-orbit free $s$-wave superconducting electrodes. We will consider the
limit of low transparency of the S-TI interfaces so that the spin-singlet
Cooper pairs tunnel from the superconducting leads into the surface states and
thus employ the Kupriyanov-Lukichev boundary conditions for the Green function
at the contacts boundary_c1 ; boundary_c2 . The leakage intensity of the
superconducting correlations is controlled by the parameter $\gamma\gg 1$,
which is the ratio of resistance per unit area of the surface of the tunneling
barrier to the resistivity of the diffusive TI surface states. To derive the
boundary conditions for the S-TI-S junction considered in this paper, one
should follow Refs. boundary_c1 ; boundary_c2  and consider the effect of
momentum-spin locking. One therefore straightforwardly arrives at the
following expression for the boundary conditions:
$2\gamma\hat{g}_{s}\mathbf{n}\cdot\hat{{\bm{\nabla}}}\hat{g}_{s}=[\hat{g}_{s},\hat{g}_{\text{SC}}]$,
where the unit vector $\mathbf{n}$ points normal to the boundary and the Green
function in the superconducting lead $\hat{g}_{\text{SC}}(\omega_{n})$ is
given by its bulk solution.

To begin, we first consider a wide S-TI-S junction with a uniform in-plane
exchange field: $\mathbf{h}(\mathbf{R})=(h_{x},h_{y},0)$, where the junction
width is larger than its length $W_{F}\gg d_{F}$, as shown in Fig. 1. Deriving
the Green function and plugging into the expression for the charge current
density Eq. (9), using $I=\int_{-W_{F}/2}^{+W_{F}/2}dyJ_{x}(y)$, we obtain:

$\displaystyle I=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}{\cal
N}\sin\left(\phi+\frac{2h_{y}d_{F}}{\alpha}\right),$ $\displaystyle{\cal
N}=T\sum_{n}\frac{|\Delta|^{2}}{\omega_{n}^{2}+|\Delta|^{2}}\frac{\mathrm{csch}(k_{n}d_{F})}{k_{n}d_{F}},$
(10)

where $R_{N}={d_{F}}/{\sigma_{N}W_{F}}$ , $|\Delta|$ and $\phi$ denote the
superconducting gap and phase difference across the contact, and
$k_{n}=\sqrt{2|\omega_{n}|/D+(2h_{x}/\alpha)^{2}}$. As seen, the $h_{x}$
component of the magnetization, along the current flow direction, plays a
depairing role and supresses the supercurrent monotonically. The transverse
component $h_{y}$, however, causes a shift in the superconducting phase
difference. Similar effects were also predicted in the fully ballistic S-TI-S
counterparts via the Bogoliubov de Gennes approach tanaka_sfs . Therefore, one
concludes that the nonmagnetic disorders in a uniformly magnetized S-TI-S
heterostructure are unable to alter the phase shift in the current phase
relation and only modify the amplitude of the critical current through the
surface states of TI, the same as temperature.

We now turn to a 2D S-TI-S junction with helical magnetization depicted in
Fig. 1. The helical pattern is given by $\mathbf{h}(\mathbf{R})=h_{0}(\cos
Qy,\sin Qy,0)$, where $Q=q/W_{F}$ and $q$ determines the actual pattern of the
magnetization. The helical magnetization with a period of $\approx 10$nm was
already observed experimentally in manganese on a tungsten substrate through
spin-polarized tunneling experiments chrl_mg_exp . In order to solve Eq. (8)
we consider limits $W_{F}\gg d_{F}$ for $q\ll 1$ and $W_{F}\gg qd_{F}$ for
$q\gg 1$. We substitute the Green function, obtained for the configuration
shown in Fig. 1, into Eq. (9) and find the supercurrent:

$\displaystyle
I=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}\int_{-W_{F}/2}^{+W_{F}/2}\frac{dy}{W_{F}}{\cal
N}_{y}\sin\bigg{(}\phi+\frac{2h_{0}d_{F}}{\alpha}\sin Qy\bigg{)},$
$\displaystyle{\cal
N}_{y}=T\sum_{n}\frac{|\Delta|^{2}}{\omega_{n}^{2}+|\Delta|^{2}}\frac{\mathrm{csch}(k_{n,y}d_{F})}{k_{n,y}d_{F}},$
(11)

where $k_{n,y}=\sqrt{2|\omega_{n}|/D+(2h_{0}/\alpha)^{2}\cos^{2}Qy}$. As seen,
the integrand in Eq. (II) is a highly nonlinear function of $y$. We have
numerically integrated Eq. (II) and plotted the critical supercurrent as a
function of $q$ and $d_{F}$ in Fig. 2. To plot the currents, we have
normalized lengths by a length scale $\xi_{S}=\sqrt{D/|\Delta|}$, energies by
the superconducting gap $|\Delta|$, defined the Thouless energy
$\varepsilon_{T}=D/d_{F}^{2}$, and normalized the critical current $I^{c}$ by
its maximum value $I^{c}_{max}$ in the interval calculated. We consider a low
temperature regime and set a fixed normalized temperature at $T/T_{c}=0.01$,
where $T_{c}$ is the superconducting critical temperature. Fig. 2a illustrates
the critical supercurrent profile as a function of $q$ at
$h_{0}=5.0\pi|\Delta|$ for different values of $d_{F}$. At small values of
junction length, $d_{F}\approx 2.0\xi_{S}$, the supercurrent changes sign at a
single value of $0.4\pi<q<0.5\pi$. At larger values of junction length,
however, we see that the supercurrent undergoes multiple reversals and the
Josephson ground state oscillates between $0$ and $\pi$ phase differences.
Fig. 2b shows the critical current as a function the junction length for
different $q$. These patterns of critical current are reminiscent of those
found in Ref. eschrig,  for a 3D metallic superconductor-ferromagnet-
superconductor junction. One should note, however, that the nature of the
$0$-$\pi$ transitions explored in this paper is essentially different from
those of Ref. eschrig, . Here, the supercurrent reversals appear solely due to
the helical magnetization which is locked to the quasiparticle momentum.
Whereas in a metallic system the exchange field plays a dephasing role on the
Cooper pairs and causes oscillations in the Cooper pair amplitude which yields
supercurrent reversal. To gain more insights, we may simplify Eq. (II) by
considering a slow rotating magnetization so that $q<1$ and neglect the $x$
component of magnetization. In this limit, we find a simple expression for the
supercurrent which clearly illustrates the damped oscillatory behaviour of
critical current as a function of $qh_{0}d_{F}/\alpha$:

$\displaystyle I=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}{\cal
N}_{0}\frac{\sin(qh_{0}d_{F}/\alpha)}{qh_{0}d_{F}/\alpha}\sin\phi.$ (12)

The normalized critical supercurrent $I^{c}$ as a function of $h_{0}$ for
various values of $d_{F}$ is shown in Fig. 3. It is clear that the critical
supercurrent changes sign upon varying the magnetization intensity $h_{0}$.
Considering Eq. (12), the supercurrent $0$-$\pi$ crossovers appear at
$qh_{0}d_{F}/\alpha=n\pi$ ($n=\pm 1,\pm 2,...$). Taking $\alpha\approx 10^{7}$
cm/s and $h_{0}\approx 5$ meV, the first $0$-$\pi$ transition happens in a
junction of length $d_{F}=200$ nm for $q\approx 0.2\pi$.

In addition to the current density the absolute value of the Cooper pair wave
function brgt_rmp3  can be expressed in the slow rotation mode by:

$\displaystyle
U(y)=\frac{T}{4\gamma}\sum_{n}\frac{|\Delta|\text{csch}(k_{n}d_{F}/2)}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}\cos\bigg{(}\frac{\phi}{2}+\frac{h_{0}d_{F}}{\alpha}Qy\bigg{)},~{}~{}~{}~{}$
(13)

which can provide information about the proximity vortices. It is evident that
both the current density and $U(y)$ vanish at locations $y/W_{F}=\frac{\phi\pm
n\pi}{2qh_{0}d_{F}}\alpha$, $n=\pm 1,\pm 3,...$ provided that
$-1/2<y/W_{F}<+1/2$. These points correspond to the normal cores of proximity
vortices in TI nphys_2015_m . It is straightforward to show that the proximity
vortices are present in the case of helical magnetization and $q$ can control
the number of vortices.

Figure 3:  (Color online). Normalized critical supercurrent through the TI
surface channel with a slow rotating magnetization, $q=0.3\pi$, as a function
of magnetization intensity $h_{0}$ for different values for the junction
thickness $d_{F}=2.0\xi_{S},2.5\xi_{S},3.0\xi_{S},3.5\xi_{S}$.

##  III Derivations of main results

###  III.1 Eilenberger Equation

The Hamiltonian describing the Rashba type surface states of a three-
dimensional topological insulator (TI) with the proximity induced
ferromagnetism with an in-plane exchange field,
$\mathbf{h}(\mathbf{r})=(h_{x}(\mathbf{r}),h_{y}(\mathbf{r}),0)$, reads:

$\hat{H}(\mathbf{r})=-i\alpha({\bm{\nabla}}\times\mathbf{e}_{z})\cdot\hat{\bm{\sigma}}+\mathbf{h}(\mathbf{r})\cdot\hat{{\bm{\sigma}}}-\mu$
(14)

where $\alpha$ is the Fermi velocity, characterizing the surface states,
$\mathbf{e}_{z}$ is a unit vector normal to the surface of TI with chemical
potential $\mu$, $\hat{\bm{\sigma}}$ is a vector comprised of the Pauli
matrices, and the hat symbol ($\hat{1}$) denotes $2\times 2$ matrices. We here
consider the Rashba type throughout our calculations. However, through the
available symmetries in the formulations, our results can be easily extended
to the Dresselhaus type, $-i\beta_{D}{\bm{\nabla}}\cdot\hat{{\bm{\sigma}}}$,
by simple prescriptions given at the end of our Sec. III.2. Therefore, without
lossing generality, we derive the Eilenberger and Usadel equations for the
Rashba type surface states.

The electron Green function is defined by:

$\displaystyle G_{ab}(\tau_{1}-\tau_{2};\mathbf{r}_{1},\mathbf{r}_{2})$
$\displaystyle=$ $\displaystyle-i\langle
T_{\tau}\Psi_{a}(\tau_{1},\mathbf{r}_{2})\Psi_{b}^{{\dagger}}(\tau_{2},\mathbf{r}_{2})\rangle~{}~{}~{}~{}~{}~{}~{}$
(15a)
$\displaystyle\bar{G}_{ab}(\tau_{1}-\tau_{2};\mathbf{r}_{1},\mathbf{r}_{2})$
$\displaystyle=$ $\displaystyle-i\langle
T_{\tau}\Psi^{{\dagger}}_{a}(\tau_{1},\mathbf{r}_{1})\Psi_{b}(\tau_{2},\mathbf{r}_{2})\rangle~{}~{}~{}~{}~{}~{}~{}$
(15b) $\displaystyle F_{ab}(\tau_{1}-\tau_{2};\mathbf{r}_{1},\mathbf{r}_{2})$
$\displaystyle=$ $\displaystyle+i\langle
T_{\tau}\Psi_{a}(\tau_{1},\mathbf{r}_{1})\Psi_{b}(\tau_{2},\mathbf{r}_{2})\rangle~{}~{}~{}~{}~{}~{}~{}$
(15c) $\displaystyle
F^{{\dagger}}_{ab}(\tau_{1}-\tau_{2};\mathbf{r}_{1},\mathbf{r}_{2})$
$\displaystyle=$ $\displaystyle+i\langle
T_{\tau}\Psi^{{\dagger}}_{a}(\tau_{1},\mathbf{r}_{1})\Psi^{{\dagger}}_{b}(\tau_{2},\mathbf{r}_{2})\rangle,~{}~{}~{}~{}~{}~{}~{}$
(15d)

where $a\equiv\uparrow,\downarrow$ and $b\equiv\uparrow,\downarrow$ define the
electron spin projections, and $\tau_{1},\tau_{2}$ are the imaginary times at
$\mathbf{r}_{1},\mathbf{r}_{2}$, respectively. To simplify our derivations, we
have introduced a unitary rotation in the particle-hole and spin spaces:

$\displaystyle\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=\left(\begin{array}[]{cc}-\hat{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&-\hat{F}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})i\hat{\sigma}^{y}\\\
-i\hat{\sigma}^{y}\hat{F}^{{\dagger}}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&\hat{\sigma}^{y}\hat{\bar{G}}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\hat{\sigma}^{y}\end{array}\right),$

here we performed a Fourier transformation
$\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=T\sum_{n}e^{-i\omega_{n}\tau}\check{G}(\tau;\mathbf{r},\mathbf{r}^{\prime})$
to Matsubara frequency $\omega_{n}=\pi T(2n+1)$ where $T$ is the temperature
and $n\in\mathbb{Z}$. Check symbol $\check{1}$ defines $4\times 4$ matrices in
the particle-hole and spin spaces. The off diagonal components of the Green
function (III.1) matches to those of superconducting bulk solutions at the
boundaries and describe the penetration of the Cooper pairs into the TI
surface states.

Considering a Gaussian distribution for the nonmagnetic impurity scattering
potential $V(\mathbf{r})$, one finds

$\langle
V(\mathbf{r})V(\mathbf{r}^{\prime})\rangle=\frac{1}{\pi\nu\tau}\delta(\mathbf{r}-\mathbf{r}^{\prime}),\;\;\;\nu=\frac{\mu}{2\pi\alpha^{2}}$
(19)

where $\tau$ is the mean free time of particles in the disordered system and
$\nu$ is the density of states per spin at the Fermi level of the normal state
of TI. In this particle-hole and spin rotated system, we arrive at the
following equation for the particle Green function averaged over the
scattering potential:

$\displaystyle\left(\begin{array}[]{cc}-i\omega_{n}+\hat{H}(\mathbf{r})&-\epsilon_{g}\\\
\epsilon_{g}^{*}&i\omega_{n}+\hat{\sigma}^{y}\hat{H}^{*}(\mathbf{r})\hat{\sigma}^{y}\end{array}\right)\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})$
(22)
$\displaystyle=\delta(\mathbf{r}-\mathbf{r}^{\prime})+\frac{1}{\pi\nu\tau}\check{G}(\omega_{n};\mathbf{r},\mathbf{r})\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime}),$
(23)

the second term on the right hand side of Eq. (22) is defined by the self-
energy part of the averaged Green function. Here we allow for the proximity
induced minigap $\epsilon_{g}$ in the surface states of the TI.

We first subtract from Eq. (22) its conjugated equation and perform a Fourier
transformation with respect to the relative space arguments:

$\check{G}(\omega_{n};\mathbf{R}+\frac{\delta\mathbf{r}}{2},\mathbf{R}-\frac{\delta\mathbf{r}}{2})=\int\frac{d\mathbf{p}}{(2\pi)^{2}}\check{G}(\omega_{n};\mathbf{R},\mathbf{p})e^{i\mathbf{p}\cdot\delta\mathbf{r}},$
(24)

in which we have defined $\mathbf{R}=(\mathbf{r}+\mathbf{r}^{\prime})/2$,
$\delta\mathbf{r}=\mathbf{r}-\mathbf{r}^{\prime}$, and denoted the momentum
vector of the quasiparticles by $\mathbf{p}$. We also define a new parameter
$\xi_{p}=\alpha p-\mu$ and by making use of the fact that the Green function
peaks at the Fermi surface, we obtain the quasiclassical Green function

$\displaystyle\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{F})=\int\frac{d\xi_{p}}{\pi
i}\check{G}(\omega_{n};\mathbf{R},\mathbf{p}).$ (25)

Using this assumption, we obtain the Eilenberger equation, eiln :

$\displaystyle\frac{\alpha}{2}\\{\hat{{\bm{\eta}}},{\bm{\nabla}}\check{g}\\}=\bigg{[}\check{g},\omega_{n}\hat{\tau}^{z}+i\mathbf{h}\cdot\hat{{\bm{\sigma}}}\hat{\tau}^{z}+i\mu\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F}+i\check{\epsilon}_{g}+\frac{\langle\check{g}\rangle}{\tau}\bigg{]},$

where $\hat{\bm{{\bm{\eta}}}}=(-\hat{\sigma}^{y},\hat{\sigma}^{x})$,
$\check{\epsilon}_{g}=\hat{\sigma}^{0}(-\epsilon_{g}\hat{\tau}^{+}+\epsilon_{g}^{*}\hat{\tau}^{-})/2$,
and $\hat{\tau}^{\pm}=\hat{\tau}^{x}\pm i\hat{\tau}^{y}$. This equation
governs moving quasiparticles in ballistic (where $1/\tau\rightarrow 0$)
systems and those with weak nonmagnetic impurities with a finite $\tau$.

###  III.2 Usadel Equation

The Eilenberger equation can be simplified in systems with strong disorders
such that $\mu>1/\tau>|\omega_{n}|,|h|,|\epsilon_{g}|$. In this case, the
quasiparticles follow diffusive trajectories which is the so-called diffusive
regime usadel . In the diffusive regime of the surface states, one can
integrate the quasiclassical Green function, Eq. (25), over all possible
directions of quasiparticles’ momentum:

$\displaystyle\langle\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{F})\rangle\equiv\int_{0}^{2\pi}\frac{d\Omega_{\mathbf{n}_{F}}}{2\pi}\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{F}),\;\mathbf{n}_{F}=\frac{\mathbf{p}_{F}}{|\mathbf{p}_{F}|},~{}~{}~{}~{}~{}\;\;$
(27)

where $\Omega_{\mathbf{n}_{F}}$ is the polar angle of the vector
$\mathbf{n}_{F}$ in the plane of TI surface. We note that Eq. (III.1) is
similar to the kinetic equation for the Keldysh component of the Green
function of TI in the normal state Schwab .

To find a solution to Eq. (III.1), we expand the Green function through the
Pauli matrices:

$\check{g}=\frac{\hat{g}^{\prime}\hat{1}+\hat{\mathbf{g}}^{\prime}\cdot\hat{{\bm{\eta}}}+\hat{g}^{\prime}_{z}\hat{\sigma}^{z}}{2},$
(28)

where
$\hat{\mathbf{g}}^{\prime}=(\hat{g}^{\prime}_{x},\hat{g}^{\prime}_{y},0)$. We
assume that the Fermi energy in the conduction band of the TI is much larger
than all other energy scales available in the system. We immediately find that
leading contributions to $\check{g}$ commute with $\hat{\eta}\mathbf{n}_{F}$.
Thus, $\hat{g}^{\prime}_{z}$ component is smaller than $\hat{g}^{\prime}$ and
$\hat{g}^{\prime}_{x,y}$ by a factor of
$\sim\mathrm{max}(\frac{1}{\tau},|\mathbf{h}|)/\mu\ll 1$ and can be neglected.
Also, the spin structure of $\check{g}$ is proportional to the conduction band
projector $(1+\hat{{\bm{\eta}}}\mathbf{n}_{F})/2$ in the limit of
$|\mathbf{h}|\ll\mu$. We then set
$\hat{\mathbf{g}}^{\prime}=\mathbf{n}_{F}\hat{g}^{\prime}$ and propose a
solution to Eq. (III.1) in the form of a direct product of two $2\times 2$
matrices:

$\displaystyle\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{F})=\hat{g}^{\prime}(\omega_{n};\mathbf{R},\mathbf{n}_{F})\frac{1+\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F}}{2},$
(29)

where we define

$\displaystyle\hat{g}^{\prime}=\left(\begin{array}[]{cc}-g&-if\\\
if^{*}&\bar{g}\end{array}\right).$ (32)

Notice that the conduction band projector satisfies:
$(1+\hat{{\bm{\eta}}}\mathbf{n}_{F})^{2}/4=(1+\hat{{\bm{\eta}}}\mathbf{n}_{F})/2$
and serves as the source of superconducting triplet correlations in the
system. Substituting expression (29) into Eq. (22), we arrive at the following
equation:

$\displaystyle\alpha(\hat{{\bm{\eta}}}+\mathbf{n}_{F})\cdot{\bm{\nabla}}\hat{g}^{\prime}$
$\displaystyle=$
$\displaystyle\bigg{[}\hat{g}^{\prime}(1+\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F}),\omega_{n}\hat{\tau}^{z}+i\mathbf{h}\cdot\hat{{\bm{\sigma}}}\hat{\tau}^{z}$
(33) $\displaystyle+$ $\displaystyle
i\check{\Delta}+\frac{\langle\hat{g}^{\prime}(1+\hat{{\bm{\eta}}}\cdot\mathbf{n}_{F})\rangle}{2\tau}\bigg{]}.~{}~{}~{}~{}$

In the limit of diffusive motion of quasiparticles, the Green function can be
expanded through the first two terms of 2D harmonics, namely:

$\hat{g}^{\prime}(\omega_{n};\mathbf{R},\mathbf{n}_{F})=\hat{g}_{s}(\omega_{n};\mathbf{R})+\mathbf{n}_{F}\cdot\hat{\mathbf{g}}_{a}(\omega_{n};\mathbf{R}),$
(34)

where the zero harmonic in the expansion is isotropic and its amplitude is
much larger than the first harmonic:
$\hat{g}_{s}\gg\mathbf{n}_{F}\cdot\hat{\mathbf{g}}_{a}$.

We now substitute the expanded Green function, Eq. (34), into Eq. (33) and
perform an integration over momentum directions. By taking a spin trace we
finally arrive at:

$\displaystyle\alpha\hat{{\bm{\nabla}}}\cdot\hat{\mathbf{g}}_{a}=2[\hat{g}_{s},\omega_{n}\hat{\tau}^{z}+i\check{\epsilon}_{g}].$
(35)

Here, we have defined a new operator $\hat{{\bm{\nabla}}}$ so that

$\hat{{\bm{\nabla}}}X={\bm{\nabla}}X+\frac{i}{\alpha}(h_{x}\mathbf{e}_{y}-h_{y}\mathbf{e}_{x})[\hat{\tau}^{z},X],$
(36)

where $\mathbf{e}_{x,y}$ are unit vectors in the $x,y$ directions. Multiplying
Eq. (33) by $\mathbf{n}_{F}$ and integrating it over momentum directions, we
find the following expression for the first harmonic term of the Green
function expansion $\hat{\mathbf{g}}_{a}$:

$\displaystyle\hat{\mathbf{g}}_{a}=-2\alpha\tau\hat{g}_{s}\hat{{\bm{\nabla}}}\hat{g}_{s}.$
(37)

To obtain the Usadel equation, it suffices we substitute Eq. (37) into Eq.
(35) which yields:

$\displaystyle
D\hat{{\bm{\nabla}}}\cdot(\hat{g}_{s}\hat{{\bm{\nabla}}}\hat{g}_{s})=[\omega_{n}\hat{\tau}^{z}+i\check{\epsilon}_{g},\hat{g}_{s}],$
(38)

where the diffusion coefficient is denoted by $D=\alpha^{2}\tau$. In the case
$\check{\epsilon}_{g}=0$, Eq. 38 leads to Eq. 8. The singlet part of the Green
function, $\hat{g}_{s}$, satisfies the Usadel equation Eq. (38). In contrast,
due to fast spin relaxation time which is given by the mean free time in the
TI, the motion of the spin part, $\hat{\mathbf{g}}_{a}$, is not diffusive, and
satisfies Eq. (37). To obtain the Eilenberger and Usadel equations for surface
states with the Dresselhaus type of spin-orbit coupling, one simply should use
${\bm{\eta}}=(\sigma_{x},\sigma_{y})$ in the Eilenberger equation, Eq.
(III.1), and change $\alpha$ to $\beta_{D}$. In the Usadel equation, Eq. (38),
however one can perform the replacement below:

$\hat{{\bm{\nabla}}}X={\bm{\nabla}}X+\frac{i}{\beta_{D}}(h_{x}\mathbf{e}_{x}+h_{y}\mathbf{e}_{y})[\hat{\tau}^{z},X].$
(39)

###  III.3 Boundary Conditions and the Current Density

We consider low transparency limit of the interface (tunneling barrier)
between the TI and superconductor boundary_c1 ; boundary_c2 . We neglect the
inverse proximity effect so that the Green function in the superconductor at
the interface is given by its bulk solution. Assuming that the electron
tunneling across the TI-S interface is spin-conserving, while it is not
momentum conserving, we eventually arrive at the following expression

$2\gamma\hat{g}_{s}\mathbf{n}\cdot\hat{{\bm{\nabla}}}\hat{g}_{s}=[\hat{g}_{s},\hat{g}_{\text{SC}}],$
(40)

where $\mathbf{n}$ is the unit vector normal to the boundary and

$\hat{g}_{\text{SC}}(\omega_{n})=\frac{1}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}[\omega_{n}\hat{\tau}^{z}+\frac{i}{2}(-\Delta\hat{\tau}^{+}+\Delta^{*}\hat{\tau}^{-})],$
(41)

is the Green function in the superconductor, in which $|\Delta|$ is the
superconducting gap. The low transparency of the interface leads to large
parameter $\gamma\gg 1$, which is the ratio of resistance per unit area of the
surface of the tunneling barrier to the resistivity of the TI.

Let us now present the current density in the TI. The current density flowing
across the surface of TI is given by:

$\mathbf{J}=\frac{e\alpha}{2}T\sum_{n}\lim_{\mathbf{r}^{\prime}\rightarrow\mathbf{r}}\mathrm{Tr}[\hat{\tau}^{z}\hat{{\bm{\eta}}}\check{G}(\omega_{n},\mathbf{r},\mathbf{r}^{\prime})],$
(42)

where $e>0$ is the absolute value of electron charge. Rewriting the current
density through the quasiclassical Green function, we find:

$\mathbf{J}=i\pi\frac{e\alpha}{4}T\sum_{n}\mathrm{Tr}[\hat{\tau}^{z}\hat{\mathbf{g}}_{a}].$
(43)

Substituting expression (37) into (42), we obtain the current density in the
diffusive limit, Eq. (9).

###  III.4 Derivation of Supercurrent through Topological Insulator with
Uniform and Helical Magnetizations

In the low proximity limit, we can expand the Green function around the bulk
solution which yields:

$\displaystyle\hat{g}_{s}(\omega_{n};\mathbf{R})=\left(\begin{array}[]{cc}\mathrm{sign}(\omega_{n})&-i{f_{+}}(\omega_{n};\mathbf{R})\\\
i{f_{-}}(\omega_{n};\mathbf{R})&-\mathrm{sign}(\omega_{n})\end{array}\right).$
(46)

The low proximity limit is experimentally relevant and can be easily achieved
in temperatures near the superconducting critical temperature or low
transparent SC-TI contacts, for instance. To derive the current density, and
consequently the total supercurrent through a helical magnetization, we use
the Usadel equation, (38). We will consider homogeneous and helical
magnetization as shown in Fig. 1 of main text, which can be described by:

$\displaystyle\mathbf{h}=(h_{x},h_{y},0),$ (47a)
$\displaystyle\mathbf{h}(y)=h_{0}(\cos Qy,\sin Qy,0),$ (47b)

where $Q=q/W_{F}$. To avoid complication, we assume that the inverse proximity
effect is neglegible and the magnetization is restricted within $0<x<d_{F}$.
Moreover, we consider a situation where the Josephson penetration length
$\lambda_{J}$ is larger than the junction width and ignore the effect of
magnetic field induced by the supercurrent itself zagoskin ; birge .
Otherwise, one should solve the Usadel equation together with the Maxwell
equations self-consistently. We assume that the junction width ($W_{F}$) is
larger than the length ($d_{F}$), $W_{F}\gg d_{F}$ and arrive at the following
differential equations for the anomalous components of the Green function:

$\displaystyle\left(\partial_{x}\mp\frac{2ih_{y}(y)}{\alpha}\right)^{2}f_{\pm}-\frac{4h^{2}_{x}(y)}{\alpha^{2}}f_{\pm}=\frac{2|\omega_{n}|}{D}f_{\pm}.$
(48)

The corresponding boundary conditions in this case result in:

$\displaystyle\gamma[\partial_{x}\mp
2ih_{y}(y)/\alpha]f_{\mp}|_{x=-d_{F}/2}=\frac{|\Delta|e^{\mp
i\phi/2}}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}},~{}~{}~{}~{}~{}~{}$ (49a)
$\displaystyle\gamma[\partial_{x}\mp
2ih_{y}(y)/\alpha]f_{\mp}|_{x=+d_{F}/2}=-\frac{|\Delta|e^{\pm
i\phi/2}}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}},~{}~{}~{}~{}~{}~{}$ (49b)

where $\phi$ is the superconducting phase difference across the junction and
the transparency of SC-TI contacts can be controlled through parameter
$\gamma$. To derive boundary conditions (49), we assume that
$\Delta\gg|\omega_{n}|f_{+}$, which is justified by the SC-TI interface with
low transparency. The Usadel equation and associated boundary conditions yeild
the following solutions:

$\displaystyle f_{+}$ $\displaystyle=$ $\displaystyle\frac{-|\Delta|}{\gamma
k_{n,y}\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}\bigg{[}\frac{\mathrm{ch}(k_{n,y}(x-d_{F}/2))}{\mathrm{sh}(k_{n,y}d_{F})}e^{i\phi/2+i\frac{2h_{y}(y)}{\alpha}(x+d_{F}/2)}+\frac{\mathrm{ch}(k_{n,y}(x+d_{F}/2))}{\mathrm{sh}(k_{n,y}d_{F})}e^{-i\phi/2+i\frac{2h_{y}(y)}{\alpha}(x-d_{F}/2)}\bigg{]},~{}~{}~{}~{}~{}~{}$
(50a) $\displaystyle f_{-}$ $\displaystyle=$
$\displaystyle\frac{-|\Delta|}{\gamma
k_{n,y}\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}\bigg{[}\frac{\mathrm{ch}(k_{n,y}(x-d_{F}/2))}{\mathrm{sh}(k_{n,y}d_{F})}e^{-i\phi/2-i\frac{2h_{y}(y)}{\alpha}(x+d_{F}/2)}+\frac{\mathrm{ch}(k_{n,y}(x+d_{F}/2))}{\mathrm{sh}(k_{n,y}d_{F})}e^{i\phi/2-i\frac{2h_{y}(y)}{\alpha}(x-d_{F}/2)}\bigg{]},~{}~{}~{}~{}~{}~{}$
(50b)

where we have introduced a wave-vector

$k_{n,y}=\sqrt{2|\omega_{n}|/D+(2h_{x}(y)/\alpha)^{2}}.$ (51)

The supercurrent density in the low proximity limit we consider here involves
the anomalous components of Green function:

$\displaystyle J_{x}(y)$ $\displaystyle=$ $\displaystyle-\frac{\pi
i\sigma_{N}}{4e}T\sum_{n}\Big{\\{}f_{+}(\partial_{x}+2ih_{y}(y)/\alpha)f_{-}$
(52) $\displaystyle-$ $\displaystyle
f_{-}(\partial_{x}-2ih_{y}(y)/\alpha)f_{+}\Big{\\}}.$

To find the total current flow across the junction in the $x$ direction, one
needs to integrate the current density, given by expression (52), over the
junction width, $W_{F}$, along the $y$ axis:
$I=\int_{-W_{F}/2}^{+W_{F}/2}dyJ_{x}(y)$. By plugging the solutions found,
(50), into (52), we obtain the total charge supercurrent, flowing through the
junction;

$\displaystyle
I=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}\int_{-W_{F}/2}^{+W_{F}/2}\frac{dy}{W_{F}}{\cal
N}_{y}\sin\bigg{(}\phi+\frac{2h_{y}(y)d_{F}}{\alpha}\bigg{)},$
$\displaystyle{\cal
N}_{y}=T\sum_{n}\frac{|\Delta|^{2}}{\omega_{n}^{2}+|\Delta|^{2}}\frac{\mathrm{csch}(k_{n,y}d_{F})}{k_{n,y}d_{F}}.$
(53)

It is instructive to consider several limiting cases of expression
(LABEL:cur_sup). In the limit of homogeneous exchange field we immediately
obtain Eq. (II). Also, by considering the slow rotation limit, where $q<1$, of
the helical magnetization given by Eq. (47) we find Eq. (12). As seen, the
critical Josephson current in this slow rotating regime, $q<1$, shows similar
features, as a function of $qh_{0}d_{F}/\alpha$, to those of conventional
metallic wide junctions subject to an external magnetic field.

Figure 4:  (Color online). Current density $J(y)$ and the absolute value of
the Cooper pair wave function $U(y)$, normalized by their maximum values, as a
function of $y$ coordinate along the junction width at the middle of junction
$x=0$ (see Fig. 1 of main text) for three different magnetization rotation
parameter $q=0.2\pi,0.4\pi$, and $0.6\pi$. The superconducting phase
difference is set fixed at $\phi=\pi/2$ and the magnetization intensity is
equal to $h_{0}=6.0\pi\Delta$.

###  III.5 Proximity Vortices

We have full numerically solved the Usadel equation in the presence of the
helical magnetization given by Eq. (47). Substituting the solutions into the
current density expression, Eq. (52), and the absolute value of the Cooper
pair wave function below

$U(x,y)\equiv-\frac{T}{8}\sum_{n}\big{\\{}f_{+}+f_{-}\big{\\}},$ (54)

we determine the spatial maps of the supercurrent and Cooper pair wave
function inside the surface states of TI. Our results demonstrate that the
current density and Cooper pair wave function are both zero at specific
locations along the junction width at $x=0$ which indicates the nucleation of
proximity vortices in this class of Josephson junctions. The existence of such
proximity vortices were quite recently confirmed experimentally in 3D metallic
SNS junctions nphys_2015_m . To gain better insights, we investigate the
proximity vortices through the analytical expressions derived from the Usadel
equation. Let us first investigate the slow rotating magnetization case which
results in clearer conclusions through the analytics and then discuss the
helical magnetization which demands numerics. Considering a slow rotating
magnetization $q<1$ and setting $x=0$ in the solutions (50) we arrive at the
following current density and Cooper pair wave function:

$\displaystyle J_{x}(y)=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}{\cal
N}_{0}\sin\bigg{(}\phi+\frac{2h_{0}d_{F}}{\alpha}Qy\bigg{)},~{}~{}~{}$ (55)

and

$\displaystyle
U(y)=\frac{T}{4\gamma}\sum_{n}\frac{|\Delta|\text{csch}(k_{n}d_{F}/2)}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}\cos\bigg{(}\frac{\phi}{2}+\frac{h_{0}d_{F}}{\alpha}Qy\bigg{)}.~{}~{}~{}~{}~{}~{}~{}$
(56)

The current density vanishes at $y/W_{F}=\frac{\phi\pm
n\pi}{2qh_{0}d_{F}}\alpha$ where $n=0,\pm 1,\pm 2,...$ while the zeros of
$U(y)$ appear at $y/W_{F}=\frac{\phi\pm m\pi}{2qh_{0}d_{F}}\alpha$ where
$m=\pm 1,\pm 3,...$ . The extra zeros in the current density can be undrestood
by noting the circulating form of quasiparticles’ paths which cancel each
other at $n=0,\pm 2,\pm 4,...$ . Therefore, at $y/W_{F}=\frac{\phi\pm
n\pi}{2qh_{0}d_{F}}\alpha$, $n=\pm 1,\pm 3,...$ provided that
$-1/2<y/W_{F}<+1/2$ the current density and $U(y)$ both vanish which
determines the location of normal core of each vortex.

Next, we consider the more complicated case where the magnetization follows
the pattern given in Eq. (47). In this case we obtain the current density and
$U(y)$ as follows:

$\displaystyle J_{x}(y)=\frac{\pi}{2e}\frac{d_{F}^{2}}{\gamma^{2}R_{N}}{\cal
N}_{y}\sin\bigg{(}\phi+\frac{2h_{0}d_{F}}{\alpha}\sin Qy\bigg{)},~{}~{}~{}$
(57)

and

$\displaystyle U(y)$ $\displaystyle=$
$\displaystyle\frac{T}{4\gamma}\sum_{n}\frac{|\Delta|\text{csch}(k_{n,y}d_{F}/2)}{\sqrt{\omega_{n}^{2}+|\Delta|^{2}}}\cos\bigg{(}\frac{\phi}{2}+\frac{h_{0}d_{F}}{\alpha}\sin
Qy\bigg{)},$ $\displaystyle k_{n,y}$ $\displaystyle=$
$\displaystyle\sqrt{2|\omega_{n}|/D+(2h_{0}\cos Qy/\alpha)^{2}}.$ (58)

To determine the location of proximity vortices, we have now normalized the
current density and $U(y)$ by their maximum values: $J_{max}$ and $U_{max}$.
Fig. 4 exhibits the normalized current density and $U(y)$ at the middle of
junction $x=0$ as a function of location along the junction width $y$ for
three different rotation degrees: $q=0.2\pi,0.4\pi,$ and $0.6\pi$. The
junction length is assumed $d_{F}=0.35\xi_{S}$, $h_{0}=6.0\pi\Delta$, and the
superconducting phase difference $\phi=\pi/2$. As seen, the rotation degree of
magnetization can change the sign of current density and induces proximity
vortices in addition to the other parameters $h_{0}$, $d_{F}$, and $\phi$
which can alter the patterns. The current density shows extra zeros compared
to $U(y)$ due to the cancelation process described above. The vertical lines
indicate the location of proximity vortices’ core where $U(y)$ and current
density vanish both. Notice that a proximity vortex can move along the
junction simply by modulating the phase different $\phi$. In other words, the
proximity induced vortices do not necessarily reside along the junction width
in a symmetric fashion with respect to $y=0$.

Finally, we note that our results on the Josephson current and proximity
vortices can be also valid in the Eilenberger limit with more sophistications
in final expressions zagoskin ; ykby . These complications arise due to the
quasiballistic motion and multiple Andreev reflections at the topological
insulator - superconductor interfaces. One the other hand diffusive regime
allows for highly simplifying and transparent calculations.

##  IV Conclusions

In conclusion, we have derived the Eilenberger and Usadel equations which
describe the quasiparticles in the disordered surface states of a TI in the
presence of superconductivity and magnetism. We employed this approach to
study the supercurrent flow through a TI with proximity induced in-plane
helical magnetization. In contrast to the case of S-TI-S junction with uniform
magnetization, our results reveal that the helical magnetization can induce
multiple supercurrent reversals and proximity vortices upon varying the
junction thickness, magnetization strength, and helical magnetization
parameters.

##  V Acknowledgments

We acknowledge support from the Swiss NF and NCCR QSIT. We would like to thank
A. Yu. Zyuzin for helpful discussions and G. Sewell for his valuable
instructions in the numerical parts of this work.

## References

  * (1) O. A. Pankratov, S. V. Pakhomov, and B. A. Volkov, Solid State Comm. 61, 93 (1987).
  * (2) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (3) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (4) G. E. Volovik, The Universe in a Helium Droplet (Oxford University Press, Oxford, 2003).
  * (5) J. R. Williams, A. J. Bestwick, P. Gallagher, Seung Sae Hong, Y. Cui, Andrew S. Bleich, J. G. Analytis, I. R. Fisher, and D. Goldhaber-Gordon, Phys. Rev. Lett. 109, 056803 (2012).
  * (6) M. Veldhorst, M. Snelder, M. Hoek, T. Gang, V. K. Guduru, X. L. Wang, U. Zeitler, W. G. van der Wiel, A. A. Golubov, H. Hilgenkamp, and A. Brinkman, Nat. Mater. 11 417 (2012).
  * (7) J. B. Oostinga, L. Maier, P. Schuffelgen, D. Knott, C. Ames, C. Brune, G. Tkachov, H. Buhmann, and L. W. Molenkamp, Phys. Rev. X 3, 021007 (2013).
  * (8) L. Galletti, S. Charpentier, M. Iavarone, P. Lucignano, D. Massarotti, R. Arpaia, Y. Suzuki, K. Kadowaki, T. Bauch, A. Tagliacozzo, F. Tafuri, and F. Lombardi, Phys. Rev. B 89, 134512 (2014).
  * (9) K. C. Nowack, E. M. Spanton, M. Baenninger, M. Konig, J. R. Kirtley, B. Kalisky, C. Ames, P. Leubner, C. Brune, H. Buhmann, L. W. Molenkamp, D. Goldhaber-Gordon, and K. A. Moler, Nat. Mater. 12, 787 (2013).
  * (10) E. M. Spanton, K. C. Nowack, L. Du, G. Sullivan, R.-R. Du, and K. A. Moler, Phys. Rev. Lett. 113, 026804 (2014).
  * (11) E. Y. Ma, M. R. Calvo, J. Wang, B. Lian, M. Muhlbauer, C. Brune, Y.-T. Cui, K. Lai, W. Kundhikanjana, Y. Yang, M. Baenninger, M. Konig, C. Ames, H. Buhmann, P. Leubner, L. W. Molenkamp, S.-C. Zhang, D. Goldhaber-Gordon, M. A. Kelly, and Z.-X. Shen, Nat. Commun. 6 7252 (2015).
  * (12) A. A. Burkov and D. G. Hawthorn, Phys. Rev. Lett. 105, 066802 (2010).
  * (13) P. Schwab, R. Raimondi and C. Gorini, Europhys. Lett. 93, 67004 (2011).
  * (14) C. Kurter, A.D.K. Finck, Y. S. Hor and D. J. Van Harlingen, Nat. Comm. 6, 7130 (2015).
  * (15) A. A. Golubov, M. Y. Kupriyanov, and E. Ilichev, Rev. Mod. Phys. 76, 411 (2004).
  * (16) F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Rev. Mod. Phys. 77, 1321 (2005).
  * (17) A. I. Buzdin, Rev. Mod. Phys. 77, 935 (2005).
  * (18) I. O. Kulik, JETP 22, 841 (1966).
  * (19) L. N. Bulaevskii, V. V. Kuzii, and A. A. Sobyanin, JETP Lett. 25, 290 (1977).
  * (20) A. I. Buzdin, L. N. Bulaevskii, and S. V. Panjukov, JETP Lett. 35, 178 (1982).
  * (21) V. V. Ryazanov, V. A. Oboznov, A. Yu. Rusanov, A. V. Veretennikov, A. A. Golubov, and J. Aarts, Phys. Rev. Lett. 86, 2427 (2001).
  * (22) T. Kontos, M. Aprili, J. Lesueur, F. Genet, B. Stephanidis, and R. Boursier, Phys. Rev. Lett. 89, 137007 (2002).
  * (23) Y. Blum, A. Tsukernik, M. Karpovski, and A. Palevski, Phys. Rev. Lett. 89, 187004 (2002).
  * (24) I. V. Krive, L. Y. Gorelik, R. I. Shekhter, and M. Jonson, Fiz. Nizk. Temp. 30, 535 (2004).
  * (25) V. Braude and Yu. V. Nazarov, Phys. Rev. Lett. 98, 077003 (2007).
  * (26) A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Phys. Rev. Lett. 101, 107001 (2008).
  * (27) A. I. Buzdin, Phys. Rev. Lett. 101, 107005 (2008).
  * (28) T. Yokoyama, M. Eto, and Y. V. Nazarov. Phys. Rev. B 89, 195407 (2014).
  * (29) D. B. Szombati, S. Nadj-Perge, D. Car, S. R. Plissard, E. P. A. M. Bakkers and L. P. Kouwenhoven, Nature Phys. 3742 (2016).
  * (30) A. Yu. Zyuzin and B. Spivak, Phys. Rev. B 61, 5902, (2000).
  * (31) Y. Tanaka, T. Yokoyama, and N. Nagaosa, Phys. Rev. Lett. 103, 107002 (2009).
  * (32) G. Eilenberger, Z. Phys. 214, 195 (1968).
  * (33) K. D. Usadel, Phys. Rev. Lett. 25, 507 (1977).
  * (34) A. V. Zaitsev, Sov. Phys. JETP 59, 1015 (1984).
  * (35) M. Y. Kuprianov and V. F. Lukichev, Sov. Phys. JETP 67, 1163 (1988).
  * (36) F. S. Bergeret, I. V. Tokatly, Phys. Rev. B 89, 134517 (2014); M. Houzet, J. S. Meyer, Phys. Rev. B 92, 014509 (2015).
  * (37) M. Bode, M. Heide, K. von Bergmann, P. Ferriani, S. Heinze, G. Bihlmayer, A. Kubetzka, O. Pietzsch, S. Blugel and R. Wiesendanger, Nature 447, 190 (2007).
  * (38) T. Champel, T. Lofwander, and M. Eschrig, Phys. Rev. Lett. 100, 077003 (2008).
  * (39) D. Roditchev, C. Brun, L. Serrier-Garcia, J. C. Cuevas, V. H. L. Bessa, M. V. Milosevic, F. Debontridder, V. Stolyarov and T. Cren, Nat. Phys. Lett. 11, 332 (2015).
  * (40) V. Barzykin and A. M. Zagoskin, Superlattices Microstruct. 25, 797 (1999).
  * (41) M. S. Crosser, J. Huang, F. Pierre, P. Virtanen, T. T. Heikkila, F. K. Wilhelm, and N. O. Birge, Phys. Rev. B 77, 014528 (2008).
  * (42) S. Hart, H. Ren, M. Kosowsky, G. Ben-Shach, P. Leubner, C. Brune, H. Buhmann, L. W. Molenkamp, B. I. Halperin, and A. Yacoby, arXiv:1509.02940.