## Phase-Dependent Damping Rate of Josephson Plasmons in the Nonequilibrium Steady State

Takanobu Jujo ∗

Division of Materials Science, Graduate School of Science and Technology, Nara Institute of Science and Technology, Ikoma, Nara 630-0101, Japan

We calculate the damping rate of Josephson plasmons due to thermally excited quasiparticles in a nonequilibrium steady state. It is known that the experimentally observed dependence of the damping rate on the phase difference behaves in the opposite way to that of the theory. It is shown that this problem can be solved by calculating the damping rate by incorporating changes in the electronic state at the resonant frequency. When the phase difference is large, the broadening of the one-particle state is large due to the presence of a dc current and brings about the effective thermalization, so that the nonequilibrium correction becomes small. Conversely, when the phase difference is small, the nonequilibrium correction becomes large and causes a smaller damping rate due to a reduction in the effective temperature. This dependence of the nonequilibrium correction on the phase difference results in a damping rate that is consistent with experimental results.

## 1. Introduction

The Josephson current, which depends on the phase difference between two superconductors in the tunnel junction, has a non-dissipative current and a dissipative current, 1, 2) and there is a so-called cosine term problem in the latter. The problem is that there is an inconsistency between theory and experiment about the sign of the phase-dependent term of the dissipative current proportional to the cosine of the phase difference. 3, 4) Theory shows that the sign of this term is positive, 5, 6) but experiments on the damping rate of the Josephson plasmon show that the sign of this term is negative. 7-9)

In recent years, this term has been studied in relation to the energy relaxation in superconducting qubits. 10,11) The experiment carried out at very low temperatures showed that in the limit of small excitation frequencies, the sign of the cosine term is

∗ E-mail address: jujo@ms.naist.jp

positive, 12) as in theory. As in this study experiments on superconducting qubits are mainly aimed at extremely low temperatures where thermally excited quasiparticles can be neglected.

On the other hand, in previous experiments in which the negative sign of the cosine term was observed, thermally excited quasiparticles are thought to be responsible for the damping rate of Josephson plasmons. This is because these experiments were conducted at relatively high temperatures ( T &gt; 0 5 . T c with T and T c being the temperature and the transition temperature, respectively), and the theoretical calculations were also performed for the case where the thermal excitation gives a finite value of the damping rate. For this reason, the previous cosine term problem at relatively high temperatures, rather than at low temperatures where superconducting qubits work, is considered to be unsolved.

Up until now, theoretical calculations for the cosine term problem have focused only on the sign of the cosine term. 13-16) Experiments have shown that as the phase difference becomes smaller, the damping rate of Josephson plasmons decreases. Therefore, rather than extracting only the cosine term and examining its sign, it is necessary to examine the dependence of the entire term in the dissipative current on the phase difference.

In this study, we consider that at the resonance frequency of Josephson plasmons, the one-particle state of electrons is in a nonequilibrium steady state due to the external field at this frequency. We calculate the damping rate based on this electronic state. A finite phase difference occurs under the dc current which changes the resonance frequency. The presence of this phase difference causes a broadening of the peak around the gap edge in the density of state. Therefore the nonequilibrium correction differs depending on the magnitude of the phase difference, and the resulting variation of the effective temperature (this is known as the Eliashberg effect ´ 17) ) leads to the phase-dependent damping rate due to quasiparticle excitations.

The structure of this paper is as follows. In Sect. 2, we give a formulation based on the quasiclassical approximation in the nonequilibrium steady state, and derive the expression of the damping rate from the current. Section 3 gives the results of the numerical calculations based on this formulation. First, we show the calculated results of the damping rate without correction terms, and show that it gives the same behavior as in the previous studies. Next, we give the results of the damping rate in a nonequilibrium steady state, and show that the nonequilibrium correction plays a major role in the dependence of the damping rate on the phase difference. Section 4 provides a discussion

of these results.

## 2. Formulation

## 2.1 The expression of the current

We give the expression of the current from which the damping rate is derived. This quantity is experimentally obtained from the impedance measurements, 7) and the expression of the impedance is given as follows. (We set /planckover2pi1 = c = 1 with c the velocity of light.) By introducing the capacitance C = /epsilon1S/ 4 πd ( d and S are the thickness and the area of the junction, respectively), we obtain Cdω A 2 ω + J ω + J ext ω = 0 from the Maxwell equation, 18) in which J ω and J ext ω are the induced and external currents, respectively 19) ( A ω is the vector potential with its frequency ω ). With use of J ω = -CdA ω ( ω 2 φ -iωγ ω ), which is obtained below (Sect. 2.3), the impedance ( Z ω ) is written as Z ω = -iωdA /J ω ext ω = ( iω/C / ω ) ( 2 -ω 2 φ + iωγ ω ). Here, ω 2 φ = ω 2 J cos( φ L -φ R ) with ω J = √ 2 eI /C 0 and I 0 = ( π ∆ 2 / eR N )tanh(∆ 2 / T ) being the plasma frequency 20) and the critical current 21) in the Josephson junction, respectively. ∆ e iφ L and ∆ e iφ R are the superconducting order parameters of two superconductors consisting the Josephson junction (the magnitude of the gap ∆ is assumed to be equal for two superconductors), and R N is the tunnel resistance in the normal state. γ ω is the damping rate and its expression is given in Sect. 2.3.

The expression of the current J ω is obtained by using the quasiclassical Green functions. The following relation holds from the kinetic equations for quasiclassical Green functions: 22)

$$\int d \epsilon \text{Tr} [ \hat { \tau } _ { 3 } \omega \hat { g } ^ { L K } _ { \epsilon + \omega / 2, \epsilon - \omega / 2 } ] & = \iint d \epsilon d \epsilon _ { 1 } \text{Tr} [ \hat { \tau } _ { 3 } \hat { \Sigma } ^ { L + } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } \hat { g } ^ { L K } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } + \hat { \tau } _ { 3 } \hat { \Sigma } ^ { L K } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } \hat { g } ^ { L - } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } \\ & \quad - \hat { g } ^ { L + } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } \hat { \Sigma } ^ { L K } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } \hat { \tau } _ { 3 } - \hat { g } ^ { L K } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } \hat { \Sigma } ^ { L - } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } \hat { \tau } _ { 3 } ]. \\. \tau \cdot \tilde { \lambda } \cdot \tilde { \lambda } \tau \quad. \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cdot \tilde { \lambda } \cd.$$

ˆ g L +( -) and ˆ g LK are the retarded (advanced) and Keldysh quasiclassical Green functions, respectively, and Σ ˆ L +( -) and Σ ˆ LK are the self-energies. The superscripts L and R are indices representing two superconductors consisting the Josephson junction. ˆ τ 3 = ( 1 0 0 -1 ), ˆ means a matrix in the Nambu space, and Tr( ) indicates taking the trace. This leads · · to the following relation between the current ( J t ) and the temporal variation of the

charge ( Q t ):

$$\text{charge} \, ( Q _ { t } ) \colon \\ \frac { d Q _ { t } ^ { L } } { d t } = - e \rho _ { 0 } \int \frac { d \omega } { 2 \pi } \omega \int \frac { d \epsilon } { 4 } \text{Try} [ \hat { r } _ { 3 } \hat { g } _ { \epsilon + \omega / 2, \epsilon - \omega / 2 } ^ { L K } ] e ^ { - i \omega t } \\ = - e \rho _ { 0 } \int \frac { d \omega } { 2 \pi } \iint \frac { d \epsilon d \epsilon _ { 1 } } { 4 } \text{Try} [ \hat { r } _ { 3 } \hat { \Sigma } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } ^ { L + } \hat { g } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } ^ { L K } + \hat { r } _ { 3 } \hat { \Sigma } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } ^ { L K } \hat { g } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } ^ { L - } \\ - \hat { g } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } ^ { L + } \hat { \Sigma } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } ^ { L K } \hat { r } _ { 3 } - \hat { g } _ { \epsilon + \omega / 2, \epsilon _ { 1 } } ^ { L K } \hat { \Sigma } _ { \epsilon _ { 1 }, \epsilon - \omega / 2 } ^ { L - } \hat { r } _ { 3 } ] e ^ { - i \omega t } = - J _ { t }. \\ \text{(on is the density of states which is assumed to be equal for two superconductors. and}$$

( ρ 0 is the density of states which is assumed to be equal for two superconductors, and the volumes of superconductors are set to be 1 and not written explicitly below.)

We take account of the tunneling effect as the self-energy correction 23) as follows.

