# Josephson quantum mechanics at odd parity

Manuel Houzet, Julia S. Meyer  Univ. Grenoble Alpes, CEA, Grenoble INP, IRIG,
Pheliqs, F-38000 Grenoble, France  Yuli V. Nazarov  Kavli Institute of
Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

## Supplemental Material

In this Supplemental Material, we present the details of the derivation of our
results.

##  I Odd parity ground state: no quantum fluctuations

Here we recall the derivation of Andreev bound states in the tunneling limit
closely following Appendix D of Ref. [1]. For the moment, we disregard the
fluctuations of the phase treating it as a number.

We describe a Josephson junction with $N_{\rm ch}$ tunneling channels by the
Hamiltonian

$H=\sum_{km\sigma}\varepsilon_{k}\alpha_{km\sigma}^{\dagger}\alpha_{km\sigma}+\sum_{km\sigma}\varepsilon_{k}\gamma_{km\sigma}^{\dagger}\gamma_{km\sigma}+H_{T}.$
(1)

Here $\alpha_{km\sigma}$ ($\gamma_{km\sigma}$) is a fermionic annihilation
operator of a Bogoliubov quasiparticle in the left (right) lead, with orbital
label $k$, spin $\sigma$, and excitation energy $\varepsilon_{k}$, which
participates in the $m^{\rm th}$ tunneling channel of the junction ($1\leq
m\leq N_{\rm ch}$, where $N_{\rm ch}$ is the number of channels). The
quasiparticle energies are $\varepsilon_{k}=\sqrt{\xi^{2}_{k}+\Delta^{2}}$,
where $\Delta$ is the superconducting gap and $\xi_{k}$ is the electron energy
of the state $k$ measured from the Fermi level in the absence of
superconductivity.

The tunneling is diagonal in channels, the corresponding term reads ($L$ being
the infinite normalization length of the channel)

$H_{T}=e^{i\varphi/2}\sum_{kk^{\prime}m\sigma}\frac{t_{m}}{L}a_{km\sigma}^{\dagger}c_{k^{\prime}m\sigma}+{\rm
H.c.}.$ (2)

It is expressed in terms of the annihilation operator of an electron in the
right lead, $a_{km\sigma}=u_{k}\alpha_{km\sigma}+\sigma
v_{k}\alpha^{\dagger}_{km-\sigma}$, with coherence factors
$u_{k},v_{k}=\sqrt{(1\pm\xi_{k}/\varepsilon_{k})/2}$, and the annihilation
operator in the left lead $c_{km\sigma}$ expressed in terms of
$\gamma_{km\sigma}$ in a similar way. The tunnel matrix element $t_{m}$ is
real by virtue of time reversibility and defines the transmission coefficient
$T_{m}$ of the corresponding channel, $T_{m}=(2\pi\nu t_{m})^{2}\ll 1$, $\nu$
being the density of states per spin and channel.

Following Appendix D of Ref. [1], we note that an effective low-energy
description of the odd parity sector involves the states with energies close
to $\Delta$, so that $u_{k}\approx v_{k}\approx 1/\sqrt{2}$. With this, the
tunneling Hamiltonian becomes

$H_{T}=\sum_{m,\sigma}i\frac{t_{m}}{L}\sin\frac{\varphi}{2}\sum_{kk^{\prime}}\left(\alpha^{\dagger}_{km\sigma}\gamma_{k^{\prime}m\sigma}-\gamma^{\dagger}_{km\sigma}\alpha_{k^{\prime}m\sigma}\right).$
(3)

It involves only the quasiparticle transfers and no terms
creating/annihilating a pair of quasiparticles.

To get further insight into the specific structure of the tunneling
Hamiltonian, let us introduce two linear combinations of the left and right
quasiparticle operators,

$\beta_{kms\sigma}=(\alpha_{km\sigma}-is\gamma_{km\sigma})/\sqrt{2}$ (4)

with the index $s=\pm 1$ that we call superposition index, or chirality. The
operators create a quasiparticle in an equal-weight superposition of the
states in the left and right lead.

Switching to these new operators, we observe that $H_{T}$ conserves the
superposition index, and the full Hamiltonian describing the quasiparticles
close to the gap edge reads

$H_{{\rm
full}}=\sum_{kms\sigma}\varepsilon_{k}\beta_{kms\sigma}^{\dagger}\beta_{kms\sigma}-\frac{t_{m}}{L}\sin\frac{\varphi}{2}\sum_{kk^{\prime}ms\sigma}s\beta_{kms\sigma}^{\dagger}\beta_{k^{\prime}ms\sigma}$
(5)

with $\varepsilon_{k}\approx\Delta+\xi_{k}^{2}/2\Delta$.

Let us consider a single quasiparticle with given $\sigma$ and $s$ in a given
channel $m$. The Hamiltonian in the space of the possible states $k$ of this
quasiparticle reads

$H_{{\rm qp}}=\Delta+\sum_{k}\frac{\xi_{k}^{2}}{2\Delta}|k\rangle\langle
k|-s\sin\frac{\varphi}{2}\frac{t_{m}}{L}\sum_{k,k^{\prime}}|k\rangle\langle
k^{\prime}|.$ (6)

The eigenstates with $E>\Delta$ correspond to delocalized quasiparticles and
are of no interest to us. We concentrate on possible bound states of the
quasiparticle with $E=\Delta-\Omega$, $\Omega$ being the positively defined
binding energy of the quasiparticle. The Schrödinger equation for the wave
function $\Psi=\sum_{k}\psi_{k}|k\rangle$ reads

$0=\left(\Omega+\frac{\xi_{k}^{2}}{2\Delta}\right)\psi_{k}-s\sin\frac{\varphi}{2}\frac{t_{m}}{L}\sum_{k^{\prime}}\psi_{k}^{\prime}$
(7)

and can be easily solved for $\Phi=\sum_{k}\psi_{k}$,

$\Phi=s\sin\frac{\varphi}{2}\frac{t_{m}}{L}\sum_{k}\frac{\Phi}{\Omega+{\xi_{k}^{2}}/{2\Delta}}.$
(8)

We replace the sum over $k$ by integration over energies, $\sum_{k}\to
L\nu\int d\xi$. The resulting integral over $\xi$ converges at the energies of
the order of $\Omega$. This gives

$\sqrt{\Omega}=\sqrt{2E^{(m)}_{J}}s\sin\frac{\varphi}{2},$ (9)

where $E^{(m)}_{J}\equiv\Delta(\pi\nu t_{m})^{2}\equiv\Delta T_{m}/4$ is the
contribution of the level $m$ to the total Josephson energy $E^{*}_{J}$ in the
even ground state, $E^{*}_{J}=\sum_{m}E^{(m)}_{J}$. For simplicity, we refer
to the contribution of the lowest level, that is, the level with the largest
$t_{m}$, just as $E_{J}$. Therefore, $E^{*}_{J}\geq E_{J}$ in general and
$E^{*}_{J}=E_{J}$ for a single-channel situation where the contribution of
other channels can be neglected.

##  II Odd parity ground state: Hamiltonian approach

Now let us take into account quantum fluctuations of $\varphi$. To this end,
we promote the phase to an operator. Let us consider first the even parity
state. In this case, the Hamiltonian comprises the Hamiltonians of the
environment and the junction,

$H_{{\rm e}}=H_{{\rm env}}-E^{*}_{J}\cos\hat{\varphi}.$ (10)

We use a bosonic description of the linear environment,

$H_{\rm
env}=\sum_{q}\omega_{q}b_{q}^{\dagger}b_{q},\qquad\hat{\varphi}=\varphi+\sum_{q}\lambda_{q}(b_{q}+b^{\dagger}_{q}).$
(11)

