# Coupling of intrinsic Josephson oscillations in layered superconductors by charge fluctuations

## 1 INTRODUCTION

The superconducting properties of the highly anisotropic cuprate–superconductors Tl<sub>2</sub>Ba<sub>2</sub>Ca<sub>2</sub>Cu<sub>3</sub>O<sub>10+δ</sub> (TBCCO) and Bi<sub>2</sub>Sr<sub>2</sub>CaCu<sub>2</sub>O<sub>8+δ</sub> (BSCCO) are well described by a stack of Josephson junctions coupling the superconducting CuO<sub>2</sub> layers in $`c`$-direction. In particular the multiple branch structure observed in the current-voltage characteristics by several groups  can be explained by this model.

Due to the low value of the critical current in $`c`$-direction the system has a small Josephson plasma frequency $`\omega _p`$ and a large value of the McCumber parameter $`\beta _c`$, which causes a strong hysteretic behaviour of the current-voltage characteristics. The low value of the Josephson plasma frequency manifests itself in the transparency of the stack with respect to THz radiation in $`c`$-direction . The longitudinal and transversal plasma oscillations have also been observed directly .

The theoretical investigation of longitudinal plasma oscillations has become popular since Koyama and Tachiki  proposed a coupling of Josephson oscillations in different barriers due to charge fluctuations. In systems with weakly coupled superconducting layers the charges on different layers need not to be constant, which is in contrast to ordinary superconductors where charge neutrality can be assumed. In the theory of Koyama and Tachiki  it follows that the gauge-invariant scalar potential

$$\mu _l=\mathrm{\Phi }_l(\mathrm{}/2e)\dot{\chi }_l$$
(1)   
does not vanish. Here $`\mathrm{\Phi }_l`$ is the electric scalar potential and $`\chi _l`$ is the phase of the superconducting order parameter $`\mathrm{\Delta }_l=|\mathrm{\Delta }_l|\mathrm{exp}(i\chi _l)`$ on layer $`l`$. On the other hand the Josephson current density $`j_c\mathrm{sin}\gamma _{l,l+1}`$ between layers $`l`$ and $`l+1`$ depends on the gauge-invariant phase difference

$$\gamma _{l,l+1}(t)=\chi _l(t)\chi _{l+1}(t)\frac{2e}{\mathrm{}}_l^{l+1}𝑑zA_z(z,t)$$
(2)   
where $`A`$ is the vector potential (in our notation $`e=|e|`$). Its time derivative (the second Josephson relation)

$$\frac{\mathrm{}}{2e}\dot{\gamma }_{l,l+1}(t)=_l^{l+1}𝑑zE_z+\mu _{l+1}\mu _l$$
(3)   
then not only depends on the voltage between the layers, but also on the potential difference $`\mu _{l+1}\mu _l`$. This finally leads to a coupling between Josephson oscillations in different barriers. A non-vanishing generalized scalar potential which is related to quasi-particle charge imbalance can also be obtained without breaking charge neutrality. The importance of this effect has been stressed by Artemenko and Kobelkov  and by Ryndyk . In the layered cuprate-superconductors probably both effects are present.

In this communication we want to derive the coupling effect in a microscopic model starting from a tunneling Hamiltonian. We arrive at an expression for the current between different layers which is formally similar to that obtained by Artemenko for a different model. In leading order in the interlayer hopping $`t_{}`$ our results are also similar to the model of Koyama and Tachiki. This simplified model will be used to study the non-linear dynamics of the system and to discuss implications of the coupling on the current-voltage characteristics.

## 2 General outline of the theory

We start from a model where the current between superconducting layers across the insulating barrier is described by a time-dependent tunneling Hamiltonian

$$H_T=\underset{l,k,k^{},\sigma }{}T_{k,k^{}}c_{l+1,k^{},\sigma }^{}c_{l,k,\sigma }^{}e^{\frac{ie}{\mathrm{}}_l^{l+1}𝑑zA_z(t)}+h.c.$$
(4)   
which depends on the vector potential $`A_z(t)`$ in the barrier. Here $`T_{k,k^{}}`$ is a tunneling matrix element describing (random) hopping between neighboring layers. In order to get a Josephson current also for d-wave superconductors we have to keep some angular dependence in this matrix element.