$$\hat { \Sigma } _ { \epsilon, \epsilon ^ { \prime } } ^ { L \pm } = \pi \rho _ { 0 } \int \frac { d \epsilon _ { 1 } d \epsilon _ { 2 } } { ( 2 \pi ) ^ { 2 } } \hat { t } _ { \epsilon - \epsilon _ { 1 } } ^ { L R } \hat { g } _ { \epsilon _ { 1 }, \epsilon _ { 2 } } ^ { R \pm } \hat { t } _ { \epsilon _ { 2 } - \epsilon ^ { \prime } } ^ { R L }.$$

Here,

$$\hat { t } _ { \omega } ^ { L R } = \int d t e ^ { i \omega t } \begin{pmatrix} t ^ { \prime } e ^ { - i e d \tilde { A } _ { t } } & 0 \\ 0 & - t ^ { \prime } e ^ { i e d \tilde { A } _ { t } } \end{pmatrix}. \\ \cdot _ { i \omega t } \, \frac { 1 } { 2 } \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \, \cdot _ { i \omega t } \, \dots \end{pmatrix}$$

( ˜ A t = ∫ dω 2 π e -iωt A ω . The case of zero voltage is studied in this paper.) e ± iedA ˜ t is replaced by e ∓ iedA ˜ t in ˆ t RL ω . With use of the perturbative expansion by edA ω

$$\lambda _ { \epsilon, \epsilon } ^ { L \pm } \simeq \pi \rho _ { 0 } t ^ { 2 } \left [ \hat { \tau } _ { 3 } \hat { g } _ { \epsilon, \epsilon } ^ { R \pm } \hat { \tau } _ { 3 } + \int \frac { d \epsilon _ { 1 } } { 2 \pi } i e d \left ( \hat { \tau } _ { 3 } g _ { \epsilon, \epsilon _ { 1 } } ^ { R \pm } A _ { \epsilon _ { 1 } - \epsilon ^ { \prime } } - A _ { \epsilon _ { 1 } } \hat { g } _ { 2, \epsilon ^ { \prime } } ^ { R \pm } \hat { \tau } _ { 3 } \right ) - \int \frac { d \epsilon _ { 1 } d \epsilon _ { 1 } ^ { \prime } } { ( 2 \pi ) ^ { 2 } } ( i e d ) ^ { 2 } A _ { \epsilon _ { - \epsilon _ { 1 } } } \hat { g } _ { 2, \epsilon _ { 1 } } ^ { R \pm } A _ { \epsilon _ { 1 } - \epsilon ^ { \prime } } + \cdots \right ].$$

We use the monochromatic field such as

$$A _ { \omega } = 2 \pi \bar { A } [ \delta ( \omega - \omega _ { 0 } ) + \delta ( \omega + \omega _ { 0 } ) ],$$

/negationslash and the oscillating terms ( /epsilon1 -/epsilon1 ′ = nω 0 with n = 0) will be neglected below:

$$\hat { g } _ { \epsilon, \epsilon ^ { \prime } } ^ { L + } \simeq \hat { g } _ { \epsilon } ^ { L + } 2 \pi \delta ( \epsilon - \epsilon ^ { \prime } ).$$

This approximation holds in the case that the strength of the external field is much smaller than ω 0 : 24) πρ 0 ( t edA ′ ¯ ) 2 /lessmuch ω 0 . From Eqs. (2) and (5) up to the linear response the current J t is written as follows.

$$J _ { t } = \frac { 1 } { 4 e R _ { N } } \int \frac { d \omega } { 2 \pi } \int \frac { d \epsilon } { 4 } \{ \text{Tr} [ \hat { g } _ { \epsilon } ^ { R + } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon - \omega / 2 } ^ { L K } + \hat { g } _ { \epsilon } ^ { R K } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon - \omega / 2 } ^ { L - } - \hat { g } _ { \epsilon + \omega / 2 } ^ { L + } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon } ^ { R K } - \hat { g } _ { \epsilon + \omega / 2 } ^ { L K } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon } ^ { R - } ] 2 \pi \delta ( \omega )$$