Here $b_{q}$ is a bosonic annihilation operator of an excitation in mode $q$,
with energy $\omega_{q}$. The operator of the phase $\hat{\varphi}$ is
composed of the constant phase bias part $\varphi$ and the fluctuating part
that is a linear superposition of the bosonic creation/annihilation operators.
The coefficients $\lambda_{q}$ in this superposition are chosen to represent
the dissipative part of the environment impedance seen by the junction,

${\rm
Re}Z(\omega)=\frac{\pi\omega}{4e^{2}}\sum_{q}\lambda_{q}^{2}\delta(\omega-\omega_{q}),$
(12)

at frequency $\omega>0$.

It is convenient for us to use this Hamiltonian with the energy counted from
its ground state $|0\rangle$, and define

$H=H_{{\rm e}}-E^{(e)}_{g},\;\;\;H|0\rangle=0.$ (13)

To obtain the odd parity state, we add a quasiparticle of chirality $s$ to the
lowest Andreev state. With this, the total Hamiltonian $H_{{\rm odd}}$ is
obtained by combining the Hamiltonian of Eq. (6) and $H$, the binding energy
$\Omega$ being an eigenvalue of this Hamiltoninan,

$\displaystyle H_{{\rm odd}}$ $\displaystyle=$ $\displaystyle
H+\sum_{k}\frac{\xi_{k}^{2}}{2\Delta}|k\rangle\langle
k|-s\sin\frac{\hat{\varphi}}{2}\frac{t_{m}}{L}\sum_{k,k^{\prime}}|k\rangle\langle
k^{\prime}|,$ (14) $\displaystyle 0$ $\displaystyle=$
$\displaystyle(\Omega+H_{{\rm odd}})|\Psi\rangle.$ (15)

With all the energy shifts we made, the ground state energy at odd parity is
given by

$E_{g}^{(o)}=E^{(e)}_{g}+\Delta-\Omega.$ (16)

The wave function $|\Psi\rangle$ is in the combined space of $k$ and
environmental degrees of freedom. Further derivation essentially repeats the
steps done in the previous Section. We substitute
$|\Psi\rangle=\sum_{k}|\psi_{k}\rangle|k\rangle$, $|\psi_{k}\rangle$ being a
wave function in the space of environmental degrees of freedom, and obtain a
closed equation for $|\Phi\rangle\equiv\sum_{k}|\psi_{k}\rangle$,

$(\sqrt{\Omega+H}-\hat{A})|\Phi\rangle;\quad\quad\hat{A}\equiv
s\sqrt{2E_{J}}\sin\frac{\hat{\varphi}}{2}.$ (17)

While Eq. (17) is an eigenvalue problem for an operator, it is more
complicated than a standard Schrödinger equation that can be readily solved by
the diagonalization of the operator. For the case in hand, the diagonalization
is not enough: as a result of it, one gets a set of eigenvalues that
parametrically depend on $\Omega$. Since $H$ is a non-negatively defined
operator, and $A$ is an operator restricted by $\pm\sqrt{2E_{J}}$, at
sufficiently big values of $\Omega$, $\Omega>\sqrt{2E_{J}}$ all eigenvalues of
the set are positive. Decreasing $\Omega$, we achieve the situation where the
lowest eigenvalue of the operator is zero. This value of $\Omega$ thus
corresponds to the actual binding energy. It can happen that the lowest
eigenvalue remains positive up to $\Omega=0$. In this case, there is no bound
state (at these settings of $s$ and phase bias $\varphi$).

##  III Odd parity ground state: Path-integral approach

The path-integral approach to Josephson quantum mechanics (at even parity) has
been developed already in the 1980s (see [2] for an extensive early review).
It has a clear advantage over any Hamiltonian method since it reduces the
description to the relevant variable $\varphi(\tau)$ only. Here we consider
only the imaginary-time, zero-temperature version of this path-integral
method. Within this method, the averages of a product
$\prod_{j}\varphi(\tau_{j})$ over the paths represent the averages of the
time-ordered products of Matsubara operators
$\hat{\varphi}(\tau)\equiv\exp(-H\tau)\hat{\varphi}\exp(H\tau)$ over the
ground state. The weight of a path $\varphi(\tau)$ is given by $\exp(-{\cal
S})$, where the action ${\cal S}(\varphi(\tau))$ comprises a quadratic part
representing the linear environment and a part representing the Josephson
energy,

$\displaystyle{\cal S}={\cal S}_{\rm{env}}-E^{*}_{J}\int
d\tau\cos\varphi(\tau),$ (18) $\displaystyle{\cal
S}_{\rm{env}}=\sum_{\omega}\frac{\omega}{{8}e^{2}{\cal
Z}(\omega)}|\varphi_{\omega}|^{2}.$ (19)

Here, ${\cal Z}(\omega)$ is the impedance at imaginary frequency. For the
simplest circuit of a resistor $R$ and a capacitor $C$ in parallel, ${\cal
Z}^{-1}(\omega)\ =R^{-1}\ {\rm sgn}(\omega)+\omega C$, while the usual
impedance $Z(\omega)$ at real frequency is given by
$Z^{-1}(\omega)=R^{-1}-i\omega C$ for this circuit.

It is worth noting already in the beginning that our approach to the odd
parity state does not involve any new action. Rather, we derive and use an
implicit equation for the binding energy $\Omega$ that involves path-integral
averages over the even parity state. This distinguishes the present approach
from the standard ones.

We start with Eq. (17). Let us note that the fact that the operator in it has
a singular value implies a divergence of the expectation value of the inverse
operator over the even parity ground state,

${\rm Ex}\equiv\langle 0|\frac{1}{\sqrt{\Omega+H}-\hat{A}}|0\rangle=\infty.$
(20)

This works if $\langle 0|\Phi\rangle\neq 0$, as we can safely assume. In fact,
the actual binding energy is determined by the lowest value of $\Omega$ at
which the expectation value diverges. Indeed, Eq. (17) at $\Omega$ smaller
than the binding energy has singular eigenvalues corresponding to the
excitations of the environment on the background of the odd parity ground
state.

Let us expand the inverse operator in terms of $\hat{A}$ and concentrate on
the second term of the expansion,

${\rm Ex}^{(2)}=\langle
0|\frac{1}{\sqrt{\Omega+H}}\hat{A}\frac{1}{\sqrt{\Omega+H}}\hat{A}\frac{1}{\sqrt{\Omega+H}}|0\rangle.$
(21)

We represent the inverse square root entering the expansion in terms of an
integral over a time-like variable

$\frac{1}{\sqrt{\Omega+H}}=\int_{-\infty}^{\infty}d\tau\;G_{0}(\tau)e^{-\Omega\tau}e^{-H\tau};\qquad
G_{0}(\tau)\equiv\Theta(\tau)/\sqrt{\pi\tau},$ (22)

where $\Theta(\tau)$ is the Heaviside step function. With this, we present the
term with the expectation value of the product of Matsubara operators,

${\rm Ex}^{(2)}=\int d\tau
d\tau_{1}d\tau_{2}\;e^{-\Omega\tau}G_{0}(\tau_{1})G_{0}(\tau_{2}-\tau_{1})G_{0}(\tau-\tau_{2})\langle
0|\hat{A}(\tau_{1})\hat{A}(\tau_{2})|0\rangle.$ (23)

The product of the Matsubara operators is already time-ordered, so its
expectation value can be replaced with the average over the paths
($A(\tau)\equiv s\sqrt{2E_{J}}\sin(\varphi(\tau)/2)$):

$\langle 0|\hat{A}(\tau_{1})\hat{A}(\tau_{2})|0\rangle\to\langle
A(\tau_{1})A(\tau_{2})\rangle.$ (24)

We do such replacement for all terms of the expansion and sum them up. With
this, the expectation value in question is expressed via a propagator
$G(\tau,\tau^{\prime})$ defined by a rather standard relation,

