A theoretical model for parallel SQUID arrays with fluxoid focussing
AA theoretical model for parallel SQUID arrays with fluxoid focussing
K.-H. M¨uller ∗ and E. E. Mitchell CSIRO Manufacturing, PO Box 218, Lindfield, NSW 2070, Australia (Dated: September 14, 2020)We have developed a comprehensive theoretical model for predicting the magnetic field responseof a parallel SQUID array in the voltage state. The model predictions are compared with ourexperimental data from a parallel SQUID array made of a YBCO thin-film patterned into widetracks, busbars and leads, with eleven step-edge Josephson junctions. Our theoretical model usesthe Josephson equations for resistively shunted junctions as well as the second Ginzburg-Landauequation to derive a system of coupled first-order nonlinear differential equations to describe thetime-evolution of the Josephson junction phase differences which includes Johnson noise. Employ-ing the second London equation and Biot-Savart’s law, the supercurrent density distribution iscalculated, using the stream function approach, which leads to a 2D second-order linear Fredholmintegro-differential equation for the stream function with time-dependent boundary conditions. Thenovelty of the model is that it calculates the stream function everywhere in the thin-film structureto determine during the time-evolution the fluxoids for each SQUID array hole. Our numericalmodel calculations are compared with our experimental data and predict the bias-current versusvoltage and the voltage versus magnetic field response with unprecedented accuracy. The modelelucidates the importance to fully take Meissner shielding and current crowding into account inorder to properly describe fluxoid focussing and bias-current injection. Furthermore, our modelillustrates the failure of the simple lumped-element approach to describe a parallel SQUID arraywith a wide thin-film structure.
I. INTRODUCTION
In this paper we develop a comprehensive theoreticalmodel for a parallel SQUID array with over-dampedJosephson junctions (JJs), made from a high- T c su-perconducting thin-film with wide tracks, busbars andleads, exposed to an applied magnetic field while drivenby a bias-current into the voltage state. Opposite tothe commonly used lumped-element model, our modelcan describe for the first time the experimental data ofa parallel SQUID array with unprecedented accuracy.Highly accurate predictions are important for the opti-misation of magnetometers, low noise current amplifiersand high frequency AC magnetic field sensors.Enhanced quantum interference in a parallel SQUIDarray was first mentioned by Feynman et al. , and waspredicted to show up as a sharpening of the array’scritical current peak, seen in measurements of the array’scritical current as a function of an applied magneticfield. Such an enhancement was observed experimentallyfor the first time in the 1960s in parallel arrays of super-conducting point contacts . A three-junction parallelSQUID array was fabricated and theoretically describedin the early 1970s using a simple lumped-element model .Later a parallel SQUID array (named superconductingquantum interference grating (SQUIG)) with 11 JJswas theoretically modelled with an improved lumped-element model by Miller et al. , revealing in detail themechanism of the enhanced quantum interference effectand its dependence on the SQUID array’s screeningparameter β L . In the late 1980s and early 1990s interestin modelling of two-dimensional (2D) arrays of JJsemerged, due to the realisation that most HTS materials are granular where grain boundaries act as JJs . Incontrast to parallel SQUID arrays which only havevertical JJ connections, 2D JJ-arrays have both verticaland horizontal JJ connections which form the plaquettesof the array . In the early 2000s parallel SQUID arrayswith varying hole area sizes between junctions, namedsuperconducting quantum interference filters (SQIFs),were investigated experimentally and theoretically,again using a lumped-element model . Serial SQUIDarrays were also studied and the question of howto optimise linearity was addressed . Furthermore,the performance of parallel SQIF arrays put in series(2D SQIF arrays) have also been investigated by Kornev et al. , Taylor et al. and Mitchell et al. . In arecent review, theoretical and experimental studies ofdifferent SQUID arrays made from LTS materials andused as miniature antennas have been compared . Thesimilarities between interference patterns in parallelSQUID arrays and optical multiple slit gratings hasbeen discussed by De Luca . Last year one and twodimensional SQUID arrays have been investigated fur-ther experimentally and theoretically . In addition,a review about design tools and progress in modelling ofsuperconducting circuits was written by Fourie .The often used lumped-element model can only beapplied if a SQUID array consists of sufficiently narrowsuperonducting tracks such that Kirchhoff’s law can beapplied at well defined current vertices . However,SQUID arrays are usually made from thin-film structureswith wide tracks, busbars and leads, where the Meissnershielding from wide superconducting structures createsstrong magnetic flux-focussing and current crowding.Neither flux-focussing nor current crowding are notpart of the lumped-element model. Going beyond a r X i v : . [ c ond - m a t . s up r- c on ] S e p the lumped-element model, a two junction SQUIDarray (the normal DC SQUID) with a wide washerstructure in the zero voltage state has been investi-gated theoretically in an approximate way by Clem andBrandt . Terauchi et al. investigate the effect of widertracks on the shape of the voltage pulses of a DC SQUID.Here, in our paper, we have developed a compre-hensive theoretical model for parallel SQUID arraysin the voltage state. In particular, we consider widesuperconducting thin-film structures in the Meissnerstate and incorporate into the Josephson array equationsthe time-dependent supercurrent density distribution ofthe array, obtained from the second London equationand Biot-Savart’s law. In contrast to the lumped-element model, our model does not make use of anylumped-element inductances but instead calculates thevalues for the fluxoids of each hole in the array duringthe time-evolution of the JJ phase differences. Thestatic supercurrent density and magnetic field distribu-tions in different superconducting thin-film geometricalstructures, based on London and Maxwell equations,with and without DC current injection, but without anyJJ’s, has been studied previously .Our paper is organised as follows. In Sec. II we outlinein detail our theoretical model for parallel SQUID arrayswith wide superconducting thin-film structures. InSec. III we briefly mention our device fabrication andexperimental set up. In Sec. IV we present the resultsof our model calculations and compare them with ourYBCO thin-film array experimental data as well as withresults from a lumped-element model calculation. Wesummarise our findings in Sec. V. II. THEORECTICAL MODEL FOR PARALLELSQUID ARRAYS
The main goal of this paper is to calculate the voltageresponse of a thin-film parallel SQUID array to anapplied perpendicular magnetic field and compare ourresults with our experimental data for a YBCO thin-filmparallel SQUID array. Contrary to the commonly usedlumped-element model , we will use the second Londonequation and Biot-Savart’s law to calculate from thesupercurrent density within the array the fluxoids of theSQUID array holes during the oscillatory time-evolutionof the array .As an example, Fig.1 displays a planar thin-filmparallel SQUID array with N = 5 JJs and N − xy plane and the magnetic induction B a is applied perpendicular in z direction. The array issymmetric about both the x and y axis. FIG. 1: Example of a parallel thin-film parallel SQUIDarray with N = 5 JJs. The domain Ω (blue) is made ofa superconducting YBCO thin-film. All JJs have thesame width w J , and the hole domains Ω k are of thesame width w h and length 2 h . ∂ Ω k is the boundary ofhole k . ϕ k are the gauge invariant phase differencesacross the JJs, where k = 1 , ..., N . The JJs areconnected via wide superconducting horizontal busbarson the top and bottom, each of area 2 a ( b − h ) with wideattached superconducting leads, 2 c wide and l long, and w J wide tracks attached to the JJs. In our fabricateddevices, l is much longer than shown here. Abias-current I b is injected into the top lead and exitsfrom the bottom lead. The time-averaged voltage V ismeasured between the ends of the two leads. A. Josephson junction phase differences andfluxoids
In our case the width w J of the JJs is much less thanthe Josephson penetration depth λ J (short junction)and the applied magnetic induction B a is sufficientlysmall such that the JJ current density is nearly constantacross junction areas. In this case, the current acrossthe junctions is described by the Josephson equation I ck sin ϕ k ( t ) , where ϕ k ( t ) is the gauge-invariant phasedifference across the k th junction at time t and I ck thejunction critical current with k =1,..., N .Because in high-temperature superconducting materi-als, like YBCO thin-films, the capacitance of the fabri-cated JJs is small, one can apply the resistively shuntedjunction (RSJ) model to describe the time-dependent to-tal current I k ( t ) flowing through the junctions, i.e. I k ( t ) = V k ( t ) R k + I ck sin ϕ k ( t ) + I Noisek ( t ) . (1)Here R k is the normal resistance of the k th junction and V k ( t ) the voltage across that junction, which accordingto the Josephson equation is V k ( t ) = Φ π dϕ k ( t ) dt , (2)where Φ is the flux quantum. I Noisek ( t ) in Eq.(1) is theJohnson noise current, originating at finite temperaturefrom the resistor R k . This noise is often also calledNyquist noise or white noise.From the second Ginzburg-Landau equation it fol-lows that ϕ k +1 ( t ) − ϕ k ( t ) = 2 π Φ [ µ λ (cid:73) ∂ Ω k j ( r , t ) · d l + Φ k ( t ) ] , (3)with k =1, ..., N -1. Here µ is the permeability offree space, λ the London penetration depth, j ( r , t ) thesupercurrent density with r a spatial vector, the symbol · means scalar product and d l is an integration lineelement. The integration contour ∂ Ω k is chosen as theinner boundary contour of hole k as indicated in Fig.1, integrated in counterclockwise direction. In Eq. (3),Φ k ( t ) is the time-dependent total magnetic flux thatpenetrates the hole area Ω k (Fig. 1). The sum of thetwo terms in square brackets in Eq. (3) is called thefluxoid and is similar to London’s fluxoid though herethe fluxoid in Eq. (3) is not quantised.The time-dependent total flux Φ k ( t ) through the k th hole of the array is the sum of two parts,Φ k ( t ) = Φ ( a ) k + Φ ( J ) k ( t ) , (4)where Φ ( a ) k is the static applied flux through the holeΩ k , i.e. Φ ( a ) k = B a A k with A k the area of the hole k and B a the applied magnetic inductance, B a = µ H a where H a is the applied magnetic field. Here we restrictour investigation to the case where all hole areas A k areof the same size A h , with A h = 2 h w h (Fig. 1). Theflux Φ ( J ) k ( t ) in Eq. (4) is the flux spilled into the k th array hole and according to Biot-Savart’s law, Φ ( J ) k ( t )originates from the supercurrent density j ( r , t ) that is flowing throughout the whole superconducting array.We will show in the following how the flux Φ ( J ) k ( t ) canbe calculated and how the junction currents I k ( t ) in Eq.(1) can be expressed in terms of the differences of gaugeinvariant phase differences, ϕ k +1 ( t ) − ϕ k ( t ), leading to asystem of coupled first-order non-linear ordinary differen-tial equations for the ϕ k ( t )’s, and an integro-differentialequation for the stream function g ( x, y ) (defined below),from which the time-averaged voltage of the array as afunction of the applied magnetic field, can be determined. B. Biot-Savart’s law, London’s equation and thestream function equation for a parallel SQUID arraywith wide tracks, busbars and leads
In order to calculate the magnetic flux Φ ( J ) k ( t ) of Eq.(4), we employ Biot-Savart’s law and the second Lon-don equation. The SQUID array we wish to model ismade out of a YBCO high-temperature superconduct-ing thin-film of thickness d = 0.125 µ m and a Londonpenetration depth (77 K) of approximately λ = 0.3 µ m.Because here λ > d , the supercurrent density throughthe thickness d of the film is nearly homogeneous, inde-pendent of the z direction, and therefore Biot-Savart’slaw in 2D can be applied . In this case, the magneticfield H ( J ) ( x, y, t ) in z direction, produced by the super-current density j ( x, y, t ) flowing in the array of domainΩ (see Fig. 1), is given by H ( J ) ( x, y, t ) = (5) d π (cid:90) Ω j x ( x (cid:48) , y (cid:48) , t ) ( y − y (cid:48) ) − j y ( x (cid:48) , y (cid:48) , t ) ( x − x (cid:48) ) (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) dx (cid:48) dy (cid:48) , where j x and j y are the x and y components of the super-current density j . One can express j x and j y in Eq. (5)in terms of the stream function g ( x, y, t ) which is definedas j x = 1 d ∂g∂y and j y = − d ∂g∂x . (6)This type of stream function approach has been usedpreviously by Khapaev and Brandt .To be allowed to integrate Eq. (5) in 2D by parts,it is required to smoothen the functions in the in-tegrand of Eq. (5) to generate continuously differ-entiable functions. We achieve this by analyticallyintegrating ( y − y (cid:48) ) / (cid:112) ( c − x (cid:48) ) + ( y − y (cid:48) ) and ( x − x (cid:48) ) / (cid:112) ( c − x (cid:48) ) + ( y − y (cid:48) ) over small intervals aroundtheir singularity points. This leads to H ( J ) ( x, y, t ) = f s ( x, y, t ) (7) − π (cid:90) Ω Q F ( x, y, x (cid:48) , y (cid:48) ) g ( x (cid:48) , y (cid:48) , t ) dx (cid:48) dy (cid:48) , where the kernel Q F ( x, y, x (cid:48) , y (cid:48) ) is defined in AppendixB, and f s ( x, y, t ) (8)= 14 π (cid:73) ∂ Ω g ( x (cid:48) , y (cid:48) , t ) (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) · n dl (cid:48) . Here dl (cid:48) is the integration line element and n is a 2Dnormal vector in the xy plane which is perpendicular onthe domain boundary ∂ Ω, pointing outwards, away fromthe area Ω. The contour ∂ Ω includes the hole-boundariesof the array.Exploiting the stream function symmetry about the x axis, i.e. g ( x (cid:48) , y (cid:48) , t ) = g ( x (cid:48) , − y (cid:48) , t ), one can restrict theintegration domain Ω in Eq. (7) and the contour integra-tion domain ∂ Ω in Eq. (8) to only the upper domains Ω u and ∂ Ω u ( y (cid:48) ≥ u ∪ Ω d , withthe superscript u referring to the upper ( y (cid:48) >
0) and d tothe lower ( y (cid:48) <
0) domain (see Fig. 2). Thus one findsfor Eq. (7) H ( J ) ( x, y, t ) = f s ( x, y, t ) (9) − π (cid:90) Ω u Q ( x, y, x (cid:48) , y (cid:48) ) g ( x (cid:48) , y (cid:48) , t ) dx (cid:48) dy (cid:48) , where the integration is now only over the upper domainΩ u with a kernel Q given by Q ( x, y, x (cid:48) , y (cid:48) ) = Q F ( x, y, x (cid:48) , y (cid:48) )+ Q F ( x, y, x (cid:48) , − y (cid:48) ) , (10)and f s ( x, y, t ) = f us ( x, y, t ) + f us ( x, − y, t ) . (11)Here f us is defined by Eq. (8) but with a contourintegration only over the upper contour ∂ Ω u of thedomain Ω u (see Fig. 2), where ∂ Ω u includes the upperboundary part of the holes. In Eq. (11), the secondterm on the right hand side, due to symmetry, accountsfor the lower domain part.In order to obtain an equation to calculate the streamfunction g ( x, y, t ) we use the second London equationwhich has the form λ ∇ × j = − H , (12)where H is the total magnetic field. Again, because λ >d , and thus j ( x, y, z ) = j ( x, y ), we can employ the 2Dsecond London equation which from Eq. (12) becomesΛ ∆ xy g ( x, y, t ) = H ( x, y, t ) . (13)Here in 2D, the magnetic field and the stream functionare governed by the 2D screening length or Pearl length Λ = λ /d . The operator ∆ xy is the 2D Laplace operator, ∂ /∂x + ∂ /∂y , while H ( x, y, t ) is the total magneticfield in z direction. Because H ( x, y, t ) = H a + H ( J ) ( x, y, t ) , (14)where H a is the applied magnetic field in z direction,one finds, using Eqs. (8), (9), (11), (13) and (14), a 2Dsecond-order linear Fredholm integro-differential equa-tion for the stream function g ( x, y, t ) of the formΛ ∆ xy g ( x, y, t ) + 14 π (cid:90) Ω u Q ( x, y, x (cid:48) , y (cid:48) ) g ( x (cid:48) , y (cid:48) , t ) dx (cid:48) dy (cid:48) (15)= H a + f us ( x, y, t ) + f us ( x, − y, t ) , with f us ( x, y, t ) = (16)14 π (cid:73) ∂ Ω u g ( x (cid:48) , y (cid:48) , t ) (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) · n dl (cid:48) . C. Array boundary condition for the streamfunction intergo-differential equation
In order to solve Eq. (15) one has to define bound-ary conditions for g ( x, y, t ) and ∆ xy g ( x, y, t ) along theboundary ∂ Ω u .Figure 2 shows the names given to different sectionsalong the ∂ Ω u boundary. As the bias-current I b is in-jected downwards into the top current lead (Fig. 1),and because j y d = − ∂g/∂x , we choose for the currentinjection boundary condition g ( x, y, t ) = I b x/ (2 c ) for( x, y ) (cid:15) ∂ Ω ( T ) . Because no current is crossing the bound-aries ∂ Ω ( L ) , ∂ Ω ( R ) and ∂ Ω k (inside holes), we find byusing Eq. (6) the boundary conditions g ( x, y, t ) = − I b / x, y ) (cid:15) ∂ Ω ( L ) and g ( x, y, t ) = I b / x, y ) (cid:15) ∂ Ω ( R ) ,while for the hole-boundary conditions along ( x, y ) (cid:15) ∂ Ω k ,one has g ( x, y, t ) = ˜ g k ( t ) where˜ g k ( t ) = k (cid:88) j =1 I j ( t ) − I b , (17)Thus, the junction currents I ( t ) and I N ( t ) at theends of the array are I ( t ) = ˜ g ( t ) + I b / I N ( t ) = − ˜ g N − ( t ) + I b /
2, while between holes, thejunction currents are I k ( t ) = ˜ g k ( t ) − ˜ g k − ( t ).With the above boundary conditions one can calculate f s ( x, y, t ) (Eqs. (11) and (16)) which can be written inthe form f s ( x, y, t ) = P ( x, y ) I b + N − (cid:88) k =1 P k ( x, y ) ˜ g k ( t ) . (18)where both P ( x, y ) and P k ( x, y ) are independent of time t and P ( x, y ) := P u ( x, y ) + P u ( x, − y ) (19)and P k ( x, y ) := P uk ( x, y ) + P uk ( x, − y ) . Here P u is the part of the integral in Eq. (16) alongthe contours ∂ Ω ( R ) ∪ ∂ Ω ( T ) ∪ ∂ Ω ( L ) , while P uk is thepart of the integral along the upper half ( y ≥
0) ofthe hole-contour ∂ Ω k (see Fig. 2). We have calculated P u ( x, y ) and P uk ( x, y ) analytically using Eq. (16) andthe results are shown in Appendix B.FIG. 2: The different boundary contour sections ∂ Ω ( L ) , ∂ Ω ( T ) , ∂ Ω ( R ) along the upper domain Ω u (shadeddarker) and the hole domains Ω k and contours ∂ Ω k for k = 1 , ..., N −
1, of the holes between domain Ω u andΩ d . D. Magnetic flux in array holes
The total magnetic flux Φ k ( t ) (Eq. (4)) in hole k ,required in Eq. (3), becomes with Eq. (14),Φ k ( t ) = µ (cid:90) Ω k ( H a + H ( J ) ( x, y, t ) ) dxdy , (20) where the integration is over the hole domain Ω k (seeFig. 2).Using Eqs. (9) and (18), one derivesΦ k ( t ) µ = H a A h + L k I b + N − (cid:88) j =1 L kj ˜ g j ( t ) (21) − π (cid:90) Ω u ˜ Q k ( x (cid:48) , y (cid:48) ) g ( x (cid:48) , y (cid:48) , t ) dx (cid:48) dy (cid:48) , where µ H a A h is the applied magnetic flux that pene-trates each hole of area A h , and L kj := (cid:90) Ω k P j ( x, y ) dxdy = 2 (cid:90) Ω k P uj ( x, y ) dxdy , (22)with j = 0 , , ...N −
1. Furthermore, in Eq. (21),˜ Q k ( x (cid:48) , y (cid:48) ) is defined as˜ Q k ( x (cid:48) , y (cid:48) ) := (cid:90) Ω k Q ( x, y, x (cid:48) , y (cid:48) ) dxdy , (23)where Q ( x, y, x (cid:48) , y (cid:48) ) is given by Eq. (10) with ( x (cid:48) , y (cid:48) ) (cid:15) Ω u . L kj in Eq. (22) as well as ˜ Q k ( x (cid:48) , y (cid:48) ) in Eq. (23) canbe calculated analytically. But here, calculations weresimply performed numerically for each hole k , sincethese quantities are solely of geometrical nature andtime-independent, and thus, have to be calculated onlyonce. E. Conversion to algebraic equations andvectorization
In order to calculate numerically the stream function g ( x, y, t ) with Eq. (15) for ( x, y ) (cid:15) Ω u , one has to convertEq. (15) into an algebraic equation. To do so, we choosea sufficiently fine square grid on Ω u (Fig. 2) and discretisethe spatial vectors r = ( x, y ) to r n = ( x n , y n ) locatedat the centre of each small square grid element of size w = (∆ x ) , where ∆ x is the square grid spacing and theindex n counts the grid elements, n = 1 , ..., N g , where N g is the total number of grid elements in the domain Ω u .The integro-differential equation Eq. (15) for g ( x, y, t )then becomes (cid:88) m (cid:15) Ω u (cid:104) Λ ∆ nm + w π Q nm (cid:105) g m ( t ) (24)= H a + P n I b + N − (cid:88) k =1 P nk ˜ g k ( t ) , for all n (cid:15) Ω u . Here g n ( t ) := g ( r n , t ) , Q nm := Q ( r n , r m ) , (25) P n := P ( r n ) and P nk := P k ( r n ) . In vector notation, the relationship in Eq. (17) be-tween the boundary values ˜ g k ( t ) and the junction cur-rents I k ( t ) and bias injection current I b becomes˜ g ( t ) = T ◦ ˜ I ( t ) − I b N − , (26)where ˜ g is the N − g = (˜ g , ..., ˜ g N − ) for the holes. In Eq. (26) the symbol ◦ means multiplication of a matrix with a vector (and alsolater matrix multiplication) and the above ( N − × ( N −
1) matrix T is defined as ( T ) kj = 1 if k ≥ j and zerootherwise, and the junction current vector ˜ I ( t ) is an N − N ) vector, ˜ I ( t ) := ( I ( t ) , ..., I N − ( t )). The vector N − is of dimension N − N − := (1 , ..., g ( t ) as g ( t ) = (cid:16) Λ D + w π Q (cid:17) − ◦ (cid:20) H a N g + P I b + P ◦ T ◦ ˜ I ( t ) + [ P − P ◦ N − ] I b − Λ d ∆ ( t ) (cid:21) . (27)The stream function vector, g ( t ) := ( g ( t ) , ..., g N g ( t )),represents the stream function g ( x, y, t ) at all grid pointelements in Ω u . The matrix D is an N g × N g matrixcorresponsing to the Laplace operator in Eq. (24),where ( D ) nm := [ − δ nm + δ Nbnm ] /w with w = (∆ x ) and δ Nbn,m = 1 if r m is a nearest neighbour of r n andzero otherwise. The matrix Q is also an N g × N g matrix, defined as ( Q ) nm := Q ( r n , r m ). The symbol N g in Eq. (27) is the vector N g = (1 , , ...,
1) ofdimension N g and P := ( P , ..., P N g ), defined in Eq.(19) and Appendix B. The symbol P is an N g × ( N − P ) nk := P k ( r n ). The components of thetime-dependent N g -dimensional vector d ∆ ( t ) in Eq. (27)originate from the part of the Laplace operator whichoperates on the domain boundary ∂ Ω u (which includesthe holes) and the boundary along junctions. Most ofthe components of d ∆ ( t ) are zero, but for grid elementsadjacent to boundaries, the components are 8 / g ( L ) /w ,8 / g ( T ) /w , 8 / g ( R ) /w , and 8 / g k ( t ) /w (along ∂ Ω k )(Fig. 2). For grid elements adjacent to junctions, thecorresponding components of d ∆ ( t ) vary linear with dis-tance along the junctions. The time dependence of d ∆ ( t )therefore arises from the time-dependent hole-boundaryconditions (Eq. (17)). The above factor of 8 / g ( t ) usingEq. (27), a very large N g × N g matrix, Λ D + w π Q , hasto be inverted. It is important to note that since thismatrix is time-independent, this large matrix inversionhas to be performed only once at the beginning of acomputation.In vector notation, using Eqs. (21) and (26), the fluxvector Φ ( t ) = (Φ ( t ) , ..., Φ N − ( t )) of the magnetic fluxinside array holes, takes the form Φ ( t ) µ = H a A h N − + ( L − L ◦ N − ) I b (28) + L ◦ T ◦ ˜ I ( t ) − w π ˜ Q ◦ g ( t ) , where L := ( L , ..., L N − , ) with L k given by Eq.(22), and L is an ( N − × ( N −
1) matrix defined as( L ) kj := L kj ( j ≥
1) where L kj is again given by Eq.(22). In Eq. (28) ˜ Q is an ( N − × N g matrix givenby ( ˜ Q ) kj := ˜ Q k ( r j ) where ˜ Q k ( r j ) is defined by Eq. (23) .Furthermore, we rewrite the Ginzburg-Landau equa-tion Eq. (3) in vector notation of the form N ◦ ϕ ( t ) = 2 π Φ ( µ Λ K ( t ) + Φ ( t ) ) , (29)where N is an ( N − × N matrix defined as ( N ) kj = − k = j and ( N ) kj = 1 for k = j − ϕ ( t ) is N dimen-sional, ϕ ( t ) := ( ϕ ( t ) , ..., ϕ N ( t )). Using Eqs. (3) and(6), the components K k ( t ) of the ( N −
1) dimensionalvector K ( t ) := ( K ( t ) , ..., K N − ( t )) are K k ( t ) = d (cid:73) ∂ Ω k j ( t ) · d l = (cid:73) ∂ Ω k ( ∂g ( t ) ∂y dx − ∂g ( t ) ∂x dy ) , (30)where the contour integration is counterclockwise alongthe boundary ∂ Ω k of the hole k (see Figs. 1 or 2). F. System of coupled differential equations for theJosephson phase differences for a parallel SQUIDarray with wide thin-film structure
Using Eq. (28) and (29) and eliminating the flux vector Φ ( t ), one derives an equation for the junction currentvector ˜ I ( t ) as a function of the phase difference vector ϕ ( t ) and the stream function vector g ( t ) of the form˜ I ( t ) = ( L ◦ T ) − (cid:20) Φ π µ N ◦ ϕ ( t ) − Λ K ( t ) − H a A h N − − (cid:18) L − L ◦ N − (cid:19) I b + w π ˜ Q ◦ g ( t ) (cid:21) . (31)Note that ˜ I ( t ) is an ( N −
1) vector and its componentsdo not contain the junction current I N ( t ) across the lastJJ. But, since I N ( t ) = I b − (cid:80) N − k =1 I k ( t ), the current I N ( t ) is well defined and thus, Eq. (31), togetherwith Eqs. (1) and (2) define a complete set of coupledfirst-order differential equations for the gauge invariantphase differences ϕ k ( t ) with k = 1 , ..., N .It is convenient to define I c as the average junctioncritical current, I c = (cid:80) Nk =1 I ck / N and R as the averagejunction resistance, R = (cid:80) Nk =1 R k / N . Note that R isnot the total array resistance. Then, Eq. (1) combinedwith Eq. (2) can be put into the form d ϕ k ( τ ) dτ = ξ k (cid:18) − η k sin ϕ k ( τ ) + I k ( τ ) I c + I Noisec ( τ ) I c (cid:19) , (32) where τ is the reduced time in dimensionless units, τ = 2 π Φ R I c t , (33)and ξ k = R k /R and η k = I ck /I c , and thus ξ k ( η k ) isa measure of deviation of R k ( I ck ) from R ( I c ). Thestandard deviations for R k and I ck in YBCO thin-filmSQUID arrays can be as large as 0 . .At this point it is convenient to write Eq. (32) in vectornotation and combining it with Eq. (31) which results in ˙ ϕ ( τ ) = S ( τ ) + (cid:98) ξ ◦ ( L ◦ T ) − (34) ◦ (cid:20) Φ π µ N ◦ ϕ ( τ ) − Λ K ( τ ) + w π ˜ Q ◦ g ( τ ) − (cid:18) L − L ◦ N − (cid:19) I b − H a A h N − (cid:21) + (cid:98) ξ ◦ ˜ I Noise ( τ ) . Here the components of the vector ˙ ϕ ( τ ) are dϕ k ( τ ) /dτ and the components of the vector S ( τ )are − ξ k η k sin ϕ k ( τ ), where k = 1 , ..., N −
1. Thematrix (cid:98) ξ is an ( N − × ( N −
1) diagonal matrix withdiagonal elements ξ k /I c , where k = 1 , ..., N −
1, and˜ I Noise ( τ ) := ( I Noise ( τ ) , ..., I NoiseN − ( τ )). Please note that K ( τ ) and g ( τ ) are both functions of time τ since the supercurrent density distribution in a SQUID arrayundergoes periodic changes with time.To make the set of coupled differential equations for ϕ k ( τ ) of Eq. (34) complete, one has to add an equa-tion for dϕ N ( τ ) /dτ using Eqs. (1) and (2), which gives,because of I b = (cid:80) Nk =1 I k ( τ ), dϕ N ( τ ) dτ = ξ N (cid:32) − η N sin ϕ N ( τ ) + (cid:34) I b − N − (cid:88) k =1 I k ( τ ) (cid:35) / I c + I NoiseN ( τ ) / I c (cid:33) . (35)Equation (34) together with Eq. (35) forms a completeset of coupled first-order differential equations for allthe ϕ k ( τ )’s. To solve this set of differential equations,initial conditions for all ϕ k ( τ = 0) have to be chosen.This is done by starting with equal junction currents I k ( τ = 0) = I b /N and using Eqs. (26) - (28) and (30) to calculate the ϕ k ( τ = 0)’s from N ◦ ϕ ( τ ) of Eq. (29),and setting ϕ ( τ = 0) = 0.Please note that Eqs. (34) and (35) together with Eq.(26), (27) and (31) are the key equations of this paper.The set of equations Eqs.(34) and (35) were solved nu-merically using the 4 th order Runge-Kutta method. Aftereach time step, chosen as ∆ τ = 0 .
1, the ϕ k ’s change andthus the Josephson currents I k ( τ ) change slightly (Eq.(31)), which then slightly changes the boundary condi-tion ˜ g ( τ ) (Eq. (26)). Thus, after each time step ∆ τ anupdated stream function g ( τ + ∆ τ ) (Eq. (27)) has to becalculated, resulting in updated g and K vectors in Eq.(34).Details about how the noise currents ˜ I Noise ( τ ) in Eq.(34) and I NoiseN ( τ ) in Eq. (35) were treated numericallyare outlined in Appendix D. G. Time-averaged voltage
The time-averaged normalised voltage,
V / ( RI c ), be-tween the leads of a parallel SQUID array, is given bytime averaging the right-hand side of the Josephson equa-tion Eq. (2) which results in VR I c = lim τ →∞ τ [ ϕ k ( τ + τ ) − ϕ k ( τ ) ] , (36)where k can be any k (cid:15) { , ..., N } . One has to choose τ large enough such that numerical self-adjustmentfor the initial ϕ k ’s has occurred and τ has to be takensufficiently large.In the case where Johnson noise can be neglected, onefinds VR I c = 2 πτ p , (37)where τ p is the period of oscillations of the ϕ k ( τ )’s. In thecase of non-negligible Johnson noise, Eq. (37) cannot beused but one can reduce the statistical error in V / ( RI c )by averaging over all the N phase differences, i.e. VRI c (cid:39) τ N N (cid:88) k =1 [ ϕ k ( τ + τ ) − ϕ k ( τ ) ] , (38)where τ has to be chosen sufficiently large. H. Effective areas of SQUID array holes
The wide tracks, busbars and leads focus magnetic fluxinto the array holes. In addition, the Meissner shieldingcurrent crowding near the holes, enhance the (cid:72) j · dl term.An effective area, A effk can be defined for each array hole k via the fluxoid it contains as A effk = lim I c ,I b → µ λ (cid:72) ∂ Ω k j · d l + Φ k B a . (39) A effk can also be defined in the limit of B a → ∞ insteadof I c and I b →
0. From the definition in Eq. (39), em-ploying Eqs. (28) and (30), it follows that A effk in vector notation, A eff , normalised to the hole area A h , is A eff A h = N − + µ B a A h (cid:16) Λ K − w π ˜ Q ◦ g (cid:17) , (40)where the subscript 0 in K and g means that g in Eq.(27) and K in Eq. (30) are evaluated with the boundarycondition g ( x, y ) = 0 along ∂ Ω u , since I c = I b = 0. III. DEVICE FABRICATION ANDEXPERIMENT
Parallel SQUID arrays were fabricated lithographi-cally by growing thin-films of YBCO on 1 cm MgOsubstrates. Steps were etched into the MgO surfaceusing a well established technique based on argonmilling . During YBCO thin-film growth by e-beamevaporation, a long grain boundary forms at the top ofthe edge of the MgO step, creating a long JJ. Films werethen lithographically patterned into parallel SQUIDarrays . The width of junctions and junction tracksis w J = 2 µ m and the width of holes w h = 4 µ m withhalf-height h = 4 µ m (see Fig. 1). For our N = 11 array,which is the array discussed in detail in this paper, thethickness of the thin-film is d = 0.125 µ m, the width ofthe busbar b − h = 8 µ m and the bias-lead half-width c = 4 µ m. For the calculation the length of the bias-leadwas chosen sufficiently long as l = 24 µ m (see Fig. 1).A micrograph of the N = 11 device is shown as an insetin Fig. 3(a).To measure the V ( B a ) response, the array was placedon a measurement probe, that generated a perpendicularapplied induction B a (magnetic field H a ), and thendipped into a dewar of liquid nitrogen and zero fieldcooled down to a temperature of 77 K. To screen outthe earth’s magnetic field, the dewar was surroundedby five layers of mu-metal shielding. A bias-current I b ,with a value which optimised the SQUID array response(maximised transfer function dV /dB a ), was injected intothe top bias-lead and the voltage V at different B a wasmeasured using the standard four-terminal method. IV. RESULTS AND DISCUSSION
In the following we discuss in detail the calculatedand experimental results for a parallel SQUID arraywith N = 11 junctions with 8 µ m wide fluxoid focussingbusbars and bias-leads. A micrograph of this device isshown as an inset in Fig. 3(a). The parameters usedin the calculation are the array bias-current I b = 200 µ A (given by experiment), the average critical currentof a junction I c = 24 µ A and the London penetrationdepth λ = 0.33 µ m. Further below we will discuss howthe values for I c and λ were determined. We foundthat choosing the square grid spacing as ∆ x = 0.5 µ m,or even ∆ x = 1 µ m, was numerically sufficient andtherefore the number of grid points N g that lie in thearray domain Ω u is N g = 3104, or 776. P ha s e s j k (r ad ) Time t N = 11B a = 0 k = 1 - 11N = 11 (a) m m k = 11 P ha s e s j k (r ad ) Time t k = 1N = 11B a = 12 mT(b) FIG. 3: The inset in (a) shows a micrograph of ourparallel SQUID array with N = 11 junctions where thehorizontal line indicates the step edge across which JJshave formed. (a) displays the time-evolution of thephase differences ϕ k for k = 1 to 11, calculated usingEqs. (34) and (35) together with Eq. (27) for B a = 0.(b) displays the time-evolution of the phase differences ϕ k ( τ ) for B a = 12 µ T.Figures 3(a) and (b) show the calculated time-evolution of the N = 11 phase differences ϕ k ( τ ) (fromEqs. (34) and (35)) for perpendicular applied magneticinductions B a = 0 and B a = 12 µ T, repectively. For - - - - - - y ( m m ) g ( x , y ) ( m A ) x ( m m ) (a) B a = 0 - - - - - - y ( m m ) g ( x , y ) ( m A ) x ( m m ) (b) B a = 12 m T FIG. 4: Stream function g ( x, y ) in the upper part( y >
0) of the N = 11 parallel SQUID array at τ =9000 for (a) B a = 0 and (b) B a = 12 µ T. B a = 0, the fluxoid in each hole is close to zero, whileat B a = 12 µ T, the fluxoid in each hole is about half aflux quantum, Φ /
2. Due to thermally activated phaseslippages, caused by an effective Johnson noise strengthΓ = 0.135 (Eq. (D1)) at T = 77K, the time-evolutionof the ϕ k ( τ ) in Figs. 3(a) and (b) is somewhat erratic.For B a = 0 coherent phase slippages by 2 π occur quitesuddenly, while in the B a = 12 µ T case, 2 π incrementsappear more sinusoidal. Furthermore, in the B a = 12 µ T case the ϕ k ( τ ) are incrementally shifted upwards byabout π with increasing k .0Figures 4(a) and (b), using Eq. (27), show the streamfunction g ( x, y ) for the upper domain Ω u at time τ =9000 (Fig. 3) for B a = 0 and B a = 12 µ T, respectively.As can be seen, the stream function values along theleft and right boundaries, ∂ Ω ( L ) and ∂ Ω ( R ) , are ∓ I b / ∓ µ A, while the step-like structure of g ( x, y )along the 10 holes corresponds to the stream functionboundary values ˜ g k in the holes which vary with time τ . -30 -20 -10 0 10 20 305101520 y ( m m ) B a = 0(a) -30 -20 -10 0 10 20 305101520 y ( m m ) x ( m m) B a = 12 m T(b)
FIG. 5: Current stream lines in the upper domain Ω u ofthe N = 11 parallel SQUID array at time τ = 9000 for(a) B a = 0 and (b) B a = 12 µ T.Figures 5(a) and (b) display the current stream linesat B a = 0 and B a = 12 µ T, obtained from the contourlines of g ( x, y ) of Figs. 4(a) and (b). Only the upperdomain Ω u is shown because of the symmetry aboutthe x axis where g ( x, y ) = g ( x, − y ) and thus from Eq.(6), j x ( x, y ) = − j x ( x, − y ) and j y ( x, y ) = j y ( x, − y ).As can be seen, for B a = 0, the current fans outfrom the bias-lead to the junctions where the currentsthrough individual junctions vary with time. Currentcrowding is visible left and right of the bias-lead and inparticular at the corners between bias-lead and busbardue to Meissner shielding. For B a = 12 µ T, circulatingMeissner shielding currents are visible in the left part ofthe busbar. Strong current crowding now only occurson the right side of the bias-lead and the right cornerbetween bias-lead and busbar. In addition, strongcurrent crowding is visible at the top of holes due to theMeissner shielding current circulating in the busbar.Figures 6(a) and (b) show the calculated perpendic-ular total magnetic field H ( x, y ), Eq. (9) and (14),inside and outside of the N = 11 parallel SQUID array( y ≥
0) for B a = 0 and B a = 12 µ T at time τ = 9000.Strong magnetic field enhancements are visible alongthe edges of the bias-lead and the upper edge of thebusbar, with a particularly strong field at the corners - - - - - - y ( m m ) m H ( x , y ) ( m T ) x ( m m ) (a) B a = 0 - - - y ( m m ) m H ( x , y ) ( m T ) x ( m m ) B a = 12 m T (b) FIG. 6: Perpendicular total magnetic field H ( x, y ) inthe upper half ( y >
0) of the N = 11 parallel SQUIDarray at time τ = 9000 for (a) B a = 0 and (b) B a = 12 µ T.between bias-lead and busbar due to strong currentcrowding. The magnetic fields around and inside of theholes look complicated and change with time as thejunction currents oscillate with time. At B a = 12 µ T theapplied magnetic induction adds to the total field andthe Meissner shielding currents induced in the busbarsand bias-leads generate additional magnetic fields alongedges. The lowest magnetic field, in the centre region ofthe busbar (not visible here), is about 3.5 µ T /µ andthus the maximum busbar-shielding is about 70%.1 -200 -150 -100 -50 0 50 100 150 2000.00.10.20.30.40.50.6 V / R I c B a (mT) l = 0.33 m mI c = 24 m AT = 77 KN = 11I b = 200 m ACalc. (a) -200 -150 -100 -50 0 50 100 150 2000.000.020.040.060.080.10 V o l t age V ( m V ) B a ( m T) I b = 200 m AN = 11Exp.T = 77 K (b)
FIG. 7: Responses V versus B a at T = 77 K for (a)calculation and (b) experiment.Figure 7(a), which is the most important result ofour paper, shows the calculated time-averaged voltage V versus the applied perpendicular magnetic induction B a , from -200 µ T to +200 µ T, for the N = 11 parallelSQUID array at a bias current I b = 200 µ A, assumingzero I ck , R k spreads. The calculation includes Johnsonnoise at T = 77 K. The average junction critical currentdensity I c and the London penetration depth λ , neededas input parameters in our calculation, were chosen togive the best agreement with our experimental datadisplayed in Fig. 7(b). These parameters are I c = 24 µ Aand λ = 0.33 µ m. The temperature dependence of theLondon penetration depth λ ( T ) for YBCO thin-films hasbeen widely investigated , and Chen et al. found for high quality YBCO thin-films, using AC-susceptibilitymeasurements, λ (77 K ) ≈ . µ m, which is similar toour value. As can be seen, our calculation in Fig. 7(a)agrees very well with our experimental data shown inFig. 7(b). The calculation reproduces accurately theexperimental ratio of maximal to minimal voltage as wellas the overall experimental envelope modulation. Also,the dips appear at the correct B a values. In addition,the shoulder-peak that initially appears near the secondside minima and then propagates outwards with increas-ing B a is closely reproduced. In the experiment thebias-current I b = 200 µ A was chosen to maximise thetransfer function dV /dB a around the centre dip. Fromthe I c value it follows that I b = 0 . N I c . In caseswhere I b < N I c , it is important to take the effects ofJohnson noise in the calculation fully into account, as wehave done here. Comparing the experimental maximumvoltage with the calculated maximum V / ( RI c ) valuein Fig. 7(a), one obtains an average single junctionresistance R = 6.2 Ω. Some of the tiny spikes around theupper parts of the calculated curve are due to numericalinaccuracies of a finite temperature calculation due to alimitation in available computation time. A e ff k / A h Hole kN = 11l = 0.33 mm
FIG. 8: Effective area enhancement factor A effk /A h ofEq. (40) versus the hole index k for N = 11. Theaverage enhancement factor is 2.72. The lower part ofthe bars is the µ λ (cid:72) j · d l fluxoid contribution to theeffective area (Eq. (39)) while the upper part is the Φ k fluxoid contribution.Our calculation reveals that the appearance of a broadenvelope modulation in the V versus B a response shownin Figs. 7(a) and (b) is due to inhomogeneous fluxoidfocussing where the effective areas A effk of the two holes( k = 5 and k = 6) closest to the centre are larger than2the effective areas of the holes at the ends. The valuesof A effk , with k = 1 , ...,
11, calculated from Eq. (40),are displayed in Fig. 8. The lower part of the bars isthe µ λ (cid:72) j · d l fluxoid contribution to the effectivearea (Eq. (39)) while the upper part is the Φ k fluxoidcontribution (which contains the applied flux). Theaverage effective area enhancement is A effav /A h = 2.72.Using this enhancement factor, one obtains for the firstside minimum position B a, := Φ o /A effav = 23.78 µ T,in close agreement with both Figs. 7(a) and (b). Theeffective area difference of 17% between the end andthe centre holes is responsible for the strong envelopemodulation seen in Figs. 7(a) and (b). The complicatedinterference pattern seen in Fig. 7(a) is very sensitive tothe actual form of the effective area A effk distribution(Fig. 8). In a simple lumped-element simulation an en-velope modulation can also be produced by varying thegeometric area sizes of the holes . The effect of varyingthe geometrical areas in parallel SQUID arrays on the V versus B a response has been investigated in detailby Oppenl¨ander et al. who simulated the behaviourof SQIF’s using a lumped-element approach. It is im-portant to note that these lumped-element simulationscannot properly account for the effect of flux focussingand Meissner shielding-currents flowing in wide tracks,busbars and leads, and in particular are not suitable tocalculate how the injected current from a wide bias-leadfans out into the junctions currents. Lumped-elementsimulations are thus unable to accurately describe the V versus B a response of a parallel SQUID array with widesuperconducting thin-film busbars and leads. without Johnsonnoise V / R I c B a (mT)Calc. N = 11I b = 200 mAl = 0.33 mmI c = 24 mAwith Johnson noise FIG. 9: Comparison of calculated response V versus B a with Johnson noise (dashed green curve, Fig. 7(a)) andwithout Johnson noise (solid red curve).Figure 9 demonstrates the importance of taking John- son noise of the junction resistors into account, whencalculating the V versus B a response. The solid redcurve in Fig. 9 was obtained without Johnson noise andshows sharper dips and narrow regions of zero voltagein contrast to the dashed green curve which is identicalto Fig. 7(a) which includes Johnson noise. This kind ofresponse difference is well known from early simulationsfor a SQUID with N = 2 . We found that numericalcalculations for the N = 11 parallel SQUID array thatinclude Johnson noise are computationally about 30times more demanding than calculations without John-son noise. At T = 77 K, particularly at small averagevoltages V , one has to calculate the time-evolution ofthe junction phase differences ϕ k ( τ ) over a very longtime period τ in order to obtain a statistically accuratetime-averaged voltage. The accuracy of calculationsshown in Fig. 7(a) is about 2%. -200 -150 -100 -50 0 50 100 150 2000.00.10.20.30.40.50.6 s Ic = 0.2 V / R I c B a (mT)Calc. N = 11 I b = 200 mAl = 0.33 mmI c = 24 mAT = 77 K FIG. 10: Effect of I ck , R k spreads with σ Ic = 0.2 (fullred curve) on V versus B a . The green dotted curve isfor no spread and is identical to Fig. 7(a). The insetshows the Gaussian random I ck /I c and R k /R used.Up to this point our calculations have not consideredspreads in critical junction currents I ck and junctionresistances R k . The effect of I ck , R k spreads with astandard deviation of σ Ic = 0.2 is shown in Fig. 10.Here we assume that I k and R k are anti-correlatedaccording to the empirical law R k I ck ∝ J / c where J c is the junction critical current density which is assumedto be a constant. The Gaussian random I ck /I c and R k /R values that were used are shown in the inset ofFig. 10. As can be seen, compared to the case of nospread (dotted green curve; identical to Fig. 7(a)), thespread (full red curve) breaks the symmetry about the V axis and V ( B a ) (cid:54) = V ( − B a ). The envelope modulation3is only mildly affected and the previously mentionedshoulder peaks stay at their positions. In contrast,the experimental data in Fig. 7(b) shows symmetricbehaviour which is due to our experimental procedurewhich averages the voltages V for direct and reversedbias-current I b , in order to cancel any apparatus voltageoffset. Thus the measurement procedure symmetrisedthe V vs B a curve in Fig. 7(b). The fact that the calcu-lation in Fig. 7(a) agrees so well with the symmetrisedexperimental data in Fig. 7(b) indicates that the I ck , R k spreads of our experimental N = 11 parallel SQUIDarray must be smaller than 0.2. -0.04 -0.02 0.00 0.02 0.04-0.2-0.10.00.10.2 R = 6.5 W T = 77 K l = 0.33 m mI c = 24 m A B i a s c u rr en t I b ( m A ) Voltage V (mV)N = 11B a = 0 Exp.Calc.
FIG. 11: Comparison between the experimental andcalculated bias-current I b versus the time-averagedvoltage V across the N = 11 parallel SQUID array for B a =0. In the calculation the resistance of each JJ isassumed to be R = 6.5 Ω.The time-averaged voltage V at B a = 0 in Figs. 7(a)and (b) is dependent on the bias current I b . Figure11 displays the calculated and experimental I b versus V ( B a = 0) dependence. The calculation included John-son noise at T = 77 K and used the same parametersas in Fig. 7(a), i.e. I c = 24 µ A and λ = 0.33 µ m,again assuming zero I ck , R k spreads. Figure 11 showsthat a junction resistance R = 6.5 Ω fits quite well theexperimental data and the thermal rounding seen in theexperimental I b versus V curve is well reproduced. TheJJ resistance value R = 6.5 Ω is close to R = 6.2 Ωwhich was extracted above from Figs. 7(a) and (b). Asseen in Fig. 11, the experimental curve has a voltageoffset error of 5 µ V found to be caused by a systematicerror from our measuring apparatus.From the supercurrents that are flowing parallel andvery close to the junctions we can calculate the fluxoidsin the junctions. We find that for the largest magnetic induction B a = 190 µ T, the value for the fluxoids injunctions is about 0.06 Φ . Thus, the critical current ofJJs is not affected by the applied magnetic field whichis sufficiently small so that our initial assumption of anearly constant critical current density across junctionareas is well justified. -200 -150 -100 -50 0 50 100 150 2000.00.10.20.30.40.50.60.70.8 V / R I c B a (mT)Calc. N = 11 b L = 0.59I c = 24 mAT = 77 K I b = 0.758 N I c l = 0.33 m m FIG. 12: V versus B a using the lumped-element modelwith Johnson noise. The dashed red curve correspondsto homogeneous bias-current injection and the full greenline represents central bias-current injection. TheSQUID-loop screening parameter used in thecalculation is β L = 0.59 .It is interesting to compare the results from ourcomprehensive model with the simpler lumped-elementmodel . In the lumped-element approach one alsosolves a system of coupled differential equations forthe junction phase differences, but these equationscontain partial inductances for all the tracks alongwhich currents are flowing. To obtain the set of coupledequations for the phase differences one has to utiliseKirchhoff’s law for all current vertices. Thus, in alumped-element model, currents fanning out from thebias-current injection lead into wide busbars cannot bemodelled properly. Therefore, in the lumped-elementapproach, one is often faced with the problem of howto best simulate the injection of the bias-current . Itis clear that the lumped-element model only works wellif the device is made up of very narrow tracks. Figure12 shows the calculated result using the lumped-elementmodel with Johnson noise for an N = 11 parallelSQUID array assuming very narrow tracks. In order toallow a fair comparison with our more elaborate model,we artificially included a constant fluxoid focussingfactor of 2.72 (Fig. 8) and chose for the screeningparameter β L = 2 I c L s / Φ = 0.59 which was obtained4by extracting the average self-inductance L s per arrayhole from our comprehensive model. The average L s was determined by forcing a loop current to flow aroundindividual holes and then evaluating the correspondinghole fluxoid. This showed that the self-inductances ofthe outside loops are about 2.3% larger than the innerones. The dotted red curve in Fig. 12 shows the resultfor homogeneous bias-current injection where all thetop array current vertices receive the same bias-current I b /N . In contrast, the solid green curve in Fig. 12 showsthe calculated result for central injection, where the fullbias-current I b is injected into the top centre currentvertex, i.e. above the central ( k = 6) junction. It canbe seen that the lumped-element model very stronglydepends on the bias-current injection scheme chosenand completely fails to reproduce our experimentaldata displayed in Fig. 7(b). Our comprehensive model,due to its greater details, takes a factor of about 10 more computation time than the computationally fastlumped-element model.We have also fabricated and measured 8 other parallelSQUID arrays with numbers of JJs ranging from N = 4to N = 81 and we intend to compare these experimen-tal data with our comprehensive model in an upcomingpaper. V. SUMMARY AND CONCLUSION
In this paper we have presented a comprehensivetheoretical model that allows us to calculate with highaccuracy the magnetic field response of an N = 11parallel SQUID array with wide thin-film geometricstructures, operated in the voltage state. The modelcalculates the fluxoids for each SQUID array holeduring the time-evolution of the phase differences ofthe JJs. This was achieved by solving numericallythe second-order linear Fredholm integro-differentialequation for the stream function, derived from thesecond London equation and Biot-Savart’s law, withboundary conditions that are updated for each timestep. The fact that the London penetration depth ofYBCO thin-films at 77 K is greater than the thicknessof our thin-film, allows us to solve the intero-differentialequation in 2D. The Josephson equations and the secondGinzburg-Landau equation for the phase differenceslead to a system of coupled first-order nonlinear differ-ential equations which depend on the stream functionwhich describes the time-varying supercurrent densitywithin the thin-film array structure. The equationsalso take into account the Johnson noise from the JJs.Compared to the much simpler and far less accuratelumped-element model approach, our comprehensivemodel, while computationally much more demanding,leads to highly accurate predictions.We have tested the predictive power of our model by comparing our model results with our experimental datafor an N = 11 parallel SQUID array with wide thin-filmstructures. The theoretical model requires only two pa-rameters, i.e. the junction critical current density I c andthe London penetration depth λ . Both parameters wereadjusted to give the best allover agreement with the ex-perimental V ( B a ) curve of the array. The model predictswith unprecedentedly high accuracy the V ( B a ) curve ofthe array over a wide applied magnetic field range andalso describes well the experimentally observed thermalrounding of the I b versus V ( B a = 0) curve. The modelreveals that the observed envelope modulation of the ex-perimental V ( B a ) curve is due to the non-equal effec-tive hole areas of the array, where the centre holes havea larger effective area. Spreads in I ck and R k lead to B a -asymmetry of the V ( B a ) curve. For wide thin-filmgeometric structures, the lumped-element model fails tocorrectly predict the V ( B a ) array response, because itcannot describe fluxoid focussing and in particular failsto appropriately handle the bias-current injection. Acknowledgments
The authors would like to thank Colin Pegrum, J¨ornBeyer, Marc Gali Labarias, Philip Fairman, Chris Lewisand John Clarke for stimulating discussions.
VI. APPENDIXESAppendix A: Kernel Q F The integral Eq. (5) is a convergent improper integralwhere the integrand is singular for ( x, y ) = ( x (cid:48) , y (cid:48) ). In or-der to apply the method of integration by parts, one cansmoothen the functions ( y − y (cid:48) ) / (cid:112) ( c − x (cid:48) ) + ( y − y (cid:48) ) and ( x − x (cid:48) ) / (cid:112) ( c − x (cid:48) ) + ( y − y (cid:48) ) in Eq. (5) by ana-lytically integrating them over the small area of a squaregrid element (∆ x ) size, so that these functions becomecontinuously differentiable functions. This leads in Eq.(7) to a kernel of the form Q F ( x, y, x (cid:48) , y (cid:48) ) = (A1)1(∆ x ) (cid:34) (cid:104) (cid:112) ¯ x + ¯ y ¯ x ¯ y (cid:105) x (cid:48) − x +∆ x/ x (cid:48) − x − ∆ x/ (¯ x ) (cid:35) y (cid:48) − y +∆ x/ y (cid:48) − y − ∆ x/ (¯ y )with the definition (cid:34) (cid:104) f ( x, y ) (cid:105) s s ( x ) (cid:35) s s ( y ) = (A2) f ( s , s ) − f ( s , s ) − f ( s , s ) + f ( s , s ) . inorder to avoid the singularity in Eq. (5). Brandt usesthe fact that 1 / (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) in Eq. (5) canbe interpreted as 4 π times the magnetic field in the xy plane of a point dipole of unit strength, positioned at( x (cid:48) , y (cid:48) ) and oriented in z direction. Since the magneticflux through the infinite xy plane is zero, this leads to anadditional equation that Brandt uses to eliminate anyunphysical singularity. We found that this method is lessaccurate than our method described above. Appendix B: Analytical expressions for thefunctions P u ( x, y ) and P uk ( x, y ) Using Eqs. (18) and (19), the function P u ( x, y ) is ob-tained by integrating the line integral in Eq. (16) alongthe contours ∂ Ω L , ∂ Ω R and ∂ Ω T (see Fig. 2), whichresults in P u ( x, y ) = 18 π [ α L ( x, y ) + α R ( x, y ) + α T ( x, y ) ] , (B1)where, α L ( x, y ) = 1 x + c ˜ y (cid:112) ( x + c ) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) b + l − y ˜ y = b − y (B2) − y − b ˜ x (cid:112) ˜ x + ( y − b ) (cid:12)(cid:12)(cid:12)(cid:12) − c − x ˜ x = − a − x + 1 x + a ˜ y (cid:112) ( x + a ) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) b − y ˜ y = − y , with c , b , l and a defined in Fig.1, and α R ( x, y ) = 1 x − c ˜ y (cid:112) ( x − c ) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) b + l − y ˜ y = b − y (B3)+ 1 y − b ˜ x (cid:112) ˜ x + ( y − b ) (cid:12)(cid:12)(cid:12)(cid:12) a − x ˜ x = c − x + 1 x − a ˜ y (cid:112) ( x − a ) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) b − y ˜ y = − y , and α T ( x, y ) = 1 c (cid:20) x ˜ xy − ( b + l ) − ( y − ( b + l )) (cid:21) (B4) 1 (cid:112) ˜ x + ( y − ( b + l )) (cid:12)(cid:12)(cid:12)(cid:12) c − x ˜ x = − c − x . The function P uk ( x, y ), where k = 1 , ..., N −
1, is ob-tained by integrating the line integral in Eq. (16) alongthe contour ∂ Ω k (Fig. 2) for y (cid:48) ≥
0, which results in4 π P uk ( x, y ) = − ˜ x ( y − h ) (cid:112) ˜ x + ( y − h ) (cid:12)(cid:12)(cid:12)(cid:12) x k + w h − x ˜ x = x k − x (B5)+ ˜ y ( x − x k ) (cid:112) ( x − x k ) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) h − y ˜ y = − y − ˜ y ( x − ( x k + w h )) (cid:112) ( x − ( x k + w h )) + ˜ y (cid:12)(cid:12)(cid:12)(cid:12) h − y ˜ y = − y , where x k := k ( w J + w h ) − w h − a . Appendix C: Laplacian
The superconducting thin-film 2D domain Ω u (Fig. 2)is divided into small square aerial elements w = (∆ x ) where the square grid spacing ∆ x was chosen as 0 . µm or 1 . µm . The grid points r n , with n = 1 , ..., N g , lie inthe centre of these grid elements.In Eq. (24), if a point r n lies more than a distance∆ x/ ∂ Ω u (which includes theboundaries ∂ Ω k ) or from the junction boundary, thenthe Laplacian on the square grid operates such that (cid:80) m ∆ nm g m = [ g ( x n + ∆ x, y n ) + g ( x n , y n + ∆ x ) − g ( x n , y n ) + g ( x n − ∆ x, y n ) + g ( x n , y n − ∆ x ) ] /w wherethe sum over m includes the four nearest neighbour sitesof r n .In Eq. (24), if a point r n lies only a distance ∆ x/ ∂ Ω u or a junction, one has to use theLaplacian for a non-equidistant grid. For example if r n is on the right side of a boundary line that runs alongthe y direction, then6 ∂ g/∂x (cid:12)(cid:12)(cid:12) ( x n ,y n ) ≈ g ( x n − h , y n ) 2 /h h + h − g ( x n , y n ) 2 h h + g ( x n + h , y n ) 2 /h h + h , (C1)where in this case h = ∆ x/ h = ∆ x , and g ( x n − h , y n ) is the stream function value on the boundary line.Equivalent equations for the Laplacian are used for other ∂ Ω u boundary lines and along junctions. Special atten-tion has to be given to points r n near corners. Appendix D: Johnson current noise
We assume that the uncorrelated Johnson noise fromthe normal resistances of the junctions are the domi-nant noise sources, compared to junction shot noise orthermal fluctuations in critical currents . When solv-ing Eqs. (34) and (35) numerically, the normalised noisecurrents I Noisek /I c , with k = 1 , ..., N , become sequences of random numbers corresponding to successive averagesover small time steps ∆ τ of the continuous noise currents I Noisek ( τ ) /I c . Thus, each noise current is an independentGaussian random variable with mean square deviation2Γ k / ∆ τ and zero average , where Γ k is the effectivenoise strength Γ k = 2 πk B TI ck Φ , (D1)and k B is the Boltzmann constant and T the tem-perature the device is held at (here T = 77 K). Whenapplying the above concept to a single resistively shuntedJJ, our numerical results agree with the Fokker-Planckcalculations of Ambegaokar and Halperin . ∗ Electronic address: [email protected] R. Feynman, R. Leighton, and M. Sands,
The FeynmanLectures on Physics, Vol. 3 (New York: Addison-Wesley,1966). J. Zimmermann and A. Silver, Phys. Rev. , 367 (1966). A. De Waele, W. H. Kraan, and R. De Bruyn Ouboter,Physica , 302 (1968). D. L. Stuehm and C. W. Wilmsen, Appl. Phys. Lett. ,456 (1972). J. Miller, G. Gunaratne, J. Huang, and T. Golding, Appl.Phys. Lett. , 3330 (1991). K.-H. M¨uller, Physica C , 717 (1989). R. Newrock, C. J. Lobb, U. Geigenm¨uller, and M. Octavio,Solid State Physics , 263 (2000). J. Oppenl¨ander, C. H¨aussler, and N. Schopohl, Phys. Rev.B , 024511 (2000). J. Oppenl¨ander, C. H¨aussler, and N. Schopohl, Physica C , 119 (2002). C. H¨aussler, J. Oppenl¨ander, and N. Schophol, J. Appl.Phys. , 1875 (2001). P. Longhini, S. Berggren, A. Palacios, A. L. de Escobar,and V. In, AIP Conference Proceedings , 254 (2011). V. K. Kornev, I. I. Soloviev, N. V. Klenov, and O. A.Mukhanov, IEEE Trans. Appl. Supercond. , 741 (2009). V. K. Kornev, I. I. Soloviev, N. V. Klenov, and O. A.Mukhanov, Supercond. Sci. Technol. , 114011 (2009). P. Longhini, S. Berggren, A. Palacios, V. In, and A. L.de Escobar, IEEE Trans. Appl. Supercond. , 391 (2011). V. K. Kornev, I. I. Soloviev, N. V. Klenov, and O. A.Mukhanov, IEEE Trans. Appl. Supercond. , 394 (2011). B. J. Taylor, S. A. E. Berggren, M. C. O’Brian, M. C.deAndrade, B. A. Higa, and A. M. Leese de Escobar, Su-percond. Sci. Technol. , 084003 (2016). E. E. Mitchell, K. E. Hannam, J. Lazar, C. J. Lewis, A. Grancea, S. T. Keenan, S. K. H. Lam, and C. P. Foley,Supercond. Sci. Technol. , 06LT01 (2016). V. K. Kornev, N. V. Kolotinskiy, A. V. Sharafiev, I. I.Soloviev, and O. A. Mukhanov, Supercond. Sci. Technol. , 103001 (2017). R. De Luca, J. Mod. Phys. , 526 (2015). E. Y. Cho, Y. W. Zhou, M. M. Khapaev, and S. A. Cybart,IEEE Trans. Appl. Supercond. , 1601304 (2019). X. Chen, L. Chen, Y. Wang, L. Wu, X. Liu, L. Ma, andZ. Wang, Scientific Reports , 9930 (2019). E. E. Mitchell, K.-H. M¨uller, W. E. Purches, S. T. Keenan,C. J. Lewis, and C. P. Foley, Supercond. Sci. Technol. ,124002 (2019). Fourie, IEEE Trans. Appl. Supercond. , 1300412 (2018). J. R. Clem and E. H. Brandt, Phys. Rev. B , 174511(2005). N. Terauchi, S. Noguchi, and H. Igarashi, IEEE Trans.Magn. , 7202804 (2015). E. H. Brandt, Phys. Rev. B , 024529 (2005). E. H. Brandt, Physica C , 327 (2007). A. Y. Khapaev, M. M. Kidiyarova-Shevchenko, P. Mag-nelind, and M. Y. Kupriyanov, IEEE, Trans. Appl. Super-cond. , 1090 (2001). M. M. Khapaev, M. Y. Kupriyanov, E. Goldobin, andM. Siegle, Supercond. Sci. Technol. , 24 (2003). S. K. Tolpygo and M. Gurvitch, Appl. Phys. Lett. , 3914(1996). B. D. Josephson,
Superconductivity (Dekker, New York,1969). M. Tinkham,
Introduction to Superconductivity (Dover,New York, 2004). F. London,
Superfluids , vol. 1 (Dover, New York, 1961). M. M. Khapaev, Differential Equations , 1019 (2005). M. M. Khapaev and M. Y. Kupriyanov, J. Phys: Conf. Ser. , 012041 (2010). J. Pearl, Appl. Phys. Lett. , 65 (1964). S. K. Lam, J. Lazar, J. Du, and C. P. Foley, Supercond.Sci. Technol. , 055011 (2014). C. P. Foley, E. E. Mitchell, S. K. H. Lam, B. Sankrithyan,Y. M. Wilson, D. L. Tilbrook, and S. J. Morris, IEEETrans. Appl. Supercond. , 4281 (1999). E. E. Mitchell and C. P. Foley, Supercond. Sci. Technol. , 065007 (2010). E. F. Talantsev and J. L. Tallon, Nat. Commun. ,doi: 10.1038/ncomms8820 (2015). D.-X. Chen, C. Navau, N. Del-Valle, and A. Sanchez, Phys-ica C , 9 (2014). C. D. Tesche and J. Clarke, J. Low Temp. Phys. , 301(1977). R. Gross, L. Alff, A. Beck, O. M. Froehlich, D. Koelle, andA. Marx, IEEE Trans. Appl. Superrcond. , 2929 (1997). R. F. Voss, J. Low Temp. Phys. , 151 (1981). V. Ambegaokar and B. I. Halperin, Phys. Rev. Lett.22