$$& J _ { t } = \frac { 1 } { 4 e R _ { N } } \int \frac { d \omega } { 2 \pi } \int \frac { d \epsilon } { 4 } \{ \text{Tr} [ \hat { g } ^ { R + } _ { \epsilon } \hat { r } _ { 3 } \hat { g } ^ { L K } _ { \epsilon - \omega / 2 } + \hat { g } ^ { R K } _ { \epsilon } \hat { r } _ { 3 } \hat { g } ^ { L - } _ { \epsilon - \omega / 2 } - \hat { g } ^ { L + } _ { \epsilon + \omega / 2 } \hat { r } _ { 3 } \hat { g } ^ { R K } _ { \epsilon } - \hat { g } ^ { L K } _ { \epsilon + \omega / 2 } \hat { r } _ { 3 } \hat { g } ^ { R - } _ { \epsilon } ] 2 \pi \delta ( \omega ) \\ & - i \text{Tr} [ \hat { r } _ { 3 } \hat { g } ^ { R + } _ { \epsilon - \omega / 2 } \hat { r } _ { 3 } \hat { g } ^ { L K } _ { \epsilon - \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { R K } _ { \epsilon - \omega / 2 } \hat { r } _ { 3 } \hat { g } ^ { L - } _ { \epsilon + \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { R K } _ { \epsilon + \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { L K } _ { \epsilon + \omega / 2 } \hat { r } _ { 3 } \hat { g } ^ { R - } _ { \epsilon + \omega / 2 }$$

$$- \hat { g } _ { \epsilon + \omega / 2 } ^ { L + } \hat { g } _ { \epsilon - \omega / 2 } ^ { R K } - \hat { g } _ { \epsilon + \omega / 2 } ^ { L K } \hat { g } _ { \epsilon - \omega / 2 } ^ { R - } - \hat { g } _ { \epsilon + \omega / 2 } ^ { R + } \hat { g } _ { \epsilon - \omega / 2 } ^ { L K } - \hat { g } _ { \epsilon + \omega / 2 } ^ { R K } \hat { g } _ { \epsilon - \omega / 2 } ^ { L - } ] e d A _ { \omega } \} e ^ { - i \omega t }$$

with R N = 1 4 / πe ρ t 2 2 0 ′ 2 . The first line of Eq. (8) reduces to the supercurrent I 0 sin( φ L -φ R ).

## 2.2 The nonequilibrium steady state

In this subsection we derive the one-particle Green functions which are used to calculate the current [Eq. (8)]. These functions are modified by the tunneling effect and the external field. The latter effect makes the system in the nonequilibrium state, and the kinetic equations for quasiclassical Green's functions are written as follows. 25)

$$\hat { \tau } _ { 3 } ( \epsilon \hat { \tau$$

and

$$\hat { \tau } _ { 3 } ( \hat { \tau } _ { 0$$

$$& \quad \, \quad \, \quad \, \quad \, \quad \, \$$

with ˆ τ 0 = ( 1 0 0 1 ), ∆ ˆ L = ( 0 ∆ iφ L ∆ -iφ L 0 ) , and

$$\begin{array} {$$

[The superscript (a) in ˆ g L a ( ) means the anomalous term in Keldysh Green function.] By omitting the oscillating part as noted above in Eq. (7) and using the perturbative expansion Eq. (5) and

$$\hat { g } _ { \epsilon } ^ { L K } = T _ { \$$

Eqs. (9) and (10) are rewritten as follows.

$$\hat { \tau } _ { 3 } ( \epsilon \hat { \tau } _ { 0 } - \hat { \Delta } _ { L } - \hat { \Sigma } _ { \epsilon } ^ { L \pm } ) \hat { g } _ { \epsilon } ^ { L \pm } - \hat { g } _ { \epsilon } ^ { L \pm } ( \epsilon \hat { \tau } _ { 0 } - \hat { \Delta } _ { L } - \hat { \Sigma } _ { \epsilon } ^ { L \pm } ) \hat { \tau } _ { 3 } = 0$$

with and

τ

ˆ (

3

/epsilon1τ

with

ˆ

0

-

ˆ

∆

L

$$\hat { \Sigma } _ { \epsilon } ^ { L \pm } = \hat { \Sigma } ^ { \prime } _ { \epsilon } + \pi \rho _ { 0 } t ^ { \prime 2 } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon } ^ { R \pm } \hat { \tau } _ { 3 } + \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { s = \pm 1 } \hat { g } _ { \epsilon + s \omega _ { 0 } } ^ { R \pm },$$

-

ˆ

Σ

L

+

-

-

/epsilon1

)ˆ

g

L a

(

)

/epsilon1

-

ˆ

g

L a

(

)

/epsilon1

(

/epsilon1τ

ˆ

0

-

ˆ

∆

L

-

ˆ

Σ

L

/epsilon1

)ˆ

τ

3

+ˆ ˆ

τ z

3

L a

(

)

/epsilon1

ˆ

g

L

/epsilon1

-

ˆ

g

L

+

/epsilon1

ˆ

z

L a

(

)

/epsilon1

ˆ

τ

3

= 0 (15)

$$\hat { z } _ { \epsilon } ^ { L ( a ) } = - \hat { \Sigma } _ { \epsilon } ^ { L ( a ) } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } ( T _ { \epsilon \pm \omega _ { 0 } } ^ { h } - T _ { \epsilon } ^ { h } ) ( \hat { g } _ { \epsilon \pm \omega _ { 0 } } ^ { R + } - \hat { g } _ { \epsilon \pm \omega _ { 0 } } ^ { R - } )$$

and

$$\hat { \Sigma } _ { \epsilon } ^ { L ( a ) } = \hat { \Sigma } _ { \epsilon } ^ { L ( a ) } + \pi \rho _ { 0 } t ^ { \prime 2 } \hat { \tau } _ { 3 } \hat { g } _ { \epsilon } ^ { R ( a ) } \hat { \tau } _ { 3 } + \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } \hat { g } _ { \epsilon \pm \omega _ { 0 } } ^ { R ( a ) }.$$

Here, Σ ˆ ′ L ± /epsilon1 and Σ ˆ ′ L a ( ) /epsilon1 in the above equations represent the inelastic scattering effect which is needed to dissipate absorbed energies and maintain the nonequilibrium steady state under the external fields, and the expression of this effect will be specified below. There also exist the equations with L and R interchanged in the above equations.

/negationslash

We introduce the following phase transformation to describe the current-carrying state due to the finite phase difference ( φ := φ L -φ R = 0);

$$\hat { \varphi } = \begin{pmatrix} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ e ^ { i ( \phi _ { L } + \phi _ { R } ) / 4 } & 0 \\ 0 & e ^ { - i ( \phi _ { L } + \phi _ { R } ) / 4 } \end{pmatrix} & ( 1 8 ) \\ \intertext { r i o h t } \ h i f t \, \text{(left)} \, \text{ of } E a \, \text{ (13)} \, \text{bv} \, \hat { \lambda } \, \text{(its commo-in-a-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-o-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-k-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i- i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-is-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-j-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-l-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-s-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-i-$$

and multiply the right (left) of Eq. (13) by ˆ (its complex conjugate ϕ ˆ ). This equation ϕ ∗ is rewritten as

$$\hat { \tau } _ { 3 } ( \tilde { \epsilon } ^ { \pm } \hat { \tau } _ { 0 } - \tilde { \Delta } ^ { \pm } _ { 1, \epsilon } \hat { \tau } _ { 1 } + i \tilde { \Delta } ^ { \pm } _ { 2, \epsilon } \hat { \tau } _ { 2 } ) \hat { g } ^ { \pm } _ { \epsilon } - \hat { g } ^ { \pm } _ { \epsilon } ( \tilde { \epsilon } ^ { \pm } \hat { \tau } _ { 0 } - \tilde { \Delta } ^ { \pm } _ { 1, \epsilon } \hat { \tau } _ { 1 } + i \tilde { \Delta } ^ { \pm } _ { 2, \epsilon } \hat { \tau } _ { 2 } ) \hat { \tau } _ { 3 } = 0.$$

Here, we put

$$\hat { g } _ { \epsilon } ^ { \pm } = \begin{pmatrix} g _ { \epsilon } ^ { \pm } & f _ { \epsilon } ^ { \pm } - h _ { \epsilon } ^ { \pm } \\ f _ { \epsilon } ^ { \pm } + h _ { \epsilon } ^ { \pm } & g _ { \epsilon } ^ { \pm } \end{pmatrix} \colon = \hat { \varphi } ^ { * } \hat { g } _ { \epsilon } ^ { L \pm } \hat { \varphi },$$

$$\tilde { \epsilon } ^ { \pm } = \epsilon - \Sigma ^ { \prime \pm } _ { \ n, \epsilon } - \pi \rho _ { 0 } t ^ { \prime 2 } g _ { \epsilon } ^ { \pm } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { s = \pm 1 } g _ { \epsilon + s \omega _ { 0 } } ^ { \pm },$$

$$\tilde { \Delta } ^ { \pm } _ { 1, \epsilon } = \Delta _ { 1 } + \Sigma ^ { \prime \pm } _ { a, \epsilon } - \pi \rho _ { 0 } t ^ { \prime 2 } f ^ { \pm } _ { \epsilon } + \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { s = \pm 1 } f ^ { \pm } _ { \epsilon + s \omega _ { 0 } },$$

$$\tilde { \Delta } _ { 2, \epsilon } ^ { \pm } = \Delta _ { 2 } + \Sigma ^ { \prime \pm } _ { \, b, \epsilon } + \pi \rho _ { 0 } t ^ { \prime 2 } h _ { \epsilon } ^ { \pm } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { s = \pm 1 } h _ { \epsilon + s \omega _ { 0 } } ^ { \pm }$$

with ∆ 1 = ∆cos( φ/ 2), ∆ 2 = -i ∆sin( φ/ 2), ˆ τ 1 = ( 0 1 1 0 ), ˆ τ 2 = ( 0 -i i 0 ), and

$$\Sigma ^ { \prime \pm } _ { n, \epsilon } \hat { \tau } _ { 0 } + \Sigma ^ { \prime \pm } _ { a, \epsilon } \hat { \tau } _ { 1 } - i \Sigma ^ { \prime \pm } _ { b, \epsilon } \hat { \tau } _ { 2 } = \hat { \varphi } ^ { * } \hat { \Sigma } ^ { \prime } _ { \epsilon } \hat { \varphi }.$$

Then the solutions to Eq. (19) are written as

$$g _ { \epsilon } ^ { \pm } = \frac { - \tilde { \epsilon } ^ { \pm } } { \tilde { \zeta } _ { \epsilon } ^ { \pm } }, \quad f _ { \epsilon } ^ { \pm } = \frac { - \tilde { \Delta } _ { 1, \epsilon } ^ { \pm } } { \tilde { \zeta } _ { \epsilon } ^ { \pm } }, \quad \text{and} \quad h _ { \epsilon } ^ { \pm } = \frac { - \tilde { \Delta } _ { 2, \epsilon } ^ { \pm } } { \tilde { \zeta } _ { \epsilon } ^ { \pm } }$$

with ˜ ζ ± /epsilon1 = √ ( ˜ ∆ ± 1 ,/epsilon1 ) 2 -( ˜ ∆ ± 2 ,/epsilon1 ) 2 -(˜ /epsilon1 ± ) 2 . From ˆ ˆ ϕ g ∗ R ± /epsilon1 ˆ and ϕ ˆ Σ ϕ ∗ ˆ ′ R ± /epsilon1 ˆ, only the sign of ˆ ϕ τ 2 component is reversed, and other than that, the same results as the above equations are obtained.

Similarly Eq. (15) is rewritten as

$$& \hat { \tau } _ { 3 } ( \tilde { \epsilon } ^ { + } \hat { \tau } _ { 0 } - \tilde { \Delta } _ { 1, \epsilon } ^ { + } \hat { \tau } _ { 1 } + i \tilde { \Delta } _ { 2, \epsilon } ^ { + } \hat { \tau } _ { 2 } ) \hat { g } _ { \epsilon } ^ { ( a ) } - \hat { g } _ { \epsilon } ^ { ( a ) } ( \tilde { \epsilon } ^ { - } \hat { \tau } _ { 0 } - \tilde { \Delta } _ { 1, \epsilon } ^ { - } \hat { \tau } _ { 1 } + i \tilde { \Delta } _ { 2, \epsilon } ^ { - } \hat { \tau } _ { 1 } ) \hat { \tau } _ { 3 } \\ & + \hat { \tau } _ { 3 } ( z _ { 0, \epsilon } \hat { \tau } _ { 0 } + z _ { 1, \epsilon } \hat { \tau } _ { 1 } - i z _ { 2, \epsilon } \hat { \tau } _ { 2 } ) \hat { g } _ { \epsilon } ^ { - } - \hat { g } _ { \epsilon } ^ { + } ( z _ { 0, \epsilon } \hat { \tau } _ { 0 } + z _ { 1, \epsilon } \hat { \tau } _ { 1 } - i z _ { 2, \epsilon } \hat { \tau } _ { 2 } ) \hat { \tau } _ { 3 } = 0.$$

Here, ˆ g ( a ) /epsilon1 = ˆ ˆ ϕ g ∗ ( a L ) /epsilon1 ˆ, ϕ

$$z _ { 0, \epsilon } = - \sum _ { n, \epsilon } ^ { \prime ( a ) } - \pi \rho _ { 0 } t ^ { \prime 2 } g _ { \epsilon } ^ { ( a ) } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } g _ { \epsilon \pm \omega _ { 0 } } ^ { ( a ) } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } ( T _ { \epsilon \pm \omega _ { 0 } } ^ { h } - T _ { \epsilon } ^ { h } ) ( g _ { \epsilon \pm \omega _ { 0 } } ^ { + } - g _ { \epsilon \pm \omega _ { 0 } } ^ { - } ),$$

$$z _ { 1, \epsilon } = - \sum _ { a, \epsilon } ^ { \prime ( a ) } + \pi \rho _ { 0 } t ^ { \prime 2 } f _ { \epsilon } ^ { ( a ) } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } f _ { \epsilon \pm \omega _ { 0 } } ^ { ( a ) } - \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } ( T _ { \epsilon \pm \omega _ { 0 } } ^ { h } - T _ { \epsilon } ^ { h } ) ( f _ { \epsilon \pm \omega _ { 0 } } ^ { + } - f _ { \epsilon \pm \omega _ { 0 } } ^ { - } ),$$

and

$$z _ { 2, \epsilon } & = - \sum _ { b, \epsilon } ^ { \prime ( a ) } - \pi \rho _ { 0 } t ^ { 2 } h _ { \epsilon } ^ { ( a ) } + \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } h _ { \epsilon \pm \omega _ { 0 } } ^ { ( a ) } + \pi \rho _ { 0 } ( t ^ { \prime } e d \bar { A } ) ^ { 2 } \sum _ { \pm } ( T _ { \epsilon$$

[Σ ′ ( a ) n,/epsilon1 , Σ ′ ( a ) a,/epsilon1 , and Σ ′ ( a ) b,/epsilon1 are the components of ˆ Σ ϕ ∗ ˆ ′ L a ( ) /epsilon1 ˆ as in Eq. (24). The ˆ -component ϕ τ 3 is omitted because it is confirmed to be small by numerical calculations and its contribution to the damping rate is negligible.] By introducing the quantity

$$x _ { \epsilon } ^ { ( a ) } \coloneqq \frac { g _ { \epsilon } ^ { ( a ) } } { g _ { \epsilon } ^ { + } - g _ { \epsilon } ^ { - } } = \frac { f _ { \epsilon } ^ { ( a ) } } { f _ { \epsilon } ^ { + } - f _ { \epsilon } ^ { - } } = \frac { h _ { \epsilon } ^ { ( a ) } } { h _ { \epsilon } ^ { + } - h _ { \epsilon } ^ { - } }, \quad \quad \quad$$

in which the equality of these three quantities follows from the normalization conditions for the quasiclassical Green functions, Eq. (26) is reduced to

$$& x _ { \epsilon } ^ { ( a ) } [ ( - \Sigma _ { n, \epsilon } ^ { \prime + } + \Sigma _ { n, \epsilon } ^ { \prime - } ) ( g _ { \epsilon } ^ { + } - g _ { \epsilon } ^ { - } ) - ( \Sigma _ { a, \epsilon } ^ { \prime + } - \Sigma _ { a, \epsilon } ^ { \prime - } ) ( f _ { \epsilon } ^ { + } - f _ { \epsilon } ^ { - } ) + ( \Sigma _ { b, \epsilon }$$

$$\times [ ( g _ { \epsilon } ^ { + } - g _ { \epsilon } ^ { - } ) ( g _ { \epsilon \pm \omega _ { 0 } } ^ { + } - g _ { \epsilon \pm \omega _ { 0 } } ^ { - } ) + ( f _ { \epsilon } ^ { + } - f _ { \epsilon } ^ { - } ) ( f _ { \epsilon \pm \omega _ { 0 } } ^ { + } - f _ { \epsilon \pm \omega _ { 0 } } ^ { - } ) - ( h _ { \epsilon } ^ { + } - h$$

here,

$$\tilde { T } _ { \epsilon } ^ { h } \coloneqq T _ { \epsilon } ^ { h } + x _ { \epsilon } ^ { ( a ) }$$

and x ( a ) /epsilon1 has a meaning of the nonequilibrium correction to the distribution function. x ( a ) /epsilon1 is obtained by solving Eq. (31).

## 2.3 Derivation of the damping rate from the current

From Eq. (8), the rf current term is written as

$$J _ { \omega } = \frac { - i } { 4 e R _ { N } } \int \frac { d \epsilon } { 4 } \text{Tr} [ \hat { r } _ { 3 } \hat { g } ^ { R + } _ { \epsilon - \omega / 2 } \hat { r } \hat { r } _ { 3 } \hat { g } ^ { L K } _ { \epsilon - \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { R K } _ { \epsilon - \omega / 2 } \hat { r } \hat { r } _ { 3 } \hat { g } ^ { L - } _ { \epsilon - \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { L + } _ { \epsilon + \omega / 2 } \hat { r } \hat { r } _ { 3 } \hat { g } ^ { R K } _ { \epsilon + \omega / 2 } + \hat { r } _ { 3 } \hat { g } ^ { L K } _ { \epsilon + \omega / 2 } \hat { r } \hat { r } _ { 3 } \hat { g } ^ { R - } _ { \epsilon + \omega / 2 }$$

$$- \hat { g } _ { \epsilon + \omega / 2 } ^ { L + } \hat { g } _ { \epsilon - \omega / 2 } ^ { R K } - \hat { g } _ { \epsilon + \omega / 2 } ^ { L K } \hat { g } _ { \epsilon - \omega / 2 } ^ { R - } - \hat { g } _ { \epsilon + \omega / 2 } ^ { R + } \hat { g } _ { \epsilon - \omega / 2 } ^ { L K } - \hat { g } _ { \epsilon + \omega / 2 } ^ { R K } \hat { g } _ { \epsilon - \omega / 2 } ^ { L - } ] e d A _ { \omega }.$$

Using the quasiclassical Green functions in the previous subsection, this current is rewritten as J ω = -CdA ω (˜ ω 2 φ -iωγ ω ) with ˜ ω 2 φ and γ ω being given as follows.

$$\tilde { \omega } _ { \phi } ^ { 2 } = \frac { 1 } { R _ { N } C } \int d \epsilon \tilde { T } _ { \epsilon } ^ { h } \left [ - R e ( g _ { \epsilon } ^ { + } - g _ { \epsilon - \omega } ^ { + } ) \text{Im} g _ { \epsilon } ^ { + } + R e ( f _ { \epsilon } ^ { + } - f _ { \epsilon - \omega } ^ { + } ) \text{Im} f _ { \epsilon } ^ { + } + \text{Im} ( h _ { \epsilon } ^ { + } - h _ { \epsilon - \omega } ^ { + } ) \text{Re} h _ { \epsilon } ^ { + } \right ]$$

and

$$\gamma _ { \omega } = \frac { 1 } { R _ { N } C } \int d \epsilon \frac { ( \tilde { T } _ { \epsilon + \omega / 2 } ^ { h } - \tilde { T } _ { \epsilon - \omega / 2 } ^ { h } ) } { 2 \omega } ( \text{Im} g _ { \epsilon + \omega / 2 } ^ { + } \text{Im} g _ { \epsilon - \omega / 2 } ^ { + } + \text{Im} f _ { \epsilon + \omega / 2 } ^ { + } \text{Im} f _ { \epsilon - \omega / 2 } ^ { + } - \text{Re} h _ { \epsilon + \omega / 2 } ^ { + } \text{Re} h _ { \epsilon - \omega / 2 } ^ { + } ).$$

We introduce the dimensionless damping rate Γ ω such as

$$\gamma _ { \omega } = \frac { 1 } { R _ { N } C } \Gamma _ { \omega }.$$

Γ ω = Γ g ω +Γ f ω with

$$\Gamma _ { \omega } ^ { g } \coloneqq \int d \epsilon \frac { ( \tilde { T } _ { \epsilon + \omega / 2 } ^ { h } - \tilde { T } _ { \epsilon - \omega / 2 } ^ { h } ) } { 2 \omega } \text{Im} g _ { \epsilon + \omega / 2 } ^ { + } \text{Im} g _ { \epsilon - \omega / 2 } ^ { + }$$

and

$$\Gamma _ { \omega } ^ { f } \coloneqq \int d \epsilon \frac { ( \tilde { T } _ { \epsilon + \omega / 2 } ^ { h } - \tilde { T } _ { \epsilon - \omega / 2 } ^ { h } ) } { 2 \omega } ( \text{Im} f _ { \epsilon + \omega / 2 } ^ { + } \text{Im} f _ { \epsilon - \omega / 2 } ^ { + } - \text{Re} h _ { \epsilon + \omega / 2 } ^ { + } \text{Re} h _ { \epsilon - \omega / 2 } ^ { + } ). \quad ( 3 8 )$$

We call Γ g ω and Γ f ω the quasiparticle term and the interference term, 3) respectively. The finite values of these two terms are brought about by the thermal excitation of quasiparticles.

## 3. Results of the Numerical Calculations

## 3.1 The damping rate in the absence of correction terms

In this subsection we show results of the damping rate in which the corrections by the tunneling effect and the nonequilibrium distribution are not included. In this case quasiclassical Green functions, Eq. (25), are approximated as follows.

$$g _ { \epsilon } ^ { \pm } = \frac { - \epsilon } { \zeta _ { \epsilon } ^ { \pm } }, \ \ f _ { \epsilon } ^ { \pm } = \frac { - \Delta \cos ( \phi / 2 ) } { \zeta _ { \epsilon } ^ { \pm } }, \ \text{ and } \ h _ { \epsilon } ^ { \pm } = \frac { i \Delta \sin ( \phi / 2 ) } { \zeta _ { \epsilon } ^ { \pm } },$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-5457ffe6-7af6-48de-9df7-e548516a7755.jpg)

0

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-11deb7b1-f0cc-4339-9280-d35e6cc973c4.jpg)

0

Fig. 1. Dependence of the damping rate of Josephson plasmons on ω J and cos φ for calculations with use of the bare Green functions: Γ 0 ω = Γ g 0 ω +Γ f 0 ω . The points indicate the damping rate at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 . The lines indicate Γ 0 ω for several cos φ whose values are indicated at the left of lines. For simplicity, ω / J ∆ is written as 0 ω J in the figure. (a) T/T c = 0 9. (b) . T/T c = 0 5. .

here ζ ± /epsilon1 = √ ∆ 2 -( /epsilon1 ± i 0) . 2 When these quantities and x ( a ) /epsilon1 = 0 are used, Eq. (35) results in

$$\gamma _ { \omega } ^ { 0 } = \frac { 1 } { R _ { N } C } ( \Gamma _ { \omega } ^ { g 0 } + \Gamma _ { \omega } ^ { f 0 } )$$

with

$$\Gamma _ { \omega } ^ { g 0 } \coloneqq & \int d \epsilon \frac { ( T _ { \epsilon + \omega / 2 } ^ { h } - T _ { \epsilon - \omega / 2 } ^ { h } ) | \epsilon + \omega / 2 | | \epsilon - \omega / 2 | \theta ( | \epsilon + \omega / 2 | - \Delta ) \theta ( | \epsilon - \omega / 2 | - \Delta ) } { 2 \omega \sqrt { ( \epsilon + \omega / 2 ) ^ { 2 } - \Delta ^ { 2 } } \sqrt { ( \epsilon - \omega / 2 ) ^ { 2 } - \Delta ^ { 2 } } } \ ( 4 1 ) \\ \text{and}$$

and

$$\Gamma _ { \omega } ^ { f 0 } \coloneqq \cos \phi \int d \epsilon \frac { ( T _ { \epsilon + \omega / 2 } ^ { h } - T _ { \epsilon - \omega / 2 } ^ { h } ) \Delta ^ { 2 } \theta ( | \epsilon + \omega / 2 | - \Delta ) \theta ( | \epsilon - \omega / 2 | - \Delta ) } { \text{sgn} ( \epsilon + \omega / 2 ) \text{sgn} ( \epsilon - \omega / 2 ) 2 \omega \sqrt { ( \epsilon + \omega / 2 ) ^ { 2 } - \Delta ^ { 2 } } \sqrt { ( \epsilon - \omega / 2 ) ^ { 2 } - \Delta ^ { 2 } } }, \\ \quad \.$$

here, θ ( ) · means a step function. This result is similar to that obtained previously in Refs. 5 and 26.

We calculate Γ 0 ω = Γ g 0 ω + Γ f 0 ω numerically. The superconducting gap at absolute zero is set to be the unit of energy ∆ 0 = 1. Figure 1 shows the dependence of the damping rate Γ 0 ω on ω J and cos φ . The lines show the dependence of the damping rate on ω and each curve takes a different value of cos φ . The points show the damping rate corresponding to each value of ω J . Since the possible range of values of ω J (cos φ ) 1 2 / differs depending on the value of ω J , the dependence of the damping rate on cos φ also differs. As the value of ω J becomes smaller, the tendency for the damping rate to

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-1f01ead0-585a-4438-b8a5-afff5234cc2c.jpg)

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-8838e97f-21cd-423e-b0dd-96ed000fb648.jpg)

0

Fig. 2. Dependence of Γ g 0 ω and Γ f 0 ω on ω J and cos φ . The points indicate the values at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 . The upper line 'g' indicates Γ g 0 ω . The five lines at the bottom indicate Γ f 0 ω for several cos φ whose values are indicated at the left of the lines. For simplicity, ω / J ∆ is written as 0 ω J in the figure. (a) T/T c = 0 9. (b) . T/T c = 0 5. .

increase becomes stronger as a function of cos φ because Γ 0 ω increases with decreasing ω . This causes the damping rate to have a positive slope as a function of cos φ . When the value of ω J is large, ω J (cos φ ) 1 2 / takes a wide range of values, so the dependence of the damping rate on cos φ becomes small. There is no significant difference between T/T c = 0 9 and . T/T c = 0 5, but the positive slope as a function of cos . φ is greater at higher temperatures.

The above dependence of the damping rate on cos φ mainly comes from the interference term Γ f 0 ω as shown in Fig. 2. Figure 2 shows the decomposition of the damping rate into the quasiparticle term (Γ g 0 ω ) and the interference term (Γ f 0 ω ). As is clear from Eq. (41), the quasiparticle term has no explicit dependence on cos φ in the absence of correction terms. Then Γ g 0 ω takes values on the same curve for all values of ω J and cos φ , and the dependence of the quasiparticle term on cos φ simply follows its dependence on ω . The damping rate due to the quasiparticle term becomes smaller as the value of cos φ increases. On the other hand, since the interference term is itself proportional to cos φ [Eq. (42)], the dependence of Γ 0 ω on cos φ is determined mainly by the interference term rather than the quasiparticle term. Then the damping rate Γ 0 ω increases with increasing the value of cos φ . At higher temperatures, there are more thermally excited quasiparticles, and the dependence of the damping rate on ω becomes weaker. The factor of cos φ becomes more effective in the damping rate at high T , and this leads to the difference

between T/T c = 0 9 and 0 5 in Fig. 1. . .

## 3.2 The damping rate in the nonequilibrium steady state

In this subsection, we consider the damping rate when a nonequilibrium distribution function is introduced as in Eq. (32) and the one-particle state is modified by a phasedependent current as in Eq. (25). Since the thermalization is brought about by lowenergy excitations, we take account of the interaction between electrons and acoustic phonons as the self-energy due to inelastic scattering introduced in Eqs. (14) and (17). The concrete expressions of ˆ Σ ϕ ∗ ˆ ′ L + /epsilon1 ˆ and ϕ ˆ Σ ϕ ∗ ˆ ′ L a ( ) /epsilon1 ˆ are similar to those of Eqs. (59) ϕ and (60) in Ref. 27 except that ˆ g ± /epsilon1 has ˆ -components in this calculation as in Eq. (20). τ 2

Numerical calculations below are performed in the same way as in Ref. 27, but there is a difference because of the existence of the phase-dependent term. Using Eq. (34) with 1 /R N C = ω /π 2 J ∆tanh(∆ 2 / T ), the numerically calculated value of the resonance frequency √ ˜ ω 2 φ at ω = ω J (cos φ ) 1 2 / shows the same dependence on ω J and cos φ as ω J (cos φ ) 1 2 / though its absolute value is slightly shifted to the high-frequency side under the influence of the external field as in Ref. 28. The smallness of this shift is related to that in the superconducting gap due to the nonequilibrium correction, which is small as compared to the change in the absorption spectrum at frequencies below the gap edge as shown in Fig. 4 of Ref. 27. We neglect this small shift, and use √ ˜ ω 2 φ /similarequal ω φ = ω J (cos φ ) 1 2 / as the resonance frequency in the calculation of the damping rate below.

Figure 3 shows the damping rate of Josephson plasmons Γ ω at several values of the resonance frequency ω J (cos φ ) 1 2 / in the nonequilibrium steady state. [ p ′ is the coupling constant between electrons and acoustic phonons introduced in Eqs. (59) and (60) in Ref. 27. The value of v k s F ( v s is the velocity of sound) in these two equations is set to be v k s F = 8 0∆ . 0 in the numerical calculations as in Ref. 27, and this value will not be written explicitly hereafter.] Unlike the case without these corrections in Fig. 1, the dependence of the damping rate on ω differs depending on the value of ω J because the nonequilibrium correction varies depending on the resonance frequency. The graphs of Γ ω in the case of cos φ = 0 61 are shown in Fig. 3, but unlike Fig. 1, the dependence on . ω differs depending on the values of ω J due to the nonequilibrium correction. (To avoid complexity, only the case of cos φ = 0 61 is shown.) It can be seen that the damping rate . becomes smaller with increasing cos φ , and the sign of the slope of Γ ω 0 is negative as a function of cos φ in contrast to the case of Γ 0 ω 0 in Fig. 1. This is because, as shown below, for the larger values of cos φ the nonequilibrium correction becomes larger, which results

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-1d113e33-8da5-4055-9541-5132ed30ec86.jpg)

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-6e3db4da-074d-44b8-bc32-2ff11d77d759.jpg)

0

Fig. 3. Dependence of the damping rate of Josephson plasmons on ω J and cos φ . The points indicate Γ ω 0 at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 . The values of cos φ are '0 15, 0 30, 0 49, 0 61, and . . . . 0 88' for . ω / J ∆ = 0 05 and 0 1, '0 15, 0 26, 0 44, 0 61, and 0 88' for 0 . . . . . . . ω / J ∆ = 0 2, and '0 15, 0 27, 0 . . . 0 42, 0 61, and 0 88' for . . . ω / J ∆ = 0 3. The four lines with no points indicate the dependence of Γ 0 . ω on ω with cos φ = 0 61 for . ω / J ∆ = 0 05, 0 1, 0 2, and 0 3. 0 . . . . πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , . 0 πρ t 0 ′ 2 = 0 01∆ , . 0 and p ′ = 0 05∆ . For simplicity, . 0 ω / J ∆ is written as 0 ω J in the figure. (a) T/T c = 0 9. (b) . T/T c = 0 5. .

in a decrease in the damping rate. This result of the dependence on cos φ is different from previous theories 5, 6) on the damping rate, and agrees with experimental results. 7-9) The results for T/T c = 0 9 and 0 5 are similar except for the magnitude of the values, . . but there are differences in the damping rate as a function of cos φ , especially for small values of ω J . At T/T c = 0 9 the dependence of Γ . ω on cos φ is smaller than that of T/T c = 0 5. This is related to the result that the sign of the slope of Γ . ω 0 as a function of cos φ changes from negative to positive as the temperature rises, which will be shown later. This also agrees with the experimental results that the sign of the slope changes near T c .

The effect of the nonequilibrium distribution on the damping rate can be seen by considering the following quantities:

$$\tau _ { x } ( \epsilon, \omega ) = ( x _ { \epsilon + \omega / 2 } ^ { ( a ) } - x _ { \epsilon - \omega / 2 } ^ { ( a ) } ) / \omega.$$

The factor ( ˜ T h /epsilon1 + ω/ 2 -˜ T h /epsilon1 -ω/ 2 ) /ω in Eq. (35) is equal to τ h ( /epsilon1, ω ) + τ x ( /epsilon1, ω ) with τ h ( /epsilon1, ω ) = ( T h /epsilon1 + ω/ 2 -T h /epsilon1 -ω/ 2 ) /ω . Figure 4(a) shows the effect of the nonequilibrium correction x ( a ) /epsilon1 on the thermal distribution tanh( /epsilon1/ 2 T ). For both T/T c = 0 9 and 0 5, the effect . . of x ( a ) /epsilon1 becomes larger with increasing cos φ . Furthermore, the sign of x ( a ) /epsilon1 changes near the energy gap and becomes negative around energies at which the density of states is large.

τ

τ

τ

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-a6ab324f-94ad-4d26-b2d0-885a79d21f56.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-aa1164ae-e758-4850-baca-aa998b474827.jpg)

0

Fig. 4. (a) Dependence of the nonequilibrium distribution function τ h ( /epsilon1, ω 0 ) + τ x ( /epsilon1, ω 0 ) on /epsilon1 at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 at T/T c = 0 9 and 0 5 with cos . . φ = 0 15 and 0 88. The four thin . . solid lines represent τ h ( /epsilon1, ω 0 ) for the same values of T/T c and cos φ . (b) Dependence of the one-particle density of states on /epsilon1 in the presence of the supercurrent: -Im g + /epsilon1 . The numbers on the left of lines indicate the values of T/T c (cos φ ). ω / J ∆ 0 = 0 05, . πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , . 0 πρ t 0 ′ 2 = 0 01∆ , and . 0 p ′ = 0 05∆ . . 0

Therefore, x ( a ) /epsilon1 has the effect of lowering the effective temperature. Since the reduction in the effective temperature is greater at the resonance frequencies with smaller phase differences, the sign of the slope of the damping rate is negative as a function of cos φ as shown in Fig. 3. This dependence of the nonequilibrium correction on cos φ is due to the broadening of the density of states in the current-carrying state.

Figure 4(b) shows the effect of the phase-dependent current on the density of states. For both T/T c = 0 9 and 0 5, the peak of the density of states becomes broad with . . decreasing cos φ , which is similar to the results in Ref. 29. This is because, as in the superconductors including paramagnetic impurities, 30,31) the effect of time-reversal symmetry breaking due to the phase-dependent current on the density of states becomes more pronounced for the smaller values of cos φ . However, as shown below, this variation of the density of states itself does not make a qualitative difference in the dependence of the damping rate on cos φ as compared to the case without corrections.