$G({\tau,\tau^{\prime}})=G_{0}({\tau-\tau^{\prime}})+\int
d\tau_{1}\;G_{0}({\tau}-\tau_{1})A({\tau_{1}})G(\tau_{1},\tau^{\prime}).$ (25)

Let us notice an analogy of the propagator in use and that of a particle in a
disordered potential [3], $A(\tau)$ playing the role of the potential. Figure
1a gives a diagram presenting a term in the perturbation expansion of
$G(\tau,\tau^{\prime})$.

The average propagator is a function of the time difference only, and

${\rm Ex}=\int d\tau e^{-\Omega t}\bar{G}(\tau);\qquad\bar{G}(\tau)=\langle
G(\tau_{0}+\tau,\tau_{0})\rangle.$ (26)

The previous reasoning implies that the average propagator grows with
increasing time difference, $\bar{G}(\tau)\propto e^{\Omega_{{\rm b}}\tau}$.
The integral in (26) thus diverges at $\Omega<\Omega_{b}$ and converges at
$\Omega>\Omega_{{\rm b}}$. This identifies $\Omega_{{\rm b}}$ as the actual
binding energy.

The diagram in Fig. 1b represents a term in the expansion of the average
propagator. Black dots where the dashed lines come together denote the
(higher-order) correlators of $A(\tau)$.

It would be nice to have a closed expression for $\Omega_{{\rm b}}$, at least
a perturbative one. Yet this is beyond our reach. Instead, we use a common
wisdom of disorder-averaged propagators [3] and introduce a self-energy
$\Sigma(\omega)$ such that, in frequency representation,

$\bar{G}(\omega)=\frac{1}{\sqrt{i\omega}-\Sigma(\omega)}.$ (27)

The self-energy (in our scheme, its dimension is the square-root of energy)
admits a perturbative expression in terms of correlators of $A(\tau)$. The
binding energy is defined by the presence of a pole or other singularity at
imaginary $\omega$, $i\omega=\Omega$, this yields

$\sqrt{\Omega}=\Sigma(-i\Omega),$ (28)

$\Omega_{{\rm b}}$ is thus the root of this equation. In Fig. 1c we give all
diagrams contributing to $\Sigma(\omega)$ up to the fourth order. This is not
the only way to draw the expansion: for instance, one can resum the propagator
including all diagramms with $\langle A\rangle$,

Figure 1:  The propagator used for the calculation of the binding energy. a. A
term in the perturbation expansion of a non-averaged propagator. b. A term in
the perturbation expansion of the averaged propagator. c. The diagrams for
$\Sigma$ up to 4th order.

$G_{0}(\omega)\to\bar{G}_{0}\equiv\frac{1}{\sqrt{i\omega}-\langle A\rangle},$
(29)

so the expansion will start with the second order and would not include black
dots with a single line. We give the explicit expressions for this case up to
the fourth order in time representation

$\displaystyle\Sigma^{(2)}(\tau)$ $\displaystyle=$
$\displaystyle\langle\langle A(0)A(\tau)\rangle\rangle\bar{G}_{0}(\tau),$ (30)
$\displaystyle\Sigma^{(3)}(\tau)$ $\displaystyle=$ $\displaystyle\int
d\tau_{1}\langle\langle
A(0)A(\tau_{1})A(\tau)\rangle\rangle\bar{G}_{0}(\tau_{1})\bar{G}_{0}(t-\tau_{1}),$
(31) $\displaystyle\Sigma^{(4)}(\tau)$ $\displaystyle=$ $\displaystyle\int
d\tau_{1}d\tau_{2}\;\bar{G}_{0}(\tau_{1})\bar{G}_{0}(\tau_{2}-\tau_{1})\bar{G}_{0}(\tau-\tau_{2})$
(32) $\displaystyle\times\left(\langle\langle
A(0)A(\tau_{1})A(\tau_{2})A(\tau)\rangle\rangle+\langle\langle
A(0)A(\tau)\rangle\rangle\langle\langle
A(\tau_{1})A(\tau_{2})\rangle\rangle+\langle\langle
A(0)A(\tau_{2})\rangle\rangle\langle\langle
A(\tau_{1})A(\tau)\rangle\rangle\right).$

Here, the double angle brackets $\langle\langle\dots\rangle\rangle$ denote the
cumulant of all variables between the brackets.

##  IV Small impedance

###  IV.1 Warm-up: single oscillator

As a warm up, let us consider an environment consisting of a single
oscillator, such that the phase fluctuation created by it is small as compared
to $\pi$. We can proceed with the Hamiltonian method. The Hamiltonian of the
oscillator can be written as

$H_{osc}=-\frac{E_{C}}{2}\partial^{2}_{f}+\frac{E_{L}}{2}f^{2}=\hbar\omega_{0}\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right),$
(33)

$f$ being the addition to the phase, $\hat{a}$ being the annihilation operator
of the oscillation. Here, $E_{L}=(4e^{2}L)^{-1}$, $E_{C}=4e^{2}/C$. The
oscillator frequency $\omega_{0}=\sqrt{E_{C}E_{L}}$ is much smaller than
$E_{L}$, this guarantees the smallness of the fluctuation,

$\hat{f}=\sqrt{\frac{\omega_{0}}{2E_{L}}}(\hat{a}+\hat{a}^{\dagger}).$ (34)

Here $\omega_{0}/E_{L}=\sqrt{L/C}$ is an effective impedance of the
oscillator. By virtue of phase bias, $E_{L}\gg E_{J}$. As to $E_{J}$ and
$\omega_{0}$, they can be in an arbitrary relation: this complicates the
analysis.

Since the fluctuations are small, we can expand

$\displaystyle\cos(\hat{\varphi})\to\cos(\varphi)-\sin(\varphi)\hat{f}-\cos(\varphi)\frac{\hat{f}^{2}}{2},$
(35)
$\displaystyle\sin\frac{\hat{\varphi}}{2}\to\sin\frac{\varphi}{2}+\cos\frac{\varphi}{2}\frac{\hat{f}}{2}-\sin\frac{\varphi}{2}\frac{\hat{f}^{2}}{8}.$
(36)

The relevant corrections to an energy are the first-order corrections in
$\langle\langle f^{2}\rangle\rangle$ and eventually $\langle f\rangle$ if
present, and the second-order corrections in $\hat{f}$.

Let us start with the even parity correction. The first and second order
corrections read correspondingly to the two terms in the following equation:

$\delta E_{g}^{(e)}=E^{*}_{J}\cos\varphi\frac{\langle
0|\hat{f}^{2}|0\rangle}{2}-\frac{|\langle
0|E_{J}^{*}\sin\varphi\hat{f}|1\rangle|^{2}}{\omega_{0}}$ (37)

with

$\langle 0|\hat{f}^{2}|0\rangle=\frac{\omega_{0}}{2E_{L}};\qquad\langle
0|\hat{f}|1\rangle=\sqrt{\frac{\omega_{0}}{2E_{L}}};\qquad\langle\hat{f}\rangle=-(E^{*}_{J}/E_{L})\sin\varphi,$
(38)

yielding

$\delta
E_{g}^{(e)}=E^{*}_{J}\frac{\,\omega_{0}}{4E_{L}}\cos\varphi+\frac{E_{J}^{*2}}{2E_{L}}\sin^{2}\varphi.$
(39)

The first term is identified as the renormalization of $E_{J}^{*}$, while the
second is an energy induced in the inductance by the supercurrent in the
junction. It causes also a shift in $f$, $\delta
f=-(E^{*}_{J}/E_{L})\sin\varphi$ . As we see, the relative strength of the
corrections reflects the ratio $\omega_{0}/E_{J}^{*}$.

Let us turn to odd parity. The correction $\delta\Omega$ to the unperturbed
binding energy $\Omega_{0}=\sqrt{2E_{J}}\sin\frac{\varphi}{2}$ is computed
from

