# Josephson quantum mechanics at odd parity

Manuel Houzet, Julia Meyer, Yuli Nazarov

► To cite this version:

Manuel Houzet, Julia Meyer, Yuli Nazarov. Josephson quantum mechanics at odd parity. Physical Review B, 2024, 110 (2), pp.L020504. 10.1103/PhysRevB.110.L020504. hal-05003791

HAL Id: hal-05003791

https://hal.science/hal-05003791v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Josephson quantum mechanics at odd parity

Manuel Houzet, Julia S. Meyer

Univ. Grenoble Alpes, CEA, Grenoble INP, IRIG, Pheliqs, 38000 Grenoble, France

Yuli V. Nazarov

Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

(Dated: September 7, 2023)

A Josephson junction may be in a stable odd parity state when a single quasiparticle is trapped in an Andreev bound state. Embedding such junction in an electromagnetic environment gives rise to a special quantum mechanics of superconducting phase that we investigate theoretically. Our analysis covers several representative cases, from the lifting of the supercurrent quench due to quasiparticle poisoning for a low ohmic impedance of the environment, to a Schmid transition in a current-biased junction that for odd parity occurs at four times bigger critical impedance. For intermediate impedances, the supercurrent in the odd state is higher than in the even one.

The energy of a tunnel junction between two superconducting leads depends periodically on the difference of superconducting phases of the two, in short, on the phase. This is the celebrated Josephson effect [1]: the phase dependence of this energy gives rise to a persistent superconducting current between the leads. Later, it has been understood that the phase becomes a quantum-fluctuating variable if a Josephson junction is embedded in an electromagnetic circuit [2]. Earlier studies concentrated on a dissipative electromagnetic environment and were essential for establishing the modern theory of dissipative quantum mechanics [3, 4]. A highlight of this research was the prediction of the Schmid transition [5]: the vanishing of the Josephson energy at a critical value of the circuit impedance $R$, $2e^2R/\pi\hbar \equiv \alpha = 1$. While this prediction is theoretically indisputable, the controversy concerning its experimental verification [6, 7] may have been resolved recently [8]. The further development of Josephson quantum mechanics evolved from dissipative circuits to dissipationless Coulomb islands. The resulting Josephson-based superconducting qubits [9, 10] are at the frontline of modern quantum technology applications.

There is something to add to this well-established field. In fact, the Josephson energy is related to Andreev bound states (ABS) in the junction [11] and does depend on their occupation. Of the two equal-weight superpositions with respect to the right/left leads in which a quasiparticle may be in, only one gives rise to a bound state. Owing to parity conservation in superconductors [12], a state with a single quasiparticle trapped in the lowest ABS (the odd parity ground state) is stable despite having a bigger energy than the state without quasiparticles (the even parity state). Physically, the parity can only be relaxed if a stray quasiparticle from a lead comes to the junction and annihilates the trapped one. Since the concentration of the quasiparticles in the leads is vanishingly small at low temperatures, the lifetime of the odd parity ground state is macroscopically long: lifetimes of several minutes have been measured [13]. We note that a single quasiparticle trapped in a spin-degenerate ABS eventually quenches the contribution of this level to the Josephson energy: this is called the quasiparticle poisoning and has been observed in [14]. When spin-

degeneracy is lifted (in finite-length junctions with spin-orbit coupling), the stability of these odd states provided the opportunity for a new kind of qubits: Andreev spin qubits, proposed in [15, 16] and realized in [17]. In recent years, there is an outburst of studies of ABS in superconducting nanostructures, including spectroscopically resolved ABS and odd parity ground states in a junction [18].

This makes it relevant to extend the Josephson quantum me-

<figure><img src="image_2.png" /><figcaption>FIG. 1. a. The odd parity Josephson junction. A single quasiparticle is trapped in the lowest Andreev level separated by $2E_J \sin^2 \frac{\varphi}{2} \ll \Delta$ from the edge of the continuous quasiparticle spectrum at the superconducting energy gap $\Delta$. In the bound state, the quasiparticle is in a certain superposition, $s = 1$, the anti-bound state corresponding to $s = -1$ (dashed curve) belongs to the continuous spectrum. b.-c. The Josephson quantum mechanics at odd parity: the odd parity Josephson junction is embedded in a linear electromagnetic environment with frequency-dependent impedance $Z(\omega)$ that causes quantum fluctuations of the phase. b. and c. correspond to phase and current bias, respectively.</figcaption></figure>
chanics to the case of a circuit embedding a Josephson junction in the odd parity ground state. Such quantum mechanics at odd parity should be quite distinct from the conventional one. For instance, for a short single-channel junction, quasiparticle poisoning is expected to completely quench the supercurrent [11]. Thus a naive and, as we will see, wrong expectation is that the junction is not present in the circuit at all. In our pivotal study, we consider a tunnel junction where the ABS are close to the superconducting gap edge, disregard weak spin-orbit interaction, and mainly concentrate on the inoperative single-level, single-junction case, see Fig. 1.

In this Letter, we provide a general description of Josephson quantum mechanics at odd parity revealing its intriguing mathematical structure. We also present the detailed analysis for three relevant cases. For low ohmic impedance, we demonstrate the incompleteness of supercurrent quenching and reveal a supercurrent jump at zero phase. For arbitrary ohmic impedance and phase bias, we establish a slower suppression of the Josephson energy in the odd state than in the even one: the supercurrent in the odd state thus becomes higher than in the even one, both remaining finite at any $\alpha$ as already shown in the even case [19]. At sufficiently large impedance, both right/left superpositions form a bound state. While their phase-dependence is suppressed upon increasing the impedance, their average binding energy tends to a constant. In addition to this, for arbitrary ohmic impedance and current bias, we encounter a Schmid transition at a higher value of the impedance than in the even state, namely, at $\alpha = 4$. The bound states persist for both superpositions and are degenerate for $\alpha > 4$. These predictions can be tested in forthcoming experiments.