Figure 5 shows the damping rate at several values of the resonance frequency ω J (cos φ ) 1 2 / in the absence of the nonequilibrium correction: Γ ( x 0) ω = Γ g x ( 0) ω + Γ f ( x 0) ω . Here, Γ g x ( 0) ω and Γ f ( x 0) ω are those of Eqs. (37) and (38) in which x ( a ) /epsilon1 is set to be 0. Unlike the results in Fig. 3, the dependence of Γ ( x 0) ω on ω is hardly affected by the value of ω J and shows almost the same curve when the value of cos φ is the same (the case of

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-c1006bc0-ac6c-44c7-8cc9-8a869aa53922.jpg)

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-52379864-7292-4429-8454-f6979891796d.jpg)

0

Fig. 5. Dependence of the damping rate of Josephson plasmons on ω J and cos φ in the absence of the nonequilibrium correction ( x ( a ) = 0): Γ ( x 0) ω = Γ g x ( 0) ω + Γ f ( x 0) ω . The points represent Γ ( x 0) ω 0 at the resonance frequency. The values of cos φ are the same as those in Fig. 3. The lines without points indicate the dependence of Γ ( x 0) ω on ω with cos φ = 0 61 for . ω / J ∆ = 0 05, 0 1, 0 2, and 0 3. 0 . . . . πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , . 0 πρ t 0 ′ 2 = 0 01∆ , and . 0 p ′ = 0 05∆ . (a) . 0 T/T c = 0 9. (b) . T/T c = 0 5. .