$\frac{\delta\Omega}{2\sqrt{\Omega_{0}}}=\delta
A+\left(\frac{\sqrt{2E_{J}}\cos\frac{\varphi}{2}}{2}\right)^{2}\frac{\langle
0|\hat{f}|1\rangle^{2}}{\sqrt{\Omega_{0}+\omega_{0}}-\sqrt{\Omega_{0}}},$ (40)

where the first and second term give the first-order and the second-order
correction, respectively.

The correction to $A$ is given as

$\delta A=\sqrt{2E_{J}}\left(\cos\frac{\varphi}{2}\frac{\delta
f}{2}-\sin\frac{\varphi}{2}\frac{\hbar\omega_{0}}{16E_{L}}\right).$ (41)

With this,

$\delta\Omega=-E_{J}\sin^{2}\frac{\varphi}{2}\frac{\hbar\omega_{0}}{4E_{L}}-\sin^{2}\varphi\frac{E_{J}E^{*}_{J}}{E_{L}}+\frac{E_{J}^{2}\sin^{2}\varphi}{2E_{L}}\frac{\sqrt{1+\omega_{0}/\Omega_{0}}+1}{2}.$
(42)

Let us note that in the classical limit, $\omega_{0}\to 0$, the correction to
the ground state energy is reduced to the inductive energy,

$\delta E_{g}^{(o)}=\delta
E_{g}^{(e)}-\delta\Omega=-\frac{(E_{J}^{*}-E_{J})^{2}}{2E_{L}}\sin^{2}\varphi.$
(43)

This energy is zero if $E_{J}^{*}=E_{J}$ owing to the poisoning. In further
considerations, we concentrate on this case. There, the phase dependence of
$\delta E_{g}^{(o)}$ eventually defines the supercurrent, while it is a small
correction to the supercurrent otherwise.

With this, we reduce the correction to the odd ground state energy to the
following expression:

$\delta E_{g}^{(o)}=\frac{\omega^{2}_{0}}{8E_{L}}\frac{{\rm
cotan}^{2}\frac{\varphi}{2}}{\left(1+\sqrt{1+\omega_{0}/\Omega_{0}}\right)^{2}}.$
(44)

This formula does not look self-explaining and needs to be elaborated.

First of all, let us note that the corresponding current is finite in the
limit $\varphi\to 0$. In this limit, $\omega_{0}\gg\Omega_{0}\approx
E_{J}\varphi^{2}/2$ at any ratio $\omega_{0}/E_{J}$. Expanding in $\varphi$
till the first order, we find ($\varphi>0$)

$\delta
E_{g}^{(o)}\approx\frac{\omega_{0}^{2}}{8E_{L}\zeta}\left(1-\frac{\varphi}{\sqrt{\zeta}}\right);\qquad\zeta\equiv\frac{\omega_{0}}{2E_{J}}.$
(45)

So the current jumps at $\varphi=0$, the value of the half-jump is given by

$\frac{I_{\rm
hj}}{2e}=-\frac{E_{J}}{2}\sqrt{\frac{E_{J}}{E_{L}}}\sqrt{\frac{\omega_{0}}{2E_{L}}}\ll
E_{J}.$ (46)

We note that the current is negative at positive phase. In general, the
minimum of $\delta E_{g}^{(o)}$ is achieved at $\varphi=\pi$.

Let us address the limiting cases. If $\omega_{0}\gg E_{J}$, the current
reduces to

$\frac{I(\varphi)}{2e}=-\frac{E_{J}}{8}\frac{\omega_{0}}{E_{L}}\sin(\varphi).$
(47)

This looks like a renormalization of $E_{J}$ by the small oscillator’s
effective impedance ${\omega_{0}}/{E_{L}}=\sqrt{L/C}$. The current jump is
small and can be neglected at this scale, except for very small $\varphi$.

The opposite limiting case $\omega_{0}\ll E_{J}$ is trickier. Here, the
current of the order of $I_{\rm hj}$ is concentrated in a narrow interval of
phase $\simeq\sqrt{\zeta}$,

$I(\varphi)=-|I_{\rm hj}|f(\varphi/2\sqrt{\zeta});\qquad
f(x)\equiv\frac{1}{\sqrt{1+x^{2}}\left(x+\sqrt{1+x^{2}}\right)^{2}}\enspace{\rm
with}\enspace f(0)=1,\enspace f(x)\to 1/(8x^{3})\enspace{\rm at}\enspace
x\to\infty.$ (48)

We plot the phase dependence of the current for several values of $\zeta$ in
Fig. 2.

Figure 2:  Small impedance, single oscillator case. The superconducting
current in the odd-parity single-channel case versus phase for several values
of $\zeta={\omega_{0}}/{(2E_{J})}$ corresponding to Eq. (44).

###  IV.2 Ohmic impedance

Let us now turn to the more complex case of an ohmic impedance. The impedance
model in use comprises a capacitance, a resistor, and an inductance connected
in parallel. The resulting admittance at imaginary frequency reads:

$|\omega|{\cal Y}(\omega)=C\omega^{2}+\frac{\omega}{R}+1/L,$ (49)

corresponding to the real-frequency admittance $Y(\omega)=-i\omega
C+1/R+1/(-i\omega L)$. The capacitance is here to cut the ohmic part at high
frequencies $\omega_{H}=1/RC$. It is important to understand that, in order to
realize of a good phase bias at zero frequency, there must be a sufficiently
small inductance here, $E_{L}\gg E^{*}_{J}$. This cuts the ohmic part at low
frequency $\omega_{L}=R/L$.

We concentrate on the case of small impedance $Re^{2}\ll 1$. In this case, we
can assume small fluctuations, $\varphi(\tau)=\varphi+f(\tau)$,
$f(\tau)\ll\pi$, and expand similar to Eqs. (35), (36) replacing
$\hat{f}\rightarrow f(\tau)$. The spectrum of the fluctuations is given by

$\langle|f_{\omega}|^{2}\rangle=4e^{2}{\cal Z}(\omega)/|\omega|.$ (50)

Let us first compute the correction to the ground state energy of the even
parity state. We expand
$\cos(\varphi(\tau))=\cos(\varphi)-\sin(\varphi)f(\tau)-\cos(\varphi)f^{2}(\tau)/2$.
The correction consists of two terms: the first one is the first-order
correction $\simeq f^{2}$, the second one is the second-order correction in
$f$ at zero frequency – it comes from the current induced in the inductance.
It also causes a shift in $f$, $\delta f=-(E^{*}_{J}/E_{L})\sin\varphi$ . As a
result, we find

$\displaystyle\delta E_{g}^{(e)}$ $\displaystyle=$ $\displaystyle
E^{*}_{J}\cos(\varphi)\frac{\langle
f^{2}\rangle}{2}-\frac{E_{J}^{*2}}{2E_{L}}\sin^{2}\varphi$ (51)
$\displaystyle=$ $\displaystyle
E^{*}_{J}\cos(\varphi)\frac{2e^{2}}{\pi}\int_{0}^{\infty}\frac{d\omega}{\omega}Z(\omega)-\frac{E_{J}^{*2}}{2E_{L}}\sin^{2}\varphi$
$\displaystyle=$
$\displaystyle\alpha\ln\left[\frac{\omega_{H}}{\omega_{L}}\right]E^{*}_{J}\cos(\varphi)-\frac{E_{J}^{*2}}{2E_{L}}\sin^{2}\varphi$
$\displaystyle=$
$\displaystyle-2\alpha\ln\left[\frac{\omega_{H}}{\omega_{L}}\right]E^{*}_{J}\sin^{2}\frac{\varphi}{2}-\frac{E_{J}^{*2}}{2E_{L}}\sin^{2}\varphi+{\rm
const},$

