Effective Models of Superconducting Quantum Interference Devices

It is well known that the electrodynamic properties of SQUIDs (Superconducting Quantum Interference Devices) are obtained by means of the dynamics of the Josephson junctions in these superconducting system (Barone & Paterno, 1982; Likharev, 1986; Clarke & Braginsky, 2004). Due to the intrinsic macroscopic coherence of superconductors, r. f. SQUIDs have been proposed as basic units (qubits) in quantum computing (Bocko et al., 1997). In the realm of quantum computing non-dissipative quantum systems with small (or null) inductance parameter and finite capacitance of the Josephson junctions (JJs) are usually considered (Crankshaw & Orlando, 2001). The mesoscopic non-simply connected classical devices, on the other hand, are generally operated and studied in the overdamped limit with negligible capacitance of the JJs and small (or null) values of the inductance parameter. Nonetheless, r. f. SQUIDs find application in a large variety of fields, from biomedicine to aircraft maintenance (Clarke & Braginsky, 2004), justifying actual scientific interest in them.


Introduction
It is well known that the electrodynamic properties of SQUIDs (Superconducting Quantum Interference Devices) are obtained by means of the dynamics of the Josephson junctions in these superconducting system (Barone & Paternò, 1982;Likharev, 1986;Clarke & Braginsky, 2004). Due to the intrinsic macroscopic coherence of superconductors, r. f. SQUIDs have been proposed as basic units (qubits) in quantum computing (Bocko et al., 1997). In the realm of quantum computing non-dissipative quantum systems with small (or null) inductance parameter and finite capacitance of the Josephson junctions (JJs) are usually considered (Crankshaw & Orlando, 2001). The mesoscopic non-simply connected classical devices, on the other hand, are generally operated and studied in the overdamped limit with negligible capacitance of the JJs and small (or null) values of the inductance parameter. Nonetheless, r. f. SQUIDs find application in a large variety of fields, from biomedicine to aircraft maintenance (Clarke & Braginsky, 2004), justifying actual scientific interest in them.
As for d. c. SQUIDs, these systems can be analytically described by means of a single junction model (Romeo & De Luca, 2004). The elementary version of the single-junction model for a d. c. SQUID takes the inductance L of a single branch of the device to be negligible, so that β = LIJ/Φ0 ≈0, where Φ0 is the elementary flux quantum and IJ is the average value of the maximum Josephson currents of the junctions. In this way, the Josephson junction dynamics is described by means of a nonlinear first-order ordinary differential equation (ODE) written in terms of the phase variable φ, which represents the average of the two gauge-invariant superconducting phase differences, φ 1 and φ 2 , across the junctions in the d. c. SQUID. By considering a device with equal Josephsons junction in each of the two symmetric branches, the dynamical equation of the variable φ can be written as follows (Barone & Paternò, 1982): where n is an integer denoting the number of fluxons initially trapped in the superconducting interference loop, τ=2πRIJt/Φ0=t/τ φ , R being the intrinsic resistive junction parameter, ψex=Φex/Φ0 is the externally applied flux normalized to Φ0 and iB= IB/IJ, is the bias 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 field-modulated maximum current IJF (IJF=|cosπψex|) in which a normalized bias current iB/2 flows. This strict equivalence comes from the hypothesis that the total normalized flux ψ=Φ/Φ0 linked to the interferometer loop can be taken to be equal to ψex. However, being ( ) 1 2 , we may say that the above hypothesis may be stated merely by means of the following identity: β=0. Therefore, for finite values of the parameter β, Eq. (1) is not anymore valid and the device behaves as if the equivalent Josephson junction possessed a non-conventional current-phase relation (CPR). In fact, for small finite values of β, one can see that the following model may be adopted (Romeo & De Luca, 2004): where Xex= cosπψex and Yex= sinπψex. A second-order harmonic in φ thus appears in addition to the usual sinφ term. The sin2φ addendum, however, arises solely from electromagnetic 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 show the additional sin2φ term for finite values of β. In order to understand how the reduction in the dimensional order of the dynamical equations is possible, it is noted that the quantities τ φ = Φ0/2πRIJ and τ ψ = L/R, denoting the characteristic time scales of the variables φ and of the number of fluxons ψ in the superconducting SQUID loop, respectively, are intimately linked to the parameter β, since τ ψ /τ φ = 2πβ. In this way, for constant applied magnetic fields, the flux dynamics for small values of β can be considered very fast with respect to the equivalent junction dynamics given in Eq. (3). As a consequence, the superconducting phase φ can be assumed to be adiabatic and the equation of motion for ψ in terms of the quasi-static variable φ can be solved by perturbation analysis. When the information for ψ is substituted back into the effective dynamical equation for φ,

