# Fractional Shapiro steps without fractional Josephson effect

Artem V. Galaktionov  I.E. Tamm Department of Theoretical Physics, P.N.
Lebedev Physical Institute, 119991 Moscow, Russia  Andrei D. Zaikin  Institute
for Quantum Materials and Technologies, Karlsruhe Institute of Nanotechnology
(KIT), 76021, Karlsruhe, Germany  I.E. Tamm Department of Theoretical Physics,
P.N. Lebedev Physical Institute, 119991 Moscow, Russia

###### Abstract

It is widely believed that superconducting junctions involving topological
insulators and hosting Majorana-like bound states may exhibit unusual
”fractional” ($4\pi$-periodic) ac Josephson effect. Accordingly, ”fractional”
Shapiro steps on the current-voltage characteristics of such junctions are
expected to occur under external microwave radiation. Here, we microscopically
evaluate Shapiro steps in topologically trivial highly transparent
superconducting weak links. The key features recovered within our analysis –
including, e.g., the so-called ”missing” Shapiro steps – turn out to be
similar to those observed in topological Josephson junctions. Our results
demonstrate that caution is needed while interpreting experimental results for
superconducting weak links in terms of Majorana physics.

###### pacs:

##  I Introduction and motivation

Over recent years a great deal of attention was paid to experimental and
theoretical investigations of superconducting junctions involving topological
insulators which can host non-trivial gapless bound states. Such Majorana-like
states are believed to exhibit rather exotic physical properties and are also
regarded as promising in the context of topologically protected quantum
computation AK .

Theoretically a doublet of Majorana-like bound states was identified in point
contacts formed by $p$-wave superconductors Yakovenko  as well as in junctions
composed of two conventional ($s$-wave) superconductors connected to each
other via two- and three-dimensional topological insulators FK ; 3D . In all
these cases the energies of these Majorana-like modes $E^{M}_{\pm}$ depend on
the superconducting phase difference $\varphi$ across the junction as
Yakovenko ; FK ; 3D

$E^{M}_{\pm}(\varphi)=\pm\delta\cos\frac{\varphi}{2}.$ (1)

For reasons and under conditions to be specified below it is widely believed
that junctions hosting such Majorana-like modes exhibit an unusual
$4\pi$-periodic current-phase relation (CPR) Yakovenko ; FK
$I_{M}\propto\sin\frac{\varphi}{2}$ which is fundamentally different from the
standard $2\pi$-periodic Josephson CPR $I_{J}=I_{c}\sin\varphi$. If so, one is
bound to conclude that – having in mind the Josephson relation
$\dot{\varphi}=2eV$ (where $e$ defines electron charge and $V$ stands for bias
voltage) – the supercurrent oscillation frequency in such junctions equals to
$eV$, i.e to a half of the standard Josephson frequency $\omega_{J}=2eV$.

The presence of Majorana-like bound states (1) could then be detected
experimentally, e.g., by observing resonances in the form of current jumps
(the so-called Shapiro steps Tinkh ; BP ; Likh ) on the the junction $I-V$
curve under the influence of external microwave radiation. Pronounced Shapiro
steps at frequencies $\omega=\omega_{J}/2$ were indeed detected in a number of
microwave experiments performed with HgTe- and BiSb-based superconducting
junctions M16 ; topmat ; Gre .

Very generally, internal properties of any weak link between two conventional
superconductors – no matter how complicated these properties are – may
influence its critical current and the form of CPR but not its periodicity in
$\varphi$. The latter, in turn, is fundamentally determined by the periodicity
of the charge space SZ90  reciprocal to the phase one. At subgap energies and
temperatures the charge space of our problem remains effectively $2e$-periodic
which immediately enforces $2\pi$-periodic in $\varphi$ CPR.

In the light of this, $4\pi$-periodic CPR could be regarded as quite unusual
since in that case the charge transfer between superconducting condensates on
both sides of the junction should be provided by quasiparticles with charge
$e$ rather than by Cooper pairs with charge $2e$. It remains unclear to us
which physical mechanism could be responsible for such kind of charge transfer
between two conventional superconductors.

For completeness, we also remark that general effective actions SZ90 ; Z94 ;
SN ; book  describing superconducting junctions with arbitrary transmissions
do, of course, include $4\pi$-periodic terms in the Josephson phase $\varphi$.
However, such (non-local in time) terms account for dissipative currents and
have nothing to do with the supercurrent flowing across the weak link.

Below we will consider a purely ballistic superconducting weak link (or an SNS
junction) with ${\mathcal{N}}$ fully transparent conducting channels and
normal state conductance $1/R_{N}={\mathcal{N}}e^{2}/\pi$. Provided the
thickness of a normal (N) layer $d$ connecting two superconducting (S)
electrodes is much smaller than the coherence length $d\ll\xi_{0}\sim
v_{F}/\Delta$ the junction hosts a pair of subgap Andreev bound states

$E^{A}_{\pm}(\varphi)=\pm\Delta\cos\frac{\varphi}{2}$ (2)

in each of ${\mathcal{N}}$ conducting modes, see also Fig. 1. Here and below
$\Delta$ stands for the superconducting order parameter in the electrodes and
$v_{F}$ is the Fermi velocity. Since the phase dependence of the bound states
energy (2) is identical to that for Majorana-like states (1), ballistic SNS
junctions can serve as a convenient playground to test some ideas and
approaches to superconducting junctions based on topological insulators.

Figure 1: A pair of Andreev bound states $E^{A}_{\pm}(\varphi)$ in short
ballistic SNS junctions.

In equilibrium CPR of short ballistic SNS junctions takes the well-known
$2\pi$-periodic in $\varphi$ form KO

$I=I_{S}(\varphi)=\frac{\pi\Delta}{eR_{N}}\sin\frac{\varphi}{2}\tanh\left(\frac{\Delta}{2T}\cos\frac{\varphi}{2}\right).$
(3)

Exactly the same relation also holds for short symmetric SINIS junctions at
resonance for arbitrary (including very low) transmissions of SN interfaces
GZ02 . The relation (3) can also be recovered from a simple formula FT ; BH

$I=2e{\mathcal{N}}\sum_{\pm}\frac{\partial
E^{A}_{\pm}}{\partial\varphi}f_{\pm},$ (4)

where $f_{\pm}\equiv f_{F}(E^{A}_{\pm})=1/[1+\exp(E^{A}_{\pm}/T)]$ are the
Fermi filling factors FN  for both Andreev bound states. This formula
indicates that equilibrium CPR in weak links with $d\ll\xi_{0}$ can be
associated only with discrete Andreev levels (2), whereas continuous electron
spectrum does not play any role in this limit.

Making use of this observation, the same formula with $E^{A}_{\pm}\to
E^{M}_{\pm}$ was employed Yakovenko ; FK ; 3D  in order to identify the
contribution from Majorana-like states (1) to the supercurrent flowing across
the junctions involving topological insulators. In equilibrium this
contribution is again $2\pi$-periodic in $\varphi$, just like CPR in Eq. (3).

One may be tempted to extend the formula (4) in order to describe ac Josephson
effect. Following Yakovenko , by sweeping the phase $\varphi$ sufficiently
slowly, but not too slowly, one would be able to drive the system in such a
way that it has no time to relax to its ground state, thereby always occupying
one of the states (1) and leaving another one empty. Assuming $\delta<\Delta$
(which typically requires violation of time-reversal symmetry) and
$eV<\Delta-\delta$ in order to separate Majorana-like levels from the
continuum of states FK  and substituting the filling factors $f_{-}=1$ and
$f_{+}=0$ into Eq. (4) with $E^{A}_{\pm}\to E^{M}_{\pm}$ one immediately
arrives at the contribution to the supercurrent in the desired form
$I_{M}\propto\sin(\varphi(t)/2)$.

In our view, this procedure cannot yet be regarded as a rigorous derivation
and in any case it should be verified by means of an elaborate quantum kinetic
analysis. This analysis should fully account, e.g., for the effect of multiple
Andreev reflection (MAR) MAR  which is known to play a crucial role in
junctions with not very low transmissions. Fortunately in the case of
ballistic SNS junctions considered here such detailed microscopic theory is
already well established Zaitsev80 ; Zaikin83 ; Uwe ; AB1 ; AB2 ; Uwe2 . This
theory demonstrates that MAR serves as a key charge transfer mechanism across
such junctions setting in already at arbitrarily small bias voltages
$V=\dot{\varphi}/2e$.

Figure 2: Schematics of multiple Andreev reflection process in ballistic SNS
junctions.

The process of multiple Andreev reflection is schematically illustrated in
Fig. 2. After each traverse across the junction a quasiparticle (hole) gets
accelerated by $eV$, suffers Andreev reflection at one of the two NS
interfaces and eventually leaves the junction after $m\approx 2\Delta/eV$ such
traverses. As a result, the distribution functions $f^{+}(E)$ and $f^{-}(E)$
for quasiparticles moving respectively in and opposite to the current
direction deviate strongly from the equilibrium distribution function
$f_{F}(E)$. They read (see, e.g., Ref. Uwe2, )

$f^{\pm}(E)=\sum_{m=0}^{\infty}f_{F}(E\mp meV)[1-{\mathcal{A}}(E\mp
meV)]\prod_{l=0}^{m-1}{\mathcal{A}}(E\mp leV),$

where ${\mathcal{A}}(E)$ is the Andreev reflection probability. At
$T,eV\ll\Delta$ one arrives at the $I-V$ curve in the form Uwe

$\bar{I}=\frac{V}{R_{N}}+\frac{2\Delta}{eR_{N}}\operatorname{sgn}V,$ (5)

where $\operatorname{sgn}x$ equals to 1, 0 and -1 respectively for $x>0$,
$x=0$ and $x<0$. The last term in Eq. (5) represents an excess current due to
MAR and holds at arbitrary values of $d$. In the absence of inelastic
relaxation this fairly large current sets in at any non-zero applied voltage
$V$.

In the limit of short SNS junctions one also finds AB1

$I(t)=\frac{V}{R_{N}}+I_{c}|\sin(\varphi/2)|\operatorname{sgn}V,\quad\varphi=2eVt,$
(6)

where $I_{c}=\pi\Delta/(eR_{N})$ is the junction critical current at $T\to 0$.
We observe that CPR (6) remains strictly $2\pi$-periodic and no
$4\pi$-periodic supercurrent component occurs at non-zero voltages. Averaging
Eq. (6) over time, we again recover Eq. (5).

As we already discussed, in the presence of external microwave radiation with
frequency $\omega$ the junction $I-V$ curve exhibits Shapiro steps which occur
under the condition AB2 ; Uwe2 ; Cuevas  $2keV=n\omega$ or, equivalently, at

$\omega=\frac{k}{n}\omega_{J},$ (7)

where $k$ and $n$ are positive integer numbers. For clarity, in what follows
we will distinguish integer and fractional Shapiro steps corresponding to
respectively integer and non-integer values of the ratio $k/n$ in Eq. (7). The
values $k\geq 2$ (irrelevant in the tunneling limit) reflect the presence of
higher harmonics of the Josephson current due to non-sinusoidal CPR, whereas
the numbers $n\geq 2$ correspond to multi-photon processes which gain
importance at higher ac signals or smaller microwave frequencies.

Hence, fractional Shapiro steps are routinely expected even in topologically
trivial Josephson junctions and the observation of such steps (e.g., at
$\omega=\omega_{J}/2$) alone is by no means sufficient to make any definite
conclusion about the presence of a $4\pi$-periodic in $\varphi$ contribution
to CPR. In addition to fractional Shapiro steps, in experiments M16 ; topmat ;
Gre  missing integer Shapiro steps at $\omega=\omega_{J}$ were reported at low
enough frequencies whereas at higher $\omega$ these steps were clearly
detected. With the aid of a phenomenological model it was argued 4pisteps
that such observations actually support the $4\pi$-periodic Josephson effect
scenario. We will return back to this issue below in Section VIII.

The structure of the paper is as follows. In Section II we outline our general
formalism which is then employed in Section III in order to evaluate electric
current across ballistic SNS junctions for a arbitrary time dynamics of the
Josephson phase. In Section IV we reconstruct a complete description of ac
Josephson effect in SNS junctions biased by a constant external voltage. The
effect of an external ac signal on both CPR and the Josephson critical current
is discussed in Section V. Sections VI and VII are devoted to the analysis of
both integer and fractional Shapiro steps respectively in the voltage- and
current-biased regimes. Discussion of our key observations is presented in
Section VIII. Some technical details of our calculation are relegated to
Appendices A, B, C and D.

##  II Basic formalism

In what follows we will consider a short $SNS$ junction with $d\ll\xi_{0}$.
This system can be conveniently described within the standard technique
combining quasiclassical Eilenberger equations Eil ; ZRev  with Zaitsev
boundary conditions Zai  matching quasiclassical propagators on both sides of
the junction. Making use of this approach it is straightforward to derive a
general expression for the effective action Z94  that accounts for arbitrary
configurations of the Josephson phase $\varphi(t)$ and holds for any junction
transmission distribution.

As we are aiming at describing non-stationary and non-equilibrium processes it
is necessary to resort to the Keldysh technique and introduce two phase
variables $\varphi_{1}(t)$ and $\varphi_{2}(t)$ defined respectively on the
forward and backward parts of Keldysh contour. Then the Keldysh version of the
effective action reads SN ; book

$iS_{t}[\varphi]=\frac{1}{2}\sum_{n}{\rm
Tr}\,\ln\left[1+\frac{{\mathcal{T}}_{n}}{4}\left(\left\\{\check{Q}_{L}(\varphi),\check{Q}_{R}\right\\}-2\right)\right].$
(8)

