**2. Two-junction quantum interference devices**

222 Superconductors – Materials, Properties and Applications

superconducting interference loop, τ=2π*RIJt*/Φ0=*t*/

field-modulated maximum current *IJF* (*IJF*=|cosπ

parameter,

ψ

identity:

β

where *Xex*= cosπ

to the usual sin

the quantities

of motion for

variables

β

ψ

φ

show the additional sin2

φ

τφ

ψ

When the information for

Eq. (3) is finally obtained.

ψ

( ) <sup>d</sup> <sup>B</sup> 1 cos sin d 2

 *i* πψ

where *n* is an integer denoting the number of fluxons initially trapped in the

current normalized to *IJ*. In what follows we shall consider zero-field cooling conditions, thus taking *n*=0. Eq. (1) is similar to the non-linear first-order ODE describing the dynamics of the gauge-invariant superconducting phase difference across a single overdamped JJ with

flows. This strict equivalence comes from the hypothesis that the total normalized flux

 β

we may say that the above hypothesis may be stated merely by means of the following

the device behaves as if the equivalent Josephson junction possessed a non-conventional

<sup>d</sup> <sup>2</sup> <sup>B</sup> sin sin 2 d 2 *ex ex <sup>i</sup> X Y*

coupling between the externally applied flux and the system, as described by Eq. (2), when

≠0. Therefore, the non-conventional CPR of the equivalent JJ in the SQUID model cannot be considered as a strict consequence of an intrinsic non-conventional CPR of the single JJs. The Josephson junctions in the device, in fact, could behave in the most classical way, obeying strictly to the Josephson current-phase relation; the interferometer, however, would still

reduction in the dimensional order of the dynamical equations is possible, it is noted that

very fast with respect to the equivalent junction dynamics given in Eq. (3). As a

The single-junction model can be adopted also when dealing with more complex systems, as one-dimensional Josephson junction arrays. In fact, by assuming small inductance

φ

ψ

*ex*. A second-order harmonic in

φ πβ

term for finite values of

ψ

( ) 1 2 , *ex*

*ex*

τφ

*ex*=Φ*ex*/Φ0 is the externally applied flux normalized to Φ0 and *iB= IB*/*IJ*, is the bias

 φ

+ − = (1)

, *R* being the intrinsic resistive junction

*ex*|) in which a normalized bias current *iB*/2

*ex*. However, being

, one can see that the

thus appears in addition

. In this way, for

can be considered

φ,

. In order to understand how the

in the superconducting SQUID loop,

can be assumed to be adiabatic and the equation

β

can be solved by perturbation analysis.

, Eq. (1) is not anymore valid and

ψ

β

=+ − *i i* (2)

β

 φ

β

β, since τψ/τφ *=* 2πβ

is substituted back into the effective dynamical equation for

φ

++ = (3)

addendum, however, arises solely from electromagnetic

*= L*/*R*, denoting the characteristic time scales of the

φ

*n*

φ

τ

=Φ/Φ0 linked to the interferometer loop can be taken to be equal to

*=0*. Therefore, for finite values of the parameter

current-phase relation (CPR). In fact, for small finite values of

φ

τ

ψ

φ

τψ

constant applied magnetic fields, the flux dynamics for small values of

in terms of the quasi-static variable

following model may be adopted (Romeo & De Luca, 2004):

*ex* and *Yex*= sinπ

term. The sin2

φ

respectively, are intimately linked to the parameter

ψ

and of the number of fluxons

*=* Φ0/2π*RIJ* and

consequence, the superconducting phase

ψψ

Let us consider a symmetric two junction interferometer with equal junctions of negligible capacitance, as shown in fig. 1. The dynamical equations for the variables φ and ψ characterizing this system, can be written in the following form (Romeo & De Luca, 2004):

$$\begin{aligned} \frac{\mathbf{d}\,\phi}{\mathbf{d}\,\pi} + \left(-1\right)^{n} \cos\left(\pi\nu\right) \sin\phi &= \frac{i\_B}{2};\\ \pi \frac{\mathbf{d}\,\nu}{\mathbf{d}\,\pi} + \left(-1\right)^{n} \sin\left(\pi\nu\right) \cos\phi + \frac{\mathcal{W}}{2\mathcal{J}} &= \frac{\mathcal{V}\_{c\mathbf{c}}}{2\mathcal{J}}.\end{aligned} \tag{4}$$

where *n* is the number of initially trapped fluxons in the superconducting loop. Let us consider a new time variable 2 *R t L* τ θ πβ = = and write the solutions for φ a and ψ in the following form: 0, 1, <sup>0</sup> <sup>1</sup> ( , ) ( ) ( ); ( , ) ( ) ( ) β β φ β τ ≈ + φ τ βφ τ ψ β τ ≈ + ψ θ βψ θ. For simplicity, set *n*=0.

Effective Models of Superconducting Quantum Interference Devices 225

 πψ θ

β

<sup>=</sup> , so that, for negligible values of this ratio, the system behaves

ψ

, one would follow the more general perturbation analysis

. Therefore, while in the present section we shall only be

<< ), in Section 5 we shall consider both evolution time, by

φ τ

, by assuming that the transient of the flux

ψ<sup>0</sup> and

, the dynamical equations for the

ψ<sup>1</sup> :

 πβθ

 φ β

(7)

(8)

2 *RIJ*

<sup>Φ</sup> <sup>=</sup> ,

π

. We have

φ τ

φ

. In this case, therefore, we

β

<sup>=</sup> , in (Grønbech-

β, one

to acquire a

ψ

τ

τ

φ

 φ τ

() () ( ) () () () ( ) ( )

( ) ( ) ( ) ( )

1 0 0,

 πψ θ

1, 0 0, 1 0 0,

<sup>d</sup> , (a) <sup>d</sup>

<sup>d</sup> 2sin cos 2 . (b) <sup>d</sup>

In Eqs. (7a-b) and (8a-b) we may notice the appearance of two different time scales the first,

effectively as if the adiabatic time evolution of the superconducting phase difference

may first let the flux variable evolve, so that a stationary magnetic state is reached, and then solve for the superconducting phase difference time evolution of the system. This is exactly

Jensen et al., 2003) and in (Romeo & De Luca, 2004). However, if one were to acquire the regular solution for the system dynamics, even when considering the approximate solution

described above, where the ratio *r* might not a priori be considered as negligible. Furthermore, considering that this ratio is proportional to the perturbation parameter

Despite the fact that the more general approach allows extension to higher order approximations of the perturbation solutions, we wish to limit our analysis to the study of the electrodynamic properties of a two-junction or a multi-junction quantum interferometer

> φ τ

= , linked to flux motion in and out the superconducting ring, the second, <sup>0</sup>

cos cos sin sin ; (b) <sup>d</sup>

cos sin , (a) d 2

*B*

 β

 φ τ πψ θ

*ex*

pertaining to the dynamics of the superconducting phase difference value

could be studied by taking asymptotic solutions of

what is done, under the assumptions of negligible value of the ratio *r*

might wish to generalize the procedure described above to higher order in

β

ψ

τ

τ

introducing a rapidly varying external magnetic field.

By considering, for the time being, only the time scale

flux variable (Eqs. (8a-b)) give the following steady-state solution for

φ

 ψ

( ) ( ) ( )

+ =

 φ τ

( )

0

ψ θ

πβ

ψθ

+ =

+ =−

β

0 0,

*i*

+ =

 πψ θ

0,

τ

φ

β

d

d

*L R* ψτ

variable

φ

for the variables

1,

τ

already noticed that 2

φ

β

β

φ τ

πψ θ

and the following for the flux number variables:

0

ψ

1

ψ

θ

ψ

τ

τ

φ and ψ

wider range of validity of the analysis.

concerned with a single time scale, namely

with very small parameter

variable rapidly vanishes ( 1

φ

θ

**Figure 1.** Schematic representation of a two-junction quantum interferometer

This approach allows us not only to account for the regular part of the solution, as seen in (Grønbech-Jensen et al., 2003) and in (Romeo & De Luca, 2004), but also to consider its singular part. Moreover, as we shall see, the role of the two time variables will become evident in what follows, since one time scale is defined for Eq. (4a) and one for Eq. (4b). Consider then φ β( ,) τ and ψ β( ,) τ to be bounded, differentiable functions, and expand the sine and cosine functions appearing in Eqs. (4a-b) to first order in β . By then collecting all coefficients of identical power of β, we can obtain a system of equations for the functions , ( ) *<sup>k</sup>* β φ τ and ( ) ψ *k* θ , with *k* = 0, 1, describing the k-th order solutions for φ and ψ, respectively. These approximate solutions are determined according to the following sequential scheme. As a first step, we use Eq. (4b) to determine 0 ψ ( ) θ . We adopt the solution found and substitute it in Eq. (4a) to determine 0, ( ) β φ τ . The latter solution, on its turn, is substituted in Eq. (4b) to find 1 ψ ( ) θ and, finally, this solution is used in Eq. (4a) to find 1, ( ) β φ τ . Note, therefore, that for defining first order solutions, knowledge of zero-th order solutions is required. Furthermore, we assume that the initial conditions are the following:

$$
\rho(\beta,0) = \phi\_{0,\beta}(0) + \beta \phi\_{1,\beta}(0); \quad \nu(\beta,0) = \nu\_0(0) + \beta \nu\_1(0) \,. \tag{5}
$$

As for initial conditions, from Eq. (4b) we may notice that 0( ) ψ *ex* τ =ψ for β = 0 , in which case we cannot even define the time variable θ. This condition, however, is inherited by the function 0 ψ ( ) θ, since the following equalities are satisfied:

$$\left.\psi\_0\left(\frac{\pi}{\beta}\right)\right|\_{\beta=0} = \left.\psi\_0\left(\theta\right)\right|\_{\theta\to\ast\ast} = \left.\psi\_{ex}\right|\_{ex}.\tag{6}$$

Furthermore, we may also notice that ( ) ( ) <sup>0</sup> 0 0 *kk k* θ τ τ ψθ ψ ψ β <sup>=</sup> <sup>=</sup> = = , for *k* = 0, 1.

By the general procedure described above we get the following differential equations for the superconducting phase variables:

$$\begin{aligned} \frac{\mathbf{d}\,\phi\_{0,\beta}}{\mathbf{d}\,\tau} + \cos\left(\pi\nu\_{0}\left(\theta\right)\right)\sin\phi\_{0,\beta}\left(\tau\right) &= \frac{i\_{\text{B}}}{2}, \\ \frac{\mathbf{d}\,\phi\_{1,\beta}}{\mathbf{d}\,\tau} + \phi\_{1,\beta}\left(\tau\right)\cos\left(\pi\nu\_{0}\left(\theta\right)\right)\cos\phi\_{0,\beta}\left(\tau\right) &= \pi\nu\_{1}\left(\theta\right)\sin\left(\pi\nu\_{0}\left(\theta\right)\right)\sin\phi\_{0,\beta}\left(\tau\right); \text{ (b)} \end{aligned} \tag{7}$$

and the following for the flux number variables:

224 Superconductors – Materials, Properties and Applications

consider a new time variable

Consider then

ΙΒ

, ( ) *<sup>k</sup>* β φ

τ

find 1, ( ) β φ τ

following:

function 0 ψ ( ) θ

φ β( ,) τ and ψ β( ,) τ

coefficients of identical power of

turn, is substituted in Eq. (4b) to find 1

φ β

case we cannot even define the time variable

superconducting phase variables:

 and ( ) ψ *k* θ

φ β τ

where *n* is the number of initially trapped fluxons in the superconducting loop. Let us

= = and write the solutions for

 ≈ + ψ θ

*JJ2*

βψ θ

to be bounded, differentiable functions, and expand the

, we can obtain a system of equations for the functions

β φ

τ

 βψ β

ψ ( ) θ

> for β

, for *k* = 0, 1.

. This condition, however, is inherited by the

and, finally, this solution is used in Eq. (4a) to

φa and

. For simplicity, set *n*=0.

ΙΒ

. By then collecting all

. The latter solution, on its

φ and ψ,

. We adopt the

= 0 , in which

(6)

ψ

in the

*R t L*

> ψ β τ

This approach allows us not only to account for the regular part of the solution, as seen in (Grønbech-Jensen et al., 2003) and in (Romeo & De Luca, 2004), but also to consider its singular part. Moreover, as we shall see, the role of the two time variables will become evident in what follows, since one time scale is defined for Eq. (4a) and one for Eq. (4b).

*JJ1* 

, with *k* = 0, 1, describing the k-th order solutions for

. Note, therefore, that for defining first order solutions, knowledge of zero-th

 ψ

. →∞ *ex*

 ψ

0

 ψ

τ

θ

 ψ

0 *kk k*

τ

β <sup>=</sup> <sup>=</sup> = = 

= + = + . (5)

ψ *ex* τ =ψ

respectively. These approximate solutions are determined according to the following

order solutions is required. Furthermore, we assume that the initial conditions are the

0, 1, <sup>0</sup> <sup>1</sup> ( ,0) (0) (0); ( ,0) (0) (0)

θ

 ψ β

 β

( ) 0 0 0

= =

 ψθ

θ

By the general procedure described above we get the following differential equations for the

=

ψθ

β

 βφ

2

βφ τ

**Figure 1.** Schematic representation of a two-junction quantum interferometer

Ι1

sine and cosine functions appearing in Eqs. (4a-b) to first order in

solution found and substitute it in Eq. (4a) to determine 0, ( )

β

As for initial conditions, from Eq. (4b) we may notice that 0( )

, since the following equalities are satisfied:

β

Furthermore, we may also notice that ( ) ( ) <sup>0</sup>

τ ψ

 φ β

sequential scheme. As a first step, we use Eq. (4b) to determine 0

ψ ( ) θ

following form: 0, 1, <sup>0</sup> <sup>1</sup> ( , ) ( ) ( ); ( , ) ( ) ( )

β

τ

 ≈ + φ

τ θ

πβ

 β

> Ι2

$$\begin{aligned} \frac{d\,\Psi\_0}{d\,\theta} + \psi\_0(\theta) &= \psi\_{xx'} \\ \frac{d\,\Psi\_1}{d\,\theta} + \psi\_1(\theta) &= -2\sin\left(\pi\psi\_0\left(\theta\right)\right)\cos\phi\_{0,\beta}\left(2\pi\beta\theta\right). \text{ (b)} \end{aligned} \tag{8}$$

In Eqs. (7a-b) and (8a-b) we may notice the appearance of two different time scales the first, *L R* ψτ = , linked to flux motion in and out the superconducting ring, the second, <sup>0</sup> 2 *RIJ* φ τ π <sup>Φ</sup> <sup>=</sup> , pertaining to the dynamics of the superconducting phase difference value φ. We have already noticed that 2 ψ φ τ πβ τ <sup>=</sup> , so that, for negligible values of this ratio, the system behaves effectively as if the adiabatic time evolution of the superconducting phase difference variable φ could be studied by taking asymptotic solutions of ψ. In this case, therefore, we may first let the flux variable evolve, so that a stationary magnetic state is reached, and then solve for the superconducting phase difference time evolution of the system. This is exactly

what is done, under the assumptions of negligible value of the ratio *r* ψ φ τ τ<sup>=</sup> , in (Grønbech-

Jensen et al., 2003) and in (Romeo & De Luca, 2004). However, if one were to acquire the regular solution for the system dynamics, even when considering the approximate solution for the variables φ and ψ, one would follow the more general perturbation analysis described above, where the ratio *r* might not a priori be considered as negligible. Furthermore, considering that this ratio is proportional to the perturbation parameter β, one might wish to generalize the procedure described above to higher order in β to acquire a wider range of validity of the analysis.

Despite the fact that the more general approach allows extension to higher order approximations of the perturbation solutions, we wish to limit our analysis to the study of the electrodynamic properties of a two-junction or a multi-junction quantum interferometer with very small parameter β. Therefore, while in the present section we shall only be concerned with a single time scale, namely φτ, by assuming that the transient of the flux

variable rapidly vanishes ( 1 ψ φ τ τ << ), in Section 5 we shall consider both evolution time, by

introducing a rapidly varying external magnetic field.

By considering, for the time being, only the time scale φ τ , the dynamical equations for the flux variable (Eqs. (8a-b)) give the following steady-state solution for ψ <sup>0</sup> and ψ<sup>1</sup> :

$$\begin{aligned} \boldsymbol{\Psi}\_{0} &= \boldsymbol{\Psi}\_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}'} \end{aligned} \qquad \text{(a)}$$

$$\boldsymbol{\Psi}\_{1} = -2\sin\left(\pi\boldsymbol{\nu}\_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}\right)\cos\phi\_{0} - 2\pi\frac{\mathbf{d}\,\boldsymbol{\nu}\_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}}{\mathbf{d}\,\boldsymbol{\tau}}, \text{ (b)}$$

Effective Models of Superconducting Quantum Interference Devices 227

π

. Having found the time dependence of the variable 0

ψ <sup>0</sup> , 0 φ and ψ<sup>1</sup> , the

*i* = , 0.4 *<sup>B</sup>*

β

= in figs. 2a, 2b, and 2c for 1.6 *<sup>B</sup>*

= and 1.6 *<sup>B</sup>*

<sup>−</sup> <sup>=</sup> (15)

(14)

φ,

*i* = ,

*i* = and

, the circulating

*i* = and for

(16)

(17)

where ( )

γ

0

ω

function 1φ

0.4 *<sup>B</sup>*

0.6 *<sup>B</sup>*

a dotted line.

Solutions for 1

ψ

<sup>2</sup> <sup>1</sup> <sup>2</sup> cos 2 4

 = − πψ

Finally, in the case *a* = 1 , we have

where ( ) ( ) ( ) <sup>0</sup>

equation. Solutions for 0

φ

adopt a more general procedure.

superconducting ring, we have

0.1 and 0.3

*i* :

*ex* = , 0.5

ψ

equal to the pseudo-period of 0,

ψ

ψ*ex* and *<sup>B</sup>*

*i* = , and 0.6 *<sup>B</sup>*

order in the parameter

*ex*

*B*

and

χ

*a*

( ) <sup>0</sup>

2

*a*

φ

<sup>−</sup> <sup>=</sup>

0

1

τ

can be found by Eq. (10b), which is a standard first order linear differential

*i* = , respectively, along with the solution obtained by numerically

πψ

. Notice that in the case of time-dependent bias currents one should

 φ τ

( )

πψ

2

cos

π

*<sup>B</sup> ex*

<sup>2</sup> .

*ex* = , respectively. The period *T* of these curves is

which is given by the following expression in terms of

<sup>−</sup> <sup>±</sup> <sup>=</sup> + −

πψ

integrating Eqs. (4a-b). In figs. 2b-c the first order approximation of the solution is shown as

*i* = , respectively. The above analysis thus leads to a solution in a closed form, to first

. *ex <sup>S</sup>*

For an arbitrary value *n*, which represents the number of fluxons initially trapped in the

() ( ) () 1 0 2sin cos . *<sup>S</sup> ex i* = =−

 πψ

β

ω

−

*a*

.

<sup>0</sup> <sup>2</sup>

() () <sup>1</sup> ( ) 0 0 sgn 2tan , 2 2 *B i*

are shown for cos 0.3 ( ) *ex*

are shown in figs. 3a-c for cos 0.3 ( ) *ex*

As a simple application, let us calculate, to first order in the parameter

ψ τ

β φ

*T*

=

 *ex* = , 0.7 and 0.9 ψ

current *iS* in the circuit, normalized to *<sup>J</sup> I* , given by (Barone & Paternò, 1982):

*i* ψ ψ

Graphs of circulating currents are shown in figs. 4a, 4b, 4c for *n* even, 2.2 *<sup>B</sup>*

2

 <sup>−</sup> 

2

*i*

tan

*i*

φ τ

tan sgn 2 4 *a* φ

φ

β

 <sup>±</sup> = + 

0

π

<sup>1</sup> can be found by Eq. (11b) by substitution. Finally, by knowledge of

where the term <sup>d</sup> <sup>2</sup> d ψ *ex* π τ has been inserted in Eq. (9b), in order to correctly take account of first order contributions in β, and the subscript β in 0, ( ) β φ τ has been elided, as it will be done for 1, ( ) βφ τ from this point on, since these functions will not depend on βin this limit.

In order to obtain some preliminary results, we start by considering ψ *ex* as constant. By the general procedure schematized above we get the following differential equations for the superconducting phase variables:

$$\begin{aligned} \frac{d\phi\_0}{d\tau} + \cos\left(\pi\nu\_{ex}\right)\sin\phi\_0\left(\tau\right) &= \frac{i\_B}{2}, \\ \frac{d\phi\_1}{d\tau} + \phi\_1(\tau)\cos\left(\pi\nu\_{ex}\right)\cos\phi\_0\left(\tau\right) &= \pi\nu\_1(\tau)\sin\left(\pi\nu\_{ex}\right)\sin\phi\_0\left(\tau\right), \text{ (b)} \end{aligned} \tag{10}$$

and the following for the flux number variables:

$$\begin{aligned} \boldsymbol{\Psi}\_{0} &= \boldsymbol{\Psi}\_{ex'} \\ \boldsymbol{\Psi}\_{1}(\boldsymbol{\tau}) &= -2\sin\left(\pi\boldsymbol{\psi}\_{ex}\right)\cos\phi\_{0}\left(\boldsymbol{\tau}\right). \text{ (b)} \end{aligned} \tag{11}$$

According to the scheme described above, by having already set ψ ψ <sup>0</sup> *ex* = , we may now solve for ( ) <sup>0</sup> φ τ in Eq. (10a). Let us therefore briefly discuss how to obtain this solution. In the case ( ) <sup>1</sup> <sup>1</sup> 2cos *B ex i* πψ *a* = > , which characterizes the running state of the junctions, we

have

$$\phi\_0\left(\tau\right) = 2\tan^{-1}\left[a + \sqrt{1 - a^2}\tan\left(\mathcal{yr} + \tan^{-1}\tilde{\xi}\_0\right)\right]. \quad + 2\mathbf{k}\,\pi\tag{12}$$

where k is an integer, ( ) <sup>2</sup> <sup>1</sup> <sup>2</sup> cos 2 4 *B ex i* γ = − πψ and ( ) <sup>0</sup> <sup>0</sup> <sup>2</sup> 0 tan 2 1 *a a* φ ξ <sup>−</sup> <sup>=</sup> − .

On the other hand, in the case ( ) <sup>1</sup> <sup>1</sup> 2cos *B ex i* πψ *a* = < , which characterizes the superconducting state of the junctions, we have

$$\phi\_0\left(\tau\right) = 2\tan^{-1}\left[a + \sqrt{a^2 - 1}\tanh\left(\tilde{\mathcal{J}}\tau + \tanh^{-1}\left(\mathcal{Z}\_0\right)\right)\right],\tag{13}$$

### Effective Models of Superconducting Quantum Interference Devices 227

$$\text{where } \tilde{\gamma} = \frac{1}{2} \sqrt{\cos^2 \left( \pi \psi\_{ex} \right) - \frac{{i\_B}^2}{4}} \text{ and } \underset{\mathcal{X}\_0}{\text{and }} \underset{\mathcal{X}\_0}{\text{and }} \frac{\tan \left( \frac{\phi\_0 \left( 0 \right)}{2} \right) - a}{\sqrt{a^2 - 1}}.$$

Finally, in the case *a* = 1 , we have

226 Superconductors – Materials, Properties and Applications

d ψ *ex* π

τ

where the term <sup>d</sup> <sup>2</sup>

done for 1, ( ) β φ τ

solve for ( ) <sup>0</sup> φ τ

have

the case ( )

2cos *B*

*i* πψ*a*

*ex*

φ τ

where k is an integer, ( )

γ

first order contributions in

superconducting phase variables:

0

φ

τ

φ

τ

*d d d*

1

φ τ

and the following for the flux number variables:

0

=

ψ ψ

ψ

β

*ex*

In order to obtain some preliminary results, we start by considering

( ) ()

+ =

0

=

According to the scheme described above, by having already set

<sup>2</sup> <sup>1</sup> <sup>2</sup> cos 2 4 *B*

*i*

= −

On the other hand, in the case ( )

superconducting state of the junctions, we have

φ τ ψ ψ

ψ τ

<sup>1</sup> <sup>1</sup>

*ex*

+ =

 πψ

πψ

*d i*

0

 φ τ πψ τ

*ex*

 φ τ

( )

*ex*

, (a)

*ex*

has been inserted in Eq. (9b), in order to correctly take account of

τ

has been elided, as it will be

β

ψ

 φ τ

= − (11)

ψ ψ

π

*a*

.

= < , which characterizes the

in this limit.

(10)

<sup>0</sup> *ex* = , we may now

*ex* as constant. By the

 in 0, ( ) β φ

= − <sup>−</sup> (9)

τ

ψ

<sup>d</sup> 2sin cos 2 , (b) <sup>d</sup>

β

from this point on, since these functions will not depend on

general procedure schematized above we get the following differential equations for the

() ( ) () () ( ) ()

*ex ex*

cos sin , (a) <sup>2</sup>

cos cos sin sin , (b)

, (a) 2sin cos . (b)

 φ τ

in Eq. (10a). Let us therefore briefly discuss how to obtain this solution. In

= > , which characterizes the running state of the junctions, we

ξ

( ) <sup>0</sup>

 <sup>−</sup> <sup>=</sup> −

2

1

φ

0

*a*

 χ

− − = +− + + (12)

<sup>0</sup> <sup>2</sup>

<sup>1</sup> <sup>1</sup>

− − = +− + (13)

tan

 πψ

1 0 1 0

() ( ) ()

*ex*

1 0

( ) ( ) 12 1 0 0

*ex*

2cos *B*

( ) ( ) ( ) 1 2 <sup>1</sup> 0 0

2tan *a a* 1 tanh tanh ,

πψ

 2tan 1 tan tan . 2k *a a* γτ

and

*i* πψ*a*

ξ

*ex*

γτ

 πψ

*B*

 φ π

1 0

, and the subscript

 πψ

$$\phi\_0\left(\tau\right) = \text{sgn}\left(a\right) \left[ 2\tan^{-1}\left(\frac{i\_B\tau}{2} + a\rho\_0^{\left(\pm\right)}\right) - \frac{\pi}{2} \right],\tag{14}$$

where ( ) ( ) ( ) <sup>0</sup> 0 0 tan sgn 2 4 *a* φ π ω <sup>±</sup> = + . Having found the time dependence of the variable 0 φ,

ψ <sup>1</sup> can be found by Eq. (11b) by substitution. Finally, by knowledge of ψ <sup>0</sup> , 0 φ and ψ <sup>1</sup> , the function 1φ can be found by Eq. (10b), which is a standard first order linear differential equation. Solutions for 0 φ are shown for cos 0.3 ( ) *ex* πψ = in figs. 2a, 2b, and 2c for 1.6 *<sup>B</sup> i* = , 0.4 *<sup>B</sup> i* = , and 0.6 *<sup>B</sup> i* = , respectively, along with the solution obtained by numerically integrating Eqs. (4a-b). In figs. 2b-c the first order approximation of the solution is shown as a dotted line.

Solutions for 1φ are shown in figs. 3a-c for cos 0.3 ( ) *ex* πψ = and 1.6 *<sup>B</sup> i* = , 0.4 *<sup>B</sup> i* = and 0.6 *<sup>B</sup> i* = , respectively. The above analysis thus leads to a solution in a closed form, to first order in the parameter β. Notice that in the case of time-dependent bias currents one should adopt a more general procedure.

As a simple application, let us calculate, to first order in the parameter β, the circulating current *iS* in the circuit, normalized to *<sup>J</sup> I* , given by (Barone & Paternò, 1982):

$$
\dot{q}\_S = \frac{\Psi - \Psi\_{ex}}{\beta}.\tag{15}
$$

For an arbitrary value *n*, which represents the number of fluxons initially trapped in the superconducting ring, we have

$$i\_S = \psi\_1(\tau) = -2\sin\left(\pi\nu\_{ex}\right)\cos\phi\_0(\tau). \tag{16}$$

Graphs of circulating currents are shown in figs. 4a, 4b, 4c for *n* even, 2.2 *<sup>B</sup> i* = and for 0.1 and 0.3 ψ *ex* = , 0.5 ψ *ex* = , 0.7 and 0.9 ψ *ex* = , respectively. The period *T* of these curves is equal to the pseudo-period of 0,β φ which is given by the following expression in terms of ψ *ex* and *<sup>B</sup> i* :

$$T = \frac{2\pi}{\sqrt{\left(\frac{\dot{I}\_B}{2}\right)^2 - \cos^2\left(\pi\varphi\_{ex}\right)}}\,\tag{17}$$

Notice that the lowest value of the period is obtained for 0.5 ψ *ex* = and that the curves for 0.1 and 0.9 ψ *ex* = and for 0.3 and 0.7 ψ *ex* = , although having the same period, as it can be argued from Eq. (17), are not equal.

Effective Models of Superconducting Quantum Interference Devices 229

as calculated by the procedure described in the text by setting

**Figure 3.** First order correction to ( ) <sup>0</sup>

= and: a) 1.6 *<sup>B</sup>*

*β=*0.02, cos 0.3 ( ) *ex* πψ

φ τ

15

10

*φ*1

5

0

*i* = ; b) 0.4 *<sup>B</sup>*

*i* = ; c) 0.6 *<sup>B</sup>*

(a)

*τ*

0 10 20 30 40

4

*φ*1

3

2

*φ*1

1

0

(b)

*τ*

0 10 20 30 40

*i* = .

(c)

*τ*

0 10 20 30 40

**Figure 2.** Average phase difference for cos 0.3 ( ) *ex* πψ = and: a) 1.6 *<sup>B</sup> i* = ; b) 0.4 *<sup>B</sup> i* = ; c) 0.6 *<sup>B</sup> i* = . Dotted blue lines represent ( ) <sup>0</sup> φ τ for *β* = 0.02, full red lines represent φ ( ) τ as calculated to first order for *β* = 0.02, and the dashed black lines represents the numerical solution of the complete system. In a) the first order approximation of the solution is not shown, for clarity reasons.

*ex* = and for 0.3 and 0.7

argued from Eq. (17), are not equal.

ψ

> 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

> > 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

*φ*

*φ*

*φ*

**Figure 2.** Average phase difference for cos 0.3 ( ) *ex*

φ τ

blue lines represent ( ) <sup>0</sup>

πψ

*β* = 0.02, and the dashed black lines represents the numerical solution of the complete system. In a) the

(a)

*τ*

0 10 20 30 40

(b)

*τ*

0 10 20 30 40

for *β* = 0.02, full red lines represent

first order approximation of the solution is not shown, for clarity reasons.

= and: a) 1.6 *<sup>B</sup>*

(c)

*τ*

0 10 20 30 40

φ ( ) τ

*i* = ; b) 0.4 *<sup>B</sup>*

*i* = ; c) 0.6 *<sup>B</sup>*

as calculated to first order for

*i* = . Dotted

0.1 and 0.9

ψ

Notice that the lowest value of the period is obtained for 0.5

ψ

*ex* = , although having the same period, as it can be

*ex* = and that the curves for

**Figure 3.** First order correction to ( ) <sup>0</sup> φ τ as calculated by the procedure described in the text by setting *β=*0.02, cos 0.3 ( ) *ex* πψ = and: a) 1.6 *<sup>B</sup> i* = ; b) 0.4 *<sup>B</sup> i* = ; c) 0.6 *<sup>B</sup> i* = .

Effective Models of Superconducting Quantum Interference Devices 231

**Figure 5.** Circulating current *iS* in the presence of an oscillating applied magnetic field shown as a function

(c)

0 50 100 150 200 250

*τ*

(a)

0 50 100 150 200 250

*τ*

0.6 0.4 0.2 0.0 0.2 0.4 0.6

0.6 0.4 0.2 0.0 0.2 0.4 0.6

0.6 0.4 0.2 0.0 0.2 0.4 0.6

*iS*

*iS*

(b)

0 50 100 150 200 250

*τ*

( ) cos ,

= 0.09 . No initially trapped flux in the superconducting loop is present.

ωτ

(18)

of the normalized time for *A*=0, *B*=0.1, *<sup>i</sup>* 2.5 *<sup>B</sup>* <sup>=</sup> , *β*=0.01 and for normalized frequencies equal to: a)

ψ *ex* τ= + *A B*

ω

= 0.03 ; b)

ω= 0.06 ; c)

*t*

Now, since <sup>0</sup>

ω

*iS*

<sup>Φ</sup> <sup>=</sup> , we can write

2 *<sup>J</sup>*

π

*RI* τ

**Figure 4.** Circulating current *iS* as a function of the normalized time for null values of the initially trapped flux and for 2.2 *<sup>B</sup> i* = and: a) 0.1 ψ *ex* = (orange), 0.3 ψ *ex* = (cyan) ; b) 0.5 ψ *ex* = ; c) 0.7 ψ *ex* = (cyan), 0.9 ψ*ex* = (orange).

The above results have been obtained for the magnetic response of the system in the presence of a constant applied flux. In the following we shall analyze the electrodynamic response of the two junction quantum interferometer in the presence of a time-dependent external flux. For this purpose, we shall take a sinusoidal forcing term, in such a way that ( ) cos *ex* ψ *t AB t* = + ω , where *A* is the normalized d. c. component of the applied flux and *B* is the normalized amplitude of the a. c. signal.

**Figure 5.** Circulating current *iS* in the presence of an oscillating applied magnetic field shown as a function of the normalized time for *A*=0, *B*=0.1, *<sup>i</sup>* 2.5 *<sup>B</sup>* <sup>=</sup> , *β*=0.01 and for normalized frequencies equal to: a) ω = 0.03 ; b) ω = 0.06 ; c) ω= 0.09 . No initially trapped flux in the superconducting loop is present.

Now, since <sup>0</sup> 2 *<sup>J</sup> t RI* τ π <sup>Φ</sup> <sup>=</sup> , we can write

230 Superconductors – Materials, Properties and Applications

*iS*

1.5 1.0 0.5 0.0 0.5 1.0 1.5

2

1

0

*iS*

1

2

**Figure 4.** Circulating current *iS* as a function of the normalized time for null values of the initially

*ex* = (orange), 0.3

(a)

*τ*

0 5 10 15 20 25 30 35

(b)

0 5 10 15 20 25 30 35

*τ*

The above results have been obtained for the magnetic response of the system in the presence of a constant applied flux. In the following we shall analyze the electrodynamic response of the two junction quantum interferometer in the presence of a time-dependent external flux. For this purpose, we shall take a sinusoidal forcing term, in such a way that

ψ

(c)

*τ*

0 5 10 15 20 25 30 35

, where *A* is the normalized d. c. component of the applied flux and *B* is

*ex* = (cyan) ; b) 0.5

ψ

 *ex* = ; c) 0.7 ψ*ex* =

*i* = and: a) 0.1 ψ

1.5 1.0 0.5 0.0 0.5 1.0 1.5

*iS*

trapped flux and for 2.2 *<sup>B</sup>*

*ex* = (orange).

ω

the normalized amplitude of the a. c. signal.

(cyan), 0.9 ψ

( ) cos *ex*

*t AB t* = +

ψ

$$\Psi\_{ex}\left(\tau\right) = A + B\cos\tilde{\partial}\tau,\tag{18}$$

where <sup>0</sup> 2 *RIJ* ω ω π <sup>Φ</sup> <sup>=</sup> . By considering normalized frequencies of the order of <sup>1</sup> φ τ<sup>−</sup> , we Effective Models of Superconducting Quantum Interference Devices 233

 φ τ

 ωτ

*ex* =+ ≈ − *A B A B* ,

(19)

 π

(20)

(21)

 π

can be determined by solving a

*=*0.01 and 2 5. *<sup>B</sup>*

*i* = is

*i* ,

(23)

( ) ( ) ( )

*A* cos Eq. (19a) can be written in the following form:

*d*

 and *bB A* = π

the parameter *b*, so that, by setting () () () 00 1

0

*d d*

η

τ

η

*d*

As before, the above expression is equal to ( )

defined as the maximum value of the current bias *<sup>B</sup>*

ω

τ

1

*a*

found, by substituting in Eq. (19b), the solution for ( ) <sup>1</sup>

φ

τ

 π

φ τ ητ

0

*<sup>B</sup> d i a*

+ =

η τ

η τ η τ

= + *b* to be a known expression, we can then write:

τ

πψ ( )

*ex*

πψ τ

*d i*

 πψ τ

+ =

0

1 2

 φ τ

+ = −

0

φ

τ

φ

τ

π

( ) π

 ωτ

where *a A* = cos( )

() () () 00 1

frequency values

oscillating signal.

of Eq. (21a) we write:

Therefore, we have

 ητ

φ τ ητ φ τ

*d d d*

() () ( ) ( ) ( ) ( ) ( ) ( )

 π

cos sin , (a) <sup>2</sup>

cos cos sin sin 2 (b)

<sup>0</sup> cos sin , <sup>2</sup>

 ητ= + *b* , we can write:

() () () ( )

 η τ

> φ τ

cos sin cos (b)

sin , (a) <sup>2</sup>

 φ τ

 πψ τ

> π π

sin( ) . In Eq. (20) we find a perturbed solution in terms of

 ωτ

> πω

 ωτ+ + *b B* (22)

and represents the circulating current *iS* in

β

*i* which can be injected in the two

1 0 0

*ex ex*

By noticing, however, that cos( ) ( ) cos cos cos sin ( ) () πψ τ

( ) () <sup>0</sup>

01 0

+ =

Notice then that the solutions to the above equations can be found by exactly the same procedure described for the case of constant applied fields. Once the solution for ( ) <sup>0</sup> *φ τ* is

first order linear differential equation with time-dependent coefficients. Assuming thus

( ) ( ) ( ) () () 0 1 2sin cos 2 sin *<sup>S</sup> ex i* = −

 η τ η τ

ψ 1 τ

= 0.03, 0.06, 0.09 , respectively, and for *A=0*,

the circuit. In figs. 5a, 5b, and 5c the time dependence of the current *iS* for normalized

represented. In these graphs we notice that the oscillating patterns, which we have already detected in the constant applied field case, are modulated by the externally applied

Another important quantity to be measured in these systems is the critical current *<sup>c</sup>*

junction interferometer without giving rise to dissipation. By considering the stationary case

( ) ( ) ( ) <sup>0</sup> 2cos sin . *<sup>B</sup> ex i* = πψ

2 cos cos . ( ) ( ) *<sup>c</sup>*

τ

 φτ

*i* = + *πA πB ωτ* (24)

*<sup>B</sup> <sup>d</sup> <sup>i</sup> a b*

+ − =

ωτ

*B*

 φ τ

may still use the analysis described above. In Section 5 we shall relax the latter hypothesis. As also specified above, a different approach, which will be developed in Section 5, needs to be used when very rapidly varying applied fields are applied to the system. We shall assume that the normalized amplitude *B* of the oscillating signal is much less than one (*B*<<*1*). The perturbation analysis is then carried out in a way at all similar as done above.

**Figure 6.** (a) Time average value *<sup>c</sup> i* of the critical current *<sup>c</sup> i* as a function of the amplitude *B* of the oscillating magnetic flux for *A*=0.1 (blue) and *A*=0.5 (orange). (b) Time average value *<sup>c</sup> i* of the critical current *<sup>c</sup> i* as a function of *A* for *B*=0.1 (blue) and *B*=0.5 (orange). Both figures are obtained by taking null values of the initially trapped flux.

We start by setting, by Eq. (9a) and (9b), () () ( ) <sup>0</sup> cos ψ τ = =+ ψ *ex* τ *A B* ωτ and ( ) ( ) 1 0 <sup>d</sup> 2sin cos 2 d *ex ex* = − − ψ ψ πψ τ φ π τ and solve the equations for the phase differences

$$\begin{aligned} \frac{d\phi\_0}{d\tau} &+ \cos\left(\pi\nu\_{ex}(\tau)\right)\sin\phi\_0\left(\tau\right) = \frac{i\_B}{2}, & \text{(a)}\\ \frac{d\phi\_1}{d\tau} &+ \phi\_1(\tau)\cos\left(\pi\nu\_{ex}(\tau)\right)\cos\phi\_0\left(\tau\right) = -\pi\sin^2\left(\pi\nu\_{ex}(\tau)\right)\sin\left(2\phi\_0(\tau)\right) \tag{b} \end{aligned} \tag{19}$$

By noticing, however, that cos( ) ( ) cos cos cos sin ( ) () πψ τ π π ωτ π π *ex* =+ ≈ − *A B A B* , ( ) π ωτ*A* cos Eq. (19a) can be written in the following form:

$$\frac{d\phi\_0}{d\tau} + \left(a - b\cos\tilde{\alpha}\sigma\right)\sin\phi\_0\left(\tau\right) = \frac{\dot{l}\_B}{2},\tag{20}$$

where *a A* = cos( ) π and *bB A* = π π sin( ) . In Eq. (20) we find a perturbed solution in terms of the parameter *b*, so that, by setting () () () 00 1 φ τ ητ ητ= + *b* , we can write:

$$\begin{aligned} \frac{d\eta\_0}{d\tau} + a \sin \eta\_0(\tau) &= \frac{i\_B}{2}, \\ \frac{d\eta\_1}{d\tau} + a \cos \eta\_0(\tau) \eta\_1(\tau) &= \sin \eta\_0(\tau) \cos \left(\delta \sigma \tau\right) \quad \text{(b)} \end{aligned} \tag{21}$$

Notice then that the solutions to the above equations can be found by exactly the same procedure described for the case of constant applied fields. Once the solution for ( ) <sup>0</sup> *φ τ* is found, by substituting in Eq. (19b), the solution for ( ) <sup>1</sup> φ τ can be determined by solving a first order linear differential equation with time-dependent coefficients. Assuming thus () () () 00 1 φ τ ητ ητ= + *b* to be a known expression, we can then write:

$$i\_{\mathbb{S}} = -2\sin\left(\pi\nu\_{ex}\left(\tau\right)\right)\cos\left(\eta\_{0}\left(\tau\right) + b\eta\_{1}\left(\tau\right)\right) + 2\pi\tilde{o}\mathcal{B}\sin\tilde{o}\sigma\tag{22}$$

As before, the above expression is equal to ( ) ψ 1 τ and represents the circulating current *iS* in the circuit. In figs. 5a, 5b, and 5c the time dependence of the current *iS* for normalized frequency values ω = 0.03, 0.06, 0.09 , respectively, and for *A=0*, β*=*0.01 and 2 5. *<sup>B</sup> i* = is represented. In these graphs we notice that the oscillating patterns, which we have already detected in the constant applied field case, are modulated by the externally applied oscillating signal.

Another important quantity to be measured in these systems is the critical current *<sup>c</sup> i* , defined as the maximum value of the current bias *<sup>B</sup> i* which can be injected in the two junction interferometer without giving rise to dissipation. By considering the stationary case of Eq. (21a) we write:

$$i\_B = 2\cos\left(\pi\nu\_{ex}\left(\tau\right)\right)\sin\phi\_0\left(\tau\right). \tag{23}$$

Therefore, we have

232 Superconductors – Materials, Properties and Applications

0.0

0.5

1.0

*ic*

1.5

2.0

0.5

1.0

*ic*

1.5

<sup>Φ</sup> <sup>=</sup> . By considering normalized frequencies of the order of <sup>1</sup>

may still use the analysis described above. In Section 5 we shall relax the latter hypothesis. As also specified above, a different approach, which will be developed in Section 5, needs to be used when very rapidly varying applied fields are applied to the system. We shall assume that the normalized amplitude *B* of the oscillating signal is much less than one (*B*<<*1*). The perturbation analysis is then carried out in a way at all similar as

*i* of the critical current *<sup>c</sup>*

(b)

0.0 0.5 1.0 1.5 2.0 2.5 3.0

*A*

(a)

0.0 0.5 1.0 1.5 2.0 2.5 3.0

*B*

We start by setting, by Eq. (9a) and (9b), () () ( ) <sup>0</sup> cos

*i* as a function of *A* for *B*=0.1 (blue) and *B*=0.5 (orange). Both figures are obtained by taking

ψ τ

and solve the equations for the phase differences

oscillating magnetic flux for *A*=0.1 (blue) and *A*=0.5 (orange). (b) Time average value *<sup>c</sup>*

d *ex*

τ

 π ψ *i* as a function of the amplitude *B* of the

 = =+ ψ *ex* τ*A B*

*i* of the critical

ωτand

φ τ<sup>−</sup> , we

where <sup>0</sup>

done above.

ω

2 *RIJ*

**Figure 6.** (a) Time average value *<sup>c</sup>*

null values of the initially trapped flux.

<sup>d</sup> 2sin cos 2

 φ

( ) ( ) 1 0

 πψ τ

*ex* = − −

current *<sup>c</sup>*

ψ

π

 ω

$$i\_c = 2 \left| \cos \left( \pi A + \pi B \cos \left( \tilde{\omega} \tau \right) \right) \right|. \tag{24}$$

Noticing that the time-averaged value <*ic*> of the critical current does not depend on the normalized frequency, it can be calculated in terms of solely *A* and *B*, the results being shown in fig. 6a and fig. 6b for null values of the initially trapped flux. In particular, in fig. 6a <*ic*> is shown as a function of the applied magnetic field amplitude *B*, for *A=0.1* and *A=0.4*, while in fig. 6b, <*ic*> vs. *A* curves are shown for *B=0.1* and *B=0.2*. In the curves in fig. 6a we notice Fraunhofer-like oscillations, while ordinary cosinusoidal oscillations are present in fig. 6b.

Effective Models of Superconducting Quantum Interference Devices 235

0 <sup>1</sup> *<sup>J</sup> LI*

across the two junctions in the cell, so

<sup>Φ</sup> , where Φ0 is

β= <<

Φ linked to the k-th cell of the array is related to the

− + Ψ= (25)

φ*N-3* φ*N-2* φ*N-1* 

**3. Multi-junction quantum interference devices** 

inductance *L* of the horizontal upper branches to be such that

*k* Φ Ψ =

0 *k*

φ

φ*3* φ

 φπ

**Figure 7.** Schematic representation of a one-dimensional array of equal Josephson junctions.

Φ Φ , where 0 *<sup>S</sup>* is the cell area, we may set

When an external magnetic field **H** , orthogonal to the plane of the array, is applied to the system so that the normalized geometric applied flux through each cell is

*I0 I1 I2 I3* I0 I0 *IN-3 IN-2 IN-1 IN*

*IB/2 IB/2* 

*L L L L L L*

normalized magnetic flux

that:

superconducting phase differences *<sup>k</sup>*−<sup>1</sup>

where *<sup>k</sup> n* is an integer and *k N* = 1, 2,..., .

0 0

φ*1* φ*2*

0 0

μ*HS* Ψ= =

*ex*

φ*0*

Φ

*ex*

In this section we shall consider the one-dimensional Josephson junction array (1D-JJA) represented in fig. 7, consisting of *N*+1 identical overdamped junctions connected in parallel. In this figure we notice that the bias current *BI* is evenly applied to the two external branches of the array. By assuming perfectly identical overdamped Josephson junctions with resistive parameter *R* and maximum Josephson current *<sup>J</sup> I* , we take the

the elementary flux quantum. This parameter has been defined as each single loop could be compared to a two-junction quantum interference device. By fluxoid quantization, the

> and *<sup>k</sup>* φ

<sup>1</sup> 2 2, *kk k k n* <sup>−</sup>

 π

For what seen above, the electrodynamic properties of a symmetric quantum interferometer containing two identical junctions with negligible capacitance can be studied by means of a perturbation approach in the parameter β, whose value gives the strength of the electromagnetic coupling between the two junction in the system. The analysis is rather similar to what done in other works in the literature (Grønbech-Jensen et al., 2003; Romeo & De Luca, 2004). However, in the present section we have presented a rather general procedure to obtain the solution to the problem to first order in the parameter β. Considering at first transient solutions, we have noticed that the function ψ β( ) ,θ governs fluxon dynamics, where θ is the ordinary time *t*, normalized to the characteristic circuital time constant *<sup>L</sup> R* τ<sup>Ψ</sup> = . By this more general approach it becomes evident that the

characteristic time constant <sup>0</sup> 2 *RIJ* <sup>Φ</sup> <sup>=</sup> φ τ π of the dynamics of the average superconducting

phase difference φ is different from the fluxon dynamics characteristic time τ <sup>Ψ</sup> , so that the asymptotic solution for the system, proposed in the analyses carried out in (Grønbech-Jensen et al., 2003) and in (Romeo & De Luca, 2004), acquires a more precise meaning in this context. Indeed, when the parameter β is sufficiently small to allow, for finite values of τ, an asymptotic evaluation of 0 ( ) <sup>2</sup> τ πβ ψ and 1 ( ) <sup>2</sup> τ πβ ψ , the general solution given in the present work coincides with the asymptotic perturbed solution proposed in (Grønbech-Jensen et al., 2003) and in (Romeo & De Luca, 2004) in the limit of negligible junction capacitance.

The perturbation analysis has been first carried out for a constant applied magnetic flux. Successively, since it could be experimentally possible to force the system with a timedependent magnetic field, it is noted that the perturbed solution for the flux number ψ, obtained for a sinusoidal magnetic flux, needs careful evaluation. In order to exhibit experimentally detectable quantities, the circulating current *iS* is evaluated as a function of time, for different values of the frequency of the forcing field, whose a. c. component is assumed to be small. Finally, the time average <*ic*> of the critical current of the device has been studied both as a function of the d. c. component *A* and of the small amplitude *B* of the oscillating part of the applied flux. In these curves two characteristic behaviors have been detected: A Fraunhofer-like pattern in <*ic*> vs. *B* curves; independence of <*ic*> from ω . We shall see in Section 5 of the present Chapter that dependence from ω appears in a more general case, i. e., when ω≈ 1 and the coefficient *B* is not assumed to be small.

## **3. Multi-junction quantum interference devices**

234 Superconductors – Materials, Properties and Applications

perturbation approach in the parameter

*R*

θ

fluxon dynamics, where

τ

characteristic time constant <sup>0</sup>

φ

context. Indeed, when the parameter

asymptotic evaluation of 0 ( ) <sup>2</sup>

general case, i. e., when

time constant *<sup>L</sup>*

phase difference

fig. 6b.

Noticing that the time-averaged value <*ic*> of the critical current does not depend on the normalized frequency, it can be calculated in terms of solely *A* and *B*, the results being shown in fig. 6a and fig. 6b for null values of the initially trapped flux. In particular, in fig. 6a <*ic*> is shown as a function of the applied magnetic field amplitude *B*, for *A=0.1* and *A=0.4*, while in fig. 6b, <*ic*> vs. *A* curves are shown for *B=0.1* and *B=0.2*. In the curves in fig. 6a we notice Fraunhofer-like oscillations, while ordinary cosinusoidal oscillations are present in

For what seen above, the electrodynamic properties of a symmetric quantum interferometer containing two identical junctions with negligible capacitance can be studied by means of a

β

electromagnetic coupling between the two junction in the system. The analysis is rather similar to what done in other works in the literature (Grønbech-Jensen et al., 2003; Romeo & De Luca, 2004). However, in the present section we have presented a rather general procedure to obtain the solution to the problem to first order in the parameter

is different from the fluxon dynamics characteristic time

asymptotic solution for the system, proposed in the analyses carried out in (Grønbech-Jensen et al., 2003) and in (Romeo & De Luca, 2004), acquires a more precise meaning in this

> πβ

work coincides with the asymptotic perturbed solution proposed in (Grønbech-Jensen et al.,

The perturbation analysis has been first carried out for a constant applied magnetic flux. Successively, since it could be experimentally possible to force the system with a timedependent magnetic field, it is noted that the perturbed solution for the flux number

obtained for a sinusoidal magnetic flux, needs careful evaluation. In order to exhibit experimentally detectable quantities, the circulating current *iS* is evaluated as a function of time, for different values of the frequency of the forcing field, whose a. c. component is assumed to be small. Finally, the time average <*ic*> of the critical current of the device has been studied both as a function of the d. c. component *A* and of the small amplitude *B* of the oscillating part of the applied flux. In these curves two characteristic behaviors have been

≈ 1 and the coefficient *B* is not assumed to be small.

Considering at first transient solutions, we have noticed that the function

2 *RIJ* <sup>Φ</sup> <sup>=</sup> φ τ

β

 and 1 ( ) <sup>2</sup> τ

2003) and in (Romeo & De Luca, 2004) in the limit of negligible junction capacitance.

detected: A Fraunhofer-like pattern in <*ic*> vs. *B* curves; independence of <*ic*> from

shall see in Section 5 of the present Chapter that dependence from

ω

ψ

τ

ψ

πβ π

, whose value gives the strength of the

is the ordinary time *t*, normalized to the characteristic circuital

of the dynamics of the average superconducting

, the general solution given in the present

ω

is sufficiently small to allow, for finite values of

<sup>Ψ</sup> = . By this more general approach it becomes evident that the

ψ β( ) ,θ

τ

<sup>Ψ</sup> , so that the

τ, an

> ψ,

ω. We

appears in a more

β.

governs

In this section we shall consider the one-dimensional Josephson junction array (1D-JJA) represented in fig. 7, consisting of *N*+1 identical overdamped junctions connected in parallel. In this figure we notice that the bias current *BI* is evenly applied to the two external branches of the array. By assuming perfectly identical overdamped Josephson junctions with resistive parameter *R* and maximum Josephson current *<sup>J</sup> I* , we take the

inductance *L* of the horizontal upper branches to be such that 0 <sup>1</sup> *<sup>J</sup> LI* β = << <sup>Φ</sup> , where Φ0 is

the elementary flux quantum. This parameter has been defined as each single loop could be compared to a two-junction quantum interference device. By fluxoid quantization, the normalized magnetic flux 0 *k k* Φ Ψ = Φ linked to the k-th cell of the array is related to the superconducting phase differences *<sup>k</sup>*−<sup>1</sup> φ and *<sup>k</sup>* φ across the two junctions in the cell, so that:

$$
\phi\_{k-1} - \phi\_k + 2\pi\Psi\_k = 2\pi n\_{k'} \tag{25}
$$

where *<sup>k</sup> n* is an integer and *k N* = 1, 2,..., .

**Figure 7.** Schematic representation of a one-dimensional array of equal Josephson junctions.

When an external magnetic field **H** , orthogonal to the plane of the array, is applied to the system so that the normalized geometric applied flux through each cell is 0 0 0 0 *ex ex* Φ μ *HS* Ψ= = Φ Φ , where 0 *<sup>S</sup>* is the cell area, we may set

$$\Psi\_k = \beta \sum\_{m=1}^k i\_m - \beta \frac{i\_B}{2} + \Psi\_{ex'} \tag{26}$$

Effective Models of Superconducting Quantum Interference Devices 237

 β

(c)

and of the *N* partial sums *ns* .

( *m* = 0,1 ) equal to zero, we obtain, for

(0) . *n ex s n* = Ψ (34)

2

<sup>−</sup> =− − <sup>−</sup> + + (36)

*B*

*i*

(32)

(33)

(35)

<sup>1</sup>*s* , which is the quantity

ββ

1 11

− + Ψ−

sin , (a) <sup>2</sup> 2

+ −

β

β

−

. Therefore, start by considering the dynamical equations of the system as

2 sin 2 sin , (b) <sup>2</sup>

The above analysis has been carried out essentially to write the dynamical equations in

We shall now develop a reduction of these variables to one, by assuming small values of the

(0) (1). *nn n ss s* ≈ + β

By substituting the above expression in Eq. (32b) and, after having multiplied both

For the first order corrections, on the other hand, we need to solve the following set of

(1) (1) (1) (1)

2

*i*

*B*

*i*

2

*B*

1 0 1 21 (1) (1) (1) 01 1

2

−= − + − −

, *n N* = 0, 1, 2,..., . By solving for (1)

( ) (1) 1 0

*<sup>s</sup> y y N N* <sup>=</sup>

*B*

*y ys s s s*

− − −

(1) (1) (1) (1)

<sup>2</sup> . <sup>2</sup>

1 1 1 , 12 1 *N*

*k*

*k*

*B*

*i*

(1) (1)

3

(1) (1)

2

304 3 2 1

−= − + − −

*yys s s s*

*N NNN*

102 1

−= − −

*yy s s*

203 2

−= − −

*yys s*

*N NN*

−= − −

*y ys s s*

−

*N i*

....

 

Substitution of the above results into Eq. (32a) gives:

β

*n n n n ex B*

+ +−− = −

2 sin 2 sin

+ + − + =−

 φ

*ds s s s*

 φ

*N N N ex ex*

*ds s ss s i*

1 1

− Ψ −Ψ

φ

β

we can set:

( )

0 0

− Ψ

*n*

 β

*s*

*B ex*

0 0

*N*

*s*

0 1 0

> φπ

 φπ

+ =+

*d s i*

φ

*d*

φ

τ

π

π

parameter

members by

*n N* = 1, 2,..., :

equations:

β

β

where ( ) <sup>0</sup> sin 2 *<sup>n</sup> ex y n* = +Ψ φ

required in Eq. (32a), we have:

 π

*d*

τ

*d*

τ

where *n N* = − 1, 2,..., 1 in Eq. (32b).

( )

terms of the effective superconducting phase difference 0

written in Eqs. (32a-b). For small values of the parameter

, by setting the coefficients of *<sup>m</sup>*

for *k N* =1 2, ,..., , where *<sup>k</sup> k J <sup>I</sup> <sup>i</sup> <sup>I</sup>* <sup>=</sup> is the normalized current flowing in the *k*-th branch

and *<sup>B</sup> B J <sup>I</sup> <sup>i</sup> <sup>I</sup>* <sup>=</sup> . The dynamical equation for each Josephson junction in the array is written by means of the resistively shunted junction (RSJ) model (Barone & Paternò, 1982) as follows:

$$\frac{d\phi\_k}{d\tau} + \sin\phi\_k = i\_{k'} \tag{27}$$

where *k N* = 0, 1, 2,..., and 0 <sup>2</sup> *<sup>J</sup> <sup>π</sup>RI <sup>τ</sup>* <sup>=</sup> *<sup>t</sup>* <sup>Φ</sup> . Eqs. (25-27) can be used to define the dynamics of the gauge-invariant superconducting phase difference *<sup>k</sup>* φ in terms of the forcing parameters Ψ*ex* and 0 *N B B k J k <sup>I</sup> i i I* <sup>=</sup> = = . In addition, the instantaneous voltage *v*( ) τ of the system can be obtained by setting:

$$\upsilon\left(\boldsymbol{\pi}\right) = \frac{1}{N+1} \sum\_{k=0}^{N} \frac{d\phi\_k}{d\boldsymbol{\pi}}.\tag{28}$$

Define now the partial sum *ns* ( 1 ≤ ≤ *n N* ) of the normalized fluxes as follows:

$$s\_n = \sum\_{k=1}^n \Psi\_k.\tag{29}$$

By fluxoid quantization (Eq. (25)), setting all *<sup>k</sup> n* 's to zero (under the hypothesis of zero initially trapped flux in the array), we can write:

$$
\phi\_k = \phi\_0 + 2\pi s\_{k'} \tag{30}
$$

for *k N* = 1, 2,..., , so that the dynamical equations (Eq. (27)) can be rewritten as follows:

$$\begin{aligned} \frac{d\phi\_0}{d\tau} + \sin\phi\_0 &= i\_0, & \text{(a)}\\ 2\pi \frac{d\mathbf{s}\_n}{d\tau} + \sin\left(\phi\_0 + 2\pi\mathbf{s}\_n\right) &= i\_n - i\_0 + \sin\phi\_0 & \text{(b)} \end{aligned} \tag{31}$$

where *n N* = 1, 2,..., in the second equation. Expressing now, by means of Eq. (26), the currents 0 *i* and *<sup>n</sup> i* in terms of the forcing terms Ψ*ex* and *<sup>B</sup> i* and of the partial sums of the flux variables *ns* , we may finally rewrite Eqs. (31a-b) as follows:

$$\frac{d\phi\_0}{d\tau} + \sin\phi\_0 = \frac{i\_B}{2} + \frac{s\_1 - \Psi\_{ex}}{\beta},\tag{a}$$

$$2\pi \frac{ds\_n}{d\tau} + \sin\left(\phi\_0 + 2\pi s\_n\right) - \sin\phi\_0 - \frac{s\_{n+1} - 2s\_n + s\_{n-1}}{\beta} = \frac{\Psi\_{ex} - s\_1}{\beta} - \frac{i\_B}{2}, \quad \text{(b)}\tag{32}$$

$$2\pi \frac{d s\_N}{d \tau} + \sin \left(\phi\_0 + 2\pi s\_N \right) - \sin \phi\_0 + \frac{s\_N - s\_{N-1}}{\beta} = \frac{\Psi\_{ex}}{\beta} - \frac{s\_1 - \Psi\_{ex}}{\beta} \tag{c}$$

where *n N* = − 1, 2,..., 1 in Eq. (32b).

236 Superconductors – Materials, Properties and Applications

*k*

*J <sup>I</sup> <sup>i</sup>*

for *k N* =1 2, ,..., , where *<sup>k</sup>*

and *<sup>B</sup> B*

follows:

Ψ*ex* and

currents 0

*i* and *<sup>n</sup>*

obtained by setting:

*<sup>I</sup> <sup>i</sup>*

*J*

where *k N* = 0, 1, 2,..., and

1

= Ψ = − +Ψ

*k m ex*

by means of the resistively shunted junction (RSJ) model (Barone & Paternò, 1982) as

sin , *<sup>k</sup> k k*

φ

φ

0 <sup>1</sup> . <sup>1</sup> *N k*

*d*

φ

τ

*k*

*N d* <sup>=</sup> <sup>=</sup> <sup>+</sup>

> 1 .

=

By fluxoid quantization (Eq. (25)), setting all *<sup>k</sup> n* 's to zero (under the hypothesis of zero

<sup>0</sup> 2 , *k k* φφ

for *k N* = 1, 2,..., , so that the dynamical equations (Eq. (27)) can be rewritten as follows:

( )

*ds s ii*

2 sin 2 sin (b) *<sup>n</sup> n n*

where *n N* = 1, 2,..., in the second equation. Expressing now, by means of Eq. (26), the

+ + =−+

0 0

 φπ

*i* in terms of the forcing terms Ψ*ex* and *<sup>B</sup>*

 π

sin , (a)

0 0 0

*n n k k*

+ =

*<sup>d</sup> <sup>i</sup>*

*d*

0 <sup>2</sup> *<sup>J</sup> <sup>π</sup>RI <sup>τ</sup>* <sup>=</sup> *<sup>t</sup>*

gauge-invariant superconducting phase difference *<sup>k</sup>*

0

initially trapped flux in the array), we can write:

0

φ

τ

π

*d*

τ

flux variables *ns* , we may finally rewrite Eqs. (31a-b) as follows:

*d*

*<sup>d</sup> <sup>i</sup>*

+ =

φ

*N B B k J k <sup>I</sup> i i I* <sup>=</sup>

φ

τ

= = . In addition, the instantaneous voltage *v*( )

( )

τ

Define now the partial sum *ns* ( 1 ≤ ≤ *n N* ) of the normalized fluxes as follows:

*s*

*v*

*<sup>i</sup> <sup>i</sup>*

 β

*k*

*m*

β

, <sup>2</sup>

*<sup>I</sup>* <sup>=</sup> is the normalized current flowing in the *k*-th branch

<sup>Φ</sup> . Eqs. (25-27) can be used to define the dynamics of the

in terms of the forcing parameters

of the system can be

τ

= Ψ (29)

= + *s* (30)

*i* and of the partial sums of the

φ

(26)

(27)

(28)

(31)

*B*

*<sup>I</sup>* <sup>=</sup> . The dynamical equation for each Josephson junction in the array is written

The above analysis has been carried out essentially to write the dynamical equations in terms of the effective superconducting phase difference 0 φand of the *N* partial sums *ns* .

We shall now develop a reduction of these variables to one, by assuming small values of the parameter β. Therefore, start by considering the dynamical equations of the system as written in Eqs. (32a-b). For small values of the parameter βwe can set:

$$\mathbf{s}\_{\boldsymbol{n}} = \mathbf{s}\_{\boldsymbol{n}}^{(0)} + \beta \mathbf{s}\_{\boldsymbol{n}}^{(1)}.\tag{33}$$

By substituting the above expression in Eq. (32b) and, after having multiplied both members by β, by setting the coefficients of *<sup>m</sup>* β ( *m* = 0,1 ) equal to zero, we obtain, for *n N* = 1, 2,..., :

$$\mathbf{s}\_{n}^{(0)} = n\Psi\_{\rm ex}.\tag{34}$$

For the first order corrections, on the other hand, we need to solve the following set of equations:

$$\begin{cases} y\_1 - y\_0 = s\_2^{(1)} - 3s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ y\_2 - y\_0 = s\_3^{(1)} - 2s\_2^{(1)} - \frac{\dot{t}\_B}{2} \\ y\_3 - y\_0 = s\_4^{(1)} - 2s\_3^{(1)} + s\_2^{(1)} - s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ \dots \\ y\_{N-1} - y\_0 = s\_N^{(1)} - 2s\_{N-1}^{(1)} + s\_{N-2}^{(1)} - s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ y\_N - y\_0 = s\_{N-1}^{(1)} - s\_N^{(1)} - s\_1^{(1)} \end{cases} \tag{35}$$

where ( ) <sup>0</sup> sin 2 *<sup>n</sup> ex y n* = +Ψ φ π , *n N* = 0, 1, 2,..., . By solving for (1) <sup>1</sup>*s* , which is the quantity required in Eq. (32a), we have:

$$s\_1^{(1)} = -\frac{N-1}{N+1} \frac{i\_B}{2} - \frac{1}{N+1} \sum\_{k=1}^{N} (y\_k - y\_0)\_{\prime} \tag{36}$$

Substitution of the above results into Eq. (32a) gives:

$$\frac{d\phi\_0}{d\tau} + \frac{1}{\left(N+1\right)}\sum\_{n=0}^{N} y\_n = \frac{\dot{I}\_B}{\left(N+1\right)}\tag{37}$$

Effective Models of Superconducting Quantum Interference Devices 239

In Figs. 8a-b *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curves are shown. In particular, in Fig. 8a the number of JJ's in the system ( ) *N* + 1 is taken to be equal to 10, while in Fig. 8b it is set equal to 15. We notice that, when the normalized flux approaches the first zero in the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curve for ( ) *N* + = 1 10 , namely 1 / 10 Ψ = *ex* , the system goes toward a purely resistive behaviour. In Fig. 8b, the same resistive behaviour is expected for 1 / 15 Ψ = *ex* , given that the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs.

**Figure 8.** Critical current *C N* <sup>1</sup> *i B* <sup>+</sup> <sup>=</sup> as a function of the applied magnetic flux for an array of 10 (a) and

(b)

1.0 0.5 0.0 0.5 1.0

ex

(a)

ex

1.0 0.5 0.0 0.5 1.0

In Figs. 9a-b *I-V* characteristics for different externally applied flux values are shown. In Fig. 9a the number of JJ's is taken to be equal to 10, while in Fig. 9b it is set equal to 15. Starting from Fig. 9a, we notice that, as the normalized flux approaches the first zero in the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curve for ( ) *N* + = 1 10 , namely 1 / 10 Ψ = *ex* , the behaviour of the system indeed becomes purely resistive. In Fig. 8b, the same resistive behaviour is obtained at 1 / 15 Ψ = *ex* , as

Ψ*ex* curve for ( ) *N* + = 1 15 reaches its first zero exactly at 1 / 15 Ψ = *ex* .

*B*15

*B*<sup>10</sup>

where only the positive branch has been chosen.

15 (b) Josephson junctions.

predicted.

The above differential equation represents the effective model for the 1D-JJA represented in Fig. 7, containing *N* + 1 identical over-damped junctions connected in parallel. We notice that the reduced model in Eq. (37), even being of first order in β, does not explicitly contain a 1 β term, given that Eq. (32a) contains <sup>1</sup> β<sup>−</sup> terms.

We can now explicitly perform the sum in Eq. (37), so that we write:

$$\frac{d\phi\_0}{d\tau} + \frac{B\_{N+1}}{\left(N+1\right)}\sin\left(\phi\_0 + N\pi\Psi\_{ex}\right) = \frac{\dot{I}\_B}{\left(N+1\right)}\tag{38}$$

where the absolute value of ( ) ( ) ( ) <sup>1</sup> sin 1 sin *ex N ex N B* π π <sup>+</sup> + Ψ <sup>=</sup> <sup>Ψ</sup> is the normalized critical current of the

device in this approximation, as already known from the literature (Likharev, 1986). We start by finding the voltage , as defined in Eq. (28), under the assumption that Eqs. (32b-c) can be replaced by the effective model given by first-order perturbation analysis above. Therefore, we may set:

$$\upsilon\left(\pi\right) = \frac{d\phi\_0}{d\tau} = \frac{i\_B}{\left(N+1\right)} - \frac{B\_{N+1}}{\left(N+1\right)}\sin\left(\phi\_0 + N\pi\Psi\_{ex}\right). \tag{39}$$

In this way, we can find the *I-V* characteristics by simply integrating Eq. (39), recalling the well-known procedure for a single overdamped junction. Indeed, noticing that

$$
\left\langle \upsilon \right\rangle = \frac{1}{T} \int\_0^T \frac{d\phi\_0}{d\tau} d\tau = \frac{2\pi}{T}\,\tag{40}
$$

where *T* is the period for the instantaneous voltage curve and a pseudo-period for the superconducting phase difference 0 φ (i.e., ( ) () 0 0 φ τ += + *T* φ τ π2 ), we first find the expression for 0 ϕfrom Eq. (39), for *B N* <sup>1</sup> *i B* <sup>+</sup> <sup>&</sup>gt; and *k* integer, so that

$$\phi\_0(\tau) + N\pi\Psi\_{ex} = 2\tan^{-1}\left[\frac{B\_{N+1}}{i\_B} - \frac{\sqrt{i\_B^2 - B\_{N+1}^2}}{i\_B}\tan\left(\frac{1}{2\left(N+1\right)}\sqrt{i\_B^2 - B\_{N+1}^2}\tau\right)\right] + 2k\pi\left.\tag{41}$$

The pseudo-period *T* of the above solution can be found by inspection, so that:

$$T = \frac{2\pi \left(N + 1\right)}{\sqrt{i\_B^2 - B\_{N+1}^2}}.\tag{42}$$

Therefore, by Eqs. (40) and (42), the *I-V* characteristics are given by the following expression:

( )2 2 <sup>2</sup> <sup>1</sup> 1 , *B N i N vB* <sup>+</sup> =+ + (43)

where only the positive branch has been chosen.

238 Superconductors – Materials, Properties and Applications

term, given that Eq. (32a) contains <sup>1</sup>

1 β

Therefore, we may set:

expression for 0

φτ

0

φ

τ

We can now explicitly perform the sum in Eq. (37), so that we write:

0 1

that the reduced model in Eq. (37), even being of first order in

φ

τ

*N*

φ

τ

*B*

where the absolute value of ( ) ( )

τ

superconducting phase difference 0

ϕ

 π ( ) ( )

*n n d i <sup>y</sup> <sup>d</sup> N N* <sup>=</sup> + =

+ +

*B*

β

(37)

(38)

(39)

(40)

(41)

(42)

, does not explicitly contain a

is the normalized critical current of the

*ex*

 π

+

τ

2 ), we first find the

 π

0 <sup>1</sup> , 1 1

The above differential equation represents the effective model for the 1D-JJA represented in Fig. 7, containing *N* + 1 identical over-damped junctions connected in parallel. We notice

<sup>−</sup> terms.

*d B <sup>i</sup> <sup>N</sup> d N N* <sup>+</sup> + +Ψ=

π

φ

( ) ( ) ( )

*ex*

device in this approximation, as already known from the literature (Likharev, 1986). We start by finding the voltage , as defined in Eq. (28), under the assumption that Eqs. (32b-c) can be replaced by the effective model given by first-order perturbation analysis above.

( ) ( ) ( ) ( ) 0 1

In this way, we can find the *I-V* characteristics by simply integrating Eq. (39), recalling the

0 0 1 2 , *<sup>T</sup> <sup>d</sup> v d Td T* = = φ

τ

where *T* is the period for the instantaneous voltage curve and a pseudo-period for the

τ

 (i.e., ( ) () 0 0 φ τ

1 1 1 2 2

*N iB kπ i i N*

> ( ) 2 2

2 1 . *B N*

*i B* <sup>+</sup> <sup>+</sup> <sup>=</sup> −

( )2 2 <sup>2</sup> <sup>1</sup> 1 , *B N*

Therefore, by Eqs. (40) and (42), the *I-V* characteristics are given by the following expression:

*N*

π

1

<sup>1</sup> 2tan tan 2 .

<sup>−</sup>

<sup>+</sup> = = − +Ψ + +

*B N*

*d B <sup>i</sup> v N d N N*

well-known procedure for a single overdamped junction. Indeed, noticing that

φ

( ) ( )

from Eq. (39), for *B N* <sup>1</sup> *i B* <sup>+</sup> <sup>&</sup>gt; and *k* integer, so that

0 1

*N B N*

*B B B i B*

The pseudo-period *T* of the above solution can be found by inspection, so that:

*T*

− + +

2 2

*ex B N*

+ Ψ= <sup>−</sup> − + <sup>+</sup>

π

+ +

 π

*ex*

<sup>0</sup> sin , 1 1 *N B*

*ex*

<sup>0</sup> sin . 1 1

π

 += + *T* φ τ

2 1

*i N vB* <sup>+</sup> =+ + (43)

φ

β

( ) <sup>1</sup> sin 1 sin

 <sup>+</sup> + Ψ <sup>=</sup> <sup>Ψ</sup>

*N*

*N*

In Figs. 8a-b *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curves are shown. In particular, in Fig. 8a the number of JJ's in the system ( ) *N* + 1 is taken to be equal to 10, while in Fig. 8b it is set equal to 15. We notice that, when the normalized flux approaches the first zero in the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curve for ( ) *N* + = 1 10 , namely 1 / 10 Ψ = *ex* , the system goes toward a purely resistive behaviour. In Fig. 8b, the same resistive behaviour is expected for 1 / 15 Ψ = *ex* , given that the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curve for ( ) *N* + = 1 15 reaches its first zero exactly at 1 / 15 Ψ = *ex* .

**Figure 8.** Critical current *C N* <sup>1</sup> *i B* <sup>+</sup> <sup>=</sup> as a function of the applied magnetic flux for an array of 10 (a) and 15 (b) Josephson junctions.

In Figs. 9a-b *I-V* characteristics for different externally applied flux values are shown. In Fig. 9a the number of JJ's is taken to be equal to 10, while in Fig. 9b it is set equal to 15. Starting from Fig. 9a, we notice that, as the normalized flux approaches the first zero in the *<sup>N</sup>* <sup>1</sup> *B* <sup>+</sup> vs. Ψ*ex* curve for ( ) *N* + = 1 10 , namely 1 / 10 Ψ = *ex* , the behaviour of the system indeed becomes purely resistive. In Fig. 8b, the same resistive behaviour is obtained at 1 / 15 Ψ = *ex* , as predicted.

**Figure 9.** (a) *I-V* characteristics of an array of 10 Josephson junctions for the following values of the applied magnetic flux Ψ*ex* : 0.05 (blue); 0.075 (cyan); 0.10 (orange); 0.15 (red). (b) *I-V* characteristics of an array of 15 Josephson junctions for the following values of the applied magnetic flux Ψ*ex* : 0.025 (blue); 0.05 (cyan); 0.075 (orange); 0.10 (red).

For fixed bias current values, the voltage versus flux curves can be obtained by Eq. (43) and is given by the following:

$$
\langle v \rangle = \frac{\sqrt{\dot{\imath}\_B}^2 - \left[ B\_{N+1} \left( \Psi\_{ex} \right) \right]^2}{\left( N+1 \right)}.\tag{44}
$$

Effective Models of Superconducting Quantum Interference Devices 241

*i* = . On the

*i* = . Notice also that, for

*i* . In (a) 4.0, 10.0, 16.0 *<sup>B</sup>*

*i* = (from

zero-voltage states occurs, for increasing normalized flux values, exactly at 10.0 *<sup>B</sup>*

**Figure 10.** Average voltage vs. applied magnetic flux for an array of 10 (a) and 15 (b) Josephson

In conclusion, by considering the dynamical equations of one-dimensional arrays containing *N*+1 identical overdamped Josephson junctions, the system of *N*+1 nonlinear first-order ordinary differential equation equations can be broken into two coupled subsystems, one consisting of only one equation for the superconducting phase of one junction in the array (arbitrarily chosen to be the first), the second describing the time evolution of *N* opportunely

(b)

0.0 0.5 1.0 1.5 2.0

ex

(a)

ex

0.0 0.5 1.0 1.5 2.0

*i* = (from bottom to top).

junctions for three different values of the normalized bias current *<sup>B</sup>*

0.0 0.2 0.4 0.6 0.8 1.0

v

0.0

0.5

1.0

v

1.5

bottom to top). In (b) 4.0, 10.0, 15.0 *<sup>B</sup>*

defined normalized flux variables.

relatively high bias currents, the Ψ*ex* intervals in which finite-voltage states are present tend to be flatter than in the case of a low number of Josephson junctions. This feature is more evident in Fig. 10b, where the region of zero-voltage states is narrower than in Fig. 10a. Therefore, for high enough values of the the number of Josephson junction in the array and of the bias current, away from integer values of the normalized applied flux, the interval in which finite-voltage states spread can be approximated by a horizontal segment.

other hand, in Fig. 10b, the same behaviour is detected at 15.0 *<sup>B</sup>*

The above expression is similar to the homologous d. c. SQUID *v* vs. Ψ*ex* curves.

In Figs. 10a and 10b we report the *v* vs. Ψ*ex* curve for ( ) *N* + = 1 10 and for ( ) *N* + = 1 15 , respectively. The normalized bias current values are 4.0, 10.0, 16.0 *<sup>B</sup> i* = (from bottom to top) for Fig. 10a and 4.0, 10.0, 15.0 *<sup>B</sup> i* = (from bottom to top) for Fig. 10b. Notice that the zero-voltage state regions on the Ψ*ex* -axis disappear at values of the bias current greater or equal to the critical current of the system. In Fig. 10a, indeed, we see that disappearance of zero-voltage states occurs, for increasing normalized flux values, exactly at 10.0 *<sup>B</sup> i* = . On the other hand, in Fig. 10b, the same behaviour is detected at 15.0 *<sup>B</sup> i* = . Notice also that, for relatively high bias currents, the Ψ*ex* intervals in which finite-voltage states are present tend to be flatter than in the case of a low number of Josephson junctions. This feature is more evident in Fig. 10b, where the region of zero-voltage states is narrower than in Fig. 10a. Therefore, for high enough values of the the number of Josephson junction in the array and of the bias current, away from integer values of the normalized applied flux, the interval in which finite-voltage states spread can be approximated by a horizontal segment.

240 Superconductors – Materials, Properties and Applications

*iB*

5

10

*iB*

15

0.05 (cyan); 0.075 (orange); 0.10 (red).

top) for Fig. 10a and 4.0, 10.0, 15.0 *<sup>B</sup>*

is given by the following:

**Figure 9.** (a) *I-V* characteristics of an array of 10 Josephson junctions for the following values of the applied magnetic flux Ψ*ex* : 0.05 (blue); 0.075 (cyan); 0.10 (orange); 0.15 (red). (b) *I-V* characteristics of an array of 15 Josephson junctions for the following values of the applied magnetic flux Ψ*ex* : 0.025 (blue);

(b)

0.0 0.2 0.4 0.6 0.8 1.0

v

(a)

v

0.0 0.2 0.4 0.6 0.8 1.0

For fixed bias current values, the voltage versus flux curves can be obtained by Eq. (43) and

( )

*B N ex i B*

*N*

In Figs. 10a and 10b we report the *v* vs. Ψ*ex* curve for ( ) *N* + = 1 10 and for ( ) *N* + = 1 15 ,

zero-voltage state regions on the Ψ*ex* -axis disappear at values of the bias current greater or equal to the critical current of the system. In Fig. 10a, indeed, we see that disappearance of

The above expression is similar to the homologous d. c. SQUID *v* vs. Ψ*ex* curves.

*v*

respectively. The normalized bias current values are 4.0, 10.0, 16.0 *<sup>B</sup>*

2 2

( )

<sup>1</sup> . <sup>1</sup>

*i* = (from bottom to top) for Fig. 10b. Notice that the

<sup>+</sup> − Ψ <sup>=</sup> <sup>+</sup> (44)

*i* = (from bottom to

**Figure 10.** Average voltage vs. applied magnetic flux for an array of 10 (a) and 15 (b) Josephson junctions for three different values of the normalized bias current *<sup>B</sup> i* . In (a) 4.0, 10.0, 16.0 *<sup>B</sup> i* = (from bottom to top). In (b) 4.0, 10.0, 15.0 *<sup>B</sup> i* = (from bottom to top).

In conclusion, by considering the dynamical equations of one-dimensional arrays containing *N*+1 identical overdamped Josephson junctions, the system of *N*+1 nonlinear first-order ordinary differential equation equations can be broken into two coupled subsystems, one consisting of only one equation for the superconducting phase of one junction in the array (arbitrarily chosen to be the first), the second describing the time evolution of *N* opportunely defined normalized flux variables.

When a solution of the latter *N* equations is found, by means of a perturbative approach to first order in the parameter β, the dynamical properties of the system are described by a single time-evolution equation. In this way, we may affirm that, for small values of β, the system may be described by an equivalent single-junction model, where the maximum Josephson current is appropriately defined in terms of the normalized applied flux Ψ*ex* . When we compare our present analysis to equivalent studies carried out for a two-junction interferometer, we realize that the degree of approximation of the present model in β is one order less than the first-order perturbation analysis carried out for the simplest two-junction system. This is a consequence of the approach followed in the present work, where we had to appropriately define partial sums of flux variables in order to separate the dynamical equations into two subsystems. When we refer to the SQUID case, then, we might state that the present analysis corresponds exactly to the β = 0 limit. Further work is therefore needed to carry out more detailed information on the system behaviour for finite values of β.

Effective Models of Superconducting Quantum Interference Devices 243

φ*M-2* φ*M-1*  φ*M*

φ*M-3*

**Figure 11.** Schematic representation of a Josephson junction array with alternate parameters. In the branches of the parallel connections overdamped 0-junctions (green circles) and π-junctions (red

*I0 I1 I2 I3 IM-3 IM-2 IM-1 IM* 

*IB/2 IB/2* 

*L02 L12 L02 L02 L12 L02*

Therefore, π-SQUIDs can be viewed as the building block of discretized models of multifacets Josephson junctions (MJJs) (Scharinger et al., 2010), in which the critical current density alternates between two opposite values along the junction length. However, even though some characteristic features of MJJs can be qualitatively reproduced by *N*×(0-π) onedimensional arrays, one should bear in mind that the latter are, in general, less complex

*IB*

For conventional arrays of overdamped Josephson junctions we have already shown that,

flux variable can be found by perturbation analysis. In this way, the multi-junction interferometer model reduces to a single non-linear ordinary differential equation. The same perturbation approach will be proposed again in the present section to derive the equivalent single-junction model of *N*×(0-π) one-dimensional arrays of overdamped Josephson junctions. Therefore, we start by considering the model system represented in fig. 11, consisting of identical overdamped junctions connected in parallel. In this system one half of the bias current *BI* is evenly applied to the two external branches of the array. This condition is obtained, for example, by injecting the current by means of a superconducting bar of width *w* equal to the length of the array. In order to have well focused bias, the penetration length of the superconducting bar would be much smaller than its width *w*.

Consider, as a first approach to the problem, the loop areas *<sup>k</sup> S* , the junctions resistive parameters *Rk* and maximum Josephson currents *Jk I* , and the inductances *<sup>k</sup> L* in the horizontal upper branches to be all equal. In this way one may write *<sup>k</sup>* <sup>0</sup> *S S* = , *R R <sup>k</sup>* <sup>0</sup> = , *Jk J*<sup>0</sup> *I I* = , *<sup>k</sup>* <sup>0</sup> *L L* = , for all allowed values of *k*. Define, as in the previous section in this case,

β

a series solution for the magnetic

diamonds) are alternately present.

φ*1* φ*2* φ*3*

**4.1. The homogeneous case** 

for small enough values of the characteristic parameter

systems than MJJs.

φ*0*

Even though part of the present analysis reproduces known results, as, for example, the expression for the maximum Josephson current, it still represents a simple way of approaching the problem of the electrodynamic response of one-dimensional arrays of overdamped Josephson junctions by an equivalent single-junction model. In fact, by starting with the simple representation of the dynamics of the system given in Eq. (38), all the results obtained for a single Josephson junction can be reproduced for an array of *N*+1 equal overdamped Josephson junctions. In addition, in case the solutions of the normalized flux variable would be extended to second order in the parameter β, following the same analysis as in the present section, effects due to finite β values in the electrodynamic properties of the system would be detected. Finally, considering that the present analysis has been carried out in the absence of flux fluctuations, its extension to noise effects can be obtained by means of already known results obtained for a single overdamped Josephson junction (Ambegaokar & Halperin, 1969; Bishop & Trullinger, 1978). In this case, however, care must be taken in considering the stochastic terms on all branches of the array.

## **4. Parallel connections of** *N* **× (0-***π***) overdamped Josephson junctions**

In the previous section we have considered an array of *N*+1 overdamped 0-junctions, without considering the possibility of inserting π-junctions in the system. We briefly recall that π-junctions (Bulaevskii et al., 1977; Geshkenbein et al., 1987; Baselman et al., 1999; Ryazanov et al., 2001 M. Weides et al., 2006), when compared to 0-junctions, possess an intrinsic phase difference exactly equal to π. By inserting a 0-junction and a π-junction in the same superconducting loop, π-SQUIDs may be realized. These non-conventional SQUIDs can be fabricated either by exploiting the symmetry properties of d-wave superconductors (Chesca, 1999; Schultz et al., 2000) or by utilizing both s-wave and dwave superconductors (Wollman et al., 1993; Smilde et al., 2004). A π-SQUID can thus be viewed as an elementary cell of a *N*×(0-π) one-dimensional array of overdamped Josephson junctions shown in fig. 11.

**Figure 11.** Schematic representation of a Josephson junction array with alternate parameters. In the branches of the parallel connections overdamped 0-junctions (green circles) and π-junctions (red diamonds) are alternately present.

Therefore, π-SQUIDs can be viewed as the building block of discretized models of multifacets Josephson junctions (MJJs) (Scharinger et al., 2010), in which the critical current density alternates between two opposite values along the junction length. However, even though some characteristic features of MJJs can be qualitatively reproduced by *N*×(0-π) onedimensional arrays, one should bear in mind that the latter are, in general, less complex systems than MJJs.

For conventional arrays of overdamped Josephson junctions we have already shown that, for small enough values of the characteristic parameter β a series solution for the magnetic flux variable can be found by perturbation analysis. In this way, the multi-junction interferometer model reduces to a single non-linear ordinary differential equation. The same perturbation approach will be proposed again in the present section to derive the equivalent single-junction model of *N*×(0-π) one-dimensional arrays of overdamped Josephson junctions. Therefore, we start by considering the model system represented in fig. 11, consisting of identical overdamped junctions connected in parallel. In this system one half of the bias current *BI* is evenly applied to the two external branches of the array. This condition is obtained, for example, by injecting the current by means of a superconducting bar of width *w* equal to the length of the array. In order to have well focused bias, the penetration length of the superconducting bar would be much smaller than its width *w*.

### **4.1. The homogeneous case**

242 Superconductors – Materials, Properties and Applications

the present analysis corresponds exactly to the

as in the present section, effects due to finite

Josephson junctions shown in fig. 11.

β

variable would be extended to second order in the parameter

considering the stochastic terms on all branches of the array.

first order in the parameter

When a solution of the latter *N* equations is found, by means of a perturbative approach to

system may be described by an equivalent single-junction model, where the maximum Josephson current is appropriately defined in terms of the normalized applied flux Ψ*ex* . When we compare our present analysis to equivalent studies carried out for a two-junction

order less than the first-order perturbation analysis carried out for the simplest two-junction system. This is a consequence of the approach followed in the present work, where we had to appropriately define partial sums of flux variables in order to separate the dynamical equations into two subsystems. When we refer to the SQUID case, then, we might state that

β

Even though part of the present analysis reproduces known results, as, for example, the expression for the maximum Josephson current, it still represents a simple way of approaching the problem of the electrodynamic response of one-dimensional arrays of overdamped Josephson junctions by an equivalent single-junction model. In fact, by starting with the simple representation of the dynamics of the system given in Eq. (38), all the results obtained for a single Josephson junction can be reproduced for an array of *N*+1 equal overdamped Josephson junctions. In addition, in case the solutions of the normalized flux

β

**4. Parallel connections of** *N* **× (0-***π***) overdamped Josephson junctions** 

In the previous section we have considered an array of *N*+1 overdamped 0-junctions, without considering the possibility of inserting π-junctions in the system. We briefly recall that π-junctions (Bulaevskii et al., 1977; Geshkenbein et al., 1987; Baselman et al., 1999; Ryazanov et al., 2001 M. Weides et al., 2006), when compared to 0-junctions, possess an intrinsic phase difference exactly equal to π. By inserting a 0-junction and a π-junction in the same superconducting loop, π-SQUIDs may be realized. These non-conventional SQUIDs can be fabricated either by exploiting the symmetry properties of d-wave superconductors (Chesca, 1999; Schultz et al., 2000) or by utilizing both s-wave and dwave superconductors (Wollman et al., 1993; Smilde et al., 2004). A π-SQUID can thus be viewed as an elementary cell of a *N*×(0-π) one-dimensional array of overdamped

system would be detected. Finally, considering that the present analysis has been carried out in the absence of flux fluctuations, its extension to noise effects can be obtained by means of already known results obtained for a single overdamped Josephson junction (Ambegaokar & Halperin, 1969; Bishop & Trullinger, 1978). In this case, however, care must be taken in

single time-evolution equation. In this way, we may affirm that, for small values of

interferometer, we realize that the degree of approximation of the present model in

to carry out more detailed information on the system behaviour for finite values of

, the dynamical properties of the system are described by a

= 0 limit. Further work is therefore needed

β

values in the electrodynamic properties of the

β, the

βis one

β.

, following the same analysis

Consider, as a first approach to the problem, the loop areas *<sup>k</sup> S* , the junctions resistive parameters *Rk* and maximum Josephson currents *Jk I* , and the inductances *<sup>k</sup> L* in the horizontal upper branches to be all equal. In this way one may write *<sup>k</sup>* <sup>0</sup> *S S* = , *R R <sup>k</sup>* <sup>0</sup> = , *Jk J*<sup>0</sup> *I I* = , *<sup>k</sup>* <sup>0</sup> *L L* = , for all allowed values of *k*. Define, as in the previous section in this case,

0 0 0 *<sup>J</sup> L I* β <sup>=</sup> <sup>Φ</sup> . Consider, finally, as also represented in fig. 11, that 0-junctions (green circles) occupy even positions ( *k M* = − 0, 2, 4,.., 1 ) and π-junctions (diamonds) occupy odd positions ( *k M* = 1, 3,.., ). The reader will notice very many similarities with the analysis in the previous section. However, due to the higher degree of complexity of the problem, we need to proceed step by step.

By fluxoid quantization, the normalized magnetic flux 0 *k k* Φ Ψ = Φ linked to the *k*-th cell of the array is related to the superconducting phase differences *<sup>k</sup>*−<sup>1</sup> φ and *<sup>k</sup>* φ across the two junctions in the cell as follows:

$$
\mu\_{k-1} - \phi\_k + \mathfrak{A}\mathfrak{A}\mathbb{1}\_k = \mathfrak{A}\mathfrak{A} \left[ n\_k + \frac{\left(-1\right)^k - 1}{4} \right],\tag{45}
$$

Effective Models of Superconducting Quantum Interference Devices 245

= Ψ (48)

(49)

(50)

*i* and of the

(51)

(52)

1 .

( )

*k*

1 1 2 , <sup>2</sup>

φ

φ

<sup>+</sup> <sup>=</sup> in (50c). Expressing now, by means of

 β

> β

=

By fluxoid quantization, setting all *<sup>k</sup> n* 's to zero in Eq. (45) under the hypothesis of zero

− − =+ +

( ) ( )

2 sin 2 sin , 2 (b)

sin , (a)

( ) *<sup>n</sup>* <sup>−</sup> 1 (c)

*M*

sin , (a) <sup>2</sup> 2 2 sin 2 sin , (b) <sup>2</sup>

+ −

2

β

−

− + Ψ−

1 1

− Ψ −Ψ

β

2 21 22 1

− + Ψ− − =−

− −

sin , (c) <sup>2</sup>

ββ

we may assume that the solution, to first order in this perturbation

*ss s s i*

*n n n ex B*

*i* , 2*<sup>n</sup> i* , and 2 1 *<sup>n</sup> i* <sup>−</sup> in terms of the forcing terms Ψ*ex* and *<sup>B</sup>*

0 21 21 0 0

partial sums of the flux variables *<sup>k</sup> s* defined in Eq. (48), we may finally rewrite Eqs. (50a-c)

− + = −+ =

0 2 20 0

+ + = −+ =

*n n*

*ds s i i kn*

2 sin 2 sin , 2

*ds s ii k*

*n n*

2 21 2 21 1

 φ

*ds s s s*

− + − + =−

 φ

*ds s ss s i*

*n n n n ex B*

+ +−− = −

0

φ

*M M M ex ex*

2 sin 2 sin , (d)

β

*<sup>n</sup>* <sup>−</sup> <sup>=</sup> in Eqs. (51b) and (51c) and where it has been set 0*<sup>s</sup>* <sup>=</sup> <sup>0</sup> .

Considering the dynamical equations of the system as written in Eqs. (51a-d), for small

(0) (1). *nn n ss s* ≈ + β

− −

 π

*n n k k*

*s*

0

φφ

( )

*n*

*k k s*

for *k M* = 1, 2,..., , so that the dynamical equations can be rewritten as follows:

 π

initially trapped flux in the array, we can write:

0 0

 φπ

> φπ

*<sup>n</sup>* <sup>−</sup> <sup>=</sup> in (50b), and <sup>1</sup> 1, 2,..., <sup>2</sup>

( )

*B ex*

( )

−+ −

( )

0 21

*n*

−

0 0

*s*

0 0

*M*

*s*

*n*

− Ψ

*s*

 β

0

φ

τ

π

π

where <sup>1</sup> 1, 2,..., <sup>2</sup>

Eq. (46), the currents 0

*d*

φ

τ

π

π

π

as follows:

*d*

2

τ

*d*

*n*

*<sup>d</sup> <sup>i</sup>*

+ =

φ

2 1

*n*

τ

*d* −

0 1 0

> φπ

 φπ

*M*

β

parameter, can be written as follows:

 φπ

+ =+

*d s i*

φ

2 sin 2

2 1

−

*n*

τ

*d*

*d*

where now <sup>1</sup> 1, 2,..., <sup>2</sup>

values of the parameter

τ

*d ds*

τ *M*

where *<sup>k</sup> n* is an integer and *k M* = 1, 2,..., (notice that the number of loops is one less the number of junctions). If an external magnetic field **H** , orthogonal to the plane of the array, is applied to the system, the normalized geometric applied flux through each cell now is 0 0 0 0 *ex ex* Φ μ *HS* Ψ= = Φ Φ . By considering Φ*k* as the sum of the applied and induced flux, we may

write:

$$\Psi\_k = \beta \left(\sum\_{m=0}^{k-1} i\_m - \frac{i\_B}{2}\right) + \Psi\_{e\alpha'} \tag{46}$$

for *k M* = 1, 2,..., , where 0 *k k J <sup>I</sup> <sup>i</sup> <sup>I</sup>* <sup>=</sup> is the normalized current flowing in the *k*-th branch and

0 *B B J <sup>I</sup> <sup>i</sup> <sup>I</sup>* <sup>=</sup> . The dynamical equation for each Josephson junction in the array can be written by

means of the RSJ model (Barone & Paternò, 1982), so that:

$$\frac{d\phi\_k}{d\tau} + \sin\phi\_k = \dot{\mathbf{r}}\_{k'} \tag{47}$$

where *k M* = 0, 1, 2,..., and 0 0 0 <sup>2</sup> *R IJ <sup>t</sup>* π τ <sup>=</sup> <sup>Φ</sup> . Eqs. (45-47) can be used to define the dynamics of the gauge-invariant superconducting phase difference *<sup>k</sup>* φ in terms of the forcing parameters Ψ*ex* and 0 *M B B k J k <sup>I</sup> i i I* <sup>=</sup> = = . Define now the partial sum *ns* ( <sup>1</sup> ≤ ≤ *n M* ) of the normalized fluxes as follows:

Effective Models of Superconducting Quantum Interference Devices 245

$$s\_n = \sum\_{k=1}^n \Psi\_k.\tag{48}$$

By fluxoid quantization, setting all *<sup>k</sup> n* 's to zero in Eq. (45) under the hypothesis of zero initially trapped flux in the array, we can write:

$$\phi\_k = \phi\_0 + 2\pi \mathbf{s}\_k + \pi \left[ \frac{1 - \left(-1\right)^k}{2} \right],\tag{49}$$

for *k M* = 1, 2,..., , so that the dynamical equations can be rewritten as follows:

244 Superconductors – Materials, Properties and Applications

By fluxoid quantization, the normalized magnetic flux

array is related to the superconducting phase differences *<sup>k</sup>*−<sup>1</sup>

1

0 *k*

*k J <sup>I</sup> <sup>i</sup>*

means of the RSJ model (Barone & Paternò, 1982), so that:

τ

the gauge-invariant superconducting phase difference *<sup>k</sup>*

where *k M* = 0, 1, 2,..., and 0 0

0

*M*

*B B k J k <sup>I</sup> i i I* <sup>=</sup>

 φπ

φ

*kk k k n* <sup>−</sup>

<sup>=</sup> <sup>Φ</sup> . Consider, finally, as also represented in fig. 11, that 0-junctions (green circles)

occupy even positions ( *k M* = − 0, 2, 4,.., 1 ) and π-junctions (diamonds) occupy odd positions ( *k M* = 1, 3,.., ). The reader will notice very many similarities with the analysis in the previous section. However, due to the higher degree of complexity of the problem, we

> 0 *k*

φ

(46)

*<sup>I</sup>* <sup>=</sup> is the normalized current flowing in the *k*-th branch and

<sup>=</sup> <sup>Φ</sup> . Eqs. (45-47) can be used to define the dynamics of

in terms of the forcing parameters

φ

= = . Define now the partial sum *ns* ( <sup>1</sup> ≤ ≤ *n M* ) of the normalized fluxes

 and *<sup>k</sup>* φ

linked to the *k*-th cell of the

across the two

(45)

(47)

*k* Φ Ψ = Φ

( )

22 , <sup>4</sup>

Φ Φ . By considering Φ*k* as the sum of the applied and induced flux, we may

*B k m ex*

*<sup>i</sup> <sup>i</sup>* −

*<sup>I</sup>* <sup>=</sup> . The dynamical equation for each Josephson junction in the array can be written by

sin , *<sup>k</sup> k k*

φ

+ =

*<sup>d</sup> <sup>i</sup>*

*d*

0 <sup>2</sup> *R IJ <sup>t</sup>* π

φ

τ

 Ψ = − +Ψ 

, <sup>2</sup>

− − − + Ψ= +

 π

where *<sup>k</sup> n* is an integer and *k M* = 1, 2,..., (notice that the number of loops is one less the number of junctions). If an external magnetic field **H** , orthogonal to the plane of the array, is applied to the system, the normalized geometric applied flux through each cell now is

1

*k*

β

*m*

0

=

1 1

*k*

0 0 0 *<sup>J</sup> L I*

need to proceed step by step.

junctions in the cell as follows:

0 0

0 0

for *k M* = 1, 2,..., , where

μ*HS* Ψ= =

*ex*

Φ

0 *B*

*J <sup>I</sup> <sup>i</sup>*

Ψ*ex* and

as follows:

*B*

*ex*

write:

β

$$\frac{d\phi\_0}{d\tau} + \sin\phi\_0 = i\_{0\prime} \tag{a}$$

$$2\pi \frac{d s\_{2n}}{d \pi} + \sin \left(\phi\_0 + 2\pi s\_{2n}\right) = i\_{2n} - i\_0 + \sin \phi\_{0'} \tag{50} \tag{50}$$

$$2\pi \frac{d s\_{2n-1}}{d \pi} - \sin \left(\phi\_0 + 2\pi s\_{2n-1}\right) = i\_{2n-1} - i\_0 + \sin \phi\_0 \quad \left(k = 2n - 1\right) \quad \text{(c)}$$

where <sup>1</sup> 1, 2,..., <sup>2</sup> *M <sup>n</sup>* <sup>−</sup> <sup>=</sup> in (50b), and <sup>1</sup> 1, 2,..., <sup>2</sup> *M n* <sup>+</sup> <sup>=</sup> in (50c). Expressing now, by means of Eq. (46), the currents 0 *i* , 2*<sup>n</sup> i* , and 2 1 *<sup>n</sup> i* <sup>−</sup> in terms of the forcing terms Ψ*ex* and *<sup>B</sup> i* and of the partial sums of the flux variables *<sup>k</sup> s* defined in Eq. (48), we may finally rewrite Eqs. (50a-c) as follows:

$$\frac{d\phi\_0}{d\tau} + \sin\phi\_0 = \frac{i\_B}{2} + \frac{s\_1 - \Psi\_{ex}}{\beta},\tag{a}$$

$$2\pi \frac{ds\_{2n}}{d\tau} + \sin\left(\phi\_0 + 2\pi s\_n\right) - \sin\phi\_0 - \frac{s\_{2n+1} - 2s\_{2n} + s\_{2n-1}}{\beta} = \frac{\Psi\_{ex} - s\_1}{\beta} - \frac{i\_B}{2},\tag{b}$$

$$2\pi \frac{d\mathbf{s}\_{2n-1}}{d\tau} - \sin\left(\phi\_0 + 2\pi \mathbf{s}\_{2n-1}\right) - \sin\phi\_0 - \frac{\mathbf{s}\_{2n} - 2\mathbf{s}\_{2n-1} + \mathbf{s}\_{2n-2}}{\beta} = \frac{\Psi\_{ex} - \mathbf{s}\_1}{\beta} - \frac{i\_\mathbf{B}}{2}, \quad \text{(c)}\tag{51}$$

$$2\pi \frac{d\mathbf{s}\_M}{d\tau} - \sin\left(\phi\_0 + 2\pi \mathbf{s}\_M\right) - \sin\phi\_0 + \frac{\mathbf{s}\_M - \mathbf{s}\_{M-1}}{\beta} = \frac{\Psi\_{ex}}{\beta} - \frac{\mathbf{s}\_1 - \Psi\_{ex}}{\beta}, \tag{d)}$$

where now <sup>1</sup> 1, 2,..., <sup>2</sup> *M <sup>n</sup>* <sup>−</sup> <sup>=</sup> in Eqs. (51b) and (51c) and where it has been set 0*<sup>s</sup>* <sup>=</sup> <sup>0</sup> .

Considering the dynamical equations of the system as written in Eqs. (51a-d), for small values of the parameter β we may assume that the solution, to first order in this perturbation parameter, can be written as follows:

$$\mathbf{s}\_{n} = \mathbf{s}\_{n}^{(0)} + \beta \mathbf{s}\_{n}^{(1)}.\tag{52}$$

By substituting the above expression in Eqs. (51b-d) and, after having multiplied both members by β, by setting to zero the coefficients of *<sup>m</sup>* β( *m* = 0,1 ), we obtain:

$$\mathbf{s}\_n^{(0)} = n\Psi\_{ex}.\tag{53}$$

Effective Models of Superconducting Quantum Interference Devices 247

Consider, next, a non-homogeneous array with alternating 0-π Josephson junctions. We shall take the parameters of all 0-junctions equal. The parameters of π-junctions, even being equal among them, are assumed to be different from those of the 0-junctions. In this case we can omit some of the calculations, having already treated the problem in detail in the

Considering again the 1D-JJA represented in fig. 11, we now assume that the loop areas *<sup>k</sup> S* , the junctions resistive parameters *Rk* and maximum Josephson currents *Jk I* , and the inductances *<sup>k</sup> L* of the horizontal upper branches alternate in their values, as we go along the

even , (a) odd

even , (b) odd

, (c) odd

α, ε, and σ

<sup>=</sup> <sup>Φ</sup> ; therefore, we may define 1 1

even

odd

(58)

are implicitly

εαβ <sup>Φ</sup> .

(59)

(60)

0 *<sup>J</sup> L I*

1

β= =

even

 even , (d) odd *<sup>k</sup> L k*

> 0 *<sup>J</sup> L I*

Fluxoid quantization give the same relation as in Eq. (45) between the superconducting phases and the normalized magnetic flux Ψ*k* linked to the *k*-th cell of the array. However,

2

*B m ex*

2

for *k M* = 1, 2,..., . The equations of the motion for the superconducting phases can now be

2 2

+ =

*<sup>d</sup> <sup>i</sup>*

*n n*

*B m ex*

*<sup>i</sup> i k*

21 21

 α

*n n*

sin , (b)

− −

sin , (a)

σ

*<sup>i</sup> i k*

*LS k*

and 0 0

β

1

−

*k*

β

*m k k*

> βσε

2

τ

*d*

φ

φ

*n*

*<sup>d</sup> <sup>i</sup>*

+ =

εα φ

φ

2 1

where *n* runs over all allowed *k*-values also in all following equations.

*n*

τ

*d* − 0 1

= −

*m*

0

− +Ψ

=

Ψ = − +Ψ

*S k*

*SS k R k*

σ

*RS k I k <sup>I</sup> II k*

*k*

 = =

 = =

*S*

*k*

*R*

*Jk*

*L*

τ<sup>=</sup> <sup>Φ</sup> 0

*J*

0 1 0

<sup>=</sup> <sup>=</sup>

 = = σ

for all allowed values of *k*. In this way, the additional parameters

0 <sup>2</sup> *R IJ <sup>t</sup>* π

1 0

ε

α

*J J*

**4.2 The non-homogeneous case** 

previous subsection.

array. In this way we write

defined. As before, we set 0 0

Eq. (46) is modified as follows

written as follows:

For the first order corrections, on the other hand, we need to solve the following set of equations:

$$\begin{cases} -y\_1 - y\_0 = s\_2^{(1)} - 3s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ y\_2 - y\_0 = s\_3^{(1)} - 2s\_2^{(1)} - \frac{\dot{t}\_B}{2} \\ -y\_3 - y\_0 = s\_4^{(1)} - 2s\_3^{(1)} + s\_2^{(1)} - s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ \dots \\ -y\_{M-2} - y\_0 = s\_{M-1}^{(1)} - 2s\_{M-2}^{(1)} + s\_{M-1}^{(1)} - s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ y\_{M-1} - y\_0 = s\_M^{(1)} - 2s\_{M-1}^{(1)} + s\_{M-2}^{(1)} - s\_1^{(1)} - \frac{\dot{t}\_B}{2} \\ -y\_M - y\_0 = s\_{M-1}^{(1)} - s\_M^{(1)} - s\_1^{(1)} \end{cases} \tag{54}$$

where ( ) <sup>0</sup> sin 2 *<sup>n</sup> ex y n* = +Ψ φ π , *n M* = 0, 1, 2,..., . By solving for (1) <sup>1</sup>*s* , which is the required quantity in Eq. (51a), we have:

$$s\_1^{(1)} = -\frac{M-1}{M+1} \frac{i\_B}{2} - \frac{1}{M+1} \sum\_{k=1}^{M} \left[ \left(-1\right)^k y\_k - y\_0 \right] \tag{55}$$

Substitution of the above result into Eq. (51a) gives:

$$\frac{d\phi\_0}{d\tau} + \frac{1}{\left(M+1\right)}\sum\_{k=0}^{M} \left(-1\right)^k \sin\left(\phi\_0 + 2\pi k \Psi\_{ex}\right) = \frac{\dot{I}\_B}{\left(M+1\right)}.\tag{56}$$

It is now possible to explicitly calculate the finite sum in Eq. (56) to get, in terms of the number <sup>1</sup> 2 *<sup>M</sup> <sup>N</sup>* <sup>+</sup> <sup>=</sup> of the individual (0-π) cells:

$$\frac{d\phi\_0}{d\tau} - \frac{A\_{2N}}{2N} \cos\left[\phi\_0 + (2N - 1)\pi\Psi\_{ex}\right] = \frac{i\_B}{2N},\tag{57}$$

where ( ) ( ) <sup>2</sup> sin 2 cos *ex N ex N A* π π <sup>Ψ</sup> <sup>=</sup> <sup>Ψ</sup> . The above equation represents the equivalent single-junction

model for the homogeneous array consisting of *N* cells, each one containing one 0-junction and one π-junction.

### **4.2 The non-homogeneous case**

246 Superconductors – Materials, Properties and Applications

members by

equations:

β

where ( ) <sup>0</sup> sin 2 *<sup>n</sup> ex y n* = +Ψ φ

number <sup>1</sup>

2

where ( ) ( ) <sup>2</sup> sin 2 cos

*N*

and one π-junction.

*A*

quantity in Eq. (51a), we have:

 π

Substitution of the above result into Eq. (51a) gives:

0

φ

τ

*<sup>M</sup> <sup>N</sup>* <sup>+</sup> <sup>=</sup> of the individual (0-π) cells:

φ

τ

*ex*

*ex*

*N*

π

π

0 2

By substituting the above expression in Eqs. (51b-d) and, after having multiplied both

For the first order corrections, on the other hand, we need to solve the following set of

(1) (1) (1) (1)

2

*B*

*i*

2

*B*

*i*

2 0 1 2 11

−= − + − −

2

− −= − + − −

1 0 1 21 (1) (1) (1) 01 1

, *n M* = 0, 1, 2,..., . By solving for (1)

*<sup>s</sup> y y M M* <sup>=</sup>

*<sup>d</sup> <sup>i</sup> <sup>k</sup> d M M* <sup>=</sup> + − + Ψ=

*y ys s s s*

− − −−

2

*y ys s s s*

( ) (1) 1 0

− − −

*M MM M*

−

*M MM*

− −= − −

*M i*

0

*k*

*y ys s s*

*M M MM*

(1) (1) (1) (1)

1 1 1 1 , 12 1 *<sup>M</sup> <sup>k</sup> <sup>B</sup>*

*k*

( ) () ( ) ( )

φ

It is now possible to explicitly calculate the finite sum in Eq. (56) to get, in terms of the

*d A <sup>i</sup> <sup>N</sup> d N N* − + −Ψ= 

model for the homogeneous array consisting of *N* cells, each one containing one 0-junction

φ

0

+ +

<sup>0</sup> cos (2 1) , 2 2 *N B*

<sup>1</sup> 1 sin 2 . 1 1 *<sup>M</sup> <sup>k</sup> <sup>B</sup>*

 π 2

*B*

*i*

(1) (1) (1) (1)

(1) (1)

3

(1) (1)

2

304 3 2 1

*yys s s s*

2

−− = − + − −

102 1

−− = − −

*yys s*

203 2

−= − −

*yys s*

....

 

β

( *m* = 0,1 ), we obtain:

(0) . *n ex s n* = Ψ (53)

2

*B*

.

<sup>1</sup>*s* , which is the required

(54)

(56)

(57)

*i*

2

*k*

*ex*

*ex*

 π

<sup>Ψ</sup> <sup>=</sup> <sup>Ψ</sup> . The above equation represents the equivalent single-junction

<sup>−</sup> =− − − − + + (55)

*B*

*i*

, by setting to zero the coefficients of *<sup>m</sup>*

Consider, next, a non-homogeneous array with alternating 0-π Josephson junctions. We shall take the parameters of all 0-junctions equal. The parameters of π-junctions, even being equal among them, are assumed to be different from those of the 0-junctions. In this case we can omit some of the calculations, having already treated the problem in detail in the previous subsection.

Considering again the 1D-JJA represented in fig. 11, we now assume that the loop areas *<sup>k</sup> S* , the junctions resistive parameters *Rk* and maximum Josephson currents *Jk I* , and the inductances *<sup>k</sup> L* of the horizontal upper branches alternate in their values, as we go along the array. In this way we write

$$\begin{aligned} \mathcal{S}\_k &= \begin{bmatrix} S\_0 & k & \text{even} \\ S\_1 = \sigma \mathcal{S}\_0 & k & \text{odd} \end{bmatrix} & \text{(a)}\\ R\_k &= \begin{cases} R\_0 & k \text{ even} \\ R\_1 = \alpha \mathcal{S}\_0 & k \text{ odd} \end{cases}, & \text{(b)}\\ I\_{jk} &= \begin{cases} I\_{f0} & k \text{ even} \\ I\_{f1} = \varepsilon I\_{f0} & k \text{ odd} \end{cases}, & \text{(c)}\\ L\_k &= \begin{cases} L\_0 & k \text{ even} \\ L\_1 = \sigma \mathcal{S}\_0 & k \text{ odd} \end{cases}, & \text{(d)} \end{aligned} \end{aligned} \tag{58}$$

for all allowed values of *k*. In this way, the additional parameters α, ε, and σ are implicitly defined. As before, we set 0 0 0 <sup>2</sup> *R IJ <sup>t</sup>* π τ <sup>=</sup> <sup>Φ</sup> and 0 0 0 *<sup>J</sup> L I* β <sup>=</sup> <sup>Φ</sup> ; therefore, we may define 1 1 1 0 *<sup>J</sup> L I* β = = εαβ <sup>Φ</sup> .

Fluxoid quantization give the same relation as in Eq. (45) between the superconducting phases and the normalized magnetic flux Ψ*k* linked to the *k*-th cell of the array. However, Eq. (46) is modified as follows

$$\Psi\_k = \begin{cases} \beta \left( \sum\_{m=0}^{k-1} i\_m - \frac{\dot{i}\_B}{2} \right) + \Psi\_{ex} & k \text{ even} \\ \beta \sigma \varepsilon \left( \sum\_{m=0}^{k-1} i\_m - \frac{\dot{i}\_B}{2} \right) + \sigma \Psi\_{ex} & k \text{ odd} \end{cases} \tag{59}$$

for *k M* = 1, 2,..., . The equations of the motion for the superconducting phases can now be written as follows:

$$\begin{aligned} \frac{d\phi\_{2n}}{d\tau} + \sin\phi\_{2n} &= i\_{2n'} \\ \frac{d\phi\_{2n-1}}{d\tau} + \varepsilon\alpha\sin\phi\_{2n-1} &= \alpha i\_{2n-1'} \end{aligned} \tag{60}$$

where *n* runs over all allowed *k*-values also in all following equations.

By adopting a first-order perturbation analysis in the parameter β, we set:

$$\begin{aligned} s\_{2n} &= s\_{2n}^{(0)} + \beta s\_{2n}^{(1)} = n \left( \sigma + 1 \right) \Psi\_{ex} + \beta s\_{2n}^{(1)} \\ s\_{2n-1} &= s\_{2n-1}^{(0)} + \beta s\_{2n-1}^{(1)} = \left[ n \left( \sigma + 1 \right) - 1 \right] \Psi\_{ex} + \beta s\_{2n-1}^{(1)} \tag{61} \end{aligned} \tag{61}$$

Effective Models of Superconducting Quantum Interference Devices 249

= = 1, 1 and *N*=10. The factors in front

= 1 , giving rise to secondary peaks in the overall curves. The functions ( ) *N E ex* ⋅ Ψ

(a)

*ex*

1.0 0.5 0.0 0.5 1.0

of these functions are chosen in such a way that they can attain the same maximum value of

**Figure 12.** (a) Envelope function ( ) *N E* Ψ*ex* as a function of 2 ( ) *ex F* Ψ . (b) Critical current as a function

(b)

1.0 0.5 0.0 0.5 1.0

*ex*

One notices that regular primary peaks of ( ) *ex F* Ψ appear at integer and half integer values of Ψ*ex* . By multiplying the functions ( ) *ex E* Ψ and ( ) *ex F* Ψ as reported in fig. 12b, we notice that primary peaks positioned at half-integer Ψ*ex* values are left unchanged, while those

By proceeding in the same way for the non-homogeneous case, starting from Eq. (62) we find that the critical current of the non homogeneous array described in Section 3 can be

( )( ). *c ex ex iE F* =Ψ Ψ (67)

of the normalized applied flux for a homogeneous parallel array with *N*=12, *ε*=1, *σ*=1.

appearing at integer Ψ*ex* values are depressed to zero.

0

5

10

*ic*

15

20

written as follows:

ε

 σ

with σ

and 2 ( ) *ex F* Ψ are separately shown in fig. 12a for

0

5

10

*NE ex* , *F ex*

15

20

the product ( )( ) *ex ex E F* Ψ⋅Ψ .

By substituting the above expression in Eqs. (49) to obtain the superconducting phases *<sup>k</sup>* φ we then write the differential equations (60a-b) in terms of the sole phase variable 0 φ and of the 2*N* normalized flux variables *ks* . By following the same steps as in the previous section, we obtain the following equivalent single-junction model for the superconducting array:

$$\frac{d\tilde{\phi}\_0}{d\tau} + \frac{\alpha \tilde{A}\_N}{(\alpha + 1)N} \left[ \sin \tilde{\phi}\_0 - \varepsilon \sin \left( \tilde{\phi}\_0 + 2\pi \sigma \Psi\_{ex} \right) \right] = \frac{\alpha i\_B}{(\alpha + 1)N} \tag{62}$$

Where here ( ) ( ) ( ) ( ) sin 1 sin 1 *ex N ex N A* π σ π σ + Ψ <sup>=</sup> + Ψ and ( )( ) 0 0 1 1 *<sup>N</sup> ex* φφπ σ = + − +Ψ . Naturally, for α, ε,

and σall equal to one, the ordinary differential equation (62) reduces to Eq. (57).

### **4.3. Critical current**

In order to find the critical current of arrays with alternating 0-π Josephson junctions, we proceed as follows. First, consider the homogeneous array described in Section 4.1. We look for the maximum value of the bias current *<sup>B</sup> i* which can be injected in the system at zero voltage

( <sup>0</sup> <sup>0</sup> *<sup>d</sup> <sup>d</sup>* <sup>=</sup> φ τ in Eq. (57)). Therefore, by maximizing with respect to 0 φthe following expression

$$\dot{a}\_B = -A\_{2N} \cos \left[\phi\_0 + (2N - 1)\pi \Psi\_{ex}\right],\tag{63}$$

we can express the critical current of the homogeneous device as follows:

$$\left|i\_c = \left|A\_{2N}\right| = \left|\frac{\sin\left(2N\pi\Psi\_{ex}^{\nu}\right)}{\cos\left(\pi\Psi\_{ex}\right)}\right|.\tag{64}$$

In order to understand the origin of the patterns we are going to show for the non homogeneous case, let us consider the result in Eq. (64) as the product of the envelop function

$$E\left(\Psi\_{ex}\right) = \sqrt{1 + \varepsilon^2 - 2\varepsilon\cos\left(2\pi\sigma\Psi\_{ex}\right)},\tag{65}$$

( ε σ= = 1, 1 in this case) and of the rapidly oscillating function

$$F\left(\Psi\_{ex}\right) = \left|\frac{\sin\left(N\pi\left(\sigma + 1\right)\Psi\_{ex}\right)}{\sin\left(\pi\left(\sigma + 1\right)\Psi\_{ex}\right)}\right|,\tag{66}$$

with σ = 1 , giving rise to secondary peaks in the overall curves. The functions ( ) *N E ex* ⋅ Ψ and 2 ( ) *ex F* Ψ are separately shown in fig. 12a for ε σ = = 1, 1 and *N*=10. The factors in front of these functions are chosen in such a way that they can attain the same maximum value of the product ( )( ) *ex ex E F* Ψ⋅Ψ .

248 Superconductors – Materials, Properties and Applications

0

 α

α

sin 1 sin 1

+ Ψ <sup>=</sup> + Ψ

*N*

π σ

π σ

( ) ( )

φ

τ

Where here ( ) ( )

the maximum value of the bias current *<sup>B</sup>*

*N*

*A*

**4.3. Critical current** 

and σ

( <sup>0</sup> <sup>0</sup> *<sup>d</sup> <sup>d</sup>* <sup>=</sup> φ

τ

( ε

 σ

By adopting a first-order perturbation analysis in the parameter

( )

*nn n ex n*

*ss s n s* −− − <sup>−</sup>

≈ + = + − Ψ+

we then write the differential equations (60a-b) in terms of the sole phase variable 0

*d A i d N N* <sup>+</sup> − +Ψ= + +

and ( )( ) 0 0 1 1 *<sup>N</sup> ex*

all equal to one, the ordinary differential equation (62) reduces to Eq. (57).

In order to find the critical current of arrays with alternating 0-π Josephson junctions, we proceed as follows. First, consider the homogeneous array described in Section 4.1. We look for

> 2 0 cos (2 1) , *B N ex iA N* =− + − Ψ φ

> > ( ) <sup>2</sup> sin 2 . cos *ex*

In order to understand the origin of the patterns we are going to show for the non homogeneous case, let us consider the result in Eq. (64) as the product of the envelop function

> ( ) ( ) <sup>2</sup> 1 2 cos 2 , *ex ex <sup>E</sup>* Ψ=+− Ψ ε

> > ( ) ( ) ( )

π σ

π σ

*N*

( ) ( ) sin 1 , sin 1

+ Ψ

+ Ψ

 ε

(0) (1) (1) 22 2 2

βσ

*nn n ex n*

*ss sn s*

βσ

*ex*

*ex*

in Eq. (57)). Therefore, by maximizing with respect to 0

we can express the critical current of the homogeneous device as follows:

*c N*

*i A*

= = 1, 1 in this case) and of the rapidly oscillating function

*F*

*ex*

Ψ =

φ ε

≈ + = +Ψ+

( )

 β

(0) (1) (1) 21 21 21 2 1

By substituting the above expression in Eqs. (49) to obtain the superconducting phases *<sup>k</sup>*

the 2*N* normalized flux variables *ks* . By following the same steps as in the previous section, we obtain the following equivalent single-junction model for the superconducting array:

( ) ( ) ( )

 φ

φφπ

0 0 sin sin 2 , 1 1 *N B*

 πσ β, we set:

> α

*i* which can be injected in the system at zero voltage

the following expression

φ

<sup>Ψ</sup> = = Ψ (64)

 π

( )

π

*N*

π

*ex*

 πσ

*ex*

*ex*

 α

 σ= + − +Ψ . Naturally, for

(61)

φ

φ

(62)

α, ε,

(63)

(65)

(66)

and of

1 , (a)

*ex*

1 1 , (b)

 β

**Figure 12.** (a) Envelope function ( ) *N E* Ψ*ex* as a function of 2 ( ) *ex F* Ψ . (b) Critical current as a function of the normalized applied flux for a homogeneous parallel array with *N*=12, *ε*=1, *σ*=1.

One notices that regular primary peaks of ( ) *ex F* Ψ appear at integer and half integer values of Ψ*ex* . By multiplying the functions ( ) *ex E* Ψ and ( ) *ex F* Ψ as reported in fig. 12b, we notice that primary peaks positioned at half-integer Ψ*ex* values are left unchanged, while those appearing at integer Ψ*ex* values are depressed to zero.

By proceeding in the same way for the non-homogeneous case, starting from Eq. (62) we find that the critical current of the non homogeneous array described in Section 3 can be written as follows:

$$i\_c = E\left(\Psi\_{ex}\right)F\left(\Psi\_{ex}\right). \tag{67}$$

Notice that the resulting pattern does not depend on α, which can be absorbed, in the effective dynamical equation, by a rescaling of the normalized time τ. Notice also that the periodicity of the envelop function ( ) *ex E* Ψ may now not be commensurable with that of the rapidly oscillating term *AN* . Consider, in fact, the case reported in figs. 13a and 13b, where the pattern is shown, as a full line, for a non homogeneous array with *N*=12, ε = 1.3 and σ= 1.5 .

Effective Models of Superconducting Quantum Interference Devices 251

is a rational number, one has periodicity in the *<sup>c</sup>*

= 2 , as in figs. 14a-c, the

<sup>+</sup> needs to be a rational

*i*

σ

σ 1 σ

ε

fundamental periods of the functions ( ) *ex E* Ψ and ( ) *ex F* Ψ are incommensurable, so that the overall pattern is not periodic. This feature can be detected by gradually increasing the range of the pattern, as it is done in fig. 14a and fig. 14b. In this way, when looking for the

σ

vs. Ψ*ex* curves; otherwise, irregular patterns, like those shown in figs. 14a-b, are found.

**Figure 14.** Critical current as a function of Ψ*ex* for a non-homogeneous parallel array with *N*=12, *ε*=1.3,

(b)

4 2 0 2 4

*ex*

(a)

1.5 1.0 0.5 0.0 0.5 1.0 1.5

*ex*

= 2 in the following ranges: (a) [-1.75, 1.75]; (b) [-5.5, 5.5]. In (a) the envelope function ( ) *N E* Ψ*ex* is

We have seen that a reduced single-junction model can be adopted to describe the overall dynamics of a one-dimensional arrays with alternating parameters of *N*×(0-π) Josephson junctions. This effective model is very useful, since it allows to obtain the critical current vs. normalized magnetic flux curves in closed analytic form. The interference patterns are seen to be qualitatively similar to recently obtained experimental results on multifacets Josephson

= 1.3 and

when we consider an array with *N*=12,

number. Therefore, when the parameter

σ

shown as a dotted line.

conditions giving periodicity, one realizes that the ratio

*ic*

*ic*

**Figure 13.** Critical current as a function of Ψ*ex* for a non-homogeneous parallel array with *N*=12, *ε*=1.3, *σ*=1.5 in the following ranges: (a) [-1.25, 1.25]; (b) [-2.5, 2.5]. In (a) the envelope function ( ) *N E* Ψ*ex* is shown as a dotted line.

In fig. 13a the envelop curve ( ) *N E ex* ⋅ Ψ is shown as a dotted line. In the interference pattern in fig. 13a we may notice the finiteness of the small peak at 0 Ψ = *ex* due to the nonvanishing value of *E*( ) 0 1 = − ε . Also, notice the reduction in height of the primary peaks close to half-integer values as compared to those in fig. 12b. Finally, notice the appearance of two extra peaks of equal height in between two successive primary peaks. In fig. 13b, where the range of Ψ*ex* is increased, we notice that this feature repeats over a period 2 *ex* ΔΨ = , as it can be directly confirmed by inspection of Eqs. (65-67). Differently from fig. 13a and 13b, when we consider an array with *N*=12, ε = 1.3 and σ = 2 , as in figs. 14a-c, the fundamental periods of the functions ( ) *ex E* Ψ and ( ) *ex F* Ψ are incommensurable, so that the overall pattern is not periodic. This feature can be detected by gradually increasing the range of the pattern, as it is done in fig. 14a and fig. 14b. In this way, when looking for the conditions giving periodicity, one realizes that the ratio σ 1 σ <sup>+</sup> needs to be a rational number. Therefore, when the parameter σ is a rational number, one has periodicity in the *<sup>c</sup> i* vs. Ψ*ex* curves; otherwise, irregular patterns, like those shown in figs. 14a-b, are found.

250 Superconductors – Materials, Properties and Applications

rapidly oscillating term *AN*

shown as a dotted line.

vanishing value of *E*( ) 0 1 = −

ε

*ic*

σ= 1.5 .

Notice that the resulting pattern does not depend on

*ic*

effective dynamical equation, by a rescaling of the normalized time

the pattern is shown, as a full line, for a non homogeneous array with *N*=12,

periodicity of the envelop function ( ) *ex E* Ψ may now not be commensurable with that of the

**Figure 13.** Critical current as a function of Ψ*ex* for a non-homogeneous parallel array with *N*=12, *ε*=1.3, *σ*=1.5 in the following ranges: (a) [-1.25, 1.25]; (b) [-2.5, 2.5]. In (a) the envelope function ( ) *N E* Ψ*ex* is

(b)

2 1 0 1 2

*ex*

(a)

1.0 0.5 0.0 0.5 1.0

*ex*

In fig. 13a the envelop curve ( ) *N E ex* ⋅ Ψ is shown as a dotted line. In the interference pattern in fig. 13a we may notice the finiteness of the small peak at 0 Ψ = *ex* due to the non-

close to half-integer values as compared to those in fig. 12b. Finally, notice the appearance of two extra peaks of equal height in between two successive primary peaks. In fig. 13b, where the range of Ψ*ex* is increased, we notice that this feature repeats over a period 2 *ex* ΔΨ = , as it can be directly confirmed by inspection of Eqs. (65-67). Differently from fig. 13a and 13b,

. Also, notice the reduction in height of the primary peaks

α

. Consider, in fact, the case reported in figs. 13a and 13b, where

, which can be absorbed, in the

. Notice also that the

ε

= 1.3 and

τ

**Figure 14.** Critical current as a function of Ψ*ex* for a non-homogeneous parallel array with *N*=12, *ε*=1.3, σ = 2 in the following ranges: (a) [-1.75, 1.75]; (b) [-5.5, 5.5]. In (a) the envelope function ( ) *N E* Ψ*ex* is shown as a dotted line.

We have seen that a reduced single-junction model can be adopted to describe the overall dynamics of a one-dimensional arrays with alternating parameters of *N*×(0-π) Josephson junctions. This effective model is very useful, since it allows to obtain the critical current vs. normalized magnetic flux curves in closed analytic form. The interference patterns are seen to be qualitatively similar to recently obtained experimental results on multifacets Josephson junctions (Scharinger et al., 2010) in which the critical current density alternates many times between two opposite values along the junction length. Discrete Josephson junction arrays, even presenting some analogies with the latter devices, are much too simple systems to describe the complete behaviour of MJJs. As a matter of fact, shielding current effects is not taken into account by analysis carried out in this section in the lowest order approximation in β. Nevertheless, the analytic results obtained for the interference patterns shed some light on the causes of the presence or absence of periodicity and on the nature of primary and secondary peaks in the *<sup>c</sup> i* vs. Ψ*ex* curves. We finally notice that the analysis relies on the choice of an even number of Josephson junction in the array. Different patterns are expected for an odd number of alternating junctions in the array, depending also on which type (0 or π) of junctions is predominant. Further work is therefore necessary to address this problem in its fullest extent.

## **5. Quantum interferometers in the presence of rapidly varying fields**

In Section 2 we have analyzed the effective model for a two-junction quantum interferometer in the presence of an oscillating magnetic flux under the hypothesis that the frequency of oscillation is comparable with the inverse of the characteristic time evolution τφ of the superconducting time variable φ. In this way, the quasi-static approach described in Section 2 has been proven to be applicable. In the present section, on the other hand, we shall consider rapidly varying externally applied fluxes, whose frequency ω is comparable with <sup>1</sup> ψτ <sup>−</sup> , so that <sup>−</sup><sup>1</sup> >> ω τ φ.

Let us therefore consider an externally applied flux having d. c. component *A* and a. c. amplitude *B*, so that

$$
\Psi\_{ex}\left(t\right) = A + B\sin\alpha t\tag{68}
$$

Effective Models of Superconducting Quantum Interference Devices 253

). This hypothesis is confirmed by what

*τφ*. Within the former

, in which the

(70)

(71)

(72)

(73)

(74)

Δτψ

way, the effective dynamics for *φ* can be found when the solution for *Ψ* is substituted in the

Start by considering the cosine term in (69b) as a quasi-static quantity (i. e., it does not vary

already stated in the previous section; i. e., while the variables *φ* varies on a characteristic time

variable *φ* does not vary appreciably. We can thus solve (69b), by perturbation analysis, by

<sup>d</sup> 2 sin cos . <sup>d</sup> + +=

 φψ θ

() () () 0 1

() () <sup>d</sup> , <sup>d</sup> *<sup>f</sup>* + = *f g* θ

solutions of the system of ordinary differential equations (8a-b) and by considering the non-

( ) <sup>0</sup> <sup>2</sup> cos sin .

Having found the solution to (8a), we can find the solution to (8b) by the same type of

( ) <sup>1</sup>

 = −2 cos , *h* θ

( ) () () ,, , sin cos , *nk nk n k*

= −

αθω

= ++ + . Therefore, the SQUID dynamics can be described, to first

 αθ

, ( ) *nJ x* being the Bessel function of order *n*, and

(75)

βψ θ

> θ

θ

<sup>+</sup>

ωθ ω ωθ

φ

 = + ψ θ

θ

at large values of

ω

1 *B* = + *A* +

ψ

πψ

*τφ* it is then possible to choose a subinterval, of the order of

ψτ

> Δ*τΨ*= 2π*β*Δ*τφ* «Δ

() () ex

ψ θ

in (8a) and (8b). Recall that an ODE of the type

*g(x)ex*d*x*. By now considering (68), by taking the non-decaying

, we have

cosine term of (69a).

Δτφ

Δ

we again find the ODEs for

θ*)*=*e* −θ 

vanishing solution of (8a) for

with () ( )

=

*n k*

*g*

( ) , <sup>2</sup> <sup>2</sup> <sup>1</sup> *n n*

+ +

() ( ) , <sup>2</sup> *n k A nk n*

order in *β*, by the following ODE:

*n k* γ

> ωθ

*J J*

θ

β

interval

time interval

first setting *τ=*2π

By now setting

has solution *f(*

where

α θ π

appreciably over an interval of time of the order of

, the variable *Ψ* varies within a time interval

*θ* and by rewriting it as follows:

β

ψ θ

ψ

θ

ψ0 and ψ1

> ψ0

ψ

reasoning. After some rather long calculations one finds

θ

γω

 

ω

π

*n k h g* +∞ +∞

=−∞ =−∞

where − − 1 1 ≈ >> ωτ τ ψ φ is the frequency of the sinusoidal term. Let us again consider Eqs. (4ab) rewritten, for *n*=0, as follows:

$$\begin{aligned} \frac{\mathbf{d}\,\phi}{\mathbf{d}\,\pi} + \cos\pi\psi\sin\phi &= \frac{i\_\text{B}}{2}, \\ \pi\frac{\mathbf{d}\,\psi}{\mathbf{d}\,\pi} + \sin\pi\psi\cos\phi &= \frac{\mathcal{V}\_{\text{ex}} - \mathcal{V}}{2\mathcal{J}} \quad \text{(b)} \end{aligned} \tag{69}$$

where the normalization 0 <sup>2</sup> *RIJ <sup>t</sup>* = =*<sup>t</sup>* Φ φ π τ τ prescribes a *τ*-dependence of the externally

applied flux as written in Eq. (18) with a normalized frequency <sup>0</sup> 2 *RIJ* <sup>Φ</sup> = = ω ω ωτ φ π . We

therefore need to consider again all steps in Section 2, having care to integrate opportunely the right hand side term of (69b), in order to obtain a solution for *Ψ* in terms of the superconducting phase by perturbation analysis on βfor arbitrary values of *A* and *B*. In this way, the effective dynamics for *φ* can be found when the solution for *Ψ* is substituted in the cosine term of (69a).

Start by considering the cosine term in (69b) as a quasi-static quantity (i. e., it does not vary appreciably over an interval of time of the order of ψτ ). This hypothesis is confirmed by what already stated in the previous section; i. e., while the variables *φ* varies on a characteristic time interval Δτφ, the variable *Ψ* varies within a time interval Δ*τΨ*= 2π*β*Δ*τφ* «Δ*τφ*. Within the former time interval Δ*τφ* it is then possible to choose a subinterval, of the order of Δτψ, in which the variable *φ* does not vary appreciably. We can thus solve (69b), by perturbation analysis, by first setting *τ=*2πβ*θ* and by rewriting it as follows:

$$\frac{\text{d}\,\text{d}\,\text{y}}{\text{d}\,\text{d}\,\text{e}} + 2\,\beta\sin\pi\text{y}\,\text{cos}\,\phi + \text{y}\,(\theta) = \text{y}\_{\text{ex}}\,(\theta). \tag{70}$$

By now setting

252 Superconductors – Materials, Properties and Applications

in β

secondary peaks in the *<sup>c</sup>*

of the superconducting time variable

<sup>−</sup> , so that <sup>−</sup><sup>1</sup> >> ω τ φ.

> τ

b) rewritten, for *n*=0, as follows:

 φ

amplitude *B*, so that

where − − 1 1 ≈ >> ωτ

ψ

where the normalization

in its fullest extent.

with <sup>1</sup> ψτ

junctions (Scharinger et al., 2010) in which the critical current density alternates many times between two opposite values along the junction length. Discrete Josephson junction arrays, even presenting some analogies with the latter devices, are much too simple systems to describe the complete behaviour of MJJs. As a matter of fact, shielding current effects is not taken into account by analysis carried out in this section in the lowest order approximation

. Nevertheless, the analytic results obtained for the interference patterns shed some light on the causes of the presence or absence of periodicity and on the nature of primary and

choice of an even number of Josephson junction in the array. Different patterns are expected for an odd number of alternating junctions in the array, depending also on which type (0 or π) of junctions is predominant. Further work is therefore necessary to address this problem

**5. Quantum interferometers in the presence of rapidly varying fields** 

φ

shall consider rapidly varying externally applied fluxes, whose frequency

ψ

φ

τ

π

τ

superconducting phase by perturbation analysis on

ψ

τ

π

0 <sup>2</sup> *RIJ <sup>t</sup>* = =*<sup>t</sup>* Φ

In Section 2 we have analyzed the effective model for a two-junction quantum interferometer in the presence of an oscillating magnetic flux under the hypothesis that the frequency of oscillation is comparable with the inverse of the characteristic time evolution

Section 2 has been proven to be applicable. In the present section, on the other hand, we

Let us therefore consider an externally applied flux having d. c. component *A* and a. c.

( ) sin *ex*

*t AB t* = +

<sup>d</sup> cos sin , (a) d 2

 φ

<sup>−</sup> + =

+ =

 πψ

φ

τ

applied flux as written in Eq. (18) with a normalized frequency <sup>0</sup>

πψ

<sup>d</sup> sin cos (b) d 2

 φ

therefore need to consider again all steps in Section 2, having care to integrate opportunely the right hand side term of (69b), in order to obtain a solution for *Ψ* in terms of the

B

*i*

ex

β

 β

ψ ψ

ω

is the frequency of the sinusoidal term. Let us again consider Eqs. (4a-

*i* vs. Ψ*ex* curves. We finally notice that the analysis relies on the

. In this way, the quasi-static approach described in

prescribes a *τ*-dependence of the externally

ω

2 *RIJ* <sup>Φ</sup> = =

π

for arbitrary values of *A* and *B*. In this

 ω ωτ φ

ω

τφ

(68)

(69)

. We

is comparable

$$
\psi\left(\theta\right) = \psi\_0\left(\theta\right) + \beta\psi\_1\left(\theta\right) \tag{71}
$$

we again find the ODEs for ψ0 and ψ1in (8a) and (8b). Recall that an ODE of the type

$$\frac{\text{d}f}{\text{d}\,\theta} + f\left(\theta\right) = g\left(\theta\right),\tag{72}$$

has solution *f(*θ*)*=*e* −θ θ *g(x)ex*d*x*. By now considering (68), by taking the non-decaying solutions of the system of ordinary differential equations (8a-b) and by considering the nonvanishing solution of (8a) for ψ0 at large values of θ, we have

$$\Psi\_0 = A + \frac{B}{1 + \tilde{\phi}^2} (\cos \tilde{\phi}\theta + \tilde{\phi}\sin \tilde{\phi}\theta). \tag{73}$$

Having found the solution to (8a), we can find the solution to (8b) by the same type of reasoning. After some rather long calculations one finds

$$\Psi\_1 = -2h(\theta)\cos\phi\_\prime \tag{74}$$

where

$$h(\boldsymbol{\theta}) = \sum\_{n=-\infty}^{+\infty} \sum\_{k=-\infty}^{+\infty} g\_{n,k} \left[ \sin \alpha\_{n,k} \left( \boldsymbol{\theta} \right) - \tilde{\boldsymbol{\alpha}} \cos \alpha\_{n,k} \left( \boldsymbol{\theta} \right) \right],\tag{75}$$

with () ( ) ( ) , <sup>2</sup> <sup>2</sup> <sup>1</sup> *n n n k J J g n k* γ γω ω = + + , ( ) *nJ x* being the Bessel function of order *n*, and

() ( ) , <sup>2</sup> *n k A nk n*π α θ π ωθ = ++ + . Therefore, the SQUID dynamics can be described, to first order in *β*, by the following ODE:

$$\frac{\mathrm{d}\,\mathrm{d}\,\phi}{\mathrm{d}\,\mathrm{d}\,\tau} + X\_0 \left(\theta\right) \mathrm{sin}\,\phi + \pi\theta\hbar \left(\theta\right) Y\_0 \left(\theta\right) \mathrm{sin}\,\mathbf{2}\,\phi = \mathrm{i}\_{\mathrm{B}} / 2\,\mathrm{m},\tag{76}$$

Effective Models of Superconducting Quantum Interference Devices 255

= cos , , (82)

 γω

φ

.

, for instance. The critical current *<sup>c</sup>*

( ) (86)

*A* . The above expression can be

0 and getting the maximum stationary

*i* vs. *B* curves for various values of the

. In fig. 15b, on the other

( ) ,*B* . The *<sup>c</sup>*

( ) ,*B* is always increasing for

*i* in terms of the a. c. component *B* of the

ρ ω

ω

term , we find:

(84)

(83)

(85)

*i* of

ω.

> ω

> ω

*i* vs.

*n n*

+ −

*n-l*

γγω

γ

is the Kronecker delta. By inserting (81) in (80), we have

 π ρω

( ) () ( ) ( ) () ( ) 0 0 2 2 1 , 21 . *<sup>n</sup>*

= +−

() () () ( ) <sup>0</sup> ,, ,

( ) ( )

*n k n-l*

 , , 2 cos , , = π ρω

β =

*i* vs. *B* curves are shown for various values of the normalized frequency

. When the d. c. component *A* and the time-varying portion of the magnetic

ρ ω

From figs. 15a-b we notice Fraunhofer-like patterns of the critical current as shown as a

The particular shape of these patterns in figs. 15a-b depends on the value of the normalized

curves are represented in fig. 16 for *A*=0.1 and for various values of *B*. In this respect, we

flux attain fixed values, and we let the normalized frequency vary with continuity, the

1 for even

1 for odd

β

( ) ,*B* is the extra-factor modifying the usual form of a d. c. SQUID,

π

( )( ) ( ) , , <sup>1</sup> 2

0) can be expressed as follows:

− −

Having expressed the effective time-averaged terms in (77) in a closed analytic form, we can understand the effect of a high-frequency field on the electrodynamic behaviour a d. c.

ξ

<sup>1</sup> . <sup>2</sup> *nml nm l n m l*

*n BJ J J J* +∞

=

where the symbol *n m*,

with

where

the device in this case (

where the quantity

value for *<sup>B</sup>*

hand, the *<sup>c</sup>*

frequency

δ

ρ ω

θ

ξ

SQUID with extremely small value of the parameter

β =

ρ ω

*i* with respect to

function of the a. c. amplitude *B*.

ω

derived by inspection from Eq. (77), by setting

In fig. 15a-b we thus show the critical current *<sup>c</sup>*

applied flux. In particular, in fig. 15a, we report the *<sup>c</sup>*

 θ

− =

 *i AB A B <sup>C</sup>* ( ) ω

expressed, in terms of a constant applied flux *A*, as 2 cos

φ.

d. c. component *A* of the applied flux. normalized frequency

critical current of the device is the same as that of the quantity

notice that, while for small fixed *B* values the quantity

( ) ( ) <sup>0</sup> *X AB*

θ

 γω

, ,

=−∞

*nml h Y gJ J* +∞

γ

Proceeding in a similar way in finding the effective coefficient of the sin2

<sup>=</sup>

*n l*

−

*nml n l*

−+ −

where ( ) 0 0 *X* θ π = cos ψ and ( ) 0 0 *Y* θ π = sin ψ.

Equation (76) thus represents the differential equation describing the dynamics of the superconducting phase difference φ in a d. c. SQUID in the presence of a time-varying externally applied flux, whose frequency ω is considered to be comparable with τψ -1 , in such a way that 1 ω ωτψ = ≈ . For slowly varying fields (ω << 1 ) one can readily verify from (8a) and (8b) that () () ψ 0 ex θ →ψ θ and 1 ex ψ ψφ → −2sin cos , respectively. In this way, the dynamics described by the quasi-static d. c. SQUID in (3) holds.

For 1 ω ωτψ = ≈ , time evolution of the two variables, φ and ψ, still occurs with two completely different time scales. In fact, as already stated, one has τψ= 2πβτφ «τφ. Therefore, flux motion is very fast with respect to the dynamics of the phase variable φ. The only difference, here, is that the externally applied flux ( ) ψ ex θ is able to follow this fast dynamics. Having carefully solved (8a) and (8b), and having found the effective single-junction dynamical equation for a d. c. SQUID, we can determine the effective time-averaged equation, by taking the time average over the fast variable ψ. Therefore, we may write

$$\frac{\mathrm{d}\,\mathrm{d}\,\phi}{\mathrm{d}\,\mathrm{d}\,\tau} + \left<\mathrm{X}\_{0}\left(\theta\right)\right>\mathrm{sin}\,\phi + \pi\beta \left\mathrm{sin}\,2\phi = \mathrm{i}\_{\mathrm{B}}\left/2\,\right.\tag{77}$$

where the symbol <*x*> stands for the time average of the variable *x*. Equation (77) can thus be considered an effective single-junction model for a d. c. SQUID in the presence of a rapidly varying magnetic field ( 1 *<sup>L</sup>* ω ω = ≈ *t* ). The average values ( ) *X*<sup>0</sup> θ and () () <sup>0</sup> *h Y* θ θ can be calculated as follows. First of all, set

$$X\_0\left(\theta\right) = \text{Re}\left[\exp(i\pi\varphi\_0)\right].\tag{78}$$

Let us next express the exponential of a cosine and a sine terms in 0 exp( ) *i*πψ through the following Bessel function identities

$$\begin{aligned} e^{ia\cos x} &= \sum\_{n=-\infty}^{+\infty} i^n f\_n \left( a \right) e^{inx} \text{ (a)}\\ e^{ia\sin x} &= \sum\_{n=-\infty}^{+\infty} f\_n \left( a \right) e^{inx} . \end{aligned} \tag{79}$$

In this way, (78) becomes

$$\left\langle X\_{0}\left(\theta\right)\right\rangle = \text{Re}\left[\sum\_{n,k} J\_{n}\left(\mathcal{Y}\right)I\_{k}\left(\mathcal{Y}\tilde{o}\right)\left\langle e^{i\alpha\_{n,k}\left(\theta\right)}\right\rangle\right].\tag{80}$$

It is now easy to show that

$$\left\langle e^{i\alpha\_{n,k}(\theta)} \right\rangle = e^{i\pi A} i^n \delta\_{n,-k} \tag{81}$$

where the symbol *n m*, δis the Kronecker delta. By inserting (81) in (80), we have

$$\left< X\_0 \left( \theta \right) \right> = \cos \pi A \,\rho \left( \left. \tilde{\alpha} \right> , \tag{82}$$

with

254 Superconductors – Materials, Properties and Applications

 π = cos ψ

superconducting phase difference

and (8b) that () () ψ

the externally applied flux ( )

varying magnetic field ( 1 *<sup>L</sup>*

calculated as follows. First of all, set

following Bessel function identities

In this way, (78) becomes

It is now easy to show that

externally applied flux, whose frequency

 0 ex θ →ψ

ψ

φ

τ

ω ω

where ( ) 0 0 *X* θ

a way that 1 ω ωτψ

For 1 ω ωτψ

over the fast variable

φ

τ

θ

 and ( ) 0 0 *Y* θ

= ≈ . For slowly varying fields (

dynamics described by the quasi-static d. c. SQUID in (3) holds.

. Therefore, we may write

 φ πβ θ

θ

Let us next express the exponential of a cosine and a sine terms in 0 exp( ) *i*

cos

sin

θ θ

θ

= ≈ , time evolution of the two variables,

very fast with respect to the dynamics of the phase variable

different time scales. In fact, as already stated, one has

ψ ex θ  φ πβ θ

 π = sin ψ.

φ

() () () 0 0B <sup>d</sup> sin sin2 2 , <sup>d</sup> ++ = *X hY i m*

Equation (76) thus represents the differential equation describing the dynamics of the

(8a) and (8b), and having found the effective single-junction dynamical equation for a d. c. SQUID, we can determine the effective time-averaged equation, by taking the time average

> ( ) () () 0 0B <sup>d</sup> sin sin 2 2 , <sup>d</sup> ++ = *X hY i*

where the symbol <*x*> stands for the time average of the variable *x*. Equation (77) can thus be considered an effective single-junction model for a d. c. SQUID in the presence of a rapidly

= ≈ *t* ). The average values ( ) *X*<sup>0</sup>

*X i* 0 0 ( )

*ia x n inx n*

*e iJ ae*

=−∞ +∞

=−∞

+∞

*n ia x inx n*

=

=

*n*

( ) ,

α θ

*e J ae*

= Re exp( ) . {

( )

, (a)

. (b)

α θ

( )

( ) () ( ) ( ) { } , <sup>0</sup> , Re . *n k <sup>i</sup> <sup>X</sup> n k n k* <sup>=</sup> *JJ e*

> , , *n k <sup>i</sup> i An n k e ei* <sup>−</sup> =

π δ

γ

 γω

 θ

πψ

ω

 and 1 ex ψ

 θ

ω

φ and ψ

> τψ= 2πβτφ «τφ

> > φ

is able to follow this fast dynamics. Having carefully solved

 φ

θ

 ψφ

 φ

in a d. c. SQUID in the presence of a time-varying

→ −2sin cos , respectively. In this way, the

is considered to be comparable with

(76)


(77)

can be

through the

(79)

(80)

(81)

τψ

<< 1 ) one can readily verify from (8a)

, still occurs with two completely

. The only difference, here, is that

 and () () <sup>0</sup> *h Y* θ

πψ

} (78)

 θ

. Therefore, flux motion is

$$f\_{\nu}(\tilde{\alpha}; \mathcal{B}) = f\_0\left(\boldsymbol{\chi}\right)f\_0\left(\boldsymbol{\chi}\tilde{\boldsymbol{\alpha}}\right) + 2\sum\_{n=1}^{+\ast} \left(-1\right)^n f\_{2n}\left(\boldsymbol{\chi}\right)f\_{2n}\left(\boldsymbol{\chi}\tilde{\alpha}\right). \tag{83}$$

Proceeding in a similar way in finding the effective coefficient of the sin2φterm , we find:

$$\left\langle h\left(\boldsymbol{\theta}\right)\boldsymbol{\Upsilon}\_{0}\left(\boldsymbol{\theta}\right)\right\rangle = \frac{1}{2} \sum\_{n,m,l=-\infty}^{+\infty} \xi\_{n,m,l} g\_{n,m} \boldsymbol{J}\_{l}\left(\boldsymbol{\chi}\right) \boldsymbol{J}\_{n+m-l}\left(\boldsymbol{\chi}\tilde{\boldsymbol{w}}\right). \tag{84}$$

where

$$
\xi\_{n,m,l} = \begin{cases}
{(-1)}^{n-l} & \text{for } {(n \cdot l)} \text{ even} \\
{-(n+k)(-1)}^{\frac{n-l-1}{2}} & \text{for } {(n \cdot l)} \text{ odd}
\end{cases} \tag{85}
$$

Having expressed the effective time-averaged terms in (77) in a closed analytic form, we can understand the effect of a high-frequency field on the electrodynamic behaviour a d. c. SQUID with extremely small value of the parameter β, for instance. The critical current *<sup>c</sup> i* of the device in this case (β = 0) can be expressed as follows:

$$\mathrm{i}\_C\left(\tilde{\alpha}\prime, A\,, B\right) = 2\left|\cos\pi\,A\right|\left|\rho\left(\tilde{\alpha}\prime B\right)\right|\,,\tag{86}$$

where the quantity ρ ω( ) ,*B* is the extra-factor modifying the usual form of a d. c. SQUID, expressed, in terms of a constant applied flux *A*, as 2 cosπ *A* . The above expression can be derived by inspection from Eq. (77), by setting β = 0 and getting the maximum stationary value for *<sup>B</sup> i* with respect to φ.

In fig. 15a-b we thus show the critical current *<sup>c</sup> i* in terms of the a. c. component *B* of the applied flux. In particular, in fig. 15a, we report the *<sup>c</sup> i* vs. *B* curves for various values of the d. c. component *A* of the applied flux. normalized frequency ω . In fig. 15b, on the other hand, the *<sup>c</sup> i* vs. *B* curves are shown for various values of the normalized frequency ω . From figs. 15a-b we notice Fraunhofer-like patterns of the critical current as shown as a function of the a. c. amplitude *B*.

The particular shape of these patterns in figs. 15a-b depends on the value of the normalized frequency ω . When the d. c. component *A* and the time-varying portion of the magnetic flux attain fixed values, and we let the normalized frequency vary with continuity, the ω critical current of the device is the same as that of the quantity ρ ω( ) ,*B* . The *<sup>c</sup> i* vs. ω curves are represented in fig. 16 for *A*=0.1 and for various values of *B*. In this respect, we notice that, while for small fixed *B* values the quantity ρ ω( ) ,*B* is always increasing for increasing values of ω (see orange line in fig. 16) for values of *B* approaching 1.0, this character is lost (see brown and cyan lines in fig. 16).

Effective Models of Superconducting Quantum Interference Devices 257

( ) ,*B* . Because of this dependence, the

*i* = , for

= 0 (orange line) decreases as we let *B* increase

can be viewed as control

*i* = and *B*=0.0 (orange);

ω= 1.0 .

ω

= 0 , in (b)

ω

= 0 (a) and

(87)

We can now calculate, for *β=0*, the flux-voltage curves (v vs. *A*) by the following well known

<sup>2</sup> <sup>d</sup> 2 2 <sup>v</sup> cos , , d 4 *B <sup>i</sup>* = =− *<sup>A</sup> <sup>B</sup>*

= 1.0 (b), and for various values of *B*. In particular, in fig. 17a we notice that the amplitude

to 0.15 first (brown line) and to 0.30 next (cyan line). The same decreasing behavior is

(a)

*A*

0 1 2 3 4

(b)

*A*

0 1 2 3 4

ω

φ

τ

amplitude of the v vs. *A* curves can be varied and *B* and

though the extra-factor

parameters. In figs. 17a-b we thus report the v vs. *A* curves for 2.5 *<sup>B</sup>*

**Figure 17.** Voltage v versus the d. c. component *A* of the applied flux for 2.5 *<sup>B</sup>*

*B*=0.15 (brown); *B*=0.30 (cyan). In (a) the normalized frequency values is

( )( )

. The important feature in this expression is that these curves

= 0.1 when *B* increases from 0.0 (orange line) to 0.15 (brown line)

ω

π ρω

ρ ω

expression (Barone & Paternò, 1982):

ω

ω

0.8

0.8

0.9

1.0

v

1.1

1.2

0.9

1.0

v

1.1

1.2

of the v vs. *A* curves obtained at *B*=0 and

for *i AB <sup>B</sup>* > 2 cos , ( )( ) π ρω

depend both on *B* and

detected in fig. 17b for

and to 0.30 (cyan line).

ω

**Figure 15.** Critical current *<sup>c</sup> i* as a function of the oscillating amplitude *B* . (a) The normalized frequency is fixed at ω = 1.0 and various values of the d. c. component *A* are represented: *A*=0.1 (orange); *A*=0.2 (brown); *A*=0.3 (cyan). (b) The d. c. component *A* is fixed at *A*=0.1 and the normalized frequency attains the following values: ω = 0.5 (orange); ω = 1.0 (brown); ω= 1.5 (cyan).

**Figure 16.** Critical current *<sup>c</sup> i* as a function of the normalized frequency ω for *A*=0.1 and for various values of the oscillating amplitude *B*: *B*=0.4 (orange); *B*=0.8 (brown); *B*=1.2 (cyan).

We can now calculate, for *β=0*, the flux-voltage curves (v vs. *A*) by the following well known expression (Barone & Paternò, 1982):

256 Superconductors – Materials, Properties and Applications

ω

character is lost (see brown and cyan lines in fig. 16).

0.0

0.5

1.0

*iC*

1.5

(see orange line in fig. 16) for values of *B* approaching 1.0, this

*i* as a function of the oscillating amplitude *B* . (a) The normalized frequency

= 1.5 (cyan).

ω

for *A*=0.1 and for various

= 1.0 and various values of the d. c. component *A* are represented: *A*=0.1 (orange); *A*=0.2

(b)

4 2 0 2 4

*B*

(a)

*B*

4 2 0 2 4

ω

(brown); *A*=0.3 (cyan). (b) The d. c. component *A* is fixed at *A*=0.1 and the normalized frequency attains

= 1.0 (brown);

*i* as a function of the normalized frequency

0 1 2 3 4

ω

increasing values of

**Figure 15.** Critical current *<sup>c</sup>*

**Figure 16.** Critical current *<sup>c</sup>*

ω

= 0.5 (orange);

0.0

0.5

1.0

*iC*

1.5

0.0

0.5

1.0

*iC*

1.5

ω

values of the oscillating amplitude *B*: *B*=0.4 (orange); *B*=0.8 (brown); *B*=1.2 (cyan).

ω

the following values:

is fixed at

$$\mathbf{v} = \left\langle \frac{\mathbf{d}\,\phi}{\mathbf{d}\,\tau} \right\rangle = \sqrt{\frac{\mathbf{i}\_B{B}^2}{4} - \cos^2\left(\pi A\right)\rho^2 \left(\tilde{\phi}, B\right)},\tag{87}$$

for *i AB <sup>B</sup>* > 2 cos , ( )( ) π ρω . The important feature in this expression is that these curves depend both on *B* and ω though the extra-factor ρ ω( ) ,*B* . Because of this dependence, the amplitude of the v vs. *A* curves can be varied and *B* and ω can be viewed as control parameters. In figs. 17a-b we thus report the v vs. *A* curves for 2.5 *<sup>B</sup> i* = , for ω = 0 (a) and ω = 1.0 (b), and for various values of *B*. In particular, in fig. 17a we notice that the amplitude of the v vs. *A* curves obtained at *B*=0 and ω = 0 (orange line) decreases as we let *B* increase to 0.15 first (brown line) and to 0.30 next (cyan line). The same decreasing behavior is detected in fig. 17b for ω = 0.1 when *B* increases from 0.0 (orange line) to 0.15 (brown line) and to 0.30 (cyan line).

**Figure 17.** Voltage v versus the d. c. component *A* of the applied flux for 2.5 *<sup>B</sup> i* = and *B*=0.0 (orange); *B*=0.15 (brown); *B*=0.30 (cyan). In (a) the normalized frequency values is ω = 0 , in (b) ω= 1.0 .

In fig. 18a-b, finally, by fixing the value of *B* first to 0.5 (a) and then to 1.0 (b), we notice that the amplitude of the v vs. *A* curves increases for increasing values of ω , as shown for ω = 0.5 (orange line), ω = 1.0 (brown line), and ω= 1.5 (cyan line).

Effective Models of Superconducting Quantum Interference Devices 259

of the applied flux is much less then <sup>1</sup>

of the dynamics of the

ψτ. The

> ψτ<sup>−</sup> .

2 *RIJ* <sup>Φ</sup> <sup>=</sup> φ τ

π

becomes evident that the characteristic time constant <sup>0</sup>

ω ωτψ= .

limit, assuming that the oscillation frequency

from the normalized frequancy

an equivalent single-junction model.

superconducting phase difference

ω

can be of interest.

a frequency

average superconducting phase difference *φ* is different from the characteristic time

perturbation analysis has been carried out for both a constant and a time-dependent applied magnetic flux. The circulating current *iS* and the time average of the critical current <*ic*> as a function of the d. c. and a. c. components of the applied flux are evaluated in the adiabatic

ω

In this limit the Fraunhofer-like pattern in <*ic*> vs. *B* curves are shown to be independent

Next, the dynamical equations of one-dimensional arrays containing *N*+1 identical overdamped Josephson junctions are considered. It has been noticed that the system of *N*+1 nonlinear first-order ordinary differential equations can be broken into two coupled subsystems, one consisting of only one equation for the superconducting phase of one junction in the array (arbitrarily chosen to be *φ0*), the second describing the time evolution of *N* opportunely defined normalized flux variables. When a solution of the latter *N* equations is found, by means of a perturbative approach to first order in the parameter *β*, the dynamical properties of the system are described by a single time-evolution equation for *φ0*. In this way, we may affirm that, for small values of *β*, the system may be described by an equivalent single-junction model, where the maximum Josephson current is appropriately defined. The analysis represents a simple way of approaching the problem of the electrodynamic response of one-dimensional arrays of overdamped Josephson junctions by

The same approach is followed for one-dimensional arrays with alternating parameters of *N*×(0-π) Josephson junction. Even in this case an effective single-junction model can be adopted to describe the overall dynamics of the system. This model is very useful, since it allows us to obtain the critical current vs. normalized magnetic flux curves in closed analytic form. The interference patterns are seen to be qualitatively similar to recently obtained experimental results on multifacets Josephson junctions (Scharinger et al., 2010) in which the critical current density alternates many times between two opposite values along the junction length. The analytic results for the interference patterns clarify the presence or absence of periodicity and the nature of primary and secondary peaks in these curves. Further investigation on the dependence of *ic* on different distributions of the JJs in the array

Finally, by allowing the magnetic flux, applied to a two-junction superconducting quantum interference device, to have an a. c. component in addition to a constant term *A*, we derive the effective reduced single-junction model describing the dynamics of the average

between this case and the one previously treated is that the alternating flux now varies with

not apply. The single-junction model for the system is again obtained by perturbation

ψτ

of the two junctions in the device. The difference

<sup>−</sup> , so that the adiabatic approach does

φ

of the same order of magnitude of <sup>1</sup>

**Figure 18.** Voltage v versus the d. c. component *A* of the applied flux for 2.5 *<sup>B</sup> i* = and ω = 0.5 (orange); ω = 1.0 (brown); ω = 1.5 (cyan). In (a) the a. c. component of the applied magnetic flux is *B*=0.50, in (b) *B*=1.0.

## **6. Conclusion**

We have studied the dynamical properties of quantum interferometers consisting of single or multiple superconducting loops, each containing two Josephson junctions.

A symmetric quantum interferometer containing two identical junctions with negligible capacitance has been considered first. The analysis of the system has been carried out by means of a perturbation approach in the parameter *β*, whose value gives the strength of the electromagnetic coupling between the two junction in the system. We have noticed that the flux-number function ψ β( ) ,θ governs fluxon dynamics, where θis the laboratory time *t*

normalized to the characteristic circuital time constant *<sup>L</sup> R* ψτ= . By this general approach it becomes evident that the characteristic time constant <sup>0</sup> 2 *RIJ* <sup>Φ</sup> <sup>=</sup> φ τ π of the dynamics of the average superconducting phase difference *φ* is different from the characteristic time ψτ . The perturbation analysis has been carried out for both a constant and a time-dependent applied magnetic flux. The circulating current *iS* and the time average of the critical current <*ic*> as a function of the d. c. and a. c. components of the applied flux are evaluated in the adiabatic limit, assuming that the oscillation frequency ω of the applied flux is much less then <sup>1</sup> ψτ <sup>−</sup> . In this limit the Fraunhofer-like pattern in <*ic*> vs. *B* curves are shown to be independent from the normalized frequancy = .

ω ωτψ

258 Superconductors – Materials, Properties and Applications

ω

v

0.95 1.00 1.05 1.10 1.15 1.20 1.25

> 1.19 1.20 1.21 1.22 1.23 1.24 1.25

v

ω

ω

*B*=1.0.

= 1.0 (brown);

**6. Conclusion** 

flux-number function

ψ β( ) ,θ

normalized to the characteristic circuital time constant *<sup>L</sup>*

ω

= 0.5 (orange line),

In fig. 18a-b, finally, by fixing the value of *B* first to 0.5 (a) and then to 1.0 (b), we notice that

ω

(a)

0 1 2 3 4

*A*

= 1.5 (cyan line).

ω

*i* = and

θ

*R* ψτ

= 1.5 (cyan). In (a) the a. c. component of the applied magnetic flux is *B*=0.50, in (b)

(b)

0 1 2 3 4

*A*

We have studied the dynamical properties of quantum interferometers consisting of single

A symmetric quantum interferometer containing two identical junctions with negligible capacitance has been considered first. The analysis of the system has been carried out by means of a perturbation approach in the parameter *β*, whose value gives the strength of the electromagnetic coupling between the two junction in the system. We have noticed that the

governs fluxon dynamics, where

ω

is the laboratory time *t*

= . By this general approach it

= 0.5 (orange);

, as shown for

the amplitude of the v vs. *A* curves increases for increasing values of

= 1.0 (brown line), and

**Figure 18.** Voltage v versus the d. c. component *A* of the applied flux for 2.5 *<sup>B</sup>*

or multiple superconducting loops, each containing two Josephson junctions.

Next, the dynamical equations of one-dimensional arrays containing *N*+1 identical overdamped Josephson junctions are considered. It has been noticed that the system of *N*+1 nonlinear first-order ordinary differential equations can be broken into two coupled subsystems, one consisting of only one equation for the superconducting phase of one junction in the array (arbitrarily chosen to be *φ0*), the second describing the time evolution of *N* opportunely defined normalized flux variables. When a solution of the latter *N* equations is found, by means of a perturbative approach to first order in the parameter *β*, the dynamical properties of the system are described by a single time-evolution equation for *φ0*. In this way, we may affirm that, for small values of *β*, the system may be described by an equivalent single-junction model, where the maximum Josephson current is appropriately defined. The analysis represents a simple way of approaching the problem of the electrodynamic response of one-dimensional arrays of overdamped Josephson junctions by an equivalent single-junction model.

The same approach is followed for one-dimensional arrays with alternating parameters of *N*×(0-π) Josephson junction. Even in this case an effective single-junction model can be adopted to describe the overall dynamics of the system. This model is very useful, since it allows us to obtain the critical current vs. normalized magnetic flux curves in closed analytic form. The interference patterns are seen to be qualitatively similar to recently obtained experimental results on multifacets Josephson junctions (Scharinger et al., 2010) in which the critical current density alternates many times between two opposite values along the junction length. The analytic results for the interference patterns clarify the presence or absence of periodicity and the nature of primary and secondary peaks in these curves. Further investigation on the dependence of *ic* on different distributions of the JJs in the array can be of interest.

Finally, by allowing the magnetic flux, applied to a two-junction superconducting quantum interference device, to have an a. c. component in addition to a constant term *A*, we derive the effective reduced single-junction model describing the dynamics of the average superconducting phase difference φ of the two junctions in the device. The difference between this case and the one previously treated is that the alternating flux now varies with a frequency ω of the same order of magnitude of <sup>1</sup> ψτ <sup>−</sup> , so that the adiabatic approach does not apply. The single-junction model for the system is again obtained by perturbation analysis to first order in the parameter β. Averaging of the rapidly varying quantities in the differential equation for φ gives the effective dynamics of the two junctions in the system. In particular, for β = 0, the critical current of the device is seen to depend on *A*, on the frequency ω and the amplitude *B* of the a. c. component of the applied magnetic flux in a closed analytic form. From the analysis of the voltage vs. applied flux curves it can be argued that the quantities ω and B can play the role of additional control parameters in the device. Further work in extending the present analysis to finite values of β is necessary. Experimental work confirming the predictions of the present analysis needs to be performed. As far as non-normalized quantities are concerned, for direct experimental confirmation of the present results, we finally notice that the junction dynamics evolves with characteristic frequencies the order of 1 THz. Therefore one needs to run the experiment with very rapidly oscillating signals (10 THz or more) in such a way that normalized frequencies of ω= 1.0 can be achieved.

Effective Models of Superconducting Quantum Interference Devices 261

Clarke, J. & Braginsky, A. I. (2004). *The SQUID Handbook,* Vol. I, Wiley-VCH, ISBN 3-527-

Crankshaw, D. S. & Orlando, T. P. (2001). Inductance Effect in the Persistent Current Qubit. *IEEE Transactions on Applied Superconductivity,* Vol.11, No.1, (March 2001), pp. 1006-

De Luca, R. (2011). Quantum interference in parallel connections of N×(0–π) overdamped Josephson junctions. *Supercon. Sci. Technol.,* Vol.24, No.6, (June 2011), pp. 065026 1-5,

Geshkenbein, V. B.; Larkin, A. I. & Barone, A. (1987). Vortices with Half Magnetic Flux Quanta in "Heavy Fermion" Superconductors . *Phys. Rev. B,* Vol.36, No.1, (July 1987),

Grønbech-Jensen, N.; Thompson, D. B.; Cirillo, M. & Cosmelli, C. (2003). Thermal escape from zero-voltage states in hysteretic superconducting interferometers. *Phys. Rev. B,* 

Likharev, K. K. (1986). *Dynamics of Josephson Junctions and Circuits,* Gordon and Breach, ISBN

Romeo, F. & De Luca, R. (2004). Effective Non-Sinusoidal Current-Phase Relations in Conventional d. c. SQUIDs. *Physics Letters A,* Vol.328, No.4-5, (August 2004), pp. 330-

Romeo, F. & De Luca, R. (2005). A reduced model for one-dimensional arrays of overdamped Josephson junctions. *Physica C,* Vol.432, No.3-4, (November 2005), pp. 159-

Ryazanov, V. V.; Oboznov, V. A.; Rusanov, A. Yu.; Veretennikov, A. V.; Golubov, A. A. & Aarts, J. (2001). Coupling of Two Superconductors Through a Ferromagnet: Evidence for a π Junction. *Phys. Rev. Lett.,* Vol.86, No.11, (March 2001), pp. 2427-2430, ISSN 0031-

Scharinger, S.; Gürlich, C.; Mints, R. G.; Weides, M.; Kohlstedt, H.; Goldobin, E.; Koelle, D. & Kleiner, R. (2010). Interference patterns of multifacet 20x(0-π) Josephson junctions with ferromagnetic barrier. *Phys. Rev. B,* Vol.81, No.17, (May 2010), pp. 174535 1-4, ISSN

Schultz, R. R.; Chesca, B.; Goetz, B.; Schneider, C. W.; Schmehl, A.; Bielefeldt, H.; Hilgenkamp, H.; Mannhart, J. & Tsuei, C. C. (2000). Design and Realization of an all dwave dc π Superconducting Quantum Interference Device. *Appl. Phys. Lett.,* Vol.76,

Smilde, H. J. H.; Ariando; Blank, D. H. A.; Hilgenkamp, H. & Rogalla, H. (2004). π SQUIDs based on Josephson Contacts Between High-Tc and Low-Tc Superconductors . *Phys. Rev.* 

Weides, M.; Kemmler, M.; Goldobin, E.; Koelle, D., Kleiner, R.; Kohlstedt, H. & Buzdin, A. (2006). High quality ferromagnetic 0 and π Josephson tunnel junctions, *Appl. Phys. Lett.*, Vol.89, No.12, (September 2006), pp. 122511 1-3, ISSN 0003-

Vol.67, No.22, (June 2003), pp. 224505 1-6, ISSN 1098-0121

2881240429, Amsterdam, The Netherlands

No.7, (February 2000), pp. 912-914, ISSN 0003-6951

*B*, Vol.70, No.2, (July 2004), pp. 024519 1-12, ISSN 1098-0121

40229-2, Weinheim, Germany

pp. 235-238, ISSN 1098-0121

334, ISSN 0375-9601

166, ISSN 0921-4534

9007

6951

1098-0121

1009, ISSN 1051-8223

ISSN 0953-2048