Let us sketch here the general derivation: all details are provided in [20]. At even parity, the Hamiltonian describing a Josephson junction embedded in a general linear environment, see Fig. 1b, reads [21]

$$H_e = H_{\text{env}} - E_J^* \cos \hat{\varphi}, \quad (1)$$

where $H_{\text{env}}$ is a Hamiltonian of non-interacting bosons, the operator of the phase drop at the junction, $\hat{\varphi}$, consists of the phase bias $\varphi$ and a linear superposition of environmental bosons, and $E_J^*$ is the even-state Josephson energy. The coefficients in the superposition are chosen such as to reproduce the frequency-dependent impedance of the environment, $Z(\omega)$. An alternative description [22] employs a path integral over a variable $\varphi(\tau)$ defined in imaginary time. The action that defines the path weight reads

$$S = \sum_{\omega} \frac{|\omega|}{8e^2 Z(i|\omega|)} |\varphi(\omega)|^2 - E_J^* \int d\tau \cos \varphi(\tau), \quad (2)$$

$\varphi(\omega)$ being the Fourier transform of $\varphi(\tau)$.

To describe the odd parity situation, we first augment the Hilbert space with the states of a single quasiparticle to reduce it at a later stage of the derivation. Without the environment, this gives the binding energy $\Omega$, measured from the edge of

the continuum, in the following form:

$$\sqrt{\Omega} = s\sqrt{2E_J} \sin \frac{\varphi}{2}. \quad (3)$$

Here, $E_J$ is the Josephson energy associated with the lowest ABS: $E_J = E_J^*$ in the single-channel case, $E_J^* > E_J$ in general, and $s = \pm 1$ characterizes the superposition between right/left leads. Equation (3) with $s = \text{sign}(\sin \frac{\varphi}{2})$ reproduces the ABS dispersion in a short junction in the tunnel limit [23]. With the environment, the above relation is modified to a singular-value equation for a wave function $|\Phi\rangle$ in the environmental degrees of freedom, that involves square-roots of Hamiltonian-like operators,

$$\left( \sqrt{\Omega + \bar{H}} - s\sqrt{2E_J} \sin \frac{\hat{\varphi}}{2} \right) |\Phi\rangle = 0, \quad (4)$$

where $H = H_e - E_g^{(e)}$ with $E_g^{(e)}$ the ground state energy in the even parity sector.