Here the summation runs over the junction conducting channels with arbitrary
transmission coefficients ${\mathcal{T}}_{n}$, $\check{Q}_{L,R}$ are $4\times
4$ Green-Keldysh matrices for the left and right superconducting electrodes.
The product of these matrices implies time convolution and curly brackets
denote anticommutation. These matrices are expressed through quasiclassical
propagators as

$\check{Q}_{R}=\check{g}=\left(\begin{array}[]{cc}\hat{g}^{R}&\hat{g}^{K}\\\
0&\hat{g}^{A}\end{array}\right)$ (9)

and

$\check{Q}_{L}(\varphi)(t,t^{\prime})=\check{\cal L}\check{\cal
M}_{+}(t)\check{\cal L}\check{g}(t,t^{\prime})\check{\cal L}\check{\cal
M}_{-}(t^{\prime})\check{\cal L}.$ (10)

The retarded, advanced and Keldysh components of the matrix (18) are in turn
$2\times 2$ matrices in the Nambu space

$\hat{g}^{R,A,K}(\epsilon)=\left(\begin{array}[]{cc}g^{R,A,K}(\epsilon)&f^{R,A,K}(\epsilon)\\\
-f^{R,A,K}(\epsilon)&-g^{R,A,K}(\epsilon)\end{array}\right),$ (11)

and we also define the matrices

$\check{\cal M}_{\pm}=\left(\begin{array}[]{cc}\exp\left[\pm
i\varphi_{1}(t)\hat{\tau}_{3}/2\right]&0\\\ 0&\exp\left[\pm
i\varphi_{2}(t)\hat{\tau}_{3}/2\right]\end{array}\right),$ (12)
$\displaystyle\check{\cal
L}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}\hat{1}&\hat{1}\\\
\hat{1}&-\hat{1}\end{array}\right),\quad\hat{\tau}_{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right).$ (17)

The retarded and advanced propagators obey the following conditions:

$\displaystyle\left(g^{R}(\epsilon)\right)^{2}-\left(f^{R}(\epsilon)\right)^{2}=1,\quad\lim_{\epsilon\rightarrow\pm\infty}g^{R}(\epsilon)=1,$
(18) $\displaystyle g^{R}(-\epsilon)=(g^{R}(\epsilon))^{*},\quad
f^{R}(-\epsilon)=-(f^{R}(\epsilon))^{*},$ (19) $\displaystyle
g^{A}(\epsilon)=-\left(g^{R}(\epsilon)\right)^{*},\quad
f^{A}(\epsilon)=-\left(f^{R}(\epsilon)\right)^{*},$ (20)

whereas the Keldysh components read

$\displaystyle g^{K}(\epsilon)=2\,{\rm
Re}\left[g^{R}(\epsilon)\right]\tanh\frac{\epsilon}{2T},$ (21)
$\displaystyle\quad f^{K}(\epsilon)=2\,{\rm
Re}\left[f^{R}(\epsilon)\right]\tanh\frac{\epsilon}{2T}.$ (22)

In the case of conventional superconducting electrodes we have

$g^{R,A}(\epsilon)=\frac{\epsilon\pm i\gamma}{\xi^{R,A}(\epsilon)},\quad
f^{R,A}(\epsilon)=\frac{\Delta}{\xi^{R,A}(\epsilon)},$ (23)

where $\xi^{R,A}(\epsilon)=\pm\sqrt{(\epsilon\pm i\gamma)^{2}-\Delta^{2}}$ and
$\gamma$ controls the strength of inelastic relaxation. In the absence of
inelastic relaxation it is necessary to keep $\gamma$ infinitesimally small.

In order to proceed let us introduce the ”classical” and ”quantum” phases in
as standard manner as $\varphi_{+}=(\varphi_{1}+\varphi_{2})/2$ and
$\varphi_{-}=\varphi_{1}-\varphi_{2}$, respectively. A general expression for
the current operator expectation value can be defined in terms of the path
integral

$\langle\hat{I}(t)\rangle=2ie\int{\cal
D}\varphi_{\pm}\frac{\delta}{\delta\varphi_{-}(t)}e^{iS_{c}[\varphi_{\pm}]+iS_{t}[\varphi_{\pm}]}.$
(24)

which should be evaluated under the conditions
$\langle\dot{\varphi}_{+}(t)\rangle=2eV(t)$ and
$\langle\dot{\varphi}_{-}(t)\rangle=0$, where $V(t)$ is the voltage across the
junction. The expression (24) can also be employed in the presence of
electron-electron interactions in which case the charging term $S_{c}$ should
be added to the effective action $S_{t}$ (8).

Perhaps we can add that the above effective action formalism is equally
convenient for the analysis of higher cumulants of the current operator. The
current-current correlation function as well as higher current correlators are
obtained by applying consecutive derivatives $\delta/\delta\varphi_{-}(t)$
under the path integral (24). In the particular case of relatively small phase
fluctuations this procedure was implemented in Ref. GZ, .

Here the ”classical” phase variable
$\varphi_{+}(t)=2e\int^{t}dt^{\prime}V(t^{\prime})$ cannot anymore be regarded
as small which makes the whole problem rather complicated to deal with. On the
other hand, in the absence of electron-electron interactions the path integral
in Eq. (24) is evaluated trivially and we obtain

$I(t)=\left.-2e\frac{\delta
S_{t}[\varphi_{\pm}]}{\delta\varphi_{-}(t)}\right|_{\varphi_{-}=0,\;\varphi_{+}(t)\equiv\varphi(t)=2e\int^{t}dt^{\prime}V(t^{\prime})}.$
(25)

From now on we will make no distinction between the phase variables
$\varphi_{+}(t)$ and $\varphi(t)$.

##  III Electric current for ballistic junctions

In the interesting for us limit ${\mathcal{T}}_{n}=1$ the action (8) becomes
simpler. With the aid of the normalization condition
$\check{Q}_{L}^{2}(t,t^{\prime})=\check{Q}_{R}^{2}(t,t^{\prime})=\delta(t-t^{\prime})$,
in this limit we have

$\displaystyle iS_{t}[\varphi]=\frac{1}{2}\sum_{n}{\rm
Tr}\,\ln\left[\frac{1}{4}\left(\check{Q}_{L}(\varphi)+\check{Q}_{R}\right)^{2}\right]$
$\displaystyle={\mathcal{N}}{\rm
Tr}\,\ln\left[\frac{1}{2}\check{I}\left(\check{Q}_{L}(\varphi)+\check{Q}_{R}\right)\right],$
(26)

where the matrix $\check{I}$ is defined as

$\check{I}=\left(\begin{array}[]{cc}\hat{\tau}_{3}&0\\\
0&-\hat{\tau}_{3}\end{array}\right).$ (27)

Expanding the argument of the logarithm in Eq. (26) up to the first order in
$\varphi_{-}(t)$ and symmetrizing the resulting expressions with respect to
the phase variables, we get

$iS_{t}={\mathcal{N}}{\rm Tr}\,\ln\left[\check{Q}_{0}+\check{Q}_{1}\right].$
(28)

Here the matrices $\check{Q}_{0}$ and $\check{Q}_{1}$ read

$\displaystyle\check{Q}_{0}=\left(\begin{array}[]{cc}\hat{a}^{R}&\hat{a}^{K}\\\
0&-\hat{a}^{A}\end{array}\right),\quad{\rm
Tr}\,\ln\left[\check{Q}_{0}\right]=0,$ (31)
$\displaystyle\check{Q}_{1}(t,t^{\prime})=\frac{\varphi_{-}(t)}{8}\left(\begin{array}[]{cc}0&-\hat{b}^{A}(t,t^{\prime})\\\
\hat{b}^{R}(t,t^{\prime})&\hat{b}^{K}(t,t^{\prime})\end{array}\right)$ (34)
$\displaystyle+\left(\begin{array}[]{cc}\hat{\tau}_{3}\hat{b}^{K}(t,t^{\prime})\hat{\tau}_{3}&\hat{\tau}_{3}\hat{b}^{R}(t,t^{\prime})\hat{\tau}_{3}\\\
-\hat{\tau}_{3}\hat{b}^{A}(t,t^{\prime})\hat{\tau}_{3}&0\end{array}\right)\frac{\varphi_{-}(t^{\prime})}{8},$
(37)

where we define

$\hat{a}^{R,A,K}(t,t^{\prime})=\left(\begin{array}[]{cc}g^{R,A,K}(t,t^{\prime})\cos\left[\frac{\varphi(t)-\varphi(t^{\prime})}{4}\right]&f^{R,A,K}(t,t^{\prime})\cos\left[\frac{\varphi(t)+\varphi(t^{\prime})}{4}\right]\\\
f^{R,A,K}(t,t^{\prime})\cos\left[\frac{\varphi(t)+\varphi(t^{\prime})}{4}\right]&g^{R,A,K}(t,t^{\prime})\cos\left[\frac{\varphi(t)-\varphi(t^{\prime})}{4}\right]\end{array}\right),$
(38)
$\hat{b}^{R,A,K}(t,t^{\prime})=\left(\begin{array}[]{cc}g^{R,A,K}(t,t^{\prime})\sin\left[\frac{\varphi(t)-\varphi(t^{\prime})}{4}\right]&f^{R,A,K}(t,t^{\prime})\sin\left[\frac{\varphi(t)+\varphi(t^{\prime})}{4}\right]\\\
f^{R,A,K}(t,t^{\prime})\sin\left[\frac{\varphi(t)+\varphi(t^{\prime})}{4}\right]&g^{R,A,K}(t,t^{\prime})\sin\left[\frac{\varphi(t)-\varphi(t^{\prime})}{4}\right]\end{array}\right).$
(39)

Combining Eq. (25) with Eqs. (28)-(37), we arrive at the general expression
for the current

$\displaystyle I(t)=\frac{ie{\mathcal{N}}}{4}\int dt^{\prime}{\rm
Tr}\left[\hat{b}^{R}(t,t^{\prime})\hat{X}^{K}(t^{\prime},t)+\hat{b}^{K}(t,t^{\prime})\hat{X}^{A}(t^{\prime},t)\right.$
$\displaystyle\left.-\hat{X}^{K}(t,t^{\prime})\hat{\tau}_{3}\hat{b}^{A}(t^{\prime},t)\hat{\tau}_{3}+\hat{\tau}_{3}\hat{X}^{R}(t,t^{\prime})\hat{\tau}_{3}\hat{b}^{K}(t^{\prime},t)\right],$
(40)

where

$\displaystyle\hat{X}^{R}=\left(\hat{a}^{R}\right)^{-1},\quad\hat{X}^{A}=-\left(\hat{a}^{A}\right)^{-1}$
(41)

and

$\displaystyle\hat{X}^{K}=-\hat{X}^{R}\circ\hat{a}^{K}\circ\hat{X}^{A}$ (42)

are the matrix elements of the inverse matrix

$\displaystyle\check{Q}_{0}^{-1}\equiv\check{X}=\left(\begin{array}[]{cc}\hat{X}^{R}&\hat{X}^{K}\\\
0&\hat{X}^{A}\end{array}\right).$ (45)

The expression (40) has an explicit causal nature and remains valid for an
arbitrary dependence of the applied voltage $V(t)$ on time. The general result
(40) coincides with those derived previously Zaitsev80 ; Zaikin83  by directly
solving the Eilenberger equations.

##  IV Time-independent bias voltage

Let us first consider the limiting case of a constant bias voltage $V$ applied
directly to the junction. Obviously, the Josephson phase then depends linearly
on time, i.e. $\varphi(t)=2eVt$. As it was demonstrated in Ref. Uwe, , in this
particular case it is possible to explicitly invert the matrices in Eq. (41)
and recover the exact expression for the current $I(t)$ across the junction
Uwe ; AB1 .

As compared to Ref. Uwe, , here we pursue a different approach outlined in
Appendix A. Introducing the notation

$a^{R}(\epsilon)=\frac{f^{R}(\epsilon)}{1+g^{R}(\epsilon)},$ (46)

which has to do with the so-called Riccati parametrization (see, e.g., Ref.
ZRev, )

$f^{R}=\frac{2a^{R}}{1-(a^{R})^{2}},\quad
g^{R}=\frac{1+(a^{R})^{2}}{1-(a^{R})^{2}},$ (47)

and employing the (corresponding to MAR) multiplicative structure of the
resulting expressions (103), from Eq. (40) we arrive at the final result

$I(t)=\sum_{l=-\infty}^{\infty}I_{l}e^{-2ielVt},$ (48)

where

$\displaystyle
I_{l=0}\equiv\bar{I}=\frac{V}{R_{N}}-\frac{1}{eR_{N}}\sum_{n=1}^{\infty}\int\limits_{-\infty}^{\infty}d\epsilon\tanh\frac{\epsilon}{2T}\left(1-\left|a^{R}(\epsilon)\right|^{2}\right)\prod_{1\leq
m\leq n}\left|a^{R}(\epsilon+meV)\right|^{2},$ (49) $\displaystyle
I_{l>0}=-\frac{1}{eR_{N}}\sum_{n=1}^{\infty}\int\limits_{-\infty}^{\infty}d\epsilon\tanh\frac{\epsilon}{2T}\left(1-\left|a^{R}(\epsilon)\right|^{2}\right)\prod_{1\leq
m\leq n}\left|a^{R}(\epsilon+meV)\right|^{2}\prod_{n+1\leq k\leq
n+2l}a^{R}(\epsilon+keV)$
$\displaystyle-\frac{1}{eR_{N}}\int\limits_{-\infty}^{\infty}d\epsilon\tanh\frac{\epsilon}{2T}\left(1-\left|a^{R}(\epsilon)\right|^{2}\right)\prod_{1\leq
k\leq 2l}a^{R}(\epsilon+keV),\quad I_{-l}=I_{l}^{*}.$ (50)