In addition to this the current is driven by the difference of scalar potentials $`\mathrm{\Phi }_l(t)`$ on the different layers. Thus the total time-dependent Hamiltonian is

$$H=\underset{l}{}(H_le\mathrm{\Phi }_l(t)N_l)+H_T(t)$$
(5)   
where $`H_l`$ is the Hamiltonian of the electrons in layer $`l`$ including superconducting interactions and pairing.

Finally we have to take into account phase fluctuations of the superconducting order parameter induced by charge fluctuations. It can be shown that the results of the calculation for physical quantities depend only on the gauge-invariant combinations $`\gamma _{l,l+1}`$ and $`\mu _l`$ of the electromagnetic potentials with the phase of the order parameter. Formally the same results are obtained if we replace in the Hamiltonian the exponential by $`\mathrm{exp}(i\gamma _{l,l+1}(t))`$, the scalar potential $`\mathrm{\Phi }_l(t)`$ by $`\mu _l(t)`$ and assume the order parameter to be real in the BCS-treatment of $`H_l`$.

The non-linear Josephson effect results from the periodic dependence on the phases $`\gamma (t)`$. In the limit of large McCumber parameter $`\beta _c`$ the phase can be written as

$$\gamma (t)=\gamma _0+\omega t+\delta \gamma (t)$$
(6)   
where $`\gamma _0`$ is the constant phase determined by the dc-current and $`\omega `$ is the Josephson frequency which is related to the dc-voltage. The part $`\delta \gamma (t)`$, which oscillates with the same frequency $`\omega `$, is small for large $`\beta _c`$.

We calculate the current response $`j_{l,l+1}`$ between neighboring superconducting layers with respect to both the tunneling Hamiltonian $`H_T`$ and the generalized scalar potential $`\mu _l`$ restricting ourselves to second order tunneling processes and linear response with respect to $`\mu _l`$.

The result can be written in the following general form:

$$j_{l,l+1}(t)=j_{l,l+1}^{qp}(t)+j_{l,l+1}^J(t)$$
(7)   
where

$$j_{l,l+1}^{qp}(t)=_{\mathrm{}}^t𝑑t_1S_0(tt_1)\frac{\dot{\gamma }(t_1)}{2}\mathrm{cos}\frac{\gamma (t)\gamma (t_1)}{2}+_{\mathrm{}}^t𝑑t_1_{\mathrm{}}^t𝑑t_2S_1(t,t_1,t_2)\mathrm{cos}\frac{\gamma (t)\gamma (t_1)}{2}(\mu _{l+1}(t_2)\mu _l(t_2))$$
(8)   
is the quasi-particle current density and

$`j_{l,l+1}^J(t)`$ $`=`$ $`j_c\mathrm{sin}\gamma (t)+{\displaystyle _{\mathrm{}}^t}𝑑t_1J_0(tt_1){\displaystyle \frac{\dot{\gamma }(t_1)}{2}}\mathrm{cos}{\displaystyle \frac{\gamma (t)+\gamma (t_1)}{2}}`$ (9)   
$`+`$ $`{\displaystyle _{\mathrm{}}^t}𝑑t_1{\displaystyle _{\mathrm{}}^t}𝑑t_2J_1(t,t_1,t_2)\mathrm{cos}{\displaystyle \frac{\gamma (t)+\gamma (t_1)}{2}}(\mu _{l+1}(t_2)\mu _l(t_2))`$   
is the Josephson current carried by the condensate ($`\gamma (t):=\gamma _{l,l+1}(t)`$). The functions $`S_0(tt_1)`$ and $`J_0(tt_1)`$ result from a folding of two normal and anomalous Green’s functions in neighboring layers (see Fig. 1), and $`j_c=J_0(t,t)=const`$ is the critical current density. In the terms $`S_1(t,t_1,t_2)`$ and $`J_1(t,t_1,t_2)`$ the additional linear dependence of the Green’s functions on $`\mu _l(t_2)`$ is considered. A similar (but more complicated expression) is obtained for the charge density response (see Fig. 2). Equations with a similar structure have been obtained by Artemenko and Kobelkov  for a different model.