cos φ = 0 61 is shown for comparison). This is almost the same as in the calculation with . the bare Green functions in Fig. 1. This shows that the negative slope of the damping rate Γ ω 0 as a function of cos φ in Fig. 3 is caused not by the current-induced change in the density of states in Fig. 4(b), but by the change in the nonequilibrium distribution in Fig. 4(a). By comparing the results between T/T c = 0 9 and 0 5, it can be seen that . . there is a difference in the slope of the damping rate Γ ( x 0) ω 0 which is similar to that in the case of Fig. 1. The dependence of the damping rate on cos φ arises from the interference term. This can be seen by decomposing Γ ( x 0) ω 0 into the quasiparticle term and the interference term.

Figure 6 shows the decomposition of the damping rate Γ ( x 0) ω 0 into the quasiparticle term (Γ g x ( 0) ω ) and the interference term (Γ f ( x 0) ω ) in the case of x ( a ) /epsilon1 = 0. The results are almost the same as those in Fig. 2, and the quasiparticle term takes values on almost the same curve regardless of the values of ω J and cos φ . (The reduction of the quasiparticle term at small resonance frequencies is due to the broadening of the density of states.) The interference terms take values on different curves depending on the values of cos φ , but its dependence on ω hardly changes depending on the values of ω J . This shows that, as mentioned above, the dependence of the damping rate on cos φ does not qualitatively change in the presence of the phase-dependent current without the nonequilibrium corrections. There are studies in which the sign of the interference term

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-7294dddc-6d81-4455-9e18-2b1d687fbfd5.jpg)

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-c1583060-2a62-49fb-a8bf-b19accca86b8.jpg)