The path integral approach is also non-trivial, bearing a similarity with the Green function treatment of a frozen disorder [24]. The key quantity is a propagator $G(\tau, \tau')$ defined in a rather standard way

$$G(\tau, \tau') = G_0(\tau - \tau') + \int d\tau_1 G_0(\tau - \tau_1) A(\tau_1) G(\tau_1, \tau'). \quad (5)$$

Here $A(\tau) \equiv s\sqrt{2E_J} \sin \frac{\varphi(\tau)}{2}$ plays the role of the disorder, $G_0(\tau) \equiv \Theta(\tau)/\sqrt{\pi\bar{\tau}}$ is the bare propagator arising from the reduction of quasiparticle continuum states. The disorder averaging should be done with the weight $e^{-S}$, that is, with respect to the even parity ground state. The averaged propagator is uniform, its Fourier component reads

$$\bar{G}(\omega) = \left( \sqrt{i\omega} - \langle A \rangle - \Sigma(\omega) \right)^{-1}, \quad (6)$$

the self-energy $\Sigma(\omega)$ being a sum of diagrams involving the correlators of $A(\tau)$ starting from the second order. Finally, the binding energy is found from

$$\sqrt{\Omega} = \langle A \rangle + \Sigma(-i\Omega). \quad (7)$$

Equations (4) and (7) demonstrate an involved structure of the resulting theory that is distinct from straightforward Hamiltonian or path-integral approaches. Nevertheless, we manage to get to experimentally verifiable predictions by using perturbation theory and renormalizations.

Let us start with the case of small ohmic impedance, $\alpha \ll 1$. For a concrete model, we add a capacitance and an inductance in parallel to the resistor $R$, $Z(\omega) = 1/(-i\omega C + 1/R + i/\omega L)$. This cuts the ohmic response both at high and low frequency, $\omega_H = 1/RC$ and $\omega_L = R/L$, respectively. The inductance providing the low cut-off is required in order to phase bias the junction, $E_J e^2 L \ll 1$. (The opposite regime may be addressed as in [19] for the even case.) We concentrate on the single-channel case of quasiparticle poisoning: the odd ground state energy $E_g^{(o)}$ does not depend on phase without
<figure><img src="image_3.png" /><figcaption>FIG. 2. a. The odd-parity supercurrent at small impedance. The curve labels are $\omega_L/E_J$, we set $\ln(\omega_H/\omega_L) = 5$. b. Bound states near zero phase for $s = \pm 1$. Here $\varphi_c = \pi\alpha\sqrt{E_J/\omega_L} \ll 1$. Dashed curves: no interaction.</figcaption></figure>

fluctuations. We aim to compute the phase-dependent correction $\delta E_g^{(o)}(\varphi)$ proportional to the fluctuations, which defines the supercurrent in the odd state.

A simple ad hoc estimation would be $\delta E_g^{(o)} \simeq \alpha E_J \cos \varphi$. While this may be a correct scale, the answer is more involved and interesting, see Fig. 2a. We note an extra dimensionless parameter $\omega_L/E_J$ that can be large or small provided $\alpha \ll 1$. We see that the current in the phase interval $\varphi \in [0, \pi]$ is negative: the minimum odd-parity Josephson junction energy corresponds to $\varphi = \pi$ rather than $\varphi = 0$. Let us note that this $\pi$-junction behaviour has a completely different origin than the one induced by magnetic correlations in the ground state of a superconductor-quantum dot-superconductor junction [25] or the one due to the continuum contribution that is left in the presence of poisoning when weak interactions and a finite length of the junction are taken into account [26]. At $\omega_L \gg E_J$,

$$ \frac{I(\varphi)}{2e} = -\frac{\alpha E_J}{2} \ln \left( \frac{\omega_H}{\omega_L} \right) \sin \varphi. \qquad (8) $$

The most interesting feature present for arbitrary ratios $\omega_L/E_J$ is the current jump at $\varphi = 0$, its half-value being

$$ \frac{I_{hj}}{2e} = -\pi\alpha E_J \sqrt{\frac{E_J}{\omega_L}}. \qquad (9) $$

At $\omega_L \ll E_J$, the supercurrent is concentrated at small $\varphi \simeq \sqrt{\omega_L/E_J}$ and reads

$$ I(\varphi) = -|I_{hj}| f(\varphi/\sqrt{2\omega_L/E_J}) \qquad (10) $$

with $f(0) = 1$ and $f(x) \to \sqrt{2}/\pi x$ at $x \to \infty$. The full expression for the monotonous function $f(x)$ is given in [20].

The current jump is associated with the fact that the perturbation theory formally ceases to hold at small $\varphi$. However, the answer beyond perturbations is really simple and shown in Fig. 2b: namely the binding energy is shifted such that the bound state reaches the continuum edge not at $\varphi = 0$, but at $\varphi = -s\varphi_c$ with $\varphi_c \equiv (|I_{hj}|/2e)/E_J$, i.e., the binding energy

<figure><img src="image_4.png" /><figcaption>FIG. 3. a. Critical currents at even and odd parity versus $\langle\langle\varphi^2\rangle\rangle$, Eq. (15). The odd parity current dominates at $\langle\langle\varphi^2\rangle\rangle > 4 \ln 2 \approx 2.8$. b. The bound regimes in the odd parity Josephson junction. A: only one superposition gives rise to a bound state ($\alpha = 0$). B: two bound states in a finite phase interval (cf. Fig. 2b). C: separatrix between B and D. D: two $4\pi$-periodic bound states are present at all phases. E: the splitting of the two bound states is much smaller than their average phase-independent energy. F: The two states $s = \pm 1$ with phase-independent energy are degenerate.</figcaption></figure>

is given by

$$ \sqrt{\Omega} = \sqrt{E_J/2} (s\varphi + \varphi_c). \qquad (11) $$

The shifts being opposite for $s = \pm 1$, this implies the presence of bound states for both superpositions in an interval $|\varphi| < \varphi_c$: this fact will become crucial for further analysis.

Let us turn to the case of an arbitrary impedance, $\alpha \simeq 1$, under conditions of phase bias. In this case, the low cut-off frequency is such that $\omega_L \gg E_J$ and does not change upon renormalization of $E_J, E_J^*$. The renormalization is thus finite at any $\alpha$: this implies that, as discussed in the even parity sector [19], no Schmid transition occurs under phase bias. While $E_J = E_J^*$ in the single-channel case, they renormalize differently. The renormalization can be computed using the relation $\langle e^{i\beta\varphi} \rangle = e^{i\beta\langle\varphi\rangle} e^{-\beta^2\langle\varphi^2\rangle/2}$, where $\langle\langle\varphi^2\rangle\rangle = \langle\varphi^2\rangle - \langle\varphi\rangle^2$, valid for Gaussian fluctuations of the phase. At even parity,

$$ \tilde{E}_J^* = E_J^* e^{-\langle\varphi^2\rangle/2} \simeq E_J^* (\omega_L/\omega_H)^\alpha. \qquad (12) $$

Here and further on the 'tilde' refers to renormalized quantities.

To understand the renormalization at odd parity, we keep terms up to the second order in the self-consistency equation (7),

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

The average $A$ is phase-dependent and strongly suppressed,

$$ \langle A \rangle = s \sqrt{2\tilde{E}_J} \sin \frac{\varphi}{2} \quad \text{with} \quad \frac{\tilde{E}_J}{E_J} = e^{-\frac{\langle\varphi^2\rangle}{4}} \simeq \left( \frac{\omega_L}{\omega_H} \right)^{\frac{\alpha}{2}}. \qquad (14) $$

This suppression is two times weaker than for even parity. The superconducting current in the odd state at $\alpha < 2$ reads

$$ \frac{I(\varphi)}{2e} = (\tilde{E}_J^* - \tilde{E}_J) \sin \varphi = E_J \left( e^{-\frac{\langle\varphi^2\rangle}{2}} - e^{-\frac{\langle\varphi^2\rangle}{4}} \right) \sin \varphi \qquad (15) $$
and is bigger than that at even parity at sufficiently large phase fluctuations, see Fig. 3a.

However, as far as the bound state spectrum is concerned, the second-order term $\Sigma^{(2)}(-i\Omega)$ can become important since it has a phase-independent part. This leads to a variety of bound regimes A-F listed in Fig. 3b. For estimates, we concentrate on the phase-independent terms in $\Sigma^{(2)}$ and, since $\Omega \ll \omega_L$, disregard the $\Omega$ dependence. This yields

$$ \Sigma^{(2)} = E_J \int_0^\infty \frac{d\tau}{\sqrt{\pi\tau}} \langle e^{i\varphi(0)/2} e^{-i\varphi(\tau)/2} \rangle. \quad (16) $$

The integrand at $\omega_H^{-1} \ll \tau \ll \omega_L^{-1}$ is proportional to $1/\tau^{1/2}(\omega_H \tau)^{\alpha/2}$. As a consequence, the integral converges at the lower cut-off if $\alpha < 1$ and at the upper cut-off if $\alpha > 1$. The estimations for $\Sigma^{(2)}$ thus read:

$$ \Sigma^{(2)} \approx \begin{cases} \tilde{E}_J / \sqrt{\omega_L}, & \alpha < 1, \\ E_J / \sqrt{\omega_H}, & \alpha > 1. \end{cases} \quad (17) $$

Comparing $\langle A \rangle$ at $\varphi = \pi$ and $\Sigma^{(2)}$, we observe that the latter dominates for $\alpha > 2[1+\ln(\omega_L/E_J)/\ln(\omega_H/\omega_L)] \equiv \alpha_c > 2$. This point (C) separates two different regimes. Now we can summarize the results.

At $\alpha < \alpha_c$, $\Sigma^{(2)}$ can be neglected in zeroth approximation. The superconducting current is given by Eq. (15). Starting from the case without fluctuations (regime A), we find that the addition of small phase-independent terms in Eq. (13) when fluctuations are weak leads to the coexistence of two bound states corresponding to the two superpositions $s = \pm 1$ in a small interval of phases around $\varphi = 0$ (regime B). At $\alpha > 1$, this interval grows with increasing $\alpha$ until an important transition (regime C) takes place at $\alpha_c$: two bound states are present at any phase. For $\alpha > \alpha_c$, both bound states are separated from the continuum by a gap (regime D). Thus the odd parity state becomes stable upon an adiabatic sweep of the phase. The bound energies are given by

$$ \Omega = \left( s \sqrt{2\tilde{E}_J} \sin \frac{\varphi}{2} + \Sigma^{(2)} \right)^2 . \quad (18) $$

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

At $\alpha$ slightly (by $\simeq 1/\ln(\omega_H/\omega_L) \ll 1$) exceeding $\alpha_c$, the binding energy $\Omega \simeq E_J^2/\omega_H \gg \tilde{E}_J$ hardly depends on the phase and $\alpha$ (regime E). The remaining phase dependence results in a strongly suppressed $4\pi$-periodic supercurrent

$$ \frac{I(\varphi)}{2e} \approx s \tilde{E}_J^{\text{eff}} \cos \frac{\varphi}{2}; \quad \tilde{E}_J^{\text{eff}} \approx \sqrt{\tilde{E}_J \Omega} \approx E_J \sqrt{\frac{\tilde{E}_J}{\omega_H}}. \quad (19) $$

Despite being suppressed, this supercurrent parametrically exceeds the one at even parity.

Let us now turn to the case of an arbitrary impedance at current bias, see Fig. 1c. In contrast with the phase bias situation,

<figure><img src="image_5.png" /><figcaption>FIG. 4. Renormalized Josephson energies $\tilde{E}_J^*$ (green) at even and $\tilde{E}_J$ (violet) at odd parity. Vertical dotted lines separate the bound regimes at odd parity indicated by capital letters. Left: phase bias, cf. Eqs. (12), (14), and (19); $\tilde{E}_J$ never vanishes. The separating regime C occurs at $\alpha = \alpha_c$. We plot $\tilde{E}_J^{\text{eff}}$ instead of $\tilde{E}_J$ at $\alpha > \alpha_c$. Right: current bias, cf. Eqs. (20), (21), and (22); the curves illustrate the suppression of $\tilde{E}_J$ as $\alpha$ increases, the Schmid transition where $\tilde{E}_J$ vanishes is at $\alpha = 1$ for even parity and at $\alpha = 4$ for odd parity. The renormalization law at odd parity changes at $\alpha = 1$. Note the different vertical scales in the left and right plot.</figcaption></figure>

there is no built-in low energy cut-off $\omega_L$: the renormalization has to be cut-off self-consistently by the renormalized Josephson energy.

Let us first reproduce the Schmid transition at even parity. The renormalized $\tilde{E}_J^*$ is given by the same Eq. (12), yet $\omega_L$ there has to be estimated as $\tilde{E}_J^*$. With this,

$$ \frac{\tilde{E}_J^*}{E_J^*} = \left( \frac{E_J^*}{\omega_H} \right)^{\frac{\alpha}{1-\alpha}}, \quad (20) $$

such that $\tilde{E}_J^*$ vanishes at the Schmid transition, $\alpha = 1$.

Let us next turn to the odd parity sector. 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 Eq. (14), we thus obtain

$$ \frac{\tilde{E}_J}{E_J} = \left( \frac{E_J}{\omega_H} \right)^{\frac{\alpha}{2-\alpha}}. \quad (21) $$

The estimation of $\Sigma^{(2)}$ with the help of Eq. (17) gives $\Sigma^{(2)} \simeq \sqrt{\tilde{E}_J}$. In contrast with the phase-biased case, the first- and second-order contributions are of the same order of magnitude, as well as all higher orders. So in the current bias case, 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$ (regimes B-C-D).

However, we still may notice and use the difference in the renormalizations of phase-dependent and phase-independent parts of $\sqrt{\Omega}$ depending on the value of $\alpha$. This becomes important at $\alpha > 1$ where, in accordance with Eq. (17), the self-energy $\Sigma^{(2)}$ does not depend on the low-energy cut-off anymore and saturates at the value $\simeq E_J/\sqrt{\omega_H}$. As to the phase-dependent part, it will further decrease with increasing
$\alpha$. This brings us to regime E: the almost degenerate bound state associated with the supercurrent described by Eq. (19). In this case, the renormalization of $E_J$ is cut off by $\tilde{E}_J^{\text{eff}}$ of Eq. (19), rather than $\tilde{E}_J$. This yields

$$ \frac{\tilde{E}_J}{E_J} = \left( \frac{E_J}{\omega_H} \right)^{\frac{3\alpha}{4-\alpha}} ; \quad \frac{\tilde{E}_J^{\text{eff}}}{E_J} = \left( \frac{E_J}{\omega_H} \right)^{\frac{\alpha+2}{4-\alpha}} . \qquad (22) $$

Therefore $\tilde{E}_J, \tilde{E}_J^{\text{eff}}$ vanish at $\alpha = 4$. This is the new Schmid transition point for half of the Cooper pair charge, indeed corresponding to $4\pi$-periodicity in phase of the supercurrent. At $\alpha > 4$, the phase-independent bound state is completely degenerate with respect to $s$ (regime F). Recalling the quasiparticle spin, we thus predict the realization of 4-fold degeneracy for the trapped quasiparticle.

In conclusion, we have formulated the Josephson quantum mechanics for a junction in the odd parity state. The nontrivial structure of the theory is encapsulated in Eqs. (4) and (7). We concentrated on the single-channel case and predicted the lifting of the supercurrent quench due to quasiparticle poisoning at small $\alpha$. The residual supercurrent is given by Eqs. (8)-(10). Furthermore, we have addressed the case of arbitrary impedance both at phase and current bias. The supercurrent at odd parity is less suppressed by quantum fluctuations and may dominate over the one at even parity. The presence of various bound regimes complicates the renormalization. For current bias, we predict a Schmid transition at $\alpha = 4$ and four-fold degenerate bound states at higher impedances.

YVN acknowledges support from the Université Grenoble Alpes for an extended stay in Grenoble during which most of the presented work was performed. MH and JSM acknowledge funding from the Plan France 2030 through the project NISQ2LSQ ANR-22-PETQ-0006.

[1] B. Josephson, Possible new effects in superconductive tunneling, Physics Letters 1, 251 (1962).

[2] A. J. Leggett, Macroscopic Quantum Systems and the Quantum Theory of Measurement, Progress of Theoretical Physics Supplement 69, 80 (1980).

[3] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).

[4] U. Weiss, *Quantum Dissipative Systems*, Series in modern condensed matter physics (World Scientific, 1999).

[5] A. Schmid, Diffusion and localization in a dissipative quantum system, Phys. Rev. Lett. 51, 1506 (1983).

[6] A. Murani, N. Bourlet, H. le Sueur, F. Portier, C. Altimiras, D. Esteve, H. Grabert, J. Stockburger, J. Ankerhold, and P. Joyez, Absence of a dissipative quantum phase transition in Josephson junctions, Phys. Rev. X 10, 021003 (2020).

[7] D. Subero, O. Maillet, D. S. Golubev, G. Thomas, J. T. Peltonen, B. Karimi, M. Marín-Suárez, A. L. Yeyati, R. Sánchez, S. Park, and J. P. Pekola, Bolometric detection of coherent Josephson coupling in a highly dissipative environment (2022), arXiv:2210.14953 [cond-mat.mes-hall].

[8] R. Kuzmin, N. Mehta, N. Grabon, R. A. Mencia, A. Bureshtein, M. Goldstein, and V. E. Manucharyan, Observation of the Schmid-Bulgadaev dissipative quantum phase transition (2023), arXiv:2304.05806 [quant-ph].

[9] Y. Makhlin, G. Schön, and A. Shnirman, Quantum-state engineering with Josephson-junction devices, Rev. Mod. Phys. 73, 357 (2001).

[10] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).