The results (48)-(50) fully conform to those of Refs. Uwe, ; AB1, . As
expected, the current $I(t)$ (48) is strictly $2\pi$-periodic in
$\varphi(t)=2eVt$ for any nonzero $V$ and no $4\pi$-periodic component occurs.

Equation (49) provides a general expression for the average current $\bar{I}$.
As we already discussed, it is characterized by the current jump at $V=0$
caused by MAR. This excess current also persists at all larger voltages.
Provided inelastic relaxation effects remain weak and can be neglected,
combining Eqs. (49) and (23) we obtain Uwe

$\bar{I}=\frac{V}{R_{N}}+\frac{2\Delta(T)}{eR_{N}}\tanh\frac{\Delta(T)}{2T}\operatorname{sgn}V,\quad
eV\ll\Delta(T)$ (51)

and Zaitsev80 ; Zaikin83 ; MAR

$\bar{I}=\frac{V}{R_{N}}+\frac{8\Delta(T)}{3eR_{N}}\tanh\frac{eV}{2T},\quad
eV\gg\Delta(T).$ (52)

Temperature effects do not eliminate the current jump although diminish its
amplitude. However, non-vanishing inelastic relaxation, if present, smears
this jump making the linear junction conductance $\bar{I}/V$ finite in the
limit $V\to 0$ Uwe .

At small voltages and temperatures $eV,T\ll\Delta$ and for weak inelastic
relaxation Eqs. (48)-(50) can be combined reducing to a simple expression for
the total current defined in Eq. (6). We will explicitly make use of this
result further below. Finally, at high voltages $eV\gg\Delta$ and low $T$ one
readily finds AB1

$I(t)\simeq\bar{I}-\frac{\pi\Delta^{2}\ln 2}{2e^{2}VR_{N}}\cos(2eVt).$ (53)

##  V Supercurrent under ac signal

As a next step we consider a somewhat different physical situation assuming
now that our SNS junction is exposed to external microwave radiation which
generates an ac voltage signal

$V(t)=V_{\rm ac}\cos(\omega t)$ (54)

across the junction. Here $V_{\rm ac}$ scales with the intensity of external
radiation and $\omega$ is the radiation frequency. The Josephson phase then
oscillates in time and takes the form

$\varphi(t)=\bar{\varphi}+2\alpha\sin(\omega t),\quad\alpha=eV_{\rm
ac}/\omega.$ (55)

In this state with zero average voltage $\bar{V}$ the junction can still carry
a non-vanishing dc supercurrent which is now affected by external radiation.
In order to evaluate this supercurrent one can employ the general formula (40)
combined with Eq. (55). However, in contrast to the case of a constant in time
bias voltage considered above, here it is not possible to exactly invert the
corresponding matrices in Eq. (40) and one should resort to certain
approximations. Physically the main complication is that in the presence of an
external ac field the quasiparticle distribution function in the contact area
is driven out of equilibrium an in a general case can only be evaluated
numerically by resolving a self-consistent quantum kinetic equation which
simultaneously accounts for Andreev reflection as well as photon absorption
and emission processes. The corresponding analysis was carried out in Ref.
Uwe2, . The influence of external radiation on both CPR and the critical
Josephson current in superconducting point contacts was previously studied in
Ref. Bergeret, .

A substantial simplification of the problem can be achieved if the
supercurrent is evaluated within the adiabatic approximation which we will
employ here. Assuming that the frequency $\omega$ is small enough to obey the
condition

$\omega\ll 2|E^{A}_{\pm}(\bar{\varphi})|,$ (56)

and neglecting (typically rather small) geometric capacitance of our junction
one finds the time-dependent current $I(t)$ which reads GZ

$I(t)\simeq I_{S}(\varphi(t))+C^{*}(\bar{\varphi})\ddot{\varphi}(t),$ (57)

where $I_{S}(\varphi)$ is defined in Eq. (3) and $C^{*}(\bar{\varphi})$
represents the renormalized junction capacitance. In the low temperature limit
it takes the form GZ

$C^{*}=\frac{\pi}{16\Delta R_{N}\cos^{4}(\bar{\varphi}/4)}$ (58)

outside an immediate vicinity of the point $\bar{\varphi}=\pi$ and formally
diverges at $\bar{\varphi}\to\pi$, thus signaling the failure of the adiabatic
approximation for $E^{A}_{\pm}(\pi)=0$.

Substituting the phase $\varphi(t)$ in the form (55) into Eq. (57), making use
of Eq. (3) and averaging the resulting expression over time, in the limit
$T\to 0$ we obtain an obvious relation (cf. also Bergeret )

$\bar{I}(\bar{\varphi})=\frac{8\Delta}{eR_{N}}\sum_{m=1}^{\infty}\frac{m(-1)^{m+1}}{4m^{2}-1}J_{0}(2m\alpha)\sin(m\bar{\varphi}),$
(59)

where $J_{0}(x)$ is the zero-order Bessel function. For $\alpha<\pi/2$ there
exists the phase interval $2\alpha-\pi\leq\bar{\varphi}\leq\pi-2\alpha$ where
the series (59) can be summed up exactly with the result

$\bar{I}(\bar{\varphi})=I_{c}J_{0}(\alpha)\sin(\bar{\varphi}/2).$ (60)

Outside this interval CPR $\bar{I}(\bar{\varphi})$ (59) can be evaluated
numerically. The corresponding dependencies are illustrated in Figs. 3-5 for
different values of $\alpha$.

Figure 3: Zero temperature current-phase relation $\bar{I}(\bar{\varphi})$
(normalized by $I_{c}\equiv I_{c}(0)$) for $\alpha=0$ (magenta), 0.6 (orange)
and 1.2 (blue).

We observe that Eq. (60) is fully confirmed under its validity conditions,
whereas beyond them CPR develops qualitatively new features including, e.g.,
negative current states within the interval $0<\bar{\varphi}<\pi$. For
$\alpha\lesssim 1.345$ CPR reaches its absolute maximum at
$\bar{\varphi}=\pi-2\alpha$ in which case the junction critical current
$I_{c}(\alpha)\equiv{\rm max}|\bar{I}(\bar{\varphi})|$ reads

$I_{c}(\alpha)=\bar{I}(\pi-2\alpha)=I_{c}J_{0}(\alpha)\cos\alpha.$ (61)

For $\alpha\gtrsim 1.345$ and $0<\bar{\varphi}<\pi$ the absolute maximum of
$|\bar{I}(\bar{\varphi})|$ is reached at negative current values (cf. Fig. 4)
and the system switches to the $\pi$-junction state which persists up to
$\alpha\simeq 2.9$ where $I_{c}(\alpha)$ achieves its second local minimum
(see Fig. 6). In this regime for $1.7\lesssim\alpha\lesssim 2.7$ the part of
CPR with $d\bar{I}/d\bar{\varphi}>0$ is well approximated by the dependence
$\bar{I}/I_{c}=-a\cos(\varphi/2)$ with $a\simeq 0.46$. At $\alpha\gtrsim 2.9$
CPR again describes a 0-junction state (cf. Fig. 5).

Figure 4: The same as in Fig. 3 for $\alpha=1.345$ (magenta), 1.5 (orange) and
2.0 (blue)(blue).

In the limit of large values of $\alpha\gg 1$ CPR can be approximately
described by keeping only the first term ($m=1$) in the series (59) which
takes the form

$\bar{I}(\bar{\varphi})\simeq\frac{8I_{c}}{3\pi^{3/2}\sqrt{\alpha}}\sin\bar{\varphi}\cos\left(2\alpha-\frac{\pi}{4}\right).$
(62)

This dependence works reasonably well except in the immediate vicinity of the
points $\alpha=3\pi/8+\pi p/2$ ($p=0,\pm 1,\pm 2,...$) where terms with $m>1$
in Eq. (59) need to be retained as well.

Figure 5: The same as in Fig. 3 for $\alpha=2.8$ (magenta), 3.1 (orange) and
3.6 (blue).

The dependence of the critical current on $\alpha$ at $T\to 0$ is displayed in
Fig. 6. We observe that – in contrast to the case of Josephson tunnel
junctions BP  – $I_{c}(\alpha)$ does not vanish at any finite value of
$\alpha$, it demonstrates oscillations which decay with increasing $\alpha$ as
$I_{c}(\alpha)\propto 1/\sqrt{\alpha}$. Local minima of $I_{c}(\alpha)$
correspond to consecutive transitions between 0- and $\pi$-junction states.
Note that such transitions take place at non-zero values of the critical
current. A similar feature was recently predicted KDZ20 ; KZ21  for
$X$-junctions driven out of equilibrium by applying a temperature gradient.
Hence, this feature is likely to be generic for junctions with non-sinusoidal
CPR. We also note that the behavior of the critical current $I_{c}(\alpha)$
displayed in Fig. 6 is fully consistent with that reported previously Bergeret
.

Figure 6: Zero temperature critical current $I_{c}(\alpha)$ normalized by
$I_{c}\equiv I_{c}(0)$. The current $I_{c}(\alpha)$ remains non-zero at any
finite value of $\alpha$. Local minima of this dependence correspond to
concequtive transitions between 0- and $\pi$-junction states.

In order to extend our analysis beyond the adiabatic limit (56) it is
necessary to include resonances associated with excitation of the system from
the lower Andreev state to the upper one by absorption of $l$ photons
($l=1,2,...$) with frequency $\omega<2\Delta$. This resonant process boosts
the (negative) contribution to the supercurrent generated by the upper Andreev
level and, hence, causes narrow dips on CPR Bergeret  at the phase values
$\bar{\varphi}=\varphi_{l}$ ($l=1,2,...$) obeying the condition
$l\omega=|E^{A}_{+}(\varphi_{l})-E^{A}_{-}(\varphi_{l})|$. Obviously, these
dips are not captured within our simple adiabatic approximation. They,
however, do not lead to any qualitative modifications of the dependence of the
critical current on $\alpha$ in Fig. 6. At even higher frequencies
$\omega\gtrsim 2\Delta$ the critical current gets further decreased due to
pairbreaking effects.

Finally, we stress that the above results apply in the interesting for us low
temperature limit $T\ll\Delta$. With increasing temperature CPR quickly
approaches the sinusoidal form and the adiabatic approximation yields
vanishing critical current at the values of $\alpha$ where the Bessel function
$J_{0}(2\alpha)$ has its zeros. On the other hand, at not too low temperatures
non-equilibrium effects caused by external radiation may also stimulate the
supercurrent Zaikin83 , i.e. under certain conditions the critical current in
SNS junctions may actually increase in the presence of an ac signal Zaikin83 ;
Bergeret .

##  VI Shapiro steps

Let us now assume that our SNS junction is biased by an external voltage
$V(t)$ which contains both constant in time and ac parts, i.e.

$V(t)=V+V_{\rm ac}\cos(\omega t+\vartheta),$ (63)

where $\vartheta$ is an arbitrary phase. Accordingly, for the Josephson phase
we now get

$\varphi(t)=2e\int^{t}dt^{\prime}V(t^{\prime})=2eVt+2\alpha\sin(\omega
t+\vartheta).$ (64)

As before, here the parameter $\alpha$ (55) effectively controls the strength
of microwave radiation effects.

The overall effect of external microwave radiation on the $I-V$ curve in fully
transparent SNS junctions was addressed in Ref. Uwe2, . The physics behind
this effect is transparent: While suffering MAR inside the junction
quasiparticles and holes may now also absorb and emit photons with frequency
$\omega$. As a result of such processes, the quasiparticle energy as well as
its distribution function may change during the MAR cycle which, in turn,
yields three major consequences Uwe2 : (i) Reduction of zero bias conductance
and excess current (particularly pronounced at smaller bias voltages
$eV\lesssim 2\Delta$), (ii) substantial modification of the subharmonic gap
structure on the $I-V$ curve and (iii) appearance of Shapiro steps at
microwave radiation frequencies (7).

Despite some previous efforts AB2 ; Uwe2 ; Cuevas  a detailed microscopic
theory of Shapiro steps in ballistic SNS junctions is not yet completed. Below
we will employ our formalism in order to achieve some progress in this
direction.

###  VI.1 Small microwave signals

We first assume that this parameter is small, $\alpha\ll 1$, and evaluate the
effect of external radiation on the $I-V$ curve perturbatively in $\alpha$.
The linear in $\alpha$ corrections $\delta\hat{X}^{R,A,K}$ to (already
evaluated in the previous section) matrices $\hat{X}^{R,A,K}$ read

$\displaystyle\delta\hat{X}^{R}=-\hat{X}^{R}\delta\hat{a}^{R}\hat{X}^{R},\;\delta\hat{X}^{A}=\hat{X}^{A}\delta\hat{a}^{A}\hat{X}^{A},$
(65)
$\displaystyle\delta\hat{X}^{K}=-\hat{X}^{R}\delta\hat{a}^{K}\hat{X}^{A}-\hat{X}^{R}\delta\hat{a}^{R}\hat{X}^{K}+\hat{X}^{K}\delta\hat{a}^{A}\hat{X}^{A}.$
(66)

We also expand the matrices $\hat{a}^{R,A,K}$ and $\hat{b}^{R,A,K}$ (cf. Eqs.
(38) and (39)) up to the linear in $\alpha$ terms