where we define the dimensionless impedance as $\alpha\equiv 2Re^{2}/\pi$. The
resulting answer is somewhat similar to Eq. (39): the first term describes a
small renormalization of $E^{*}_{J}$, while the second term gives the
inductive energy.

Let us address the corrections to the ground state energy at odd parity. For
concreteness, we assume $s=1$, $\sin(\varphi/2)>0$. We expand $A$ in small
deviations:

$A=\sqrt{2E_{J}}\left(\sin\frac{\varphi}{2}-\cos\frac{\varphi}{2}\frac{f(\tau)}{2}-\sin\frac{\varphi}{2}\frac{f^{2}(\tau)}{8}\right).$
(52)

Since we work in the lowest order in fluctuations, we need the self-energy
part in the second order only, and the self-consistency equation reads

$\sqrt{\Omega}=\langle
A\rangle+\Sigma^{(2)}(-i\Omega);\qquad\Sigma^{(2)}(-i\Omega)=\int\frac{d\omega}{2\pi}\frac{\langle\langle|A|^{2}_{\omega}\rangle\rangle}{\sqrt{\Omega+i\omega}-\langle
A\rangle}.$ (53)

We expand it to the first order in the correction $\delta\Omega$,
$\Omega_{0}\equiv 2E_{J}\sin^{2}(\varphi/2)$ being the binding energy in
zeroth order. This yields

$\frac{\delta\Omega}{2\sqrt{\Omega}_{0}}=\delta\langle
A\rangle+\frac{E_{J}}{2}\cos^{2}\frac{\varphi}{2}\int\frac{d\omega}{2\pi}\frac{4e^{2}Z(\omega)}{|\omega|\left(\sqrt{\Omega_{0}+i\omega}-\sqrt{\Omega_{0}}\right)};\qquad\delta\langle
A\rangle=-\frac{\sqrt{2E_{J}}}{8}\sin\frac{\varphi}{2}\langle\langle
f^{2}\rangle\rangle+\frac{\sqrt{2E_{J}}}{2}\cos\frac{\varphi}{2}\delta f.$
(54)

It is instructive to identify two parts in this correction. One is
proportional to $\delta\langle A\rangle$, and we will call it coherent. The
second part is proportional to $\Sigma^{(2)}(-i\Omega)$ and is called
incoherent. The coherent correction presents a renormalization of $E_{J}$ and
a part of the inductive energy,

$\delta\Omega_{\rm{coh}}=-\alpha\ln\left[\frac{\omega_{H}}{\omega_{L}}\right]E_{J}\sin^{2}\frac{\varphi}{2}-\frac{E_{J}E^{*}_{J}}{E_{L}}\sin^{2}\varphi.$
(55)

If we compare this with the renormalization of $E^{*}_{J}$ in the even parity
state, we note that the renormalization in the odd parity sector is two times
smaller: a fact that we will recall when addressing the arbitrary impedance
case.

The analysis of the incoherent correction is complicated by the fact that it
depends both on $\omega_{L}$ and $\Omega_{0}$. The upper cut-off is not
relevant since the integral defining the correction converges at the maximum
of these two frequencies. If the impedance is not small, $\alpha\simeq 1$,
$\omega_{L}$ surely exceeds $\Omega_{0}$ as by definition of phase bias
$E_{L}\gg E_{J}$. If the impedance is small,
$\omega_{L}/\Omega_{0}\simeq\alpha/(E_{J}/E_{L})$ is the ratio of two small
numbers and can be large or small. We give the general formulas below.

For our impedance model,

$\Sigma^{(2)}(-i\Omega)=\frac{E_{J}}{2}\cos^{2}\frac{\varphi}{2}\int\frac{d\omega}{2\pi}\frac{4e^{2}R}{\left(|\omega|+\omega_{L}\right)\left(\sqrt{\Omega+i\omega}-\sqrt{\Omega_{0}}\right)}.$
(56)

This is a tricky integral: the denominator contains a delta-functional
contribution that is easy to miss at $\Omega\to\Omega_{0}$. To see this, let
us transform the denominator assuming $\Omega\to\Omega_{0}$,
$\omega\ll\Omega_{0}$

$\frac{1}{\sqrt{\Omega+i\omega}-\sqrt{\Omega_{0}}}=\frac{\sqrt{\Omega+i\omega}+\sqrt{\Omega_{0}}}{\Omega-\Omega_{0}+i\omega}\approx\frac{2\sqrt{\Omega_{0}}}{\Omega-\Omega_{0}+i\omega}\to
2\pi\sqrt{\Omega}_{0}{\rm sgn}(\delta\Omega)\delta(\omega).$ (57)

With this,

$\Sigma^{(2)}(-i\Omega)={E_{J}}\cos^{2}\frac{\varphi}{2}\sqrt{\Omega_{0}}2e^{2}L+\frac{E_{J}}{2}\cos^{2}\frac{\varphi}{2}\;{\rm
p.v.}\int\frac{d\omega}{2\pi}\frac{4e^{2}R}{\left(|\omega|+\omega_{L}\right)\left(\sqrt{\Omega_{0}+i\omega}-\sqrt{\Omega_{0}}\right)},$
(58)

where p.v. indicates the principal-part integration. This gives the following
correction:

$\delta\Omega_{\rm{incoh}}=\frac{E^{2}_{J}}{2E_{L}}\sin^{2}\varphi+\alpha
E_{J}\cos^{2}\frac{\varphi}{2}X(x);\qquad
x\equiv\frac{\omega_{L}}{2E_{J}\sin^{2}\frac{\varphi}{2}}.$ (59)

$X(x)$ is an elementary but rather complex function of $x$. Its asymptotes are

$X(x)\approx\ln\left(\frac{4e}{x}\right)\quad{\rm at}\quad x\to 0;\qquad
X(x)\approx\pi\sqrt{\frac{2}{x}}\quad{\rm at}\quad x\to\infty.$ (60)

A compact expression is as follows $(z>1)$:

$X((z^{2}-z^{-2})/2)=4z\left(\frac{{\rm arccoth}(z)}{z^{2}+1}+\frac{{\rm
arccot}(1/z)}{z^{2}-1}-\frac{\pi z}{2(z^{4}-1)}\right).$ (61)

The correction to the ground state energy at odd parity defines the
superconducting current that survives poisoning. We concentrate on the single-
level case, $E_{J}=E^{*}_{J}$. In this case, the parts of theinductive energy
in Eqs. (51), (55), (59) cancel each other as they did in the previous
subsection, and the answer reads

$\delta E_{g}^{(o)}=\delta E_{g}^{(e)}-\delta\Omega=-\alpha E_{J}\left({\rm
ln}\left[\frac{\omega_{H}}{\omega_{L}}\right]\sin^{2}\frac{\varphi}{2}+\cos^{2}\frac{\varphi}{2}X\left(\frac{\omega_{L}}{2E_{J}\sin^{2}\frac{\varphi}{2}}\right)\right).$
(62)

Let us note that the first term gives a negative supercurrent. If it dominates
(this is the case if the logarithmic factor is really bigger than 1), the
Josephson junction becomes a $\pi$ junction. As for the second term, the
energy correction is negative and reaches 0 both at $\varphi=0$ and
$\varphi=\pi$. So if the second term is sufficiently big in comparison with
the first one, there are two equivalent energy minima: one at
$\varphi_{0},0<\varphi_{0}<\pi$, and the symmetric one at $2\pi-\varphi_{0}$.

We observe the current jump at $\varphi\to 0$. In this case, $x\to\infty$
irrespective the ratio $\omega_{L}/E_{J}$ and the value of half-jump is given
by

$\frac{I_{{\rm hj}}}{2e}=-\pi\alpha
E_{J}\sqrt{\frac{E_{J}}{\omega_{L}}}=-\frac{E_{J}}{2}\sqrt{\frac{E_{J}}{E_{L}}}\sqrt{\frac{\omega_{L}}{E_{L}}}.$
(63)