Eq. (3) is finally obtained.
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 values in the N current loops of a one-dimensional array containing N+1 identical overdamped Josephson junctions, the dynamical equations for the gauge-invariant superconducting phase differences can be reduced to a single non linear differential equation (Romeo & De Luca, 2005). The resulting time-evolution equation is seen to be similar to the single-junction dynamical equation with an appropriately defined currentphase relation. As specified before, the critical current, the I-V characteristics and the fluxvoltage curves of the array can be determined analytically by means of the effective model. Furthermore, a one dimensional array of N cells of 0-and π-junctions in parallel can be considered (De Luca, 2011). In this case, by assuming that junctions parameters and effective loop areas alternate as one moves along the longitudinal direction of the array, going from 0-to π-junctions, an effective single junction model for the system can be derived. It can be shown that, by this model, interference patterns of the critical current as a function of the applied magnetic flux can be analytically found and compared with existing experiments (Scharinger et al., 2010).
Finally, a single-junction model for a d.c. SQUID is derived when we consider the effect of rapidly varying applied fields whose frequency ω is comparable with 1 ψ τ − . By letting the applied magnetic flux have, in addition to a constant term A, an a. c. component, we can find, by similar reasoning as in the case of a constant applied magnetic flux, the effective reduced single-junction model for the system. In particular, for β = 0, the critical current of the device is seen to depend on A, and 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.
The work will thus be organized as follows. In Section 2 the derivation of an effective singlejunction model for a symmetric d. c. SQUID containing two equal junctions will be briefly reviewed. In Section 3 the extension to this model to Josephson junction arrays with equal junctions in all branches will be considered. In Section 4 the case of the alternate presence of 0-and π-junctions in the array is considered, the system being similar to multifacets Josepshon junctions. In Section 5 the effective single-junction model for a d. c. SQUID in the presence of rapidly varying field is derived. Finally, in Section 6 conclusions are drawn and further investigations are suggested.

Two-junction quantum interference devices
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):  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 , ( ) k β φ τ 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: 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: Furthermore, we may also notice that , for k = 0, 1.
By the general procedure described above we get the following differential equations for the superconducting phase variables: 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, pertaining to the dynamics of the superconducting phase difference value φ. We have already noticed that 2 ψ φ τ πβ τ = , 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 ψ φ τ τ = , 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 0 ψ and 1 ψ : where the term d 2 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: and the following for the flux number variables: where k is an integer, On the other hand, in the case Finally, in the case 1 a = , we have where ( ) , 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 J I , given by (Barone & Paternò, 1982): For an arbitrary value n, which represents the number of fluxons initially trapped in the superconducting ring, we have 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.  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 ( ) where A is the normalized d. c. component of the applied flux and B is the normalized amplitude of the a. c. signal. . By considering normalized frequencies of the order of 1 φ τ − , we 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.
By noticing, however, that (19a) can be written in the following form: In Eq. (20) we find a perturbed solution in terms of the parameter b, so that, by setting ( ) ( ) ( ) 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 φ τ η τ η τ to be a known expression, we can then write: 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 . B 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 c i , defined as the maximum value of the current bias B 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: 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.
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 β.  , 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.