$\displaystyle\delta\hat{a}^{R}(t,t^{\prime})=-\frac{\alpha}{2}\left(\begin{array}[]{cc}g^{R}(t,t^{\prime})\sin\frac{eV(t-t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)-\sin(\omega
t^{\prime}+\vartheta)\right)&f^{R}(t,t^{\prime})\sin\frac{eV(t+t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)+\sin(\omega t^{\prime}+\vartheta)\right)\\\
f^{R}(t,t^{\prime})\sin\frac{eV(t+t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)+\sin(\omega
t^{\prime}+\vartheta)\right)&g^{R}(t,t^{\prime})\sin\frac{eV(t-t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)-\sin(\omega t^{\prime}+\vartheta)\right)\end{array}\right),$ (69)
$\displaystyle\delta\hat{b}^{R}=\frac{\alpha}{2}\left(\begin{array}[]{cc}g^{R}(t,t^{\prime})\cos\frac{eV(t-t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)-\sin(\omega
t^{\prime}+\vartheta)\right)&f^{R}(t,t^{\prime})\cos\frac{eV(t+t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)+\sin(\omega t^{\prime}+\vartheta)\right)\\\
f^{R}(t,t^{\prime})\cos\frac{eV(t+t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)+\sin(\omega
t^{\prime}+\vartheta)\right)&g^{R}(t,t^{\prime})\cos\frac{eV(t-t^{\prime})}{2}\left(\sin(\omega
t+\vartheta)-\sin(\omega t^{\prime}+\vartheta)\right)\end{array}\right)$ (72)

and similarly for $a^{A,K}$ and $b^{A,K}$. Combining all these expressions
with Eq. (40) we arrive at the linear in $\alpha$ correction to the $I-V$
curve of our voltage-biased SNS junction due to the presence of external
radiation. It is easy to observe that this correction differs from zero only
under the condition $\omega=k\omega_{J}$ with $k=1,2,...$ giving rise to
integer $\vartheta$-dependent Shapiro steps $\delta I_{k}(\vartheta)$.

Explicit analytic expressions for $\delta I_{k}$ are specified in Appendix B,
cf. Eq. (105). They take the form

$\delta I_{k}(\vartheta)=a_{k}\cos\vartheta+b_{k}\sin\vartheta,$ (73)

implying that the overall magnitude of the $k$-th Shapiro step is $\delta
I_{k}=2\sqrt{a_{k}^{2}+b_{k}^{2}}$. Employing Eq. (105) one can evaluate
$\delta I_{k}$ numerically for different values of radiation frequency
$\omega$, temperature $T$ and effective inelastic relaxation rate $\gamma$.

Figure 7: The magnitude of integer Shapiro steps $\delta I_{k}$ (normalized by
$I_{0}=\alpha\Delta/(2eR_{N})$) evaluated from Eq. (105) for
$\omega/\Delta=0.5$, $\gamma/\Delta=0.01$ and different values of
$T/\Delta=0.1$ (diamonds), 1 (squares) and 3 (circles). Inset: The same for
$\omega/\Delta=0.5$, $T/\Delta=0.1$ and $\gamma/\Delta=0.2$. Figure 8:  The
same as in Fig. 7 for $T/\Delta=0.1$, $\gamma/\Delta=0.01$ and different
values of $\omega/\Delta=0.1$ (squares), 0.3 (diamonds) and 1.2 (circles). The
steps indicated by diamonds are shifted by 0.4 upwards to avoid overlapping.

Figure 7 demonstrates the dependence of integer Shapiro steps $\delta I_{k}$
on $k$ at three different values of $T$ and the external radiation frequency
$\omega=0.5\Delta$ in the limit of very weak inelastic relaxation
$\gamma=0.01\Delta$ in which case the phase coherence remains well preserved.
We observe that $\delta I_{k}$ demonstrates decaying oscillations with
increasing $k$. In addition, there is an overall decrease of $\delta I_{k}$
with increasing $T$ at any given $k$ due to temperature smearing. Note that
such oscillations disappear completely even at low values of $T$ as soon as
inelastic relaxation becomes more pronounced. This effect is illustrated in
the inset of Fig. 7 demonstrating $\delta I_{k}$ decaying monotonously with
increasing $k$.

Figure 8 displays the dependence of the magnitude of Shapiro steps on $k$ at
different values of $\omega$ for $\gamma=0.01\Delta$ and $T=0.1\Delta$. We
again observe decaying oscillations of $\delta I_{k}$ with increasing $k$.

###  VI.2 Beyond perturbation theory

The above perturbative in $\alpha$ procedure is sufficient to evaluate all
integer Shapiro steps with $\omega=k\omega_{J}$, whereas in order to recover
fractional Shapiro steps (7) with $n>1$ corresponding to multi-photon
processes it is necessary to proceed to higher orders in $\alpha$.
Accordingly, for $\alpha\ll 1$ such fractional steps remain parametrically
smaller than the integer ones and, hence, can be neglected in this limit.
Note, however, that for small values of the bias voltage $V$ the parameter
$\alpha=eV_{\rm ac}/\omega\sim V_{\rm ac}/V$ is in general not at all small,
thus making any perturbation theory in $\alpha$ insufficient. The task at hand
is to go beyond such perturbation theory and evaluate both integer and
fractional Shapiro steps for arbitrary values of $\alpha$.

In order to accomplish this goal one again has to compute the inverse matrices
$\hat{X}^{R,A}=(\hat{a}^{R,A})^{-1}$ and substitute the result into Eq. (40)
for the current. The procedure in the presence of an ac signal (63) is
outlined in Appendix C. Here, however, we take a somewhat different route. We
will first give a simple estimate for fractional Shapiro steps which turns out
to be sufficient for $\alpha\lesssim 1$. Then we will specifically address the
interesting for us low bias regime in which case one can develop a more
accurate calculation that will be applicable for all values of $\alpha$.

Our simple estimate for the magnitude of (fractional) Shapiro steps $\delta
I_{\frac{1}{n}}$ which occur at $\omega=\omega_{J}/n$ is outlined in Appendix
D. It yields