0

Fig. 6. Dependence of the quasiparticle term (Γ g x ( 0) ω ) and the interference term (Γ f ( x 0) ω ) in the absence of the nonequilibrium correction ( x ( a ) /epsilon1 = 0) on ω J and cos φ . The points indicate the values of these terms at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 . The line 'g' indicates the dependence of Γ g x ( 0) ω on ω with ω / J ∆ = 0 1 and cos 0 . φ = 0 61. The five lines at the bottom show the dependence of . Γ f ( x 0) ω on ω with ω / J ∆ = 0 1 and cos 0 . φ (its value is shown to the left of lines). πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , . 0 πρ t 0 ′ 2 = 0 01∆ , and . 0 p ′ = 0 05∆ . For simplicity, . 0 ω / J ∆ is written as 0 ω J in the figure. (a) T/T c = 0 9. . (b) T/T c = 0 5. .

becomes negative when a phenomenological parameter that smooths out the Riedel peak is introduced, 14,32) but this calculation showed that the sign of the interference term remains positive in the current-carrying state.

When the nonequilibrium correction is taken into account, the above behaviors of the quasiparticle and the interference terms change qualitatively. Figure 7 shows the decomposition of the damping rate Γ ω in Fig. 3 into the quasiparticle term and the interference term [Eqs. (37) and (38), respectively]. Compared to the results in Figs. 2 and 6, the quasiparticle term has a greater dependence on cos φ , and the damping rate is suppressed with increasing cos φ due to the effective temperature lowering by the nonequilibrium correction. In the same way this correction suppresses also the interference term as a function of cos φ as compared to the results in Figs. 2 and 6. These behaviors bring about a negative slope in the damping rate Γ ω as a function of cos φ as shown in Fig. 3.