Note the similarity with Eq. (46) if associating $\omega_{L}$ and
$\omega_{0}/\sqrt{2}$.

Let us look at the limiting cases. If $\omega_{L}\gg E_{J}$, the argument of
$X$ is always bigger than 1 and the limiting form of the correction reads

$\delta E_{g}^{(o)}=\alpha E_{J}\left(-{\rm
ln}\left[\frac{\omega_{H}}{\omega_{L}}\right]\sin^{2}\frac{\varphi}{2}-\pi\sqrt{\frac{E_{J}}{\omega_{L}}}\cos\frac{\varphi}{2}\sin\varphi\right).$
(64)

The first term dominates over the whole range of $\varphi$ except $\varphi\to
0$, the second one gives rise to a jump at $\varphi\to 0$.

In the opposite case, $\omega_{L}\ll E_{J}$, the argument of $X$ is swept from
large to small values within a narrow interval of the phase,
$\varphi\simeq\sqrt{\omega_{L}/E_{J}}$. The current is concentrated in this
interval,

$I(\varphi)=-|I_{\rm hj}|f(\varphi/\sqrt{2\omega_{L}/E_{J}});\qquad
f(x)\equiv-\frac{\sqrt{2}}{\pi}X^{\prime}(x^{-2})/x^{3}\enspace{\rm
with}\enspace f(0)=1,\enspace f(x)\to\sqrt{2}/(\pi x)\enspace{\rm at}\enspace
x\to\infty.$ (65)

Note again some similarities with Eq. (48): The current scales with $I_{\rm
hj}$ and is concentrated at small $\varphi$, though the scaling functions are
different.

Figure 3:  Small ohmic impedance. Left: The superconducting current in the
odd-parity single-channel case versus phase. Right: corresponding energies. We
set $\ln(\omega_{H}/\omega_{L})=5$, the curve labels correspond to various
values of $\omega_{L}/E_{J}$.

We plot the resulting current in Fig. 3 for $\omega_{L}/E_{J}=10,2,0.4,0.16$.
For this example, we set $\ln(\omega_{H}/\omega_{L})=5$. In this case, the
energy minimum is at $\varphi=\pi$ for all $\omega_{L}/E_{J}$ in use.

###  IV.3 Beyond perturbation theory

It follows from the previous analysis that $\delta\Omega\propto\varphi$ at
$\varphi\to 0$ while $\Omega_{0}\propto\varphi^{2}$. Therefore, $\delta\Omega$
becomes comparable with $\Omega_{0}$ at sufficiently small $\varphi$ and the
perturbation theory should not work. However, there is a simple workaround
common for both small-impedance models we consider.

To start with, let us understand when it is plausible to disregard the higher-
order diagrams in the self-energy. The second-order $\Sigma^{(2)}$ can be
estimated as $E_{J}\langle f^{2}\rangle/\sqrt{{\rm max}(\Omega,\omega)}$,
where $\omega$ here is either $\omega_{0}$ or $\omega_{L}$ from previous
subsections. The fourth-order $\Sigma^{(4)}$ contains 3 $G_{0}$ and two
A-correlators, therefore it can be estimated as $E_{J}^{2}(\langle
f^{2}\rangle)^{2}/(\sqrt{{\rm max}(\Omega,\omega)})^{3}$. They differ by a
factor $E_{J}\langle f^{2}\rangle/{\rm max}(\Omega,\omega)$ that has to be
small to neglect higher orders. (The estimation of the relative value of
$\Sigma^{(3)}$ and all odd order diagramms is even smaller since we consider
$\varphi\ll\pi$.)

Let us now assume $\Omega\ll\omega$ at $\varphi=0$. Then $\Sigma^{(2)}$ does
not depend on $\varphi$, and $\Omega$ can be estimated as
$\Omega\simeq(\Sigma^{(2)})^{2}\simeq E_{J}^{2}(\langle
f^{2}\rangle)^{2}/\omega$. Since $\langle f^{2}\rangle\simeq\omega/E_{L}$, one
obtains $\Omega/\omega\simeq(E_{J}/E_{L})^{2}\ll 1$ under the assumption of
phase bias. The higher orders can be neglected if $E_{J}\langle
f^{2}\rangle/\omega\simeq E_{J}/E_{L}\ll 1$, that is, under the same
assumption.

The phase dependence is incorporated through $\langle
A\rangle=\sqrt{E_{J}/2}\,s\varphi$. The resulting expression for $\Omega$ can
be presented in the following form:

$\sqrt{\Omega}=\sqrt{E_{J}/2}\left(s\varphi+\varphi_{c}\right)$ (66)

with $\varphi_{c}\simeq\sqrt{E_{J}/E_{L}}\sqrt{\omega/E_{L}}$. The value of
$\varphi_{c}$ can be also obtained from Eqs. (46), (63), $\varphi_{c}=(|I_{\rm
hj}|/2e)/E_{J}$.

We note that in the narrow interval of the phase
$-\varphi_{c}<\varphi<\varphi_{c}$ the bound state is present for both values
of $s$ (see Fig. 4): a qualitative change as compared to the situation without
interaction.

Figure 4:  Bound states in the narrow interval of the phases for $s=\pm 1$
(blue and red). Thin curves: no interaction.

##  V Arbitrary impedance: phase bias

In this Section, we address larger impedances $\alpha\simeq 1$. It is
essential that we keep the conditions of phase bias, $E_{L}\gg
E^{*}_{J},E_{J}$. In this case, the low cut-off frequency $\omega_{L}$ is much
bigger than $E^{*}_{J},E_{J}$ and does not change upon renormalization of
$E_{J}$. This renormalization is thus always finite. This implies that there
is no Schmid transition at phase bias, contrary to what happens at current
bias. This situation will be addressed in the next Section.

For the consideration, we will need the averages of the phase exponents over
the fluctuations induced by the environment: let us cite these results here
($\beta$ is arbitrary at the moment, we will need $\beta=1,1/2$):

$\displaystyle\langle e^{-i\beta\varphi(\tau)}\rangle$ $\displaystyle=$
$\displaystyle e^{-i\beta(\varphi+\langle
f\rangle)}e^{-\frac{\beta^{2}}{2}\langle\langle
f^{2}\rangle\rangle};\qquad\langle\langle
f^{2}\rangle\rangle=\int\frac{d\omega}{2\pi}\frac{4e^{2}{\cal
Z}(\omega)}{|\omega|};\qquad\langle
e^{-i\beta(\varphi(0)-\varphi(\tau))}\rangle=e^{-\frac{\beta^{2}}{2}U(\tau)};$
(67) $\displaystyle U(\tau)$ $\displaystyle\equiv$
$\displaystyle\langle(f(0)-f(\tau))^{2}\rangle=4\int\frac{d\omega}{2\pi}\sin^{2}(\omega\tau/2)\frac{4e^{2}{\cal
Z}(\omega)}{|\omega|}\quad{\rm with}\quad U(0)=0,\enspace
U(\tau\to\infty)=2\langle f^{2}\rangle;$ (68) $\displaystyle\langle
e^{-i\beta(\varphi(0)+\varphi(\tau))}\rangle$ $\displaystyle=$ $\displaystyle
e^{-i2\beta(\varphi+\langle
f\rangle)}e^{\frac{\beta^{2}}{2}U(\tau)-2\beta^{2}\langle\langle
f^{2}\rangle\rangle}.$ (69)

In the leading approximation, $\langle\langle
f^{2}\rangle\rangle=2\alpha\ln(\omega_{H}/\omega_{L})$ and
$U(\tau)=4\alpha\ln(\omega_{H}\tau)$ at
$\omega_{H}^{-1}\ll\tau\ll\omega_{L}^{-1}$. One needs to use these formulas
with care since the ratio $\omega_{H}/\omega_{L}\sim(E_{C}/E_{L})\alpha^{-2}$
depends on $\alpha$ as well.