Multi-junction quantum interference devices
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 B I 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 J I , we take the inductance L of the horizontal upper branches to be such that 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 linked to the k-th cell of the array is related to the superconducting phase differences 1 k − φ and k φ across the two junctions in the cell, so that: where k n is an integer and 1, 2,..., k N = . 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 , where 0 S is the cell area, we may set is the normalized current flowing in the k-th branch and B B J 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: sin , Define now the partial sum n s ( 1 n N ≤ ≤ ) of the normalized fluxes as follows: By fluxoid quantization (Eq. (25)), setting all k n 's to zero (under the hypothesis of zero initially trapped flux in the array), we can write: , so that the dynamical equations (Eq. (27)) can be rewritten as follows: 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 n s .
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: (1) . n n n s s s ≈ +β By substituting the above expression in Eq. (32b) and, after having multiplied both members by β, by setting the coefficients of m β ( 0, 1 m = ) equal to zero, we obtain, for 1, 2,..., n N = : . n e x s n = Ψ For the first order corrections, on the other hand, we need to solve the following set of equations: (1) (1) 1 0 2 1 The above differential equation represents the effective model for the 1D-JJA represented in Fig. 7, containing 1 N + 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 1 β − terms.
We can now explicitly perform the sum in Eq. (37), so that we write: where the absolute value of ( ) 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: 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 0 0 where T is the period for the instantaneous voltage curve and a pseudo-period for the superconducting phase difference 0 φ (i.e., The pseudo-period T of the above solution can be found by inspection, so that: Therefore, by Eqs. (40) and (42), the I-V characteristics are given by the following expression: ( ) where only the positive branch has been chosen.
In Figs. 8a-b 1 N B + vs. ex Ψ curves are shown. In particular, in Fig. 8a the number of JJ's in the system ( ) 1 N + 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 1 N B + vs.   For fixed bias current values, the voltage versus flux curves can be obtained by Eq. (43) and is given by the following: 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. equal to the critical current of the system. In Fig. 10a, indeed, we see that disappearance of 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. 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 β.
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.

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., 2001M. 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. 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 B I 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.

The homogeneous case
Consider, as a first approach to the problem, the loop areas k S , the junctions resistive   ). 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 ( ) where k n is an integer and 1, 2,..., k M = (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 The dynamical equation for each Josephson junction in the array can be written by means of the RSJ model (Barone & Paternò, 1982), so that: for 1, 2,..., k M = , so that the dynamical equations can be rewritten as follows: where now 1 1, 2,..., 2 in Eqs. (51b) and (51c) and where it has been set 0 0 s = .
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: (1) . n n n s s s ≈ +β (52) By substituting the above expression in Eqs. (51b-d) and, after having multiplied both members by β, by setting to zero the coefficients of m β ( 0, 1 m = ), we obtain: (0) . n e x s n = Ψ (53) For the first order corrections, on the other hand, we need to solve the following set of equations: (1) (1) 1 0 2 1 It is now possible to explicitly calculate the finite sum in Eq. (56) to get, in terms of the of the individual (0-π) cells:

The non-homogeneous case
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 k S , the junctions resistive parameters k R and maximum Josephson currents Jk I , and the inductances k L of the horizontal upper branches alternate in their values, as we go along the array. In this way we write for all allowed values of k. In this way, the additional parameters α, ε, and σ are implicitly defined. As before, we set for 1, 2,..., k M = . The equations of the motion for the superconducting phases can now be written as follows: 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: By substituting the above expression in Eqs. (49) to obtain the superconducting phases k φ we then write the differential equations (60a-b) in terms of the sole phase variable 0 φ and of the 2N normalized flux variables k s . By following the same steps as in the previous section, we obtain the following equivalent single-junction model for the superconducting array: Naturally, for α, ε, and σ all equal to one, the ordinary differential equation (62) reduces to Eq. (57).

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 B i which can be injected in the system at zero voltage (57)). Therefore, by maximizing with respect to 0 φ the following expression 2 0 cos (2 1) , we can express the critical current of the homogeneous device as follows: 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)   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: ( ) ( ).   fig. 13a we may notice the finiteness of the small peak at 0 ex Ψ = due to the nonvanishing value of ( ) 0 1 E ε = − . 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 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 σ σ + needs to be a rational number. Therefore, when the parameter σ is a rational number, one has periodicity in the c i vs. ex Ψ curves; otherwise, irregular patterns, like those shown in figs. 14a-b, are found. 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  (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 c 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.

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 1 ψ τ − , so that where the normalization ω ω τ π . 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 ψ τ , 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 where the symbol <x> stands for the time average of the variable x. Equation (77) 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 c i of the device in this case (β = 0) can be expressed as follows: 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 B i with respect to φ.
In fig. 15a-b we thus show the critical current c i in terms of the a. c. component B of the applied flux. In particular, in fig. 15a, we report the c 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 c 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    We can now calculate, for β=0, the flux-voltage curves (v vs. A) by the following well known expression (Barone & Paternò, 1982):

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 ( ) 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 1 ψ τ − .
In this limit the Fraunhofer-like pattern in <ic> vs. B curves are shown to be independent from the normalized frequancy ψ ω ωτ =  .
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 1 ψ τ − , 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.