In Fig. 3, as a function of cos φ , the negative slope at T/T c = 0 9 is less steep than . that at T/T c = 0 5. This is attributed to the large inelastic scattering . effect at high temperatures, which results in relatively small nonequilibrium correction at high T due to the strong thermalization as shown in Fig. 4(a). In Fig. 8, this effect of the inelastic

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-6c509e65-d466-4963-9b6f-e834460cc1f2.jpg)

0

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-e77988a2-ad7b-48b6-b6b0-13164bece058.jpg)

0

Fig. 7. Dependence of the quasiparticle term (Γ g ω ) and the interference term (Γ f ω ) on ω J and cos φ . The points represent the values of these terms at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 , and the numbers indicate the values of ω / J ∆ with ( ) and ( 0 g f ) meaning Γ g ω and Γ f ω , respectively. The values of cos φ for each ω J are the same as those in Fig. 3. The two solid lines show the dependence of Γ g ω (top) and Γ f ω (bottom) on ω in the case of ω / J ∆ = 0 1 and cos 0 . φ = 0 61. . πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , . 0 πρ t 0 ′ 2 = 0 01∆ , and . 0 p ′ = 0 05∆ . (a) . 0 T/T c = 0 9. (b) . T/T c = 0 5. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01042v1-figure-fff041f7-9a3b-4436-b4ec-ea650db6250c.jpg)

0

Fig. 8. Dependence of the damping rate of the Josephson plasmons on cos φ at the resonance frequency ω 0 = ω J (cos φ ) 1 / 2 . The points indicate Γ ω 0 with cos φ = 0 15, 0 30, . . 0 49, . 0 61, . and 0 88. . The numbers on the left of lines indicate the values of T/T c ( p / ′ ∆ ). 0 ω / J ∆ = 0 05, 0 . πρ 0 ( t edA ′ ¯ ) 2 = 0 0005∆ , and . 0 πρ t 0 ′ 2 = 0 01∆ , . 0

scattering is shown in a different way by changing the values of the coupling constant ( p ′ ) between electrons and acoustic phonons. Figure 8 shows the dependence of the damping rate Γ ω on cos φ for different values of T/T c and p ′ . The sign of the slope of the damping rate as a function of cos φ changes from negative to positive as temperature

increases for a fixed p ′ . This change in sign occurs similarly when increasing p ′ at a fixed temperature. This shows that the inelastic scattering induces a change in the sign of the slope at high temperatures.

## 4. Summary and Discussion