For a systematic evaluation one has to insert the ansatz for $`\gamma _{l,l+1}(t)`$ into these expressions and has to separate different harmonics in the Josephson frequency $`\omega `$. A considerable simplification is obtained if one keeps the non-linear effects only in the sin-term of the Josephson current and linearizes the other terms with respect to $`\gamma `$. Then one obtains for the current density across the barrier with thickness b

$$j_{l,l+1}=j_c\mathrm{sin}\gamma _{l,l+1}+\sigma _0\frac{\mathrm{}}{2eb}\dot{\gamma }_{l,l+1}+\sigma _1(\mu _l\mu _{l+1})/b$$
(10)   
and for the density response

$$\delta \rho _l=\chi _{\rho \rho }^{(0)}\mu _l+\chi _{\rho \rho }^{(2)}(\mu _{l+1}+\mu _{l1}2\mu _l)+\frac{\mathrm{}}{2eb}\sigma _1(\gamma _{l1,l}\gamma _{l,l+1}).$$
(11)   
In the normal state $`\sigma _0=\sigma _1`$ and $`j_{l,l+1}`$ depends only on the voltage across the barrier. In the superconducting state $`\gamma _{l,l+1}`$ and $`\mu _l`$ have separate physical meaning and we need the second equation for $`\delta \rho _l`$ to determine $`\mu _l`$ as a function of the voltage. The density reponse function $`\chi _{\rho \rho }^{(0)}`$ is finite in the superconducting state. At low temperatures it is only weakly frequency dependent. Approximately it is given by $`\chi _{\rho \rho }^{(0)}2e^2N_2(0)`$, where $`N_2(0)`$ is the two-dimensional density of states of the electron gas in the CuO<sub>2</sub>-layers. The conductivities $`\sigma _{0,1}`$ as well as $`\chi _{\rho \rho }^{(2)}`$ describe the charge exchange with the neighboring layers and are proportional to $`t_{}^2`$. Adding finally the displacement current we obtain a relation between the electronic current density across the barrier with the external current density,

$$j=j_{ł,l+1}+ϵϵ_0\dot{E}_{l,l+1}.$$
(12)   
In a first step we eliminate the scalar potential difference $`\mu _l`$ in favour of the gauge-invariant phase difference $`\gamma _{l,l+1}`$ and the electric field by using the Josephson relation

$$\frac{\mathrm{}}{2e}\dot{\gamma }_{l,l+1}(t)=bE_{l,l+1}+\mu _{l+1}\mu _l.$$
(13)   
In the next step we express the charge fluctuations with help of the Maxwell equation

$$\delta \rho _l=ϵ_0ϵ(E_{l,l+1}E_{l1,l}).$$
(14)   
We then finally arrive at the following differential equation for the phase:

$$\frac{j}{j_c}=\left(1\alpha \mathrm{\Delta }^{(2)}\right)\mathrm{sin}\gamma _{l,l+1}+\frac{1}{\omega _c}\left(1\eta \mathrm{\Delta }^{(2)}\right)\dot{\gamma }_{l,l+1}+\frac{1}{\omega _p^2}\left(1\zeta \mathrm{\Delta }^{(2)}\right)\ddot{\gamma }_{l,l+1}$$
(15)   
where $`\omega _p^2=j_c2eb/(ϵϵ_0\mathrm{})`$, $`1/\omega _c=\sigma _0/(ϵϵ_0\omega _p^2)`$ The dimensionless quantities $`\alpha ,\eta ,\zeta `$ describe the coupling of the phase-difference in different layers via the derivative operator $`\mathrm{\Delta }^{(2)}`$, which is defined as $`\mathrm{\Delta }^{(2)}f_l=f_{l+1}+f_{l1}2f_l`$. In particular