$\frac{\delta
I_{\frac{1}{n}}}{I_{c}}\sim\left\\{\begin{array}[]{l}J^{2}_{n/2}\left(\alpha/2\right)\;{\rm
for\,even}\;n,\\\
|J_{(n+1)/2}\left(\alpha/2\right)J_{(n-1)/2}\left(\alpha/2\right)|\;{\rm
for\,odd}\;n,\end{array}\right.$ (74)

where, as above, $J_{n}(x)$ are Bessel functions.

All these Shapiro steps correspond to $k=1$. The magnitude of other fractional
Shapiro steps $\delta I_{\frac{k}{n}}$ decays with growing $k$ similarly to
integer Shapiro steps considered above.

It is obvious from the above estimate that for $\alpha\sim 1$ the magnitude of
fractional Shapiro steps is in general of the same order as that for integer
Shapiro steps. Furthermore, the former can easily exceed the latter for some
values of $\alpha$. For instance, the ratio between the magnitudes of Shapiro
steps $\delta I_{1}$ and $\delta I_{\frac{1}{2}}$ which occur respectively at
$\omega=\omega_{J}$ and $\omega=\omega_{J}/2$ reads

$\frac{\delta I_{1}}{\delta
I_{\frac{1}{2}}}\sim\frac{|J_{0}(\alpha)|}{|J_{1}(\alpha)|}.$ (75)

Clearly, for $\alpha\ll 1$ the integer Shapiro step always dominates $\delta
I_{1}\gg\delta I_{\frac{1}{2}}$, whereas for $\alpha\sim 1$ we already have
$\delta I_{1}\sim\delta I_{\frac{1}{2}}$ or possibly even $\delta
I_{1}\ll\delta I_{\frac{1}{2}}$ since Bessel functions of different order have
different zeroes. Hence, for small enough voltages $V$ the fractional
(”$4\pi$”) Shapiro step can dominate of the integer (”$2\pi$”) one even though
the so-called ”$4\pi$-Josephson effect” is totally absent. In fact, this
observation should not be taken as any surprise because it is well known that
in the case of tunnel junctions one has Tinkh ; BP ; Likh  $\delta
I_{\frac{1}{n}}\propto J_{n}(2\alpha)$.

At sufficiently small values of the applied voltage $eV\ll\Delta$ it is also
possible to evaluate the magnitude of fractional Shapiro steps more accurately
by extending Eq. (6) to the time-dependent bias voltage in the form (63). The
accuracy of this approximation will be discussed below in the next section.
Substituting the phase $\varphi(t)$ (64) into Eq. (6) and averaging the
resulting current $I(t)$ over time we obtain

$\displaystyle\bar{I}=\frac{V}{R_{N}}+I_{c}\left\langle\left|\sin\left[z+\vartheta+\alpha\sin\left(\frac{\omega}{eV}z\right)\right]\right|\right.$
$\displaystyle\left.\times{\rm
sgn}\left[1+\frac{\alpha\omega}{eV}\cos\left(\frac{\omega}{eV}z\right)\right]\right\rangle_{z}$
(76)

Here we are interested in the $\vartheta$-dependent terms emerging from Eq.
(76) provided $\omega$ and $eV$ are commensurate. In this case the current
periodically depends on the phase $\vartheta$ reaching both maximum and
minimum values $I_{\rm max}$ and $I_{\rm min}$ within the period $\pi$. The
magnitude of a Shapiro step is then defined as $\delta I=I_{\rm max}-I_{\rm
min}$. Three different Shapiro steps $\delta I_{1}$, $\delta I_{\frac{1}{2}}$
and $\delta I_{\frac{1}{3}}$ corresponding respectively to $\omega=2eV$,
$\omega=eV$ and $\omega=2eV/3$ are displayed in Fig. 9 as functions of the
parameter $\alpha$.

Figure 9: Shapiro steps $\delta I_{1}(\alpha)$ (green), $\delta
I_{\frac{1}{2}}(\alpha)$ (orange) and $\delta I_{\frac{1}{3}}(\alpha)$ (blue)
derived from Eq. (76) respectively at $\omega=\omega_{J}$,
$\omega=\omega_{J}/2$ and $\omega=\omega_{J}/3$. For $\alpha>0$ these three
curves reach their first minima respectively at $\alpha\simeq$ 0.862, 1.407
and 1.943 where they take very small but positive values min$(\delta
I_{1})\simeq 0.006$, min$(\delta I_{\frac{1}{2}})\simeq 0.0053$ and
min$(\delta I_{\frac{1}{3}})\simeq 0.005$. Note that the above numbers as well
as the sharp peak feature observed in all these curves may to a certain extent
depend on the approximation employed here.

In accordance with our perturbative results we observe that at small values of
$\alpha$ the principal Shapiro step $\delta I_{1}\propto\alpha$ dominates over
the fractional ones $\delta I_{\frac{1}{n}}\propto\alpha^{n}$. On the other
hand, the peaks of all three curves at $\alpha\sim 1$ are comparable in
magnitude and each of the Shapiro steps can be significantly larger than the
other ones due to different positions of these peaks. For instance, at low
enough voltages $V$ within the interval $0.6\lesssim\alpha\lesssim 1.3$ the
fractional Shapiro step $\delta I_{\frac{1}{2}}$ dominates over $\delta I_{1}$
which almost vanishes in the middle part of this interval. On the other hand,
at higher voltages the condition $\omega=2eV$ is fulfilled al lower values of
$\alpha$ where the integer step $\delta I_{1}(\alpha)$ already becomes
significant. This behavior is qualitatively consistent with the one observed
in experiments M16 ; topmat ; Gre  but, of course, it by no means can serve as
an evidence for $4\pi$-periodic CPR in our junctions.

Comparing our simple estimate (74) with more accurate results following from
Eq. (76) we conclude that – in the agreement with our expectations (see
Appendix D) – the former works reasonably well for $\alpha\lesssim 1$ whereas
for larger values of $\alpha$ the fractional Shapiro steps are better
approximated, e.g., by the formula $\delta I_{\frac{1}{n}}\sim
I_{c}J_{n/2}^{2}(3\alpha)$, cf. Fig. 9.

##  VII Current-biased regime

Let us now consider a somewhat different physical situation assuming that our
SNS junction is biased by an external current in the form

$I(t)=I+I_{\rm ac}\cos(\omega t+\vartheta).$ (77)

Such current biased regime appears to be realized in a number of recent
experiments M16 ; topmat ; Gre . It is obvious that in this case the voltage
across the junction cannot anymore be constant and should vary in time
$V=V(t)$ even in the absence of an ac signal $I_{\rm ac}\to 0$.

Within the adiabatic approximation, i.e. provided the voltage $V(t)$ changes
in time slowly enough, it is still possible to employ the same CPR (6) and
(53) initially derived for the time independent voltage bias. This
approximation should work provided $V(t)$ does not change significantly during
the whole MAR cycle, i.e. during the period of time

$\delta t_{\rm mar}\sim\frac{2\Delta}{eV(t)}\frac{d}{v_{F}}$ (78)

it takes for a quasiparticle to cross the normal layer between two
superconductors $\sim 2\Delta/(eV(t))$ times and escape into one of the
superconducting terminals, see Fig. 2.

In fact, we have already made use of this approximation in the previous
section evaluating the $I-V$ curve and Shapiro steps in the presence of an ac
voltage signal (63) simply by substituting the Josephson phase in the form
(64) into CPR (6). In that case it suffices to require $\delta t_{\rm mar}\ll
1/\omega$ or, equivalently, $d\ll\xi_{0}eV(t)/\omega$. Here, we will also
assume that $V(t)$ changes slowly enough and again stick to the same adiabatic
approximation which we verify in the end of our calculation. Within this
approximation and provided $eV(t)$ remains much smaller than $\Delta$ the
phase and voltage dynamics is controlled by the following equation

$I+I_{\rm ac}\cos(\omega
t+\vartheta)=\dot{\varphi}/(2eR_{N})+I_{c}|\sin(\varphi/2)|\operatorname{sgn}\dot{\varphi},$
(79)

which describes the current balance in our SNS junction.

In the absence of an ac current component, i.e. in the limit $I_{\rm ac}\to
0$, the solution of this equation can be written in the form

$F(\varphi,a)=eI_{c}R_{N}t,\quad a=I/I_{c}>1,$ (80)

where

$\displaystyle
F(\varphi,a)=\frac{2}{\sqrt{a^{2}-1}}\left[\arctan\left(\frac{a\tan\left(\frac{\varphi}{4}\right)-{\rm
sgn}\left[\tan\left(\frac{\varphi}{4}\right)\right]}{\sqrt{a^{2}-1}}\right)\right.$
$\displaystyle\left.+\pi\chi_{1}+\left(2\arctan\frac{1}{\sqrt{a^{2}-1}}-\pi\right)\chi_{2}\right],$
(81)

where we defined $\chi_{n}={\rm floor}\left[\varphi/(2\pi n)\right]$ ($n=1,2$)
as integer parts of $\varphi/(2\pi n)$. This solution allows to determine the
oscillating in time voltage across the junction $V(t)=\dot{\varphi}/2e$ which
reads

$V(t)=\frac{I_{c}R_{N}(a^{2}-1)}{a+\sin\left[\sqrt{a^{2}-1}\left|s\right|-\arctan(1/\sqrt{a^{2}-1})\right]},$
(82)

where the parameter

$\displaystyle
s=eI_{c}R_{N}t+\frac{2}{\sqrt{a^{2}-1}}\arctan\frac{1}{\sqrt{a^{2}-1}}$ (83)

is defined within the interval

$\displaystyle-\frac{\pi b}{\sqrt{a^{2}-1}}<s<\frac{\pi
b}{\sqrt{a^{2}-1}},\quad b=\frac{2}{\pi}\arctan\sqrt{\frac{a+1}{a-1}}.$ (84)

The function $V(t)$ is continued periodically outside this interval with the
period $2\pi b/(eI_{c}R_{N}\sqrt{a^{2}-1})$.

The average voltage $\bar{V}$ is determined either by averaging of Eq. (82)
over this period or simply by combining the asymptotic form $F(\varphi\gg
1,a)\simeq\varphi b/\sqrt{a^{2}-1}$ with Eq. (80). Either of these two ways
yields the following $I-V$ curve for our SNS junction (cf. also AB2 )

$\displaystyle\bar{V}=\frac{\pi}{4}I_{c}R_{N}\frac{\sqrt{(I/I_{c})^{2}-1}}{\arctan\sqrt{\frac{I+I_{c}}{I-I_{c}}}}.$
(85)

Bearing in mind the relation $I_{c}R_{N}=\pi\Delta/e$ we note that the
condition $e\bar{V}\ll\Delta$ implies that Eqs. (82) and (85) remain valid as
long as $I-I_{c}\ll I_{c}$. In this limit from Eq. (85) we obtain

$\displaystyle\bar{V}\simeq\frac{\pi\Delta}{2e}\sqrt{(I/I_{c})^{2}-1},\quad
I-I_{c}\ll I_{c}.$ (86)

In addition, with the aid of Eqs. (78) and (82) we conclude that our adiabatic
approximation holds for sufficiently short SNS junctions with $d\ll\xi_{0}$.

Let us now include an ac current bias into our consideration. Provided the
amplitude of an ac signal remains small as compared to the critical current
$I_{\rm ac}\ll I_{c}$ one can solve Eq. (79) perturbatively in $I_{\rm ac}$.
In the first order in this parameter we obtain

$\delta\varphi(t)=2eR_{N}V(t)\int^{t}dt^{\prime}\frac{\delta I+I_{\rm
ac}\cos(\omega t^{\prime}+\vartheta)}{V(t^{\prime})},$ (87)

where $V(t)$ is defined in (82) and $\delta I$ denotes the correction to the
current flowing across the junction. The time derivative of the first order
correction to the phase (87) determines an extra voltage value $\delta
V=\delta\dot{\varphi}\propto I_{\rm ac}$ generated by an ac current signal. It
is necessary to verify that the time average $\langle\delta
V\rangle_{t}=\delta\varphi(t)/t|_{t\to\infty}$ equals to zero AL ; Th . At all
frequencies except for particular values $\omega=2ek\bar{V}$ this condition is
justified for $\delta I=0$ in which case we are back to the $I-V$ curve (85)
derived in the absence of an ac signal. For $\omega=2ek\bar{V}$ we arrive at
integer Shapiro steps in the form

$\displaystyle\delta I_{k}(\vartheta)=-\frac{I_{c}I_{\rm ac}b\sin(\pi
b)\cos\vartheta}{\pi
I\left(b^{2}-k^{2}\right)}\left[1+\frac{\sqrt{a^{2}-1}}{\pi
a^{2}b}\right]^{-1}.$ (88)

In the interesting for us limit $I-I_{c}\ll I_{c}$ Eq. (88) yields $\delta
I_{k=1}\sim I_{\rm ac}$ and $\delta I_{k\geq 2}\sim I_{\rm
ac}\sqrt{a-1}/k^{2}$, i.e. the amplitudes of all Shapiro steps with $k\geq 2$
contain extra small factor $\sqrt{a-1}\ll 1$ and, on top of that, decay
quickly with increasing $k$.

As before, fractional Shapiro steps emerge only in the higher orders in the ac
signal and are mostly pronounced for sufficiently large values of $I_{\rm
ac}$. In the latter limit Eq. (79) can easily be resolved perturbatively in
$I_{c}$. In the lowest order one can neglect the last term in Eq. (79) and
then easily integrate this equation with the result
$\varphi(t)=2eIR_{N}t+2\alpha\sin(\omega t+\vartheta)$. Substituting this
expression for $\varphi(t)$ into the last term in Eq. (79) and averaging the
result over time we again arrive at Eq. (76) with

$\bar{I}\to I,\quad V\to\bar{V},\quad\alpha=\frac{eI_{\rm ac}R_{N}}{\omega}.$
(89)

As we already know, this equation describes fractional Shapiro steps $\delta
I_{\frac{1}{n}}$ which occur at frequencies $\omega=2e\bar{V}/n$. For small
$\alpha$ one again has $\delta I_{\frac{1}{n}}\propto\alpha^{n}$, whereas for
bigger values of $\alpha$ the fractional Shapiro steps $\delta
I_{\frac{1}{n}}(\alpha)$ are displayed in Fig. 9. Note that the interesting
for us parameter range $\alpha\sim 1$ (or, equivalently, $I_{\rm ac}\sim
I_{c}$) is now close to the border of applicability of our approach and,
hence, for such values of $\alpha$ one could in principle expect some
corrections to Eq. (76) which, however, cannot alter any of our conclusions.
In particular, we may conclude that fractional Shapiro steps can dominate over
integer ones at sufficiently small frequencies $\omega$ of an external ac
signal (i.e. for $\alpha\gtrsim 1$) also in the current-biased regime.

##  VIII Discussion

In this work we constructed a microscopic theory describing Shapiro steps in
topologically trivial ballistic SNS junctions subject to external microwave
radiation of arbitrary intensity. Our analysis was mainly focused on the most
interesting limit of subgap voltages and temperatures revealing a number of
features similar to those observed in recent experiments M16 ; topmat ; Gre
with superconducting junctions involving topological insulators. The main such
features are: i) the presence of non-vanishing ”$n=0$ Shapiro step”, i.e. the
critical current $I_{c}(\alpha)$ that does not vanish at any finite $\alpha$,
ii) the presence of Shapiro steps on the $I-V$ curve at all voltages and
frequencies obeying the condition (7), including the ”fractional” one with
$\omega=\omega_{J}/2$ sometimes interpreted as an evidence for the so-called
$4\pi$-periodic Josephson effect and iii) the possibility for ”missing”
Shapiro step with $\omega=\omega_{J}$ (along with other Shapiro steps) at
certain radiation intensities and frequencies.

The feature i) is a direct consequence of a non-sinusoidal CPR inherent to all
types of highly transparent superconducting weak links at low temperatures
$T\ll\Delta$. Hence, it is by no means surprising that this feature has been
observed in experiments Gre  where junctions with high quality interfaces have
been employed.

The feature ii) has to do with multi-photon processes. They play little role
at low intensities of microwave radiation (i.e. for $\alpha\ll 1$) in which
case integer Shapiro steps $\delta I_{k}(\alpha)$ (with $k\geq 1$) are mainly
pronounced, see Figs. 7 and 8. However, at higher ac signals or smaller
microwave frequencies fractional Shapiro steps $\delta
I_{\frac{1}{n}}(\alpha)$ (with $n\geq 2$) gain importance reaching almost the
same magnitude as integer ones for $\alpha\gtrsim 1$, cf., e.g., Fig. 9.

Finally, the feature iii), i.e. ”missing” Shapiro steps, can be explained with
the aid of our results displayed in Fig. 9. Obviously, at around
$\alpha\approx 1$ the step $\delta I_{1}(\alpha)$ almost vanishes being much
smaller than, e.g., $\delta I_{\frac{1}{2}}(\alpha)$ and $\delta
I_{\frac{1}{3}}(\alpha)$. Keeping the amplitude of an ac signal unchanged and
increasing the radiation frequency $\omega$ one decreases $\alpha$ eventually
reaching the regime where $\delta I_{1}(\alpha)$ becomes of order or even
bigger than fractional Shapiro steps. Then in the experiment one would observe
either ”missing” Shapiro step $\delta I_{1}(\alpha)$ or, to the contrary,
fully restored one respectively at the lower and higher radiation frequencies
$\omega$. This is exactly what has been observed in the experiments M16 ;
topmat ; Gre .

Thus, our analysis demonstrates that all experimentally detected features i),
ii) and iii) can be reproduced within a microscopic model describing
topologically trivial SNS junctions with $2\pi$-periodic CPR, cf., e.g., the
results presented in Figs. 6 and 9 of this work and those in Fig. 2 of Ref.
Gre, . None of these features actually requires $4\pi$-periodic CPR which
never pops up in our calculation. In addition, we point out that missing
integer Shapiro steps at $\omega=\omega_{J}$ – along with well pronounced
fractional ones – were recently observed in topologically trivial Josephson
junctions based on InAs quantum wells MS . It is also well known that Shapiro
steps with odd $n=1,3,...$ can be significantly reduced (as compared to ones
with even $n=2,4,...$) or even vanish completely due to size effects FN1 .
Hence, caution is needed while unambiguously interpreting experimental results
for superconducting weak links in terms of Majorana physics.

We believe that all the results obtained here can be directly applied also to
superconducting junctions hosting Majorana-like bound states (1) in the case
of preserved time-reversal symmetry $\delta\approx\Delta$. No $4\pi$-periodic
CPR is expected for any voltage bias $V$ in that case. The situation with
$\delta<\Delta$ may still deserve a separate microscopic analysis.

Note, however, that also in the latter case the treatment of ac Josephson
effect exploiting Eq. (4) generally requires a more rigorous justification. In
the limit of short junctions this equation is verified only in equilibrium in
which case it can also be extended to more complicated situations, e.g., to
hybrid superconducting structures involving triplet pairing GKZ ; KGZ . The
situation out of equilibrium could be more tricky since, e.g., the effect of
MAR needs to be included into consideration. At not too small junction
transmissions MAR sets in already at low voltages and, as we have
demonstrated, becomes an essential ingredient of the whole physical picture.
Multiple Andreev reflection also plays an important role in experiments M16 ;
topmat ; Gre  with junctions involving topological insulators FN2  and, hence,
it needs to be properly accounted for in any theoretical analysis of ac
Josephson effect in such structures.

Acknowledgements

We acknowledge useful discussions with M.S. Kalenkov and A.G. Semenov.

##  Appendix A Matrix inversion

To begin with, let us introduce convenient notations for time dependent
quasiclassical propagators employed in our calculation. These propagators in
general depend on two time variables. They can be expressed in terms of the
following expansion over the voltage harmonics

$x(t,t^{\prime})=\sum_{n=-\infty}^{\infty}\int\frac{d\epsilon}{2\pi}x(\epsilon,n)e^{-i\epsilon(t-t^{\prime})}e^{-ineV(t+t^{\prime})/2},$
(90)

which defines $x(\epsilon,n)$. We will also make use of the notation $z(t,n)$
which essentially means that

$z(t,t^{\prime})=z(t-t^{\prime})\exp(-ineV(t+t^{\prime})/2).$

Below we will employ the following property: Provided $z(t,t^{\prime})$ can be
represented as a convolution

$z(t,t^{\prime})=f(t,m)\circ g(t,n),$

then $z(\epsilon,m+n)$ has only $m+n$ components with

$z(\epsilon,m+n)=f\left(\epsilon+\frac{neV}{2},m\right)g\left(\epsilon-\frac{meV}{2},n\right).$
(91)

In order to proceed let us express the matrix $\hat{X}^{R}(\epsilon,n)$ in the
form

$\hat{X}^{R}(\epsilon,n)=\left(\begin{array}[]{cc}\beta^{R}(\epsilon,n)&\gamma^{R}(\epsilon,n)\\\
\gamma^{R}(\epsilon,n)&\beta^{R}(\epsilon,n)\end{array}\right).$ (92)