[11] C. W. J. Beenakker, Universal limit of critical-current fluctuations in mesoscopic Josephson junctions, Phys. Rev. Lett. 67, 3836 (1991).

[12] D. V. Averin and Y. V. Nazarov, Single-electron charging of a superconducting island, Phys. Rev. Lett. 69, 1993 (1992).

[13] D. J. van Woerkom, A. Geresdi, and L. P. Kouwenhoven, One minute parity lifetime of a NbTiN Cooper-pair transistor, Nature Physics 11, 547 (2015).

[14] M. Zgirski, L. Bretheau, Q. Le Masne, H. Pothier, D. Esteve, and C. Urbina, Evidence for long-lived quasiparticles trapped in superconducting point contacts, Phys. Rev. Lett. 106, 257003 (2011).

[15] N. M. Chtchelkatchev and Y. V. Nazarov, Andreev quantum dots for spin manipulation, Phys. Rev. Lett. 90, 226806 (2003).

[16] C. Padurariu and Y. V. Nazarov, Theoretical proposal for superconducting spin qubits, Phys. Rev. B 81, 144519 (2010).

[17] M. Hays, V. Fatemi, D. Bouman, J. Cerrillo, S. Diamond, K. Serniak, T. Connolly, P. Krogstrup, J. Nygård, A. L. Yeyati, A. Geresdi, and M. H. Devoret, Coherent manipulation of an Andreev spin qubit, Science 373, 430 (2021).