$$\alpha =ϵϵ_0/(b\chi _{\rho \rho }^0)+O(t_{}^2),$$
(16)   
$$\eta =ϵϵ_0/(b\chi _{\rho \rho }^0)(12\sigma _1/\sigma _0)+O(t_{}^2).$$
(17)   
$`\zeta `$ is proportional to $`t_{}^2`$. For $`\omega \mathrm{\Delta }`$ the quantity $`\alpha `$ is only weakly frequency dependent (see Fig. 3). Therefore it will be approximated by its value at $`\omega =0`$, $`\alpha (0)=ϵϵ_0/(2e^2bN_2(0))`$. If we neglect $`\eta `$ and $`\zeta `$ we arrive back at the theory of Koyama and Tachiki . In fact at $`\omega T\mathrm{\Delta }`$ we find $`\sigma _1\sigma _0/2`$ for d-wave superconductors. Thus neglecting $`\eta `$ and $`\zeta `$ seems to be a good approximation for small values of $`\omega `$. For strong coupling in $`c`$-direction one ends up with a situation, where charge fluctuations on the layers are suppressed, $`\delta \rho _l`$ = 0, but a finite $`\mu _l`$ is generated by charge-imbalance of quasi-particles . A study of these effects in this model will be done in the future.

## 3 Influence on the current-voltage characteristics

In our further discussion of the experimental implications of this effect on the current-voltage characteristics we restrict ourselves to the approximation (15) with $`\alpha =const,\eta =\zeta =0`$ supplemented by

$$\frac{\mathrm{}}{2e}\dot{\gamma }_{l,l+1}=V_{l,l+1}\alpha (V_{l+1,l+2}+V_{l1,l}2V_{l,l+1})$$
(18)   
where we have defined the voltage $`V_{l,l+1}=_l^{l+1}E_z𝑑z=bE_{l,l+1}`$. In particular, we obtain a relation for the dc-voltages, if we replace $`\dot{\gamma }_{ł,l+1}`$ by its time-average $`<\dot{\gamma }_{l,l+1}>`$. A junction is called to be in the resistive state if $`<\dot{\gamma }_{l,l+1}(t)>0`$. In the case of one junction in the resistive state there is a finite voltage-drop also in the neighboring junctions due to the coupling $`\alpha `$. This is shown in Fig. 4. Note that the total voltage is given by

$$V=\underset{l}{}V_{l,l+1}=\frac{\mathrm{}}{2e}\underset{l}{}\dot{\gamma }_{l,l+1}$$
(19)   
as in the absence of the coupling.

We want to study in particular the influence of the dynamical coupling between different barriers in the resistive state. This will be done both by numerical simulations and analytical calculations using a Green’s function technique similar to that used by Takeno in a different context. Details of the calculations will be published elsewhere.

Let us begin with the discussion of one barrier in the resistive state. It may be helpful to visualize the dynamics of such a system by considering the mechanical analog of the RSJ model: The dynamics of each phase difference can be described by a pendulum with a constant torque being proportional to the bias current $`I`$ which is the same in all barriers (Fig. 5). In the absence of the coupling $`\alpha `$ the pendulum is either rotating, this corresponds to the resistive state of the barrier, or has a constant phase, this corresponds to the superconducting state. In the presence of the coupling which produces an additional torque, one still can distinguish between these two types of motion: a rotating state and a non-rotating vibrating state. In the rotating state there is a running phase, $`\gamma (t)=\gamma _0+\omega t+\delta \gamma (t)`$ with a finite value of the Josephson frequency $`\omega =<\dot{\gamma }>`$. In the non-rotating state $`<\dot{\gamma }>=0`$, but there are still oscillations. Such localized solutions are known in non-linear dynamics as roto-breathers . The running phase in the resistive barrier causes finite phase oscillations in the neighboring barriers. The amplitude of these oscillation depends on the ratio of the rotation frequency, i.e. the Josephson frequency $`\omega `$, and the eigenfrequencies of the oscillations of the coupled system at small amplitudes. The latter are determined by the Josephson plasma frequency $`\omega _p`$ and have a bandwidth which is proportional to the coupling $`\alpha `$. In Fig. 6 we show the result of an analytical calculation of the oscillation amplitude for different barriers as a function of the distance from the resistive barrier and as a function of the Josephson frequency $`\omega `$. In this example the plasma frequency is $`\omega _p/(2\pi )=0.6`$ THz and $`\alpha =0.2`$. If the rotation frequency is high, $`\omega \omega _p`$, which is usually the case for intrinsic Josephson systems, the oscillation amplitude in the neighboring barriers is small and falls off exponentially with the distance from the resistive barrier. In the case of a Josephson frequency $`\omega `$ within the plasma band, the phase-rotation in one barrier leads to long-range plasma oscillations in the neighborhood. These are the longitudinal plasma-oscillations considered by several authors which in principle can also be excited by longitudinal electric fields in the purely superconducting state.