Note that $\beta^{R}(\epsilon,n)$ differs from zero for even $n$, while
$\gamma^{R}(\epsilon,n)$ is non-zero for odd $n$. Then the matrix equation

$(\hat{a}^{R}\circ\hat{X}^{R})(\epsilon,n)=\delta_{n,0}$ (93)

can be identically rewritten in terms of the following equations

$\displaystyle\frac{1}{2}\left[g^{R}\left(\epsilon+(n+1)\frac{eV}{2}\right)+g^{R}\left(\epsilon+(n-1)\frac{eV}{2}\right)\right]\beta^{R}(\epsilon,n)+$
(94)
$\displaystyle\frac{1}{2}f^{R}\left(\epsilon+(n-1)\frac{eV}{2}\right)\gamma^{R}\left(\epsilon-\frac{eV}{2},n-1\right)+\frac{1}{2}f^{R}\left(\epsilon+(n+1)\frac{eV}{2}\right)\gamma^{R}\left(\epsilon+\frac{eV}{2},n+1\right)=\delta_{n,0}$

and

$\displaystyle\frac{1}{2}\left[g^{R}\left(\epsilon+(n+1)\frac{eV}{2}\right)+g^{R}\left(\epsilon+(n-1)\frac{eV}{2}\right)\right]\gamma^{R}(\epsilon,n)+$
(95)
$\displaystyle\frac{1}{2}f^{R}\left(\epsilon+(n-1)\frac{eV}{2}\right)\beta^{R}\left(\epsilon-\frac{eV}{2},n-1\right)+\frac{1}{2}f^{R}\left(\epsilon+(n+1)\frac{eV}{2}\right)\beta^{R}\left(\epsilon+\frac{eV}{2},n+1\right)=0.$

It is convenient to introduce the variables $\zeta^{R}(\epsilon,n)$ and
$\tilde{\zeta}^{R}(\epsilon,n)$ which read

$\zeta^{R}(\epsilon,n)=\left\\{\begin{array}[]{c}\beta^{R}(\epsilon,n)\quad{\rm
for\;even}\;n,\\\ \gamma^{R}(\epsilon,n)\quad{\rm
for\;odd}\;n,\end{array}\right.\quad\tilde{\zeta}^{R}(\epsilon,n)=\zeta^{R}\left(\epsilon+\frac{neV}{2},n\right).$
(96)

Then Eqs. (94), (95) can be rewritten as

$\displaystyle\frac{1}{2}\left[g^{R}\left(\epsilon+(2n+1)\frac{eV}{2}\right)+g^{R}\left(\epsilon+(2n-1)\frac{eV}{2}\right)\right]\tilde{\zeta}^{R}(\epsilon,n)+$
(97)
$\displaystyle\frac{1}{2}f^{R}\left(\epsilon+(2n-1)\frac{eV}{2}\right)\tilde{\zeta}^{R}(\epsilon,n-1)+\frac{1}{2}f^{R}\left(\epsilon+(2n+1)\frac{eV}{2}\right)\tilde{\zeta}^{R}(\epsilon,n+1)=\delta_{n,0}.$

This equation has essentially the structure
$\hat{T}\hat{\zeta}^{R}=\delta_{n,0}$, where $\hat{\zeta}^{R}$ is the column
composed of $\tilde{\zeta}^{R}(\epsilon,n)$ and $\hat{T}$ is the symmetric
tridiagonal matrix with the elements

$\displaystyle
T_{n,n}=\frac{1}{2}\left[g^{R}\left(\epsilon+\left(n+\frac{1}{2}\right)eV\right)+g^{R}\left(\epsilon+\left(n-\frac{1}{2}\right)eV\right)\right],$
(98) $\displaystyle
T_{n,n+1}=T_{n+1,n}=\frac{1}{2}f^{R}\left(\epsilon+\left(n+\frac{1}{2}\right)eV\right).$

The equations allowing to determine the matrix $\hat{X}^{A}$ are derived in a
similar manner. We introduce

$\hat{X}^{A}(\epsilon,n)=\left(\begin{array}[]{cc}\beta^{A}(\epsilon,n)&\gamma^{A}(\epsilon,n)\\\
\gamma^{A}(\epsilon,n)&\beta^{A}(\epsilon,n)\end{array}\right),\quad\zeta^{A}(\epsilon,n)=\left\\{\begin{array}[]{c}\beta^{A}(\epsilon,n)\quad{\rm
for\;even}\;n,\\\ \gamma^{A}(\epsilon,n)\quad{\rm
for\;odd}\;n,\end{array}\right.\quad\tilde{\zeta}^{A}(\epsilon,n)=\zeta^{A}\left(\epsilon-\frac{neV}{2},n\right).$
(99)

Making use of Eq. (20) we can rewrite the equation
$-(\hat{X}^{A}\circ\hat{a}^{A})(\epsilon,n)=\delta_{n,0}$ in the form

$\displaystyle\frac{1}{2}\tilde{\zeta}^{A}(\epsilon,n)\left[g^{R}\left(\epsilon-\left(n-\frac{1}{2}\right)eV\right)+g^{R}\left(\epsilon-\left(n+\frac{1}{2}\right)eV\right)\right]^{*}+$
(100)
$\displaystyle\frac{1}{2}\tilde{\zeta}^{A}(\epsilon,n-1)\left[f^{R}\left(\epsilon-\left(n-\frac{1}{2}\right)eV\right)\right]^{*}+\frac{1}{2}\tilde{\zeta}^{A}(\epsilon,n+1)\left[f^{R}\left(\epsilon-\left(n+\frac{1}{2}\right)eV\right)\right]^{*}=\delta_{n,0}.$

Equation (100) can again be expressed as $\hat{T}\hat{\zeta}^{A}=\delta_{n,0}$
with

$\displaystyle
T_{n,n}=\frac{1}{2}\left[g^{R}\left(\epsilon-\left(n+\frac{1}{2}\right)eV\right)+g^{R}\left(\epsilon-\left(n-\frac{1}{2}\right)eV\right)\right]^{*}$
(101) $\displaystyle
T_{n,n+1}=T_{n+1,n}=\frac{1}{2}\left[f^{R}\left(\epsilon-\left(n+\frac{1}{2}\right)eV\right)\right]^{*}.$

The inversion procedure for tridiagonal and block tridiagonal matrices is
outlined in Ref. M, . Employing the theorem 2.3 from M  one can explicitly
evaluate the inverse of tridiagonal matrix $\hat{T}$, defined by Eq. (98). To
be more specific, let us choose $\hat{T}_{nk}$ with $-N\leq n,k\leq N$, i.e.
we are now dealing with a $(2N+1)\times(2N+1)$ square matrix where $N$ is
large. Then, employing the second Eq. (18), we set

$g^{R}\left(\epsilon+eV\left(-N-\frac{1}{2}\right)\right)=g^{R}\left(\epsilon+eV\left(N+\frac{1}{2}\right)\right)=1.$
(102)

All the other functions involved in the matrix elements are defined in Eqs.
(98). Making use of the theorem 2.3 from M  together with the first Eq. (18),
we get

$\displaystyle\left(\hat{T}^{-1}\right)_{0,0}=1,\quad\left(\hat{T}^{-1}\right)_{n,0}=\prod_{n\leq
k\leq-1}\left(-\frac{f^{R}\left(\epsilon+eV\left(k+\frac{1}{2}\right)\right)}{1+g^{R}\left(\epsilon+eV\left(k+\frac{1}{2}\right)\right)}\right),\;{\rm
if}\;n<0$ $\displaystyle\left(\hat{T}^{-1}\right)_{n,0}=\prod_{0\leq k\leq
n-1}\left(-\frac{f^{R}\left(\epsilon+eV\left(k+\frac{1}{2}\right)\right)}{1+g^{R}\left(\epsilon+eV\left(k+\frac{1}{2}\right)\right)}\right),\;{\rm
if}\;n>0.$ (103)

Introducing the notation (46), employing the multiplicative structure of Eqs.
(103) together with the relationship

$\displaystyle\left(g^{R}(\epsilon+neV)-g^{R}(\epsilon+(n+1)eV)\right)\tilde{\zeta}^{R}\left(\epsilon+\frac{eV}{2},n\right)+f^{R}(\epsilon+neV)\tilde{\zeta}^{R}\left(\epsilon+\frac{eV}{2},n-1\right)$
$\displaystyle-f^{R}(\epsilon+(n+1)eV)\tilde{\zeta}^{R}\left(\epsilon+\frac{eV}{2},n+1\right)=-2\,{\rm
sgn}\,n\;\tilde{\zeta}^{R}\left(\epsilon+\frac{eV}{2},n\right),$ (104)

from Eq. (40) we recover Eqs. (48)-(50).

##  Appendix B Shapiro steps at small microwave signals

Employing the perturbative in $\alpha$ procedure outlined in Sec. VIA and
making use of Eqs. (18), (19), we arrive at general expressions for Shapiro
current steps $\delta I_{k}$ at $\omega=2keV$ which take the form

$\displaystyle\delta I_{k}=\,{\rm
Re}\left\\{\frac{\pi\alpha}{16eR_{N}}\sum_{n,p,q}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\left({\rm
sign}\,p-{\rm
sign}\,n\right)\tilde{\zeta}^{R}(\epsilon,\,q)\tilde{\zeta}^{A}(\epsilon,\,p)\tilde{F}(\epsilon,\,q)\right.$
$\displaystyle\times\tilde{\zeta}^{R}\left(\epsilon-(p+n)eV,\,n\right)Z^{\prime}\left(\epsilon-(p-q+n)\frac{eV}{2},\,p+q+n\right)$
(105)
$\displaystyle+\frac{\pi\alpha}{16eR_{N}}\sum_{n,m}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\tilde{\zeta}^{R}(\epsilon-
neV,\,n)Z^{\prime}\left(\epsilon-\frac{(n-m)eV}{2},\,n+m\right)\widetilde{\widetilde{F}}(\epsilon,m)\tilde{\zeta}^{R}\left(\epsilon,\,m\right)$
$\displaystyle+\frac{\pi\alpha}{16eR_{N}}\sum_{n,m}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\,{\rm
sign}\,n\,\tilde{\zeta}^{R}\left(\epsilon-\frac{(n+m)eV}{2},\,n\right)Z^{\prime\prime}\left(\epsilon,\,n+m\right)\tilde{\zeta}^{A}\left(\epsilon+\frac{(n+m)eV}{2},\,m\right)$
$\displaystyle-\frac{\pi\alpha}{4eR_{N}}\sum_{n}\sum_{p+q=n}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}Y^{\prime}\left(\epsilon-\frac{(p-q)eV}{2},\,n\right)\tilde{\zeta}^{R}\left(\epsilon,\,q\right)\tilde{\zeta}^{A}\left(\epsilon,\,p\right)\tilde{F}\left(\epsilon,\,q\right)$
$\displaystyle\left.+\frac{\pi\alpha}{4eR_{N}}\sum_{n}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\,Y^{\prime\prime}(\epsilon,n)\tilde{\zeta}^{R}\left(\epsilon-\frac{neV}{2},\,n\right)\right\\}.$

Here we defined

$\displaystyle
Z^{\prime}(\epsilon,\,n+m)=\left[g^{R}\left(\epsilon+\left(k-\frac{1}{2}\right)eV\right)+g^{R}\left(\epsilon+\left(\frac{1}{2}-k\right)eV\right)-g^{R}\left(\epsilon+\left(k+\frac{1}{2}\right)eV\right)\right.$
$\displaystyle\left.-g^{R}\left(\epsilon-\left(k+\frac{1}{2}\right)eV\right)\right]e^{i\vartheta}\delta\left(2k-n-m\right)+\left[f^{R}\left(\epsilon+keV\right)+f^{R}\left(\epsilon-
keV\right)\right]e^{i\vartheta}$
$\displaystyle\times\left[\delta\left(2k-1-n-m\right)-\delta\left(2k+1-n-m\right)\right]$
(106)

and

$\displaystyle
Z^{\prime\prime}(\epsilon,\,n+m)=\left[g^{K}\left(\epsilon+\left(k-\frac{1}{2}\right)eV\right)+g^{K}\left(\epsilon+\left(\frac{1}{2}-k\right)eV\right)-g^{K}\left(\epsilon+\left(k+\frac{1}{2}\right)eV\right)\right.$
$\displaystyle\left.-g^{K}\left(\epsilon-\left(k+\frac{1}{2}\right)eV\right)\right]e^{i\vartheta}\delta\left(2k-n-m\right)+\left[f^{K}\left(\epsilon+keV\right)+f^{K}\left(\epsilon-
keV\right)\right]e^{i\vartheta}$
$\displaystyle\times\left[\delta\left(2k-1-n-m\right)-\delta\left(2k+1-n-m\right)\right].$
(107)

Similarly we defined

$\displaystyle
Y^{\prime}(\epsilon,n)=\left[g^{R}\left(\epsilon+\left(k+\frac{1}{2}\right)eV\right)+g^{R}\left(\epsilon+\left(k-\frac{1}{2}\right)eV\right)-g^{R}\left(\epsilon-\left(k+\frac{1}{2}\right)eV\right)\right.$
$\displaystyle\left.-g^{R}\left(\epsilon-\left(k-\frac{1}{2}\right)eV\right)\right]e^{i\vartheta}\delta(2k-n)+\left[f^{R}(\epsilon+keV)+f^{R}(\epsilon-
keV)\right]e^{i\vartheta}$
$\displaystyle\times\left[\delta(2k+1-n)+\delta(2k-1-n)\right]$ (108)

and

$\displaystyle
Y^{\prime\prime}(\epsilon,n)=\left[g^{K}\left(\epsilon+\left(k+\frac{1}{2}\right)eV\right)+g^{K}\left(\epsilon+\left(k-\frac{1}{2}\right)eV\right)-g^{K}\left(\epsilon-\left(k+\frac{1}{2}\right)eV\right)\right.$
$\displaystyle\left.-g^{K}\left(\epsilon-\left(k-\frac{1}{2}\right)eV\right)\right]e^{i\vartheta}\delta(2k-n)-\left[f^{K}(\epsilon+keV)+f^{K}(\epsilon-
keV)\right]e^{i\vartheta}$
$\displaystyle\times\left[\delta(2k+1-n)+\delta(2k-1-n)\right].$ (109)

In addition, we employed the definitions

$\displaystyle\tilde{F}(x,q)=F_{1}\left(x-\frac{eV}{2}\right)+F_{2}\left(x+\frac{eV}{2}\right),\;{\rm
if}\;q>0,$ (110)
$\displaystyle\tilde{F}(x,q)=F_{1}\left(x-\frac{eV}{2}\right)+F_{1}\left(x+\frac{eV}{2}\right),\;{\rm
if}\;q=0,$
$\displaystyle\tilde{F}(x,q)=F_{1}\left(x+\frac{eV}{2}\right)+F_{2}\left(x-\frac{eV}{2}\right),\;{\rm
if}\;q<0$

and

$\displaystyle\widetilde{\widetilde{F}}(x,m)=F_{2}\left(x+\frac{eV}{2}\right)-F_{1}\left(x-\frac{eV}{2}\right),\;{\rm
if}\;m>0,$ (111)
$\displaystyle\widetilde{\widetilde{F}}(x,m)=F_{1}\left(x+\frac{eV}{2}\right)-F_{1}\left(x-\frac{eV}{2}\right),\;{\rm
if}\;m=0,$
$\displaystyle\widetilde{\widetilde{F}}(x,m)=F_{1}\left(x+\frac{eV}{2}\right)-F_{2}\left(x-\frac{eV}{2}\right),\;{\rm
if}\;m<0,$

where

$\displaystyle
F_{1}(x)=\tanh\left(\frac{x}{2T}\right)\frac{1-\left|a^{R}\right|^{2}(x)}{1-\left(a^{R*}\right)^{2}(x)},$
(112) $\displaystyle
F_{2}(x)=-\tanh\left(\frac{x}{2T}\right)\frac{a^{R*}(x)}{a^{R}(x)}\frac{1-\left|a^{R}\right|^{2}(x)}{1-\left(a^{R*}\right)^{2}(x)}=-\frac{a^{R*}(x)}{a^{R}(x)}F_{1}(x).$

The functions $\tilde{\zeta}^{R,A}$ were already introduced in Appendix A.
They read

$\tilde{\zeta}^{R}\left(\epsilon+\frac{eV}{2},l\right)=\left\\{\begin{array}[]{l}(-1)^{l}\prod_{1\leq
k\leq l}a^{R}(\epsilon+eVk),\quad{\rm if}\;l>0,\\\ 1,\quad{\rm if}\;l=0,\\\
(-1)^{l}\prod_{l+1\leq k\leq 0}a^{R}(\epsilon+eVk),\quad{\rm
if}\;l<0,\end{array}\right.$ (113)
$\tilde{\zeta}^{A}\left(\epsilon-\frac{eV}{2},l\right)=\left\\{\begin{array}[]{l}(-1)^{l}\prod_{1\leq
k\leq l}a^{R*}(\epsilon-eVk),\quad{\rm if}\;l>0,\\\ 1,\quad{\rm if}\;l=0,\\\
(-1)^{l}\prod_{l+1\leq k\leq 0}a^{R*}(\epsilon-eVk),\quad{\rm
if}\;l<0.\end{array}\right.$ (114)

##  Appendix C Arbitrary microwave signals

Provided the phase $\varphi(t)$ is defined in Eq. (64), we may write

$\displaystyle\cos\frac{\varphi(t)-\varphi(t^{\prime})}{4}=\frac{1}{2}e^{ieV(t-t^{\prime})/2}\left[\sum_{k,l=-\infty}^{\infty}J_{k+l}\left(\frac{\alpha}{2}\right)J_{k-l}\left(\frac{\alpha}{2}\right)e^{ik\omega(t-t^{\prime})+il\omega(t+t^{\prime})}+\right.$
(115)
$\displaystyle\left.+\sum_{k,l=-\infty}^{\infty}J_{k+l+1}\left(\frac{\alpha}{2}\right)J_{k-l}\left(\frac{\alpha}{2}\right)e^{i\left(k+\frac{1}{2}\right)\omega(t-t^{\prime})+i\left(l+\frac{1}{2}\right)\omega(t+t^{\prime})}\right]+c.c.$

Here the parameter $\alpha$ is not necessarily small. We also have

$\displaystyle\cos\frac{\varphi(t)+\varphi(t^{\prime})}{4}=\frac{1}{2}e^{i\vartheta+ieV(t+t^{\prime})/2}\left[\sum_{k,l=-\infty}^{\infty}J_{k+l}\left(\frac{\alpha}{2}\right)J_{k-l}\left(\frac{\alpha}{2}\right)e^{ik\omega(t+t^{\prime})+il\omega(t-t^{\prime})}+\right.$
(116)
$\displaystyle\left.+\sum_{k,l=-\infty}^{\infty}J_{k+l+1}\left(\frac{\alpha}{2}\right)J_{k-l}\left(\frac{\alpha}{2}\right)e^{i\left(k+\frac{1}{2}\right)\omega(t+t^{\prime})+i\left(l+\frac{1}{2}\right)\omega(t-t^{\prime})}\right]+c.c.$

Similarly to Eq. (90) let us define

$x(t,t^{\prime})=\sum_{n,m=-\infty}^{\infty}\int\frac{d\epsilon}{2\pi}x(\epsilon,n,m)e^{-i\epsilon(t-t^{\prime})}e^{-ineV(t+t^{\prime})/2}e^{-im\omega(t+t^{\prime})/2}.$
(117)

For the convolution $z(t,t^{\prime})=x(t,n,m)\circ y(t,k,l)$ we then have

$z(\epsilon,n+k,m+l)=x\left(\epsilon+\frac{keV}{2}+\frac{l\omega}{2},n,m\right)y\left(\epsilon-\frac{neV}{2}-\frac{m\omega}{2},k,l\right).$
(118)

The relationship defining the inverse matrix $\hat{X}^{R}$ can be expressed in
the form

$\left(\hat{a}^{R}\circ\hat{X}^{R}\right)\left(\epsilon,n,m\right)=\delta_{n,0}\delta_{m,0},$
(119)

where for even and odd values of $m$ we have respectively

$\displaystyle\hat{a}^{R}(\epsilon,n,m)=$ (120)
$\displaystyle\frac{1}{2}\delta_{n,0}\sum_{k}J_{k+\frac{m}{2}}\left(\frac{\alpha}{2}\right)J_{k-\frac{m}{2}}\left(\frac{\alpha}{2}\right)\left[g^{R}\left(\epsilon+\frac{eV}{2}-k\omega\right)+g^{R}\left(\epsilon-\frac{eV}{2}-k\omega\right)\right]\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)+$ (123)
$\displaystyle\frac{1}{2}\left(\delta_{n,1}e^{-i\vartheta}+\delta_{n,-1}e^{i\vartheta}\right)\sum_{l}J_{\frac{m}{2}+l}\left(\frac{\alpha}{2}\right)J_{\frac{m}{2}-l}\left(\frac{\alpha}{2}\right)f^{R}(\epsilon-l\omega)\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right)$ (126)

and

$\displaystyle\hat{a}^{R}(\epsilon,n,m)=\frac{1}{2}\delta_{n,0}\sum_{k}J_{k+\frac{1}{2}+\frac{m}{2}}\left(\frac{\alpha}{2}\right)J_{k+\frac{1}{2}-\frac{m}{2}}\left(\frac{\alpha}{2}\right)\times$
(127)
$\displaystyle\times\left[-g^{R}\left(\epsilon+\frac{eV}{2}-\left(k+\frac{1}{2}\right)\omega\right)+g^{R}\left(\epsilon-\frac{eV}{2}-\left(k+\frac{1}{2}\right)\omega\right)\right]\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)+$ (130)
$\displaystyle\frac{1}{2}\left(\delta_{n,1}e^{-i\vartheta}-\delta_{n,-1}e^{i\vartheta}\right)\sum_{l}J_{\frac{m}{2}+l+\frac{1}{2}}\left(\frac{\alpha}{2}\right)J_{\frac{m}{2}-l-\frac{1}{2}}\left(\frac{\alpha}{2}\right)f^{R}\left(\epsilon-\left(l+\frac{1}{2}\right)\omega\right)\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right).$ (133)