Let us address the even parity sector first. The phase-dependent part of the
energy is given by the renormalized (and strongly suppressed) Josephson
energy,

$E_{g}^{(e),1}=-\tilde{E}^{*}_{J}\cos\varphi\enspace{\rm
with}\enspace\frac{\tilde{E}^{*}_{J}}{E^{*}_{J}}=e^{-\langle\langle
f^{2}\rangle\rangle/2}\simeq\left[\frac{\omega_{L}}{\omega_{H}}\right]^{\alpha}.$
(70)

There is also a second-order contribution that is not exponentially
suppressed.

$E_{g}^{(e),2}=-E^{*}_{J}\int d\tau\left[e^{-U(\tau)/2}-e^{-\langle\langle
f^{2}\rangle\rangle}\right].$ (71)

For strong suppression ($\tilde{E}_{J}^{*}/E_{J}^{*}\ll 1$), the integrand at
$\omega_{H}^{-1}\ll\tau\ll\omega_{L}^{-1}$ is $(\omega_{H}\tau)^{-2\alpha}$.
Thus, if $\alpha>1/2$, the integral converges at the upper cut-off
$\omega_{H}\tau\simeq 1$ whereas it converges at the lower cut-off otherwise.
So we estimate (skipping the prefactors):

$\displaystyle E_{g}^{(e),2}\simeq\frac{E^{*2}_{J}}{\omega_{H}}\ \ \ {\rm if}\
\alpha>1/2;\qquad
E_{g}^{(e),2}\simeq\frac{E^{*2}_{J}}{\omega_{L}}\left[\frac{\omega_{L}}{\omega_{H}}\right]^{2\alpha}\simeq\frac{\tilde{E}^{*2}_{J}}{\omega_{L}}\
\ \ {\rm if}\ \alpha<1/2.$ (72)

The correction at $\alpha<1/2$ is related to the one discussed in Ref. [4]. We
conclude that the second-order contribution dominates at $\alpha>1/2$.

We turn to odd parity and keep terms up to the second order in the self-
consistency equation (cf. Eq. (53)),

$\sqrt{\Omega}=\langle A\rangle+\Sigma^{(2)}(-i\Omega).$ (73)

The average $A$ is strongly suppressed,

$\langle A\rangle=\sqrt{2\tilde{E}_{J}}s\sin\frac{\varphi}{2}\enspace{\rm
with}\enspace\frac{\tilde{E}_{J}}{E_{J}}=e^{-{\langle
f^{2}\rangle}/{4}}\simeq\left[\frac{\omega_{L}}{\omega_{H}}\right]^{\alpha/2}.$
(74)

As was noted when considering the small impedance limit, the suppression of
the odd-parity Josephson coupling is two times weaker than the one of the
even-parity Josephson coupling (see Eq. (70)).

At $\alpha<1$, $\Sigma^{(2)}$ can be neglected in zeroth approximation. The
superconducting current is given by the renormalized expression

$\frac{I(\varphi)}{2e}=(\tilde{E}^{*}_{J}-\tilde{E}_{J})\sin\varphi=\left(E^{*}_{J}e^{-\frac{\langle
f^{2}\rangle}{2}}-E_{J}e^{-\frac{\langle f^{2}\rangle}{4}}\right)\sin\varphi.$
(75)

Since $E_{J}$ is less suppressed than $E_{J}^{*}$, there is a chance that the
supercurrent in the odd state is bigger in magnitude than that in the even
state. The current at phases $\varphi\in[0,\pi]$ is negative in this case. For
the single-channel case $E^{*}_{J}=E_{J}$ the current is always negative in
this phase interval.

However, the second order term $\Sigma^{(2)}(-i\Omega)$ can become important
since it has a phase-independent part that controls the position of the bound
state with respect to the continuum edge. This leads to a variety of bound
regimes A-F listed in Fig. 5. To estimate $\Sigma^{(2)}$, we can concentrate
on the phase-independent terms that are not exponentially suppressed at large
$\alpha$ and, since $\Omega\ll\omega_{L}$, disregard their $\Omega$
dependence. This yields

$\Sigma^{(2)}={E_{J}}\int_{0}^{\infty}\frac{d\tau}{\sqrt{\pi\tau}}\langle\langle
e^{i\varphi(0)/2}e^{-i\varphi(\tau)/2}\rangle\rangle.$ (76)

The integrand at $\omega_{H}^{-1}\ll\tau\ll\omega_{L}^{-1}$ is
$(\omega_{H}\tau)^{-\alpha/2}\tau^{-1/2}$. The integral thus converges at the
upper cut-off if $\alpha>1$ and at the lower cut-off if $\alpha<1$. The
estimations for $\Sigma^{(2)}$ then read

$\displaystyle\Sigma^{(2)}\simeq\frac{E_{J}}{\sqrt{\omega_{H}}}\ \ \ {\rm if}\
\alpha>1;\;$ (77)
$\displaystyle\Sigma^{(2)}\simeq\frac{E_{J}}{\sqrt{\omega_{L}}}\left[\frac{\omega_{L}}{\omega_{H}}\right]^{\alpha/2}\simeq\frac{\tilde{E}_{J}}{\sqrt{\omega_{L}}}\
\ \ {\rm if}\ \alpha<1.$ (78)

If $\langle A\rangle$ dominates but the finite $\Sigma^{(2)}$ is taken into
account, there is a small current jump at $\varphi\to 0$, namely $|I_{\rm
hj}|/(2e)\simeq\tilde{E_{J}}\sqrt{\tilde{E}_{J}/\omega_{L}}$ corresponding to
the level structure described by Eq. (66) with
$\varphi_{c}=\Sigma^{(2)}/\sqrt{\tilde{E}_{J}/2}\simeq\sqrt{\tilde{E}_{J}/\omega_{L}}$
(regime B). Eventually, $\varphi_{c}$ and the relative value of the jump
$|I_{\rm hj}|/(2e\tilde{E}_{J})$ decrease with increasing $\alpha$ at
$1/\ln(\omega_{H}/\omega_{L})<\alpha<1$.

We compare $\langle A\rangle$ and $\Sigma^{(2)}$ to establish that the latter
dominates at
$\alpha>2(1+\ln(\omega_{L}/E_{J})/\ln(\omega_{H}/\omega_{L}))\equiv\alpha_{c}>2$.
At $\alpha_{c}>\alpha>1$ we are still in the regime B, with the only
difference that $\Sigma^{(2)}$ saturates at the value $\simeq
E_{J}/\sqrt{\omega_{H}}$ and therefore $\varphi_{c}$ and the relative current
jump increase with increasing $\alpha$, namely $\varphi_{c}\simeq
E_{J}/\sqrt{\omega_{H}\tilde{E_{J}}}$.

At $\alpha\approx\alpha_{c}$, where $\sqrt{2\tilde{E}_{J}}=\Sigma^{(2)}$, an
important transition (regime C) takes place: the bound state is present at any
phase for both $s=\pm 1$ (regime D). As a consequence, the odd parity state
becomes stable upon an adiabatic sweep of the phase. The bound state energies
are given by

$\Omega=\left(s\sqrt{2\tilde{E}_{J}}\sin\frac{\varphi}{2}+\Sigma^{(2)}\right)^{2}.$
(79)

The resulting superconducting current at a given $s$ becomes $4\pi$ periodic,
a phenomenon similar to that signifying the presence of Majorana modes [5].
The $2\pi$ periodicity is restored upon relaxation to the lowest energy state
within the odd sector.

At $\alpha>\alpha_{c}$, the bound state energy $\Omega\simeq
E_{J}^{2}/\omega_{H}$ hardly depends on the phase and $\alpha$. The remaining
phase dependence results in a strongly suppressed supercurrent (regime E).