[18] R. Aguado, A perspective on semiconductor-based superconducting qubits, Applied Physics Letters 117, 240501 (2020).

[19] 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).

[20] See Supplemental Material, which includes Refs. [28, 29], at [URL will be inserted by publisher] for the details of the derivations.

[21] A. O. Caldeira and A. J. Leggett, Influence of dissipation on quantum tunneling in macroscopic systems, Phys. Rev. Lett. 46, 211 (1981).

[22] 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).

[23] 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).

[24] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinskii, Methods of quantum field theory in statistical physics (Dover, New York, NY, 1975).

[25] L. N. Bulaevskii, V. V. Kuzii, and A. A. Sobyanin, Superconducting system with weak coupling to the current in the ground state, JETP Lett. 25, 290 (1977).

[26] P. D. Kurilovich, V. D. Kurilovich, V. Fatemi, M. H. Devoret, and L. I. Glazman, Microwave response of an Andreev bound state, Phys. Rev. B 104, 174517 (2021).

[27] 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).

[28] A. Larkin and Y. Ovchinnikov, Density of states in inhomogeneous superconductors, Sov. Phys. JETP 34, 1144 (1972).

[29] 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).
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.

### 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_{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. \quad (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^{th}$ tunneling channel of the junction ($1 \le m \le N_{ch}$, where $N_{ch}$ is the number of channels). The quasiparticle energies are $\varepsilon_k = \sqrt{\xi_k^2 + \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'm\sigma} \frac{t_m}{L} a_{km\sigma}^{\dagger} c_{k'm\sigma} + \text{H.c.} \quad (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_{km-\sigma}^{\dagger}$, 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'} \left( \alpha_{km\sigma}^{\dagger} \gamma_{k'm\sigma} - \gamma_{km\sigma}^{\dagger} \alpha_{k'm\sigma} \right). \quad (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} \quad (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_{\text{full}} = \sum_{kms\sigma} \varepsilon_k \beta_{kms\sigma}^{\dagger} \beta_{kms\sigma} - \frac{t_m}{L} \sin \frac{\varphi}{2} \sum_{kk'ms\sigma} s \beta_{kms\sigma}^{\dagger} \beta_{k'ms\sigma} \quad (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_{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'} |k\rangle\langle k'|. \qquad (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'} \psi_{k'}' \qquad (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}. \qquad (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_J^{(m)}} s \sin \frac{\varphi}{2}, \qquad (9)$$