As before, we define

$\hat{X}^{R}(\epsilon,n,m)=\left(\begin{array}[]{cc}\beta^{R}(\epsilon,n,m)&\gamma^{R}(\epsilon,n,m)\\\
\gamma^{R}(\epsilon,n,m)&\beta^{R}(\epsilon,n,m)\end{array}\right),$ (134)

where $\beta^{R}$ differs from zero for even $n$ whereas $\gamma^{R}$ is non-
zero for odd $n$. We again introduce

$\zeta^{R}(\epsilon,n,m)=\left\\{\begin{array}[]{c}\beta^{R}(\epsilon,n,m)\quad{\rm
for\;even}\;n,\\\ \gamma^{R}(\epsilon,n,m)\quad{\rm
for\;odd}\;n,\end{array}\right.\quad\tilde{\zeta}^{R}(\epsilon,n,l)=\zeta^{R}\left(\epsilon+\frac{neV}{2}+\frac{l\omega}{2},n,l\right)$
(135)

and compose the column

$\hat{\zeta}^{R}(\epsilon)=\left(\begin{array}[]{c}\vdots\\\
\tilde{\zeta}^{R}(\epsilon,n-1,...)\\\ \tilde{\zeta}^{R}(\epsilon,n,...)\\\
\tilde{\zeta}^{R}(\epsilon,n+1,...)\\\ \vdots\end{array}\right),$ (136)

where dots in $\tilde{\zeta}^{R}(\epsilon,n,...)$ indicate that the index $l$
is running in its range.

The system of equations (119) can then be expressed in the form

$\hat{T}\hat{\zeta}^{R}=\delta_{k,0}\delta_{n,0},$ (137)

where the matrix $\hat{T}$ has the block tridiagonal structure

$\hat{T}=\left(\begin{array}[]{ccccc}\hat{D}_{1}&-e^{i\vartheta}\hat{A}_{2}^{T}&\phantom{d}&\phantom{d}&\phantom{d}\\\
-e^{-i\vartheta}\hat{A}_{2}&\hat{D}_{2}&-e^{i\vartheta}\hat{A}_{3}^{T}&\phantom{d}&\phantom{d}\\\
\phantom{d}&\ddots&\ddots&\ddots&\phantom{d}\\\
\phantom{d}&\phantom{d}&-e^{-i\vartheta}\hat{A}_{N-1}&\hat{D}_{N-1}&-e^{i\vartheta}\hat{A}_{N}^{T}\\\
\phantom{d}&\phantom{d}&\phantom{d}&-e^{-i\vartheta}\hat{A}_{N}&\hat{D}_{N}\end{array}\right)$

similar to that discussed in Ref. M,  and the indices $k$ and $n$ again
correspond respectively to frequency and voltage harmonics. The factors
$e^{\pm i\vartheta}$ factors can be eliminated by means of a unitary
transformation

$\hat{T}=\hat{U}^{\dagger}\hat{T}_{0}\hat{U},\quad\hat{T}_{0}=\left(\begin{array}[]{ccccc}\hat{D}_{1}&-\hat{A}_{2}&\phantom{d}&\phantom{d}&\phantom{d}\\\
-\hat{A}_{2}&\hat{D}_{2}&-\hat{A}_{3}&\phantom{d}&\phantom{d}\\\
\phantom{d}&\ddots&\ddots&\ddots&\phantom{d}\\\
\phantom{d}&\phantom{d}&-\hat{A}_{N-1}&\hat{D}_{N-1}&-\hat{A}_{N}\\\
\phantom{d}&\phantom{d}&\phantom{d}&-\hat{A}_{N}&\hat{D}_{N}\end{array}\right)$
(138)

with $\hat{U}$ being a diagonal unitary matrix with elements
$U_{nn}=e^{i\phi_{n}}$. Equation (138) holds for
$\phi_{n+1}-\phi_{n}=\vartheta$. Provided the difference $k-l$ is even, for
the elements $D^{\left[n\right]}_{kl}(\epsilon)$ of the matrix $\hat{D}_{n}$
we get

$\displaystyle
D^{\left[n\right]}_{kl}(\epsilon)=\frac{1}{2}\sum_{s=-\infty}^{\infty}J_{s+\frac{k-l}{2}}\left(\frac{\alpha}{2}\right)J_{s+\frac{l-k}{2}}\left(\frac{\alpha}{2}\right)\left[g^{R}\left(\epsilon+\left(n+\frac{1}{2}\right)eV+\left(\frac{l+k}{2}-s\right)\omega\right)+\right.$
$\displaystyle+\left.g^{R}\left(\epsilon+\left(n-\frac{1}{2}\right)eV+\left(\frac{l+k}{2}-s\right)\omega\right)\right],$
(139)

whereas for odd values of $k-l$ we obtain

$\displaystyle
D^{\left[n\right]}_{kl}(\epsilon)=\frac{1}{2}\sum_{s=-\infty}^{\infty}J_{s+\frac{1+k-l}{2}}\left(\frac{\alpha}{2}\right)J_{s+\frac{1+l-k}{2}}\left(\frac{\alpha}{2}\right)\left[g^{R}\left(\epsilon+\left(n-\frac{1}{2}\right)eV+\left(\frac{l+k-1}{2}-s\right)\omega\right)\right.$
$\displaystyle\left.-g^{R}\left(\epsilon+\left(n+\frac{1}{2}\right)eV+\left(\frac{l+k-1}{2}-s\right)\omega\right)\right].$
(140)

The above equations indicate that the symmetry relation
$D^{\left[n\right]}_{kl}(\epsilon)=D^{\left[n\right]}_{lk}(\epsilon)$ is
obeyed.

Likewise, the elements $A^{\left[n\right]}_{kl}(\epsilon)$ of the matrix
$\hat{A}_{n}$ respectively for even and odd values of the difference $k-l$
read

$A^{\left[n\right]}_{kl}(\epsilon)=-\frac{1}{2}\sum_{s}J_{\frac{k-l}{2}+s}\left(\frac{\alpha}{2}\right)J_{\frac{k-l}{2}-s}\left(\frac{\alpha}{2}\right)f^{R}\left(\epsilon+\left(n-\frac{1}{2}\right)eV+\left(\frac{l+k}{2}-s\right)\omega\right)$
(141)

and

$A^{\left[n\right]}_{kl}(\epsilon)=-\frac{1}{2}\sum_{s}J_{\frac{k-l+1}{2}+s}\left(\frac{\alpha}{2}\right)J_{\frac{k-l-1}{2}-s}\left(\frac{\alpha}{2}\right)f^{R}\left(\epsilon+\left(n-\frac{1}{2}\right)eV+\left(\frac{l+k-1}{2}-s\right)\omega\right).$
(142)

These equations imply that the matrix $A^{\left[n\right]}_{kl}$ is symmetric
for even $k-l$ and antisymmetric for odd $k-l$.

##  Appendix D Estimate for fractional Shapiro steps

The matrix $\check{Q}_{0}^{-1}\equiv\check{X}$ defined in Eq. (45) can be
identically rewritten as

$\check{X}=\check{X}_{0}+\sum_{m=1}^{\infty}(-1)^{m}\check{X}_{0}\left(\delta\check{a}\check{X}_{0}\right)^{m}.$
(143)

Here $\check{X}_{0}$ is the inverse matrix $\check{Q}_{0}^{-1}$ evaluated for
$\check{a}_{0}=\delta(t-t^{\prime})$ and
$\delta\check{a}=\check{a}-\check{a}_{0}$.

In order to illustrate the idea of our estimate let us consider one of the
terms in Eq. (40), e.g., the term containing the combination
$\tau_{3}\hat{X}^{R}\tau_{3}\hat{b}^{K}$. One of the contributions to the
current (proportional to $\delta\check{a}$) generated by this term takes the
form (cf. the first term in the square brackets in Eq. (116))

$\displaystyle\frac{ie{\mathcal{N}}}{4}\int
dt^{\prime}f^{K}(t-t^{\prime})f^{R}(t^{\prime}-t)e^{ieV(t+t^{\prime})+i\vartheta}$
$\displaystyle\times\sum_{p,l=-\infty}^{\infty}J_{p+l}\left(\frac{\alpha}{2}\right)J_{p-l}\left(\frac{\alpha}{2}\right)e^{ip\omega(t+t^{\prime})+il\omega(t^{\prime}-t)}$
(144)

Under the condition $eV=-p\omega$ the dependence of the combination (144) on
$t+t^{\prime}$ drops out and it depends only on the time difference
$t-t^{\prime}$. Performing the Fourier transformation we receive the
contributions containing the integrals

$\sim\sum_{l}\int\frac{d\epsilon}{2\pi}f^{K}(\epsilon)f^{R}(\epsilon+l\omega)$
(145)

It is easy to see that the main contribution to this sum is provided by the
term with $l=0$. In order to demonstrate that let us consider the combination

$\displaystyle\frac{ie{\mathcal{N}}}{4}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}f^{K}(\epsilon)\left(f^{R}(\epsilon)+f^{A}(\epsilon)\right)$
(146) $\displaystyle=$
$\displaystyle\frac{ie{\mathcal{N}\Delta^{2}}}{4}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\tanh\frac{\epsilon}{2T}\left(\mathcal{R}_{+}^{-2}(\epsilon)-\mathcal{R}_{-}^{-2}(\epsilon)\right)$