$I(\varphi)/e\sim
E_{J}\sqrt{\frac{\tilde{E}_{J}}{\omega_{H}}}s\cos\frac{\varphi}{2}.$ (80)

Despite being suppressed, this supercurrent is parametrically bigger than the
one $\propto\tilde{E}_{J}^{*}$ in the even parity state.

The resulting supercurrent in this regime is $4\pi$ periodic provided $s$ is
conserved during the measurement time. The relaxation to the energetically
favourable value of $s$ would restore $2\pi$ periodicity. More involved
research is required to estimate a typical relaxation time in this situation:
let us explain why. The terms in the Hamiltonian that break the conservation
of $s$ requires an asymmetry of the quasiparticle wave function in the
right/left lead. Such asymmetry is not manifested under standard assumptions
of semiclassical theory of superconductivity, where the difference between the
leads effectively averages out. The asymmetry may arise owing to disorder-
induced fluctuations of the superconducting order parameter [6] or mesoscopic
fluctuations of Andreev scattering [7]. This provides a small factor
associated with the ratio of the electronic wavelength and the spread of an
Andreev state over the leads. In addition to this symmetry violation, the
relaxation requires an inelastic processes to bridge the energy difference
between the split states with different $s$. This inelastic process may arise
from the electromagnetic environment. With this, a typical relaxation time may
be even longer than the life time of the odd parity state under consideration.

Figure 5:  The various bound regimes that can occur in the odd-parity
Josephson junction. Qualitatively, they can be described with Eq. (79). A:
only one superposition gives rise to a bound states (realized at no
interaction, $\alpha=0$). B: two bound states in a finite phase interval near
$\varphi=0,2\pi$. (cf. Fig 4) C: transition between B and D. D: both bound
states are present at all phases, the resulting superconducting current is
$4\pi$ periodic. E: the splitting of two bound states is much smaller than
their average phase-independent energy. F: Two bound states are degenerate,
and bear no phase-dependence.

##  VI Arbitrary impedance: current bias

In contrast with the phase bias situation, there is no built-in low-energy
cut-off at current bias: the lower cut-off is determined by the renormalized
Josephson energy.

Let us remind how this works by addressing the even parity sector. The
renormalized $\tilde{E}^{*}_{J}$ is given by the same formula as in the
previous section, see Eq. (70), where $\omega_{L}$ is estimated as
$\tilde{E}^{*}_{J}$. (This estimation is fine in the region of interest
$\alpha\simeq 1$. At $\alpha\ll 1$, $\tilde{E}^{*}_{J}\alpha$ would be a more
accurate approximation. Yet taking into account the difference between these
two estimations exceeds the accuracy of the approach.) With this,

$\frac{\tilde{E}^{*}_{J}}{E^{*}_{J}}=\left[\frac{E^{*}_{J}}{\omega_{H}}\right]^{\frac{\alpha}{1-\alpha}},$
(81)

that is, $\tilde{E}_{J}^{*}$ vanishes at the Schmid transition, $\alpha=1$.

Let us now turn to the odd parity sector and assume $E^{*}_{J}=E_{J}$. (The
analysis of a many-channel situation $E^{*}_{J}\gg E_{J}$ leads to a more
involved situation where both $\tilde{E}^{*}_{J}$ and $\tilde{E}_{J}$ may play
a role of the low-energy cut-off. This analysis is beyond the scope of the
present paper.)

To start with, let us concentrate on the interval $\alpha<1$. In this case,
the lower cut-off can be unambiguously identified as $\tilde{E}_{J}$. Applying
the results of the previous section (see Eq. (74)), we obtain

$\frac{\tilde{E}_{J}}{E_{J}}\simeq\left[\frac{E_{J}}{\omega_{H}}\right]^{\alpha/2}\quad\to\quad\frac{\tilde{E}_{J}}{E_{J}}=\left[\frac{E_{J}}{\omega_{H}}\right]^{\frac{\alpha}{2-\alpha}}.$
(82)

The estimation of $\Sigma^{(2)}$ with the help of Eq. (78) gives
$\Sigma^{(2)}\simeq\sqrt{\tilde{E}_{J}}$. We see that, by contrast to the
situation of phase bias, the first- and second-order contributions are of the
same order of magnitude, as well as all higher orders. Therefore the accuracy
of the method does not allow to predict the phase dependence of the energy,
nor if bound states persist for both values of $s$, as it was the case under
phase bias (regimes B-C-D). However, we still note and use the difference in
the renormalizations for phase-dependent ($\tilde{E}_{J}$) and phase-
independent parts of $\sqrt{\Omega}$.

This becomes important at $\alpha>1$, where in accordance with Eq. (77)
$\Sigma^{(2)}$ does not depend on the low cut-off anymore and saturates at the
value $\simeq E_{J}/\sqrt{\omega_{H}}$. As to the phase-dependent part, it
further decreases with increasing $\alpha$. This brings us to regime E: almost
degenerate bound states described by Eq. (80). In this case, the lower cut-off
is $\simeq E_{J}\sqrt{\frac{\tilde{E}_{J}}{\omega_{H}}}$ rather than
$\simeq\tilde{E}_{J}$, which yields

$\frac{\tilde{E}_{J}}{E_{J}}=\left[\frac{E_{J}}{\omega_{H}}\right]^{\frac{3\alpha}{4-\alpha}}.$
(83)

We see that $\tilde{E}_{J}$ vanishes at $\alpha=4$. This is the Schmid
transition point for a half of the Cooper pair charge corresponding to a
$4\pi$ periodicity in phase.

The bound state is completely degenerate with respect to $s$ at $\alpha>4$
(regime F). If we recall the spin, we observe the realization of a 4-fold
degeneracy.

The results of the two last Sections are summarized in Fig. 6.

Figure 6:  The summary of the results at arbitrary impedance. Renormalized
Josephson energies $\tilde{E}_{J}$ at both parities. Capital letters indicate
the bound regimes at odd parity, dotted vertical lines separate the regimes.
Left: phase bias, $\tilde{E}_{J}$ never vanishes. The separating regime $C$
occurs at $\alpha=\alpha_{c}$. Right: current bias. Schmid transition is at
$\alpha=1$ for even parity and at $\alpha=4$ for odd one. The renormalization
law at odd parity changes at $\alpha=1$.

## References

  * Catelani  _et al._ [2012] G. Catelani, S. E. Nigg, S. M. Girvin, R. J. Schoelkopf, and L. I. Glazman, Decoherence of superconducting qubits caused by quasiparticle tunneling, Phys. Rev. B 86, 184514 (2012).
  * Schön and Zaikin [1990] G. Schön and A. Zaikin, Quantum coherent effects, phase transitions, and the dissipative dynamics of ultra small tunnel junctions, Physics Reports 198, 237 (1990).
  * Abrikosov  _et al._ [1975] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinskii,  _Methods of quantum field theory in statistical physics_ (Dover, New York, NY, 1975).
  * Hekking and Glazman [1997] F. W. J. Hekking and L. I. Glazman, Quantum fluctuations in the equilibrium state of a thin superconducting loop, Phys. Rev. B 55, 6551 (1997).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Superconducting proximity effect and Majorana fermions at the surface of a topological insulator, Phys. Rev. Lett. 100, 096407 (2008).
  * Larkin and Ovchinnikov [1972] A. Larkin and Y. Ovchinnikov, Density of states in inhomogeneous superconductors, Sov. Phys. JETP 34, 1144 (1972).
  * Feigel’man and Skvortsov [2012] M. V. Feigel’man and M. A. Skvortsov, Universal broadening of the Bardeen-Cooper-Schrieffer coherence peak of disordered superconducting films, Phys. Rev. Lett. 109, 147002 (2012).