where $E_J^{(m)} \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_J^{(m)}$. 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^* \ge E_J$ in general and $E_J^* = E_J$ for a single-channel situation where the contribution of other channels can be neglected.

## 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_e = H_{env} - E_J^* \cos \hat{\varphi}. \qquad (10)$$

We use a bosonic description of the linear environment,

$$H_{env} = \sum_q \omega_q b_q^\dagger b_q, \qquad \hat{\varphi} = \varphi + \sum_q \lambda_q (b_q + b_q^\dagger). \qquad (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,

$$\mathrm{Re}Z(\omega) = \frac{\pi\omega}{4e^2} \sum_q \lambda_q^2 \delta(\omega - \omega_q), \qquad (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_e - E_g^{(e)}, \quad H|0\rangle = 0. \qquad (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_{odd}$ is obtained by combining the Hamiltonian of Eq. (6) and $H$, the binding energy $\Omega$ being an eigenvalue of this Hamiltonian,

$$H_{odd} = H + \sum_k \frac{\xi_k^2}{2\Delta} |k\rangle\langle k| - s \sin \frac{\varphi}{2} \frac{t_m}{L} \sum_{k,k'} |k\rangle\langle k'|, \qquad (14)$$

$$0 = (\Omega + H_{odd})|\Psi\rangle. \qquad (15)$$
With all the energy shifts we made, the ground state energy at odd parity is given by

$$E_g^{(o)} = E_g^{(e)} + \Delta - \Omega. \qquad (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 + \hat{H}} - \hat{A})|\Phi\rangle; \qquad \hat{A} \equiv s\sqrt{2E_J} \sin \frac{\hat{\varphi}}{2}. \qquad (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 $\hat{H}$ is a non-negatively defined operator, and $\hat{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$).

## 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(-S)$, where the action $S(\varphi(\tau))$ comprises a quadratic part representing the linear environment and a part representing the Josephson energy,

$$S = S_{\text{env}} - E_J^* \int d\tau \cos \varphi(\tau), \qquad (18)$$

$$S_{\text{env}} = \sum_{\omega} \frac{\omega}{8e^2 Z(\omega)} |\varphi_{\omega}|^2. \qquad (19)$$

Here, $Z(\omega)$ is the impedance at imaginary frequency. For the simplest circuit of a resistor $R$ and a capacitor $C$ in parallel, $Z^{-1}(\omega) = R^{-1} \text{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,

$$\mathrm{Ex} \equiv \langle 0 | \frac{1}{\sqrt{\Omega + \hat{H}} - \hat{A}} | 0 \rangle = \infty. \qquad (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,

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

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

$$\frac{1}{\sqrt{\Omega + \hat{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}, \qquad (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,

$$ \mathrm{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. \quad (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 \rightarrow \langle A(\tau_1) A(\tau_2) \rangle. \quad (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')$ defined by a rather standard relation,

$$ G(\tau, \tau') = G_0(\tau - \tau') + \int d\tau_1 G_0(\tau - \tau_1) A(\tau_1) G(\tau_1, \tau'). \quad (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')$.

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

$$ \mathrm{Ex} = \int d\tau e^{-\Omega t} \bar{G}(\tau); \quad \bar{G}(\tau) = \langle G(\tau_0 + \tau, \tau_0) \rangle. \quad (26) $$

The previous reasoning implies that the average propagator grows with increasing time difference, $\bar{G}(\tau) \propto e^{\Omega_b \tau}$. The integral in (26) thus diverges at $\Omega < \Omega_b$ and converges at $\Omega > \Omega_b$. This identifies $\Omega_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_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)}}. \quad (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), \quad (28) $$

$\Omega_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$,

$$ G_0(\omega) \rightarrow \bar{G}_0 \equiv \frac{1}{\sqrt{i\omega - \langle A \rangle}}, \quad (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

$$ \Sigma^{(2)}(\tau) = \langle\langle A(0)A(\tau) \rangle\rangle \bar{G}_0(\tau), \quad (30) $$

$$ \Sigma^{(3)}(\tau) = \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), \quad (31) $$

$$ \Sigma^{(4)}(\tau) = \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) \times (\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). \quad (32) $$

Here, the double angle brackets $\langle\langle \dots \rangle\rangle$ denote the cumulant of all variables between the brackets.
<figure><img src="image_6.png" /><figcaption>FIG. 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 Σ up to 4th order.</figcaption></figure>

## SMALL IMPEDANCE

### 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_f^2 + \frac{E_L}{2} f^2 = \hbar \omega_0 \left( \hat{a}^\dagger \hat{a} + \frac{1}{2} \right), \quad (33)$$

$f$ being the addition to the phase, $\hat{a}$ being the annihilation operator of the oscillation. Here, $E_L = (4e^2L)^{-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). \quad (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

$$ \cos(\hat{\varphi}) \rightarrow \cos(\varphi) - \sin(\varphi)\hat{f} - \cos(\varphi)\frac{\hat{f}^2}{2}, \qquad (35) $$

$$ \sin \frac{\hat{\varphi}}{2} \rightarrow \sin \frac{\varphi}{2} + \cos \frac{\varphi}{2} \frac{\hat{f}}{2} - \sin \frac{\varphi}{2} \frac{\hat{f}^2}{8}. \qquad (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} \qquad (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, \qquad (38) $$

yielding

$$ \delta E_g^{(e)} = E_J^* \frac{\omega_0}{4E_L} \cos \varphi + \frac{E_J^{*2}}{2E_L} \sin^2 \varphi. \qquad (39) $$

The first term is identified as the renormalization of $E_J^*$, while the second is an energy induced in the inductance by the super-current 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}}, \qquad (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). \qquad (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}. \qquad (42) $$

Let us note that in the classical limit, $\omega_0 \rightarrow 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. \qquad (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_0^2}{8E_L} \frac{\cotan^2 \frac{\varphi}{2}}{\left(1 + \sqrt{1 + \omega_0 / \Omega_0}\right)^2}. \qquad (44) $$

This formula does not look self-explaining and needs to be elaborated.
<figure><img src="image_7.png" /><figcaption>FIG. 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).</figcaption></figure>

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); \quad \zeta \equiv \frac{\omega_0}{2E_J}. \qquad (45) $$

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

$$ \frac{I_{hj}}{2e} = -\frac{E_J}{2} \sqrt{\frac{E_J}{E_L}} \sqrt{\frac{\omega_0}{2E_L}} \ll E_J. \qquad (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 \omega_0}{8 E_L} \sin(\varphi). \qquad (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_{hj}$ is concentrated in a narrow interval of phase $\simeq \sqrt{\zeta}$,

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

We plot the phase dependence of the current for several values of $\zeta$ in Fig. 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|\mathcal{Y}(\omega) = C\omega^2 + \frac{\omega}{R} + 1/L, \qquad (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} \to f(\tau)$. The spectrum of the fluctuations is given by

$$ \langle |f_{\omega}|^2 \rangle = 4e^2 \mathcal{Z}(\omega) / |\omega|. \qquad (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 $\approx 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

$$ \begin{aligned} \delta E_g^{(e)} &= E_J^* \cos(\varphi) \frac{\langle f^2 \rangle}{2} - \frac{E_J^{*2}}{2E_L} \sin^2 \varphi \\ &= 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 \\ &= \alpha \ln \left[ \frac{\omega_H}{\omega_L} \right] E_J^* \cos(\varphi) - \frac{E_J^{*2}}{2E_L} \sin^2 \varphi \\ &= -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 + \text{const}, \end{aligned} \qquad (51) $$

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). \qquad (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); \quad \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}}. \qquad (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| (\sqrt{\Omega_0 + i\omega} - \sqrt{\Omega_0})}; \quad \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. \qquad (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_{\text{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. \qquad (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}{(|\omega| + \omega_L)(\sqrt{\Omega} + i\omega - \sqrt{\Omega_0})}. \qquad (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}\operatorname{sgn}(\delta\Omega)\delta(\omega). \quad (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} \text{p.v.} \int \frac{d\omega}{2\pi} \frac{4e^2 R}{(|\omega| + \omega_L)(\sqrt{\Omega_0 + i\omega} - \sqrt{\Omega_0})}, \quad (58) $$

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

$$ \delta\Omega_{\text{incoh}} = \frac{E_J^2}{2E_L} \sin^2 \varphi + \alpha E_J \cos^2 \frac{\varphi}{2} X(x); \quad x \equiv \frac{\omega_L}{2E_J \sin^2 \frac{\varphi}{2}}. \quad (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 \text{at } x \to 0; \quad X(x) \approx \pi\sqrt{\frac{2}{x}} \quad \text{at } x \to \infty. \quad (60) $$

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

$$ X((z^2 - z^{-2})/2) = 4z \left( \frac{\operatorname{arccoth}(z)}{z^2 + 1} + \frac{\operatorname{arccot}(1/z)}{z^2 - 1} - \frac{\pi z}{2(z^4 - 1)} \right). \quad (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 the inductive 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( \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). \quad (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_{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}}. \quad (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( -\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). \quad (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 \approx \sqrt{\omega_L/E_J}$. The current is concentrated in this interval,

$$ I(\varphi) = -|I_{hj}| f(\varphi/\sqrt{2\omega_L/E_J}); \quad f(x) \equiv -\frac{\sqrt{2}}{\pi} X'(x^{-2})/x^3 \quad \text{with} \quad f(0) = 1, \quad f(x) \to \sqrt{2}/(\pi x) \quad \text{at} \quad x \to \infty. \quad (65) $$

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

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.
<figure><img src="image_8.png" /><figcaption>FIG. 3: Small ohmic impedance. Left: The superconducting current in the odd-parity single-channel case versus phase. Right: corresponding energies. We set ln(ω<sub>H</sub>/ω<sub>L</sub>) = 5, the curve labels correspond to various values of ω<sub>L</sub>/E<sub>J</sub>.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 4: Bound states in the narrow interval of the phases for s = ±1 (blue and red). Thin curves: no interaction.</figcaption></figure>

Beyond perturbation theory

It follows from the previous analysis that δΩ ∝ ϕ at ϕ → 0 while Ω₀ ∝ ϕ². Therefore, δΩ becomes comparable with Ω₀ at sufficiently small ϕ 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 Σ⁽²⁾ can be estimated as E<sub>J</sub>⟨f²⟩/√max(Ω, ω), where ω here is either ω₀ or ω_L from previous subsections. The fourth-order Σ⁽⁴⁾ contains 3 G₀ and two A-correlators, therefore it can be estimated as E<sub>J</sub>²⟨(f²)²⟩/(√max(Ω, ω))³. They differ by a factor E<sub>J</sub>⟨f²⟩/max(Ω, ω) that has to be small to neglect higher orders. (The estimation of the relative value of Σ⁽³⁾ and all odd order diagramms is even smaller since we consider ϕ ≪ π.)

Let us now assume Ω ≪ ω at ϕ = 0. Then Σ⁽²⁾ does not depend on ϕ, and Ω can be estimated as Ω ≃ (Σ⁽²⁾)² ≃ E<sub>J</sub>²((f²)²)/ω. Since ⟨f²⟩ ≃ ω/E<sub>L</sub>, one obtains Ω/ω ≃ (E<sub>J</sub>/E<sub>L</sub>)² ≪ 1 under the assumption of phase bias. The higher orders can be neglected if E<sub>J</sub>⟨f²⟩/ω ≃ E<sub>J</sub>/E<sub>L</sub> ≪ 1, that is, under the same assumption.

The phase dependence is incorporated through ⟨A⟩ = √E<sub>J</sub>/2 sϕ. The resulting expression for Ω can be presented in the following form:

$$
\sqrt{\Omega} = \sqrt{E_J/2} (s\varphi + \varphi_c) \qquad (66)
$$

with $\varphi_c \approx \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_{hj}|/2e)/E_J$.

We note that in the narrow interval of the phase −ϕ_c < ϕ < ϕ_c the bound state is present for both values of s (see Fig. 4): a qualitative change as compared to the situation without interaction.
## 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$):

$$ \langle e^{-i\beta\varphi(\tau)} \rangle = e^{-i\beta(\varphi+\langle f \rangle)} e^{-\frac{\beta^2}{2}\langle f^2 \rangle}; \quad \langle \langle f^2 \rangle \rangle = \int \frac{d\omega}{2\pi} \frac{4e^2\mathcal{Z}(\omega)}{|\omega|}; \quad \langle e^{-i\beta(\varphi(0)-\varphi(\tau))} \rangle = e^{-\frac{\beta^2}{2}U(\tau)}; \quad (67) $$

$$ U(\tau) \equiv \langle (f(0) - f(\tau))^2 \rangle = 4 \int \frac{d\omega}{2\pi} \sin^2(\omega\tau/2) \frac{4e^2\mathcal{Z}(\omega)}{|\omega|} \quad \text{with} \quad U(0) = 0, U(\tau \to \infty) = 2\langle f^2 \rangle; \quad (68) $$

$$ \langle e^{-i\beta(\varphi(0)+\varphi(\tau))} \rangle = e^{-i2\beta(\varphi+\langle f \rangle)} e^{\frac{\beta^2}{2}U(\tau)-2\beta^2\langle f^2 \rangle}. \quad (69) $$

In the leading approximation, $\langle f^2 \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 \quad \text{with} \quad \frac{\tilde{E}_J^*}{E_J^*} = e^{-\langle f^2 \rangle / 2} \approx \left[ \frac{\omega_L}{\omega_H} \right]^\alpha. \quad (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 f^2 \rangle} \right]. \quad (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):

$$ E_g^{(e),2} \approx \frac{E_J^{*2}}{\omega_H} \quad \text{if } \alpha > 1/2; \quad E_g^{(e),2} \approx \frac{E_J^{*2}}{\omega_L} \left[ \frac{\omega_L}{\omega_H} \right]^{2\alpha} \approx \frac{\tilde{E}_J^{*2}}{\omega_L} \quad \text{if } \alpha < 1/2. \quad (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). \quad (73) $$

The average $A$ is strongly suppressed,

$$ \langle A \rangle = \sqrt{2\tilde{E}_J s} \sin \frac{\varphi}{2} \quad \text{with} \quad \frac{\tilde{E}_J}{E_J} = e^{-\langle f^2 \rangle / 4} \approx \left[ \frac{\omega_L}{\omega_H} \right]^{\alpha/2}. \quad (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. \quad (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. \quad (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

$$ \Sigma^{(2)} \approx \frac{E_J}{\sqrt{\omega_H}} \quad \text{if } \alpha > 1; \quad (77) $$

$$ \Sigma^{(2)} \approx \frac{E_J}{\sqrt{\omega_L}} \left[ \frac{\omega_L}{\omega_H} \right]^{\alpha/2} \approx \frac{\tilde{E}_J}{\sqrt{\omega_L}} \quad \text{if } \alpha < 1. \quad (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_{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_{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. \quad (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}. \quad (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.

## 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$.
<figure><img src="image_10.png" /><figcaption>FIG. 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.</figcaption></figure>

$\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}}, \qquad (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} \approx \left[ \frac{E_J}{\omega_H} \right]^{\alpha/2} \quad \rightarrow \quad \frac{\tilde{E}_J}{E_J} = \left[ \frac{E_J}{\omega_H} \right]^{\frac{\alpha}{2-\alpha}}. \qquad (82) $$

The estimation of $\Sigma^{(2)}$ with the help of Eq. (78) gives $\Sigma^{(2)} \approx \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{\tilde{\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 $\approx 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 $\approx E_J\sqrt{\tilde{E}_J/\omega_H}$ rather than $\approx \tilde{E}_J$, which yields

$$ \frac{\tilde{E}_J}{E_J} = \left[ \frac{E_J}{\omega_H} \right]^{\frac{3\alpha}{4-\alpha}}. \qquad (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.

<footnote>[1] 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).</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 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$.</figcaption></figure>

[2] 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).

[3] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinskii, *Methods of quantum field theory in statistical physics* (Dover, New York, NY, 1975).

[4] 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).

[5] 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).

[6] A. Larkin and Y. Ovchinnikov, Density of states in inhomogeneous superconductors, Sov. Phys. JETP **34**, 1144 (1972).

[7] 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).