In this paper, we calculated the dependence of the damping rate on the phase difference at the resonance frequency. We considered that the electronic state is in a nonequilibrium steady state due to an external field with this resonance frequency. In this state, the nonequilibrium distribution function is dependent on the phase difference, and then the damping rate changes depending on this difference. When the phase difference is large, the time-reversal symmetry breaking causes a broadening of the oneparticle density of states as in the presence of paramagnetic impurities. As compared to the case where the phase difference is small and the density of states has a sharp peak, this broadening has effects of strengthening the thermalization and reducing the nonequilibrium correction. Therefore, when the phase difference is small, the reduction of the effective temperature is larger than that in the case of a large phase difference. This difference in the effective temperature depending on the phase difference causes the dependence of the damping rate due to thermally excited quasiparticles on the phase difference.

As a result, the damping rate becomes smaller with decreasing phase difference (increasing cos φ ). On the other hand, the interference term has originally an explicit dependence on the phase difference, which causes an increase in the damping rate as the phase difference becomes smaller. Therefore, the actual dependence of the damping rate on the phase difference is determined by the competition between the nonequilibrium correction and the phase difference term in the interference term. The calculations in this paper showed that the former can be predominant over the latter and the damping rate becomes smaller with decreasing the phase difference. This is consistent with experimental results showing that the damping rate has a negative slope as a function of cos φ . As the temperature rises, the effect of the inelastic scattering by acoustic phonons increases, and the nonequilibrium correction becomes smaller due to the strong thermalization. Therefore, at high temperatures, the interference term becomes predominant in the damping rate, and it is possible that the damping rate increases with decreasing phase difference. This means that the sign of the slope of the damping rate as a function of cos φ can change at high temperatures. This behavior qualitatively agrees with the

experimental results.

Other previous studies have assumed that the dependence of the damping rate on the phase difference arises only from the interference term, and have discussed the sign of this term. 14,15) Since these studies considered the zero frequency limit, the dependence of the quasiparticle term on the frequency was not taken into account. In experiments, the frequency is finite, and ω / J ∆ is about 1 0 / 20. 7) In our study, it was shown that even at a value of ω / J ∆ = 0 05, the dependence of the damping rate on the frequency exists 0 . due to the nonequilibrium correction in the quasiparticle term.

Considering the comparison with the experiment, in our calculation γ /ω ω J = Γ ω / R ( N Cω J ) /similarequal 0 007 and 0 038 at . . T/T c = 0 5 and . T/T c = 0 9, respectively, by using . an experimental value (1 /R N C /ω ) J /similarequal 0 016. . 7) This is almost consistent with the experimental value γ /ω ω J /similarequal 0 014 at . T/T c = 0 583. Though it is difficult to control the . effect of inelastic scattering experimentally, the magnitude of electron-phonon scattering varies depending on materials. 33) Therefore, it is possible to verify this theory by examining the temperatures at which the dependence of the damping rate on the phase difference changes for various materials.

## Acknowledgment

The numerical computation in this work was carried out at the Yukawa Institute Computer Facility.

## References

- 1) B. D. Josephson, Phys. Lett. 1 , 251 (1962).
- 2) B. D. Josephson, Adv. Phys. 14 , 419 (1965).
- 3) D. N. Langenberg, Rev. Phys. Appl. 9 , 35 (1974).
- 4) A. Barone and G. Patern`, o Physics and Applications of the Josephson Effect (Wiley, New York, 1982). Chap. 2, Sect. 6.
- 5) R. E. Harris, Phys. Rev. B 10 , 84 (1974).
- 6) U. K. Poulsen, Rev. Phys. Appl. 9 , 41 (1974).
- 7) N. F. Pedersen, T. F. Finnegan, and D. N. Langenberg, Phys. Rev. B 6 , 4151 (1972).
- 8) O. H. Soerensen, J. Mygind, and N. F. Pedersen, Phys. Rev. Lett. 39 , 1018 (1977).
- 9) S. Rudner and T. Claeson, J. Appl. Phys. 50 , 7070 (1979).
- 10) G. Catelani, J. Koch, L. Frunzio, R. J. Schoelkopf, M. H. Devoret, and L. I. Glazman Phys. Rev. Lett. 106 , 077002 (2011).
- 11) J. Lepp¨kangas, M. Marthaler, and G. Sch¨n, Phys. Rev. B a o 84 , 060505(R) (2011).
- 12) I. M. Pop, K. Geerlings, G. Catelani, R. J. Schoelkopf, L. I. Glazman, and M. H. Devoret, Nature 508 , 369 (2014).
- 13) K. Hida and Y. Ono, J. Low Temp. Phys. 29 , 355 (1977).
- 14) A. B. Zorin, I. O. Kulik, K. K. Likharev, and J. R. Schrieffer, Sov. J. Low Temp. Phys. 5 , 537 (1979).
- 15) A. M. Gulyan and G. F. Zharkov, Sov. Phys. JETP, 62 , 89 (1985).
- 16) A. Levy Yeyati, A. Mart´ ın-Rodero, and J. C. Cuevas, J. Phys. Condens. Matter 8 , 449 (1996).
- 17) G. M. Eliashberg, JETP Lett. ´ 11 , 114 (1970).
- 18) A. J. Dahm, A. Denenstein, T. F. Finnegan, D. N. Langenberg, and D. J. Scalapino, Phys. Rev. Lett. 20 , 859 (1968).
- 19) D. Pines and P. Nozi`res, e The Theory of Quantum Liquids (W.A. Benjamin, New York, 1966). Chap. 4, Sect. 7.
- 20) P. W. Anderson, in Lectures on the Many-Body Problem , ed. by E. R. Caianiello (Academic Press, New York, 1964), Vol. 2.
- 21) V. Ambegaokar and A. Baratoff, Phys. Rev. Lett. 10 , 486 (1963).

- 22) J. Voutilainen, T. T. Heikkil¨, and N. B. Kopnin, Phys. Rev. B a 72 , 054505 (2005).
- 23) A. F. Volkov, Sov. Phys. JETP, 41 , 376 (1974)
- 24) A. V. Semenov, I. A. Devyatov, P. J. de Visser, and T. M. Klapwijk, Phys. Rev. Lett. 117 , 047002 (2016).
- 25) G. M. Eliashberg, Sov. Phys. JETP ´ 34 , 668 (1972).
- 26) A. I. Larkin and Yu. N. Ovchinnikov, Sov. Phys. JETP 24 , 1035 (1967).
- 27) T. Jujo, J. Phys. Soc. Jpn. 92 , 074706 (2023).
- 28) P. J. de Visser, D. J. Goldie, P. Diener, S. Withington, J. J. A. Baselmans, and T. M. Klapwijk, Phys. Rev. Lett. 112 , 047004 (2014).
- 29) E. V. Bezuglyi, A. S. Vasenko, V. S. Shumeiko, and G. Wendin, Phys. Rev. B 72 , 014501 (2005).
- 30) A. A. Abrikosov and L. P. Gor'kov, Sov. Phys. JETP 12 , 1243 (1961).
- 31) S. Skalski, O. Betbeder-Matibet, and P. R. Weiss, Phys. Rev. 136 , A1500 (1964).
- 32) F. Sheldon, S. Peotta, and M. Di Ventra, Eur. Phys. J. Appl. Phys. 81 , 10601 (2018).
- 33) S.B.Kaplan, C.C.Chi, D.N.Langenberg, J.J.Chang, S.Jafarey, and D.J.Scalapino, Phys. Rev. B 14 , 4854 (1976).
Takanobu Jujo
Phase-dependent damping rate of Josephson plasmons in the nonequilibrium steady state
We calculate the damping rate of Josephson plasmons due to thermally excited
quasiparticles in a nonequilibrium steady state. It is known that the
experimentally observed dependence of the damping rate on the phase difference
behaves in the opposite way to that of theory. It is shown that this problem
can be solved by calculating the damping rate by incorporating changes in the
electronic state at the resonant frequency. When the phase difference is large,
the broadening of the one-particle state is large due to the presence of a dc
current and brings about the effective thermalization, so that the
nonequilibrium correction becomes small. Conversely, when the phase difference
is small, the nonequilibrium correction becomes large and causes a smaller
damping rate due to a reduction in the effective temperature. This dependence
of the nonequilibrium correction on the phase difference results in a damping
rate that is consistent with experimental results.