entering into the expression for the current. Keeping track of a small
imaginary part $\pm i\gamma$ in the expressions for
${\mathcal{R}_{\pm}(\epsilon)=\sqrt{(\epsilon\pm i\gamma)^{2}-\Delta^{2}}}$
and making use of the Sokhotski-Plemelj theorem, for the combination (146) we
obtain

$\displaystyle
e{\mathcal{N}}\Delta^{2}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}\tanh\frac{\epsilon}{2T}\frac{\gamma\epsilon}{\left[(\epsilon+i\gamma)^{2}-\Delta^{2}\right]\left[(\epsilon-i\gamma)^{2}-\Delta^{2}\right]}$
(147) $\displaystyle=$
$\displaystyle\frac{e{\mathcal{N}}\Delta^{2}}{8}\int\limits_{-\infty}^{\infty}\frac{d\epsilon}{\epsilon}\tanh\frac{\epsilon}{2T}\left[\delta(\epsilon+\Delta)+\delta(\epsilon-\Delta)\right]$
$\displaystyle=$
$\displaystyle\frac{e{\mathcal{N}}\Delta}{4}\tanh\frac{\Delta}{2T}\sim I_{c}.$

Thus we observe that the main contribution to the integral with $l=0$ in Eq.
(145) comes from $\epsilon=\pm\Delta$ since the expression under this integral
strongly peaked at these two values of $\epsilon$. The terms with $l\neq 0$
contain no such peaks, hence, their contributions are smaller and can be
safely neglected.

Substituting the estimate (147) into Eq. (144) and identifying $n=-2p$ we
arrive at the result for the magnitude of fractional Shapiro steps $\delta
I_{\frac{1}{n}}$ which occur at $\omega=\omega_{J}/n$ with even $n$:

$\delta I_{\frac{1}{n}}\sim I_{c}J^{2}_{n/2}(\alpha/2).$ (148)

Repeating the whole analysis with the second term in the square brackets of
Eq. (116) we recover the analogous estimate for odd values of $n$:

$\delta I_{\frac{1}{n}}\sim I_{c}J_{(n+1)/2}(\alpha)J_{(n-1)/2}(\alpha/2).$
(149)

Finally, we note that the above simple analysis includes only linear in
$\delta\check{a}$ terms in the formal series (143). Hence, one may expect that
the estimates (148) and (149) should work reasonably well for $\alpha\lesssim
1$ and may become less accurate for bigger values of $\alpha$.

## References

  * (1) See, e.g., R. Aguado and L.P. Kouwenhoven, Physics Today 73, 44 (2020).
  * (2) H.-J. Kwon, K. Sengupta, and V.M. Yakovenko, Eur. Phys. J. B 37, 349 (2004).
  * (3) L. Fu and C.L. Kane, Phys. Rev. B 79, 161408(R) (2009).
  * (4) G. Tkachov and E.M. Hankiewicz, Phys. Rev. B 88, 075401 (2013).
  * (5) M. Tinkham, Introduction to Superconductivity, second edition, McGraw-Hill, New York, 1996.
  * (6) A. Barone and G. Paterno, Physics and Applications of the Josephson Effect, John Wiley & Sons, New York, 1982.
  * (7) K.K. Likharev, Dynamics of Josephson Junctions and Circuits, Gordon and Breach Science Publishers, New York, 1986.
  * (8) J. Wiedenmann, E. Bocquillon, R. S. Deacon, S. Hartinger, O. Herrmann, T.M. Klapwijk, L. Maier, C. Ames, C. Brüne, C. Gould, A. Oiwa, K. Ishibashi, S. Tarucha, H. Buhmann, and L.W. Molenkamp, Nat. Commun. 7, 10303 (2016).
  * (9) E. Bocquillon, J. Wiedenmann, R.S. Deacon, T.M. Klapwijk, H. Buhmann, L.W. Molenkamp, Microwave Studies of the Fractional Josephson Effect in HgTe-Based Josephson Junctions, in: Topological Matter, eds. D. Bercioux, J. Cayssol, M. Vergniory, M. Reyes Calvo, Springer Series in Solid-State Sciences, vol. 190, p. 115, Springer, Cham, 2018.
  * (10) K. Le Calvez, L. Veyrat, F. Gay, P. Plaindoux, C. Winkelmann, H. Courtois, and B. Sacépé, Commun. Phys. 2, 4 (2019).
  * (11) G. Schön and A.D. Zaikin, Phys. Rep. 198, 237 (1990).
  * (12) A.D. Zaikin, Influence of Coulomb and Proximity Effects on Electron Tunneling through Normal Metal-Superconductor Interfaces, Physica B 203, 255 (1994).
  * (13) I. Snyman and Yu.V. Nazarov, Phys. Rev. B 77, 165118 (2008).
  * (14) A.D. Zaikin and D.S. Golubev, Dissipative Quantum Mechanics of Nanostructures: Electron Transport, Fluctuations and Interactions (Jenny Stanford, Singapore, 2019).
  * (15) I.O. Kulik and A.N. Omel’yanchuk, Fiz. Nizk. Temp. 3, 945 (1977) [Sov. J. Low Temp. Phys. 3, 459 (1977)].
  * (16) A.V. Galaktionov and A.D. Zaikin, Phys. Rev. B 65, 184507 (2002).
  * (17) A. Furusaki and M. Tsukada, Physica B 165–166, 967 (1990).
  * (18) C.W.J. Beenakker and H. van Houten, Phys. Rev. Lett. 66, 3056 (1991).
  * (19) Here we only consider grand canonical ensembles where the total number of electrons is not fixed. Generalization of Eq. (4) to canonical ensembles has been worked out, e.g., by S.V. Sharov and A.D. Zaikin, Phys. Rev. B 71, 014518 (2005).
  * (20) T.M. Klapwijk, G.E.Blonder, and M.Tinkham, Physica B+C 109-110, 1657 (1982).
  * (21) A.V. Zaitzev, Zh. Eksp. Teor. Fiz. 78, 221 (1980) [Sov. Phys. JETP 51, 111 (1980)].
  * (22) A.D. Zaikin, Zh. Eksp. Teor. Fiz. 84, 1560 (1983) [Sov. Phys. JETP 57, 910 (1983)].
  * (23) U. Gunsenheimer and A.D. Zaikin, Phys. Rev. B 50, 6317 (1994).
  * (24) D. Averin and A. Bardas, Phys. Rev. Lett. 75, 1831 (1995).
  * (25) D. Averin and A. Bardas, Phys. Rev. B 53, R1705 (1996).
  * (26) U. Gunsenheimer and A.D. Zaikin, EPL 41, 195 (1998).
  * (27) J.C. Cuevas, J. Heurich, A. Martín-Rodero, A. Levy Yeyati, and G. Schön, Phys. Rev. Lett. 88, 157001 (2002).
  * (28) F. Dominguez, O. Kashuba, E. Bocquillon, J. Wiedenmann, R.S. Deacon, T.M. Klapwijk, G. Platero, L.W. Molenkamp, B. Trauzettel, and E.M. Hankiewicz, Phys. Rev. B 95, 195430 (2017).
  * (29) G. Eilenberger, Z. Phys. 214, 195 (1968).
  * (30) W. Belzig, F.K. Wilhelm, C. Bruder, G. Schön, and A.D. Zaikin, Superlatt. and Microstr. 25, 1251 (1999).
  * (31) A.V. Zaitsev, Zh. Eksp. Teor. Fiz. 86, 1742 (1984) [Sov. Phys. JETP 59, 1015 (1985)].
  * (32) A.V. Galaktionov and A.D. Zaikin, Phys. Rev. B 82, 184520 (2010).
  * (33) F.S. Bergeret, P. Virtanen, T.T. Heikkilä, and J.C. Cuevas, Phys. Rev. Lett. 105, 117001 (2010).
  * (34) M.S. Kalenkov, P.E. Dolgirev, and A.D. Zaikin, Phys. Rev. B 101, 180505(R) (2020).
  * (35) M.S. Kalenkov and A.D. Zaikin, Eur. Phys. J. Spec. Top. 230, 813 (2021).
  * (36) L.G. Aslamazov and A.I. Larkin, Pis’ma Zh. Eksp. Teor. Fiz. 9, 150 (1969) [JETP Lett. 9, 87 (1969)].
  * (37) E.D. Thompson, J. Appl. Phys. 44, 5587 (1973).
  * (38) M.C. Dartiailh, J.J. Cuozzo, B.H. Elfeky, W. Mayer, J. Yuan, K.S. Wickramasinghe, E. Rossi, and J. Shabani, Nature Comm. 12, 78 (2021).
  * (39) See, e.g., Sec. 11.2.2. of Ref. BP, .
  * (40) A.V. Galaktionov, M.S. Kalenkov, and A.D. Zaikin, Phys. Rev. B 77, 094520 (2008).
  * (41) M.S. Kalenkov, A.V. Galaktionov, and A.D. Zaikin, Phys. Rev. B 79, 014521 (2009).
  * (42) This conclusion is fully consistent, e.g., with the observation of a fairly large excess current in the experiments M16 ; topmat ; Gre .
  * (43) G. Meurant, SIAM Journal on Matrix Analysis and Applications 13, 707 (1992).