From an experimental point of view it is important to know how this effect influences the current-voltage characteristics of an intrinsic Josephson system, showing a multi-branch structure, where the n-th branch corresponds to n resitive barriers. As we have seen above, the total voltage measured for the stack of Josephson junctions is still given by the sum of the Josephson frequencies of the resistive barriers. For the total dc-voltage only the values of $`<\dot{\gamma }_{l,l+1}>`$ for the resistive barriers contribute. In the case of two or more resistive barriers in the stack it makes a difference whether these resistive barriers are next neighbors or are separated by one or more non-resistive junctions. We have checked this both by numerical simulations and analytical calculations. We find that in the case of two resistive barriers next to each other the rotation frequency of both junctions is the same, but (at the same current) is slightly higher than in the case of well-separated junctions. For two uncoupled junctions, of course, the total voltage is just the double of one resistive junction (the first branch). This is shown in Fig. 7 where we compare the second branch of the current voltage characteristics for two neighboring resistive junctions with the total voltage of two well separated resistive junctions. For realistic values of $`\alpha <1`$ the effect is very small, therefore we have exaggerated in the figure this effect. More pronounced is the difference near the plasma frequency: If we lower the voltage towards the upper edge of the plasma band, at least one of the resistive junctions returns to the superconducting state. This return point (the minimum of the curves) is different for the two situations: for coupled junctions the return-voltage and the dc-current is smaller than for uncoupled resistive junctions.

Our numerical simulations also show that in the case of two neighboring resistive junctions the two phases rotate coherently. This remains true even if the critical currents of the two junctions are slightly different. This phase-locking is very important for application of coupled Josephson oscillations for the generation or amplification of radiation.

Our experiments on the intrinsic Josephson systems BSCCO and TBCCO show a very precise additive structure of the different branches at least for the lower branches. Deviations from the additivity at higher order branches may be attributed to heating effects. This poses an upper limit to the coupling parameter of $`\alpha <1`$. An upper limit of $`\alpha `$ is also obtained from the hysteretic return point of the current voltage characteristics from the first resistive branch to the superconducting state. A lower limit of this return point is given by the upper edge of the plasma band given by $`\omega _p\sqrt{1+2\alpha }`$. From the experimental value for $`\omega _p`$ and the return voltage $`\omega _{return}`$ one can also estimate $`\alpha <1`$. Finally one may calculate $`\alpha `$ for a 2-dimensional electron gas: for $`ϵ25`$ one obtains $`\alpha =0.1`$.

In this paper we have studied the coupling of Josephson oscillations in different barriers due to charge fluctuations. Another mechanism which also leads to such a coupling is the excitation of $`c`$-axis phonons in the stack by Josephson oscillations, which we have studied recently . The phonons excitated in one resistive barrier induce field oscillations also in neighboring barriers and thus support coherence of Josephson oscillations.

## 4 Summary

Starting from a microscopic model for the Josephson effect in a stack of superconducting layers coupled by tunneling barriers we have derived a set of equations for the gauge-invariant phase difference in different barriers. A simplified model has then been used to study the non-linear dynamics of the coupled system. In particular, we have studied the difference in the current-voltage characteristics of two neighboring resistive junctions in comparison with well separated resistive junctions. From the regular additive structure of the current-voltage branches found in experiments on BSCCO and TBCCO and the return voltage one can derive an upper limit of the coupling constant of the order of $`\alpha <1`$.

## ACKNOWLEDGMENTS

This work was supported by the Bayerische Forschungsstiftung (Ch.P.), the German Science Foundation (Ch.H.), and the Humboldt-Foundation (A.S.).