Effects of strong capacitive coupling between meta-atoms in rf SQUID metamaterials

We consider, for the first time, the effects of strong capacitive and inductive coupling between radio frequency superconducting quantum interference devices (rf SQUIDs) in an overlapping metamaterial geometry when driven by rf flux at and near their self-resonant frequencies. The equations of motion for the gauge-invariant phases on the Josephson junctions in each SQUID are set up and solved. Our model accounts for the high-frequency displacement currents through capacitive overlap between the wiring of SQUID loops. We begin by modeling two overlapping SQUIDs and studying the response in both the linear and nonlinear high-frequency driving limits. By exploring a sequence of more and more complicated arrays, the formalism is eventually extended to the N×N×2 overlapping metamaterial array, where we develop an understanding of the many ( 8N2−8N+3 ) resulting resonant modes in terms of three classes of resonances. The capacitive coupling gives rise to qualitatively new self-resonant responses of rf SQUID metamaterials, and is demonstrated through analytical theory, numerical modeling, and experiment in the 10–30 GHz range on capacitively and inductively coupled rf SQUID metamaterials.


I. INTRODUCTION
A Radio Frequency (rf) Superconducting Quantum Interference Device (SQUID) is simply a superconducting loop interrupted by a single Josephson junction (JJ).This device was originally introduced to measure small dc magnetic fields and to operate as a magnetic-fluxto-frequency transducer [1][2][3][4].Later, it was recognized that rf SQUIDs, due to their low loss, self resonance, and strong interaction with electromagnetic fields, can be used as meta-atoms in a metamaterial, both in the quantum [5,6], and classical limits [7,8].The rf SQUID effectively acts as a macroscopic-quantum split-ring resonator.The rf SQUID metamaterials were initially realized by covering a plane with rf SQUIDs to act as a nonlinear and highly tunable metasurface [9][10][11].
The properties of rf SQUID metamaterials can be tuned by means of dc magnetic flux Φ dc , rf magnetic flux Φ rf , temperature, and the structure of the metamaterial, which dictates the interactions between the meta-atoms.The closed superconducting loop creates the condition for magnetic flux quantization, while the single Josephson junction in the loop introduces the Josephson inductance L JJ = Φ0 2πIc(T ) cos δ [12][13][14].Here Φ 0 = h/2e is the magnetic flux quantum, with h Planck's constant and e the electronic charge, I c is the critical current of the junction, and δ is the gauge-invariant phase difference across the junction.The self resonant frequency of the rf SQUID is dictated by the junction capacitance and shunt capacitance (which are in parallel and have a total capacitance C) along with the geometrical inductance of the loop L, and the Josephson inductance L JJ .A dc magnetic flux applied to the SQUID loop creates screening currents that adjust the gauge-invariant phase δ in the JJ, and thus tune L JJ through a wide range of positive and negative values [15,16].The result is a highly tunable self resonant frequency that can span the mm-wave, microwave, and RF frequency ranges [9][10][11].In nonhysteretic SQUIDs (β rf ≡ 2πLI c /Φ 0 < 1) the flux tuning is periodic with period Φ 0 .The application of rf flux to the SQUID modifies the impedance of the JJ through the Josephson effect nonlinearity [17], I = I c (T ) sin δ, and changes the self-resonant frequency and tunability with dc flux [11,18].The rf SQUID response is also tunable with temperature through the temperature-dependent critical current I c and inductance of the superconducting loop.
The key enabling characteristic of the rf SQUID as a meta-atom is its nonlinearity, which has been thoroughly examined and demonstrated by the two-tone intermodulation distortion (IMD) experimental results on rf SQUID metamaterials [19,20].The intrinsic nonlinearity of the Josephson effect, along with the extreme tunability of rf SQUIDs, leads to bistability [21,22] and multistability [23,24] in their response to rf and dc driving fields.This in turn leads to complex and hysteretic behavior, including the phenomenon of transparency [18].Theory predicts that, under appropriate circumstances, driven rf SQUIDs will display strange nonchaotic attractors [25] and chaos [26,27].The rf SQUID metamaterials can also act as nonlinear gain media when immersed in passing electromagnetic waves [28][29][30][31][32][33][34], which is based on the nonlinear processes enabled by the Josephson effect that transfer energy to a signal at frequency f s from a strong pump signal at frequency f p .
In early theoretical and experimental works on rf SQUID metamaterials, the rf SQUIDs were packed together side-by-side in either one [9,10,35,36] or two dimensions, [11] with substantial long-range (dipoledipole) mutual inductance of the SQUID loops due to their close lateral proximity in the plane.This coupling gives rise to remarkable collective behaviors of the metasurface, such as chimera states [37][38][39][40][41][42][43], disorderdominated states [21,44,45], and coherent modes of oscillation [46].Such states have been directly imaged by laser scanning microscopy in the superconducting state under microwave magnetic flux driving illumination [47,48].Prior work examining collections of SQUID-like entities in two dimensions, not necessarily metamaterials, include the following: superinductors made up of planar ladders of superconducting wires/loops incorporating Josephson junctions [49], and Josephson transmission lines utilizing SQUID arrays to create magneto-inductive waveguides [50,51].Another type of Josephson metamaterial recently realized utilizes a tunable plasma edge created by current-biased linear arrays of Josephson junctions embedded in a three-dimensional waveguide [52].In contrast to the present work, this metamaterial interacts mainly with high frequency electric fields, rather than magnetic fields.
There have been proposals for three-dimensional versions of quantum metamaterials [27,53], and previous experimental work on three dimensional superconducting metamaterials based on spiral resonators [54].In this work, we experimentally realized three dimensional arrays of rf SQUIDs employed as metamaterials for the first time.By stacking the SQUIDs vertically, we introduce positive mutual inductive coupling for nearest neighbors, very different from the co-planar geometry and, more importantly, add the qualitatively new aspect of strong capacitive coupling that permits high frequency currents to flow between SQUID loops.To the best of our knowledge, such coupling has not been considered in the past in any aspect of SQUID physics or technology, and can lead to dramatic new properties of coupled SQUIDs.Our three-dimensional (3D) SQUID metamaterials have fluxquantized loops mixed with non-flux-quantized loops.These latter loops are enabled by the capacitors that host displacement currents between SQUIDs.Faraday's law is applied to all of the non-SQUID loops, in addition to the flux quantization condition in the appropriate loops.
Parasitic capacitive coupling can also occur in superconducting digital electronics (SCDE), based on the propagation of ps-duration single-flux quantum voltage pulses between logic circuit elements [55].The pulses are processed by means of inductive loops, typically based on superconducting wires including Josephson junctions, essentially acting as non-resonant SQUIDs.It is wellestablished that state-of-the-art SCDE suffers from an inefficient use of space on chip in many practical computing and signal processing applications.This is due to the fact that SCDE is based on magnetic flux, as opposed to the monopole electric charge utilized in CMOS electronics, and the resulting need to create and control dipole sources, such as current loops, transformers, inductors, etc. [56].One way to mitigate this problem is to create three-dimensional circuit layouts in which logic and wiring layers are distributed in the third dimension, separated by ground planes [57].However, this threedimensional geometry can introduce new and unexpected coupling effects between circuits.Our inductively and capacitively coupled rf SQUID metamaterials can act as a surrogate test-bed to study coupling effects in future highly-integrated SCDE circuits.
Quantum computers utilize large arrays of qubits with controlled interactions (typically either inductive or capacitive) between many pairs of qubits [58][59][60][61].For charge and phase qubits, the nearest-neighbour interactions are enabled by capacitors, rather than inductors [60,[62][63][64][65].Of recent interest is the design of a tunable coupler transmon that is capacitively coupled to a pair of qubits to achieve high-fidelity two-qubit gates [66][67][68][69].Our rf SQUID metamaterials differ in that multiple coupling capacitors are included, creating a highly integrated network of both capacitive and inductive coupling between all of the SQUIDs.We note that the effects of both capacitive and inductive coupling between rf SQUIDs has been considered as a step in deriving the quantum Hamiltonian of arrays of interacting qubits [70,71].For the flux qubits, the nearest-neighbor inductive coupling can be adjusted by introducing an intermediary SQUID between the qubits to be coupled, whose properties are tuned by a local magnetic flux [72][73][74][75].Another approach is to have two qubits share a common wiring loop bond.This bond may have a variable kinetic inductance, or Josephson inductance, that depends on the sum of the currents flowing in the two qubit loops through that bond [76][77][78].Our coupling design is uniquely different in that it introduces interactions through a combination of inductive and high-frequency capacitive coupling.The possibility exists to tune the capacitive coupling through external manipulation of the dielectric material in the capacitor.
The outline of the paper is as follows.In Section II A, the theory of inductively coupled rf SQUID meta-atoms is briefly reviewed.In Section II B we then introduce the strong capacitive coupling between two SQUID loops in what we call the corner-coupled geometry and discuss its effect on the dynamics of the SQUIDs.In Section II E, a sequence of larger overlapping rf SQUID structures are studied to uncover a total of three distinct classes of resonant modes from these unique metamaterial structures.Analytical results for the resonant modes and their dispersion with dc flux are obtained in the linear response limit, and comparisons are made to full nonlinear numerical simulations.This culminates in consideration of the N×N×2 metamaterial where 2 refers to the number of overlapping layers formed by the corner-coupled geometry and the calculation of the number of resonant modes in such system.In Section III we then compare the results of theory to experimental data on contrasting samples: a single-layer 12×12×1 and an overlapping 12×12×2 Nb-based rf SQUID metamaterial all utilizing the same rf SQUID meta-atom.In Section IV, we discuss the generalization of our model to non-regular arrays of overlapping corner-coupled SQUIDs, the limitations of our analysis, and opportunities for further studies of this remarkable class of superconducting metamaterials.

II. THEORY
A. Model for a system of inductively coupled SQUIDs based on the resistively and capacitively shunted junctions First we shall review the consequences of flux quantization in a single rf SQUID loop to establish our notation and approach to setting up and solving the equations of motion for the gauge-invariant phase.The total flux Φ in the rf SQUID loop and gauge-invariant phase difference δ across the junction are related through the statement of flux quantization in a superconducting loop, which requires the order parameter to be single-valued upon going on a continuous and closed loop C through the superconducting material.Mathematically, this self-consistency condition results in [2,79], where n = 0, ±1, ±2, .... and Φ is the magnetic flux through any surface that terminates on the continuous circuit C.This can be re-written as, Φ = nΦ 0 − Φ0 2π δ, where Φ 0 = h 2e is the flux quantum.Without loss of generality, taking n = 0 in the SQUID loop, the expression is simplified to Φ = − Φ0 2π δ.We shall assume that the SQUID loop maintains quasi-static flux quantization through the microwave frequency range so that Eq. ( 1) holds for both the rf and dc flux in a single galvanically-connected rf SQUID loop.One would expect that Eq. ( 1) remains valid for all situations in which the superconducting order parameter has a well-defined phase throughout the material, which should extend to time scales as short as the order parameter relaxation time, expected to be in the ps range, except close to T c [80,81].
The total flux in the loop can be expressed as [2,3,16,79], where Φ app stands for the applied flux and the minus sign is chosen to account for the diamagnetic relationship between the applied and induced fluxes.The induced flux is Φ ind = LI for a single SQUID with a loop inductance L carrying current I.Note that this relationship can be generalized to a system of many interacting SQUIDs as ⃗ Φ ind = ← → L ⃗ I, where the inductance matrix ← → L contains self-inductances on its diagonal and mutual inductances between SQUIDs in the off-diagonal elements [46].
The current I in the junction is expressed using the resistively and capacitively shunted junction (RCSJ) model where the junction is treated as a parallel combination of three branches: an ideal Josephson junction, a capacitor C and a shunt resistor R [82].The total current from the three branches is thus: Using the second Josephson equation, one can relate the voltage drop on the junction V to the time derivative of the gauge-invariant phase as V = Φ0 2π δ, where the overdot denotes time derivative.In this case, the response current can be written as, The flux quantization condition, Eq. ( 2) can now be written for a system of inductively-coupled identical rf SQUIDs as, [46] whose normalized form reads where ⃗ Φ dc and ⃗ Φ rf are the vectors of dc and rf magnetic flux applied to each SQUID in the array, respectively, and we assume time-harmonic rf flux at a single frequency ω.Here ⃗ δ is the array of gauge-invariant phase differences on the junctions in all of the rf SQUIDs in the array, and the over-dot □ represents derivative with respect to time.The equation can be reduced into the dimensionless form with the following substitutions: ϕ dc,rf = 2πΦ dc,rf /Φ 0 , ← → κ = ← → L /L geo , where L geo is the geometric inductance for a single SQUID in the system, β rf = 2πL geo I c /Φ 0 , γ = L geo /C/R, τ = ω geo t = t/ L geo C, and Ω = ω/ω geo = ω L geo C. Note that we also introduce the geometric resonance ω geo = 1/ L geo C, which is the resonant frequency of the rf SQUID meta-atom in the absence of the Josephson effect.
Equation (4) is a system of driven second-order coupled nonlinear differential equations for the gauge-invariant phase vector as a function of time, ⃗ δ(τ ), which dictates the response of the metamaterial to external electromagnetic fields.Solving for ⃗ δ(τ ) allows one to calculate all observable properties of the system.Prior work has explored solutions to these equations for purely inductively-coupled rf SQUID metamaterials [19,20,46].The following sections focus on extending this model to capacitively-coupled rf SQUIDs.We begin with the simplest case of a pair of corner-coupled overlapping rf SQUIDs, and then consider a sequence of larger and larger arrays of overlapping corner-coupled SQUIDs, eventually addressing the problem of the N × N × 2 system.The parameters for the SQUIDs in the following calculation are given in Appendix B and Table.II.

B. Model for two corner-coupled overlapping SQUIDs
The simplest model for overlapping SQUIDs is a pair of rf SQUID loops having wiring layers overlapping each other at the corner, and the overlapping portions are separated by a thin dielectric layer, which forms two capacitors whose capacitance C ov is comparable to that of the junction, C (see Fig. 1 (b)).It should be noted that these overlapping capacitors do not include Josephson coupling between the superconducting wires, but do create a route for high-frequency displacement currents to flow between the wiring loops of neighboring SQUIDs.The capacitors shunt the SQUID loop, breaking the uniformity of highfrequency currents in the loops, which leads to different currents in the non-junction branches I a1(b1) from the currents in the junction branches I a0(b0) as illustrated in Fig. 1 (a).Note that the overlapping capacitors have no direct influence on the dc currents, which are constrained to flow only through individual SQUID loops.Following the same flux-to-current approach in Sec.II A, one can write down the flux quantization conditions for the two corner-coupled SQUIDs (called a and b) as follows: where the induced flux is expressed as contributions from different segments of the two SQUID loops on the second line.The elements of the first (second) row in the inductance matrix are determined as the partial inductance between the individual segments denoted by the subscript, and the galvanically connected SQUID loop a (b).Partial inductance is a concept that generalizes the inductance of a closed loop to that of segments in the loop [83,84].Consider a segment 1 in a closed loop c that is experiencing a time-varying magnetic field due to the current in segment Although Eq. ( 5) relates the flux with the currents, the non-junction currents I a1,b1 are not yet expressed as functions of the gauge-invariant phase differences δ a,b and their time derivatives.To obtain the equations of motion in ⃗ δ, one needs to invoke Faraday's law on the center overlapping loop from capacitor 1 to 2 to junction b and back to capacitor 1 in Fig. 1 (a), as well as current conservation on the capacitor nodes.This will allow us to solve for I a1,b1 in terms of the gauge-invariant phases.

Faraday's law applied to the center overlapping loop
There is no need to consider Faraday's law in the formulation of the equations for the non-overlapping SQUIDs, since it is implicitly incorporated in the flux quantization condition Eq. (4).In fact, applying Faraday's law to a single SQUID loop will result in the time derivative of Eq. ( 4).On the other hand, its application is necessary for a superconducting loop interrupted by capacitors, such as the small center loop carrying currents I b0 and I a1 formed by the two overlapping corner-coupled SQUIDs in Fig. 1 (a).By applying Faraday's law to this center loop (denoted cen), one finds: with Φ ind cen = M cen,a0 I a0 + L cen,b0 I b0 + L cen,a1 I a1 + M cen,b1 I b1 .V 1 , V 2 are the voltages across the capacitors at nodes 1 and 2. The sign convention of the capacitor voltages are chosen so that the positive voltage corresponds to the electric field pointing out of the plane, from the bottom Nb wiring layer to the top Nb wiring layer or from SQUID b to SQUID a. V b is the voltage across the junction of SQUID loop b, and the currents I a0,b0,a1,b1 are labeled in Fig. 1 (a).Again, the partial inductances M cen,a0 , L cen,b0 , L cen,a1 , M cen,b1 from each segment to the center overlapping loop (cen) are involved in the expression for the induced flux Φ ind cen .

Conservation of current through the overlapping capacitors
The effect of capacitive coupling on the corner-coupled SQUIDs are understood through the conservation of currents applied to the overlapping capacitors: where nodes 1 and 2 have identical capacitance C ov based on our design.The current conservation statements, Eqs.(7), reduce to the following relations: The flux equation Eq. ( 5) allows us to express the nonjunction currents I a1 and I b1 in terms of the junction currents and the gauge-invariant phases, in other words I a1 (I a0 , I b0 , δ a , δ b ), I b1 (I a0 , I b0 , δ a , δ b ) as: The definition of CD is given in Appendix A. The current conservation statement Eq.( 8) now becomes a constraint on δ a , δ b , I a0 , I b0 after substituting I a1,b1 from Eqs. (10): where κ a,b and L δa,δb are defined in Appendix A. After substituting the expressions for I a1,b1 in Eqs. ( 10) into Faraday's law Eq.( 6), the time derivatives of junction voltages in Eq. ( 9), V1 = − V2 , are determined as: The parameters κ va , κ vb , L Ia , L Ib are defined in Appendix A. One can see that Eq. ( 12) contains 4thorder derivatives of δ a,b , by noting that I a0,b0 given in the RCSJ model Eq. ( 3) brings two more time derivatives into the equation.The expression for V1 Eq. ( 12) can then be substituted into the current laws Eq. ( 7) to obtain solutions for the non-junction currents I a1,b1 (I a0,b0 , Ïa0,b0 , δa,b ).

Equation of motion for gauge invariant phase differences
The flux equations for the two SQUID loops can now be set up.Assuming that the applied dc and rf flux amplitudes are the same in both SQUID loops, and that the rf flux is sinusoidal at a single frequency ω with amplitude Φ rf , the flux equation, Eq. ( 5), becomes: where the induced flux (last two terms on the right hand side of Eq. ( 13)) is separated into two contributions: the conventional inductively-coupled SQUIDs with mutual inductance M , and the correction due to the overlapping capacitors.The term with the inductance matrix can be regarded as the limit without capacitive coupling, when C ov = 0. Consequently, the current becomes uniform inside the SQUID loops such that I a0 = I a1 , I b0 = I b1 .The 2 × 4 inductance matrix in Eq. ( 5) is then reduced to the 2 × 2 matrix above with self inductance of the SQUID loop determined as , and the mutual inductance between the two SQUID loops , which can be positive, zero, or negative depending on the overlapping area between the two SQUID loops.The last term in Eq. ( 13) involving C ov brings in qualitatively new phenomena in the high frequency response of coupled rf SQUIDs.

Linear-limit solutions
To analytically understand the two corner-coupled SQUID system, one can start by simplifying the equations in the low-driving-amplitude linear limit when |ϕ rf | ≪ 1.The full solution to δ(t) of any individual SQUID can be separated into its dc and rf components: [19].Under a weak applied rf flux, the rf response is also vanishingly small, |δ rf | ≪ 1, so that any nonlinear rf response is negligible.Therefore, the solution takes the following form δ = δ dc +δ rf exp(iωt), where ω is the driving frequency [19].The junction current from the RCSJ model under the weak rf flux approximation is , where the terms with second or higher order in δ rf are dropped.After substituting the expressions for I a0 and I b0 in the equation of motion Eq. ( 13), and converting to the dimensionless vector format, one obtains a system of algebraic equations: with where ← → κ , ← → κ loop , ← → κ δ , and ← → κ I are defined in Appendix A.Here ← → κ is the 2 × 2 conventional inductive coupling matrix for the two SQUIDs without capacitive coupling, as in the middle term on the right hand side of Eq. ( 13).
Consider the case when the applied flux ⃗ ϕ rf is zero, the resonance condition occurs when nontrivial solutions for ⃗ δ rf exist, which requires a non-invertible response tensor, or det( ← → χ ) = 0.The real parts of the solutions to this characteristic equation in Ω = ω/ω geo are the resonance frequencies, plotted in Fig. 2 as a function of dc flux applied to the SQUIDs.The characteristic equation is only sixth order in Ω, since the matrix ← → κ I in front of the (−iγΩ 3 + Ω 4 ) term in Eq. ( 15) is non-invertible with zero determinant.There are thus six roots from solving the sixth order equation, which come in pairs where the real parts are opposite to each other.Only the three positive roots are shown in Fig. 2. The tuning curves of eigenfrequencies are periodic in applied dc flux with a periodicity of Φ 0 , with each curve centered at an integer multiple of Φ 0 .However, each resonant solution curve extends beyond the dc flux range of 1Φ 0 and overlaps the adjacent curves due to the hysteresis from the SQUID loop, since β rf > 1.
The three resonances cover a much broader frequency range compared to that of a single SQUID, which is shown as a black dotted curve in Fig. 2. Mode II (the yellow curve) out of the three resonances closely follows the dc-flux tunability of a single SQUID loop resonance.To better understand the nature of the other two modes, one can examine the solutions ⃗ δ rf to the linearized equations, Eqs. ( 14) at the corresponding eigenfrequencies.In particular, the current values (I a0,b0,a1,b1 ) in the SQUID loops as a function of frequency at zero dc flux are shown in Fig. 3, expressed as the dimensionless currents ι = 2πL geo I/Φ 0 .The currents indeed undergo resonances near the three eigenfrequencies Ω ∼ 2.3, 2.55, 4.88 in Fig. 2 for zero dc flux.After comparing the current values from the three different eigenmodes, the dominant current distribution for each mode can be summarized in the schematics in Fig. 4, where the branches with strong currents are highlighted in red.Mode II at Ω ∼ 2.55 clearly stands out as the only mode where the current remains uniform inside one SQUID loop, just as in side-by-side pairs of inductively-coupled SQUIDs, which explains the match between the single SQUID eigenfrequency and the mode II in Fig. 2.
FIG. 4: The dominant rf (not dc) current distribution for the three linearized eigenmodes.The three modes are organized from left to right to match the corresponding eigenfrequencies from low to high frequency at zero dc flux in Fig. 2. The dominating junction and the "loop" in the mode are shaded in red.The capacitor nodes are shaded red when the rf current passes through the capacitors in the corresponding mode.
The frequencies for the modes can also be estimated quantitatively from the current distribution in the circuit.
For a single SQUID loop, the resonance can be predicted from the lumped element model, Ω res = ω res /ω geo = (LC) −0.5 /ω geo = This circuit model can be generalized to one SQUID loop a in a large coupled system: where the loop inductance has changed from the geometric inductance L geo to the effective inductance L a,eff (ω) = Φ ind a (ω)/I a (ω), accounting for the coupling from other SQUID loops in the large system.For instance, in the low power linear limit, a planar inductively coupled system has antiferromagnetic couplings among the SQUIDs.Thus, L eff is always lower than L geo , resulting in a slightly higher Ω res than expected for the single SQUID design parameter.This property no longer holds true in an overlapping system.In particular, for the two corner-coupled SQUIDs model, Φ ind a is given in Eq. ( 5).The resulting L eff and Ω res calculated for the loop a and b at the three resonance modes are listed in Table I.As a consequence of the nonuniform current distribution in one SQUID loop, the real part of the effective inductance can take on much wider range of values from higher than L geo to large negative values.The corresponding resonant frequencies agree with the eigenfrequencies Ω eig solved from the linearized characteristic equation shown in Fig. 2. The inductance of a circuit generally scales with its size.Therefore, the longer dominant current branches in modes I & II in Fig. 4 contribute to larger effective inductance magnitudes compared to mode III, where the short segments dominate the current distribution.A small positive effective inductance, 0 < L eff ≪ L geo , leads to a resonance frequency much higher than that of the single SQUID, √ 1 + β rf cos δ .In contrast, the larger magnitudes of the effective inductance, |L eff | ≳ L geo , in modes I & II in Fig. 4 lead to lower resonance frequencies close to the prediction for a single SQUID.The apparent difference in dc flux tunability between the highest frequency mode and the other modes can also be explained by the magnitude of effective inductance.The small effective inductance leads to a large L geo /L eff that renders the dc tuning term represented by β rf cos δ less effective in Eq. (16).

C. Full nonlinear numerical solutions to the two corner-coupled SQUIDs
Although analytical solution to the system of equations at rf flux driving levels beyond the linear limit is difficult, we can obtain the full nonlinear solution numerically.For the convenience of the numerical solver, the equations are first converted into dimensionless form as in Eq. ( 4) with the additional introduction of dimensionless currents: ι = 2πL geo I/Φ 0 .In the general practice of numerically solving a system of differential equations, the equations are first reformulated as a system of first order initial value problems.The equations of motion Eq. ( 13) consist of two flux equations, each a 4th-order differential equation for δ.However, due to the constraint in Eq. ( 11) relating δ a,b and I a0,b0 , there are only six degrees of freedom, two less than otherwise expected.This can be illustrated by manipulating the matrix expression in Eq. ( 13) as follows: L −1 δa Row 1 + L −1 δb Row 2, which is equivalent to the constraint in Eq. (11).Therefore, instead of solving the overdetermined system in Eq. ( 13) directly with eight variables, one should reduce the system to one equation of motion for one of the SQUIDs, along with the constraint Eq. ( 11), and establish the initial value problems with six variables: δ a,b , δa,b , ι a0 , ιa0 , as follows.
where ϕ app = ϕ dc + ϕ rf sin(Ωτ ) is the dimensionless applied flux.The second time derivatives in δ in Eqs.(17 c, d) are related to the currents using the RCSJ model, Eq. ( 3).The constraint Eq. ( 11) is invoked to express ι b0 in Eq. (17 d).The second time derivative of current (ϊ a0 ) in Eq. (17 f) is obtained from the equation of motion Eq. (13).All other parameters are defined in Appendix A. The resulting system of initial value problems was solved with the LSODA function from SciPy based on the FORTRAN library ODEPACK.
To compare with experimental results, one needs to convert the solutions for ⃗ δ(t) into measurable quantities.Here, the power dissipated by the resistive channel of the junction, and the resulting change in transmission magnitude through the metamaterial are calculated as [18], where the sum is over Josephson junctions in the metamaterial (here i = a, b and V i is the voltage drop on junction i), P incident is the total rf power incident on the metamaterial (which provides the rf flux bias to the SQUIDs), and S 21 is the transmission coefficient through the metamaterial at frequency f .P incident is related to the applied rf flux Φ rf as follows.For the experimental setup in Fig. 12, the sample lies in the center of a rectangular waveguide and is perpendicular to the rf magnetic field of the propagating TE 10 mode.The rf magnetic field, and thus the rf flux Φ rf at the location of the SQUIDs, can be calculated from the incident power P incident , the waveguide dimensions, and the frequency [11].The total dissipated power is summed over the individual contribution for each SQUID i.This calculation assumes that the only lossy element in the setup is the normal charge carriers tunneling in the junction, represented by the parameter R.
The resulting transmission at low applied rf flux amplitude (Φ rf ∼ 10 −3 Φ 0 ), near the linear limit, as a function of dimensionless driving frequency Ω and dc flux is plotted in Fig. 5.The dark bands represent the resonances in the rf SQUID system, where the amplitudes of the rf currents in the SQUID loops, and thus the dissipated powers in the junctions, are maximized.Unlike the case for a pair of side-by-side rf SQUIDs, which have two resonant modes, there are now three distinct modes tuned by the applied dc magnetic flux.The red lines in Fig. 5 show the dispersion of the linearized solutions from Eq. ( 14), and show good agreement with the solutions to the full nonlinear equations, Eqs.(17), in the weak-driving limit.

D. Nonlinear Properties of the Corner-Coupled SQUIDs
Here we examine the evolution of the three modes as the amplitude of the rf driving flux is increased.Figure 6 shows the evolution of the resonant modes from the linear limit Φ rf /Φ 0 ∼ 10 −3 at zero dc flux to higher rf flux amplitudes.All three modes show a suppression of their resonant frequencies in a manner similar to that observed for single rf SQUIDs [11].Note that the two lower frequency modes show substantial tuning with rf flux amplitude, but the high-frequency mode is only weakly affected due to the small effective inductance of this mode leading to a resonance which is less sensitive to the applied magnetic flux, according to Eq.( 16).Also note that all three modes achieve linear response again at high driving amplitudes, in the sense that the resonant frequencies are independent of rf flux amplitude when Φ rf ≳ Φ 0 .This behavior was observed before [11] and can be attributed to the term linear in ⃗ δ in Eq. 4. At large rf driving amplitudes, the sin ⃗ δ term is bounded between ±1, much smaller than the leading-order ⃗ δ term, which reduces the equation of motion to the form of a harmonic oscillator.Further examination of the nonlinear properties of these hysteretic SQUID metamaterials will the subject of future work.FIG.6: Nonlinear solutions to the two corner-coupled SQUIDs for Φ dc = 0.The quantity plotted is |S 21 |(dB) on a logarithmic color scale, as a function of applied rf magnetic flux amplitude in units of the flux quantum Φ 0 and dimensionless frequency Ω = ω/ω geo .The linear limit discussed in Sec.II B 5 is reproduced at low applied power on the left of the plot.Upon increasing Φ rf above Φ 0 , the high-power linear limit is achieved where the SQUID loop resonance is suppressed to the geometric frequency (Ω = 1).

E. Model for larger overlapping SQUID systems
Building upon the formalism developed for the two corner-coupled SQUIDs in sections II B-II D, we now consider the larger systems of overlapping SQUIDs as exemplified in Fig. 7.One major difference in the larger systems compared to two corner-coupled SQUIDs is the introduction of a new kind of circuit loop enclosing an area outside any galvanically connected SQUID loop.To better distinguish the different loop circuits in the large systems, and to streamline the discussion, some common vocabulary should be established.Highlighted in blue in Fig. 7 (a) is the new loop, named the "extra-SQUID" loop, while the loop formed at the corners of two overlapping SQUIDs, as studied in the two-SQUID case, is referred to as a "partial loop", colored red in Fig. 7 (a).Together with the conventional galvanicallyconnected SQUID loops, these three types of loop circuits dictate the dynamics of the gauge-invariant phases through either Faraday's law (e.g.Eq. ( 6)), or the flux quantization condition in the SQUID loop (e.g.Eq. ( 4)).Another challenge in modelling the larger system is the fact that each SQUID loop is broken into many segments by the additional capacitive nodes formed by the overlap with wiring of several neighboring SQUIDs.For instance, each SQUID in Fig. 7 (a) has 4 capacitor nodes and thus 4 segments.There are in total 16 segments, each with a different current.Attempting to follow a similar treatment to that used for the two corner-coupled SQUIDs in Sec.II B, one would need to express the 12 non-junction currents in terms of 4 junction currents, 4 gauge-invariant phase differences and their time derivatives.The number of unknown currents can be further reduced by invoking current conservation laws on the SQUID loops requiring zero net current into any SQUID.The constraints need to be applied to all SQUID loops but one, since the net current into the last SQUID is equal to the net current leaving the rest of the SQUID loops which is already zero due to the current conservation laws.The number of in-dependent non-junction currents is thus 12 − (4 − 1) = 9, which still could not be completely determined from 4 flux quantization conditions (Eq. ( 4)) from the 4 SQUID loops.We should note that applying Faraday's laws to the 5 non-SQUID loops only relates the time derivatives of the currents, not the currents themselves.
To circumvent this problem, the model should be reformulated in terms of the voltages across the capacitor nodes as the variables instead of the junction currents as commonly done in the Lagrangian formalism for a magnetic circuit.All the junction currents are expressed using the RCSJ model in Eq. ( 3) in terms of the set of ⃗ δ, and their time derivatives.The non-junction currents are obtained from current conservation laws on the capacitor nodes (e.g.Eq. 7).For the four corner-coupled SQUIDs in Fig. 7 (a), there are in total 8 capacitor nodes.Due to the current conservation laws mentioned above, the number of independent voltages is reduced to 5, corresponding to the number of non-SQUID loops.This system can then be set up and solved analytically as demonstrated in Sec.II E 1.This voltage formalism turns out to be a more general approach for modeling overlapping SQUIDs compared to the current formalism followed in Sec.II B for the two corner-coupled SQUIDs.Appendix C illustrates the application of the voltage formalism to the two corner-coupled SQUIDs, and arrives at the same eigenfrequency solutions as found in Sec.II B.

Four corner-coupled overlapping SQUIDs
With reference to the system shown in Fig. 7 (a), the dynamics of the system is now described by the 4 gaugeinvariant phase differences δ a,b,c,d and their time derivatives, as well as 5 independent capacitor nodal voltages V 1,2,3,4,5 and their time derivatives.The voltages across the rest of the capacitor nodes can be expressed in terms of V 1,2,3,4,5 through the current conservation laws inside each continuous SQUID loop.Applying Faraday's law to the non-SQUID loops and flux quantization conditions to the SQUID loops, we obtain the following system of equations of motion. where ).The top 4 rows are from the flux quantization conditions on the 4 SQUID loops, and the bottom 5 rows from application of Faraday's law to the non-SQUID loops.The non-SQUID loops are labeled by the SQUID loops involved in forming their circuit.For example, the top left partial loop in Fig. 7 (a) is labelled as ac.The voltages in the parentheses, e.g.(V 1 + V 2 + V 3 ) in the sixth row, are dependent nodal voltages obtained from applying current conservation laws to the continuous SQUID loops.The induced flux Φ ind xx is calculated from the partial inductances from the individual branches just as in Sec.II B, which is a function of junction currents ⃗ I and node voltages ⃗ V .
The low rf flux amplitude linear limit solution can then be obtained by expressing I a,b,c,d in terms of δ a,b,c,d and their time derivatives using the RCSJ model in Eq. ( 3), and replacing time derivatives with iΩ under Fourier transform.The resulting system has a size 9 × 9 with the 9 variables: δ a,b,c,d and V 1,2,3,4,5 .The eigenfrequencies for this linear system can be calculated and are plotted in Fig. 8 as a function of dc flux.As discussed in section II B 5, and illustrated in Figs. 2 and 4, the high frequency modes are dominated by the smaller loops, which contribute to smaller effective inductances, explaining their higher resonance frequency and weaker tunability, according to Eq. ( 16).The same loop-resonance correspondence can be established for the system of four corner-coupled SQUIDs containing 4 SQUID loops, 4 partial loops and 1 extra-SQUID loop, as labeled in Fig. 8.As expected, the much smaller extra-SQUID loop brings about a very high resonance at Ω ≈ 12 in Fig. 8. Similar to the two corner-coupled SQUIDs discussed in Sec.II C, the full nonlinear solution is obtained numerically.The details for setting up the numerical solver can be found in Appendix D. Fig. 9 shows the resulting transmission |S 21 | as a function of driving frequency Ω and applied dc flux Φ dc in the low rf flux amplitude linear limit, Φ rf ∼ 10 −3 Φ 0 .The transmission calculated from the numerical solution has good agreement with the eigenfrequencies from the linear limit solutions, as illustrated by the coincidence of the dark (absorbing) features with the red dotted curves.A separate calculation is performed around Ω = 12 which resolves the extra-SQUID loop modes that can not be captured in the frequency range in Fig. 9.

2 × 2 × 2 corner-coupled overlapping SQUIDs
A full numerical solution, or even an analytical solution in the linear limit, is very computationally expensive to obtain for the large N × N × 2 system studied experimentally.However, the case of N = 2 (Fig. 7 (b)) can be tackled easily and should illustrate the general properties of the model.For the 2 × 2 × 2 system, there are 8 SQUID loops, 9 partial loops, and 2 extra-SQUID loops.In the absence of any symmetries, we would thus expect a total of 19 resonant modes, comprised of 8 lower frequency modes near the single SQUID resonance, 9 partial loop modes at about twice the single-SQUID resonance frequency, and 2 extra-SQUID loop modes near Ω = 12.
Following the same treatment as in Sec.II E 1, the eigenfrequencies in the low rf flux amplitude linear limit, and the full nonlinear numerical solution, can be obtained, and are shown in Figs. 10 and 11, respectively.Indeed, 19 eigenfrequencies fall in the expected range for their corresponding loop modes as depicted in Fig. 10.We note once again the general decreasing sensitivity to dc flux with increasing frequency of the modes, consistent with Eq.( 16).
The full numerical solution in the low rf flux amplitude limit in Fig. 11 shows resonances coincident with the eigenfrequencies obtained from the linear-limit solution.In the general case of an N × N × 2 SQUID array, there are 2N 2 SQUID loops, (2N − 1) 2 partial loops, and 2(N − 1) 2 extra-SQUID loops, resulting in a total number of 8N 2 − 8N + 3 distinct loops.As illustrated above, each loop corresponds to an equation of motion in the voltage formalism.Thus, one would expect 8N 2 − 8N + 3 equations for an N × N × 2 SQUID array.There are two capacitor nodes for each partial loop.However, the current conservation law from the SQUID loops will constrain the number of independent nodal voltages to 2×# of partial loops−(# of SQUID loops−1) = 6N 2 − 8N + 3. Together with the 2N 2 gauge-invariant phase differences δ for each SQUID, the dynamics is described by a total of 8N 2 − 8N + 3 variables.The resulting (8N 2 − 8N + 3) × (8N 2 − 8N + 3) system can be solved exactly (see Appendix D) and will generate 8N 2 −8N +3 modes in the absence of any symmetry.
Equipped with the loop-mode correspondence established in the analysis for the three model systems discussed in section II, we can extrapolate to the prediction for a larger system, where many partial loop modes should be visible between Ω = 3 and 6, with very weak dc flux tunability, while the SQUID loop modes occur at lower frequency range with high dc flux tunability.The extra-SQUID loop modes around Ω = 12 is almost invariant in dc flux and is beyond the frequency range of our measurement capability.

III. EXPERIMENT A. SQUID samples
The samples in this work are fabricated by the STAR Cryoelectronics Selective Niobium Anodization Process (SNAP) with Josephson junction critical current density J c = 1µA/µm 2 [85].The Nb films have a critical temperature of T c = 9.2K, and the samples were measured at T = 4.6K.A representative three-dimensional rf SQUID metamaterial sample is shown in the inset of Fig. 12.The main concern in the design is to maintain the rf-SQUID self-resonant frequency in the 10-30 GHz range while keeping β rf = L geo /L JJ reasonably small (≲ 10) to avoid extreme hysteresis.The low resonant frequency of the rf-SQUID requires a large junction inductance, corresponding to a small critical current.The requirements of small junction size and low critical current density is an uncommon constraint for Josephson junction foundries, since most superconducting electronics applications prefer high Josephson energy (E J = Φ 0 I c /(2π)) and thus high critical current of the junction.

SNAP fabrication
SNAP fabrication starts by first depositing the Nb/Al-AlO x /Nb trilayer on the wafer with J c = 1µA/µm 2 .This trilayer is then patterned with a wet etch that defines the base layer of the SQUID loop, the vias between the top and the base wiring layers of the SQUID, as well as the anodization rails connecting all of the SQUIDs to the edge of the wafer.The anodization rails are required to supply a voltage bias to the Nb layers for the anodization process in the next step.The junctions, anodization rails, and the area surrounding the vias [86] are protected by photoresist, while the remaining exposed Nb is anodized to form the 100 nm-thick insulating dielectric layer of Nb 2 O 5 between the top and base Nb layers.Previous measurements of the dielectric constant of Nb 2 O 5 made by this process yield ϵ r = 29 [87,88].Next, the second Nb layer is deposited and patterned to form the top wiring layer of the SQUID.After the SQUIDs are defined, the anodization rails connecting the SQUIDs are finally removed with a wet etch.

SNEP fabrication
The first generation of the overlapping SQUID samples, as shown in Fig. 1, was made with a slightly more complicated process, the Selective Niobium Etch Process (SNEP) from STAR Cryoelectronics.This process begins with the same Nb/Al-AlO x /Nb trilayer deposition.However, instead of anodization, the junctions are defined by a reactive ion etch (RIE) on the top Nb layer of the trilayer.The entire trilayer is then patterned and etched to form the base wiring layer.A layer of SiO 2 with a thickness of 300 nm is then deposited with plasma enhanced chemical vapor deposition (PECVD), and this serves as the insulating dielectric between the top and base layers.The dielectric deposition is then followed by another etch to open contact vias through the SiO 2 layer.Next, the top Nb wiring layer is deposited and patterned.The process concludes with a final passivation of the wafer by depositing a layer of SiO 2 .In addition to its complexity, the unit area capacitance from the dielectric SiO 2 in SNEP is very low ≈ 0.15 fF/µm 2 , compared to that for the Nb 2 O 5 dielectric in SNAP ≈ 2.5 fF/µm 2 .Therefore, to maintain a low self-resonance for the rf SQUID, a larger capacitor pad is needed for samples made by the SNEP process, which further constrains the design.

SQUID array design
Since the fabrication processes only allow for a single trilayer for the junction definition, one cannot simply stack two independently-defined layers of 2D rf SQUID arrays in the third dimension.One of the layers must be shifted in-plane so that its junction avoids that of the other layer of SQUIDs, and this constraint creates the peculiar overlapping geometry studied here.To achieve the most symmetric configuration, the overlapping area between the SQUIDs from the two layers is designed to be roughly a quarter of the single-SQUID loop area.The shifted stack between the two 2D rf SQUID arrays results in the many overlapping capacitors C ov between SQUID loops from different layers.The design parameters for the overlapping SQUID metamaterial sample SNAP 161A is summarized in  The mutual inductances in the last two rows of the Table II represent the inductive coupling in the absence of the overlapping capacitors.However, under the strong capacitive coupling studied in this work, the rf currents can leave the SQUID loop through the capacitor nodes, breaking the uniformity of the current within one SQUID loop.The mutual inductance alone is thus no longer sufficient to correctly treat the coupling between SQUIDs.

B. Experimental Setup
The rf SQUID metamaterial is embedded in the center of a brass WR-42 waveguide, which has a cut-off frequency of approximately 14.1 GHz, and provides good transmission above 15 GHz.The sample is oriented along the direction of the wave propagation and its planar surface is perpendicular to the rf magnetic field in the first propagating TE mode of the waveguide.The dc magnetic flux is provided by a superconducting Helmholtz coil mounted on the outside of the waveguide and oriented to apply a nominally uniform field perpendicular to the SQUID array surface.The microwave transmission S 21 through the sample-loaded waveguide is measured by a microwave vector network analyzer (VNA), as shown in Fig. 12.The amplitude of the incident signal is reduced by input attenuators to carefully control the rf flux ϕ rf witnessed by the sample.The transmitted signal is amplified by a cryogenic amplifier and by a room temperature amplifier, before being measured by the VNA.The measurements are carried out at a temperature of 4.6 K.
Note that no wires are connected to the rf SQUIDs, and there are no galvanic contacts between different SQUIDs, hence all changes to the properties of the metamaterial occur by strictly 'wireless' means.FIG.12: Schematic of the rf SQUID metamaterial transmission measurement, showing the waveguide which hosts the metamaterial, the rf and dc field components, the microwave network analyzer, as well as attenuators and amplifiers.Inset: optical image of the overlapping 12×12×2 SQUID array sample made with the SNAP process.The blue shaded links between SQUIDs are the remnants of the anodization rails after wet etch.Variables that can be controlled in the experiment include Φ rf , Φ dc , rf driving frequency f , and sample temperature T .

C. Measurement results
The rf SQUID response is typically very small in the transmission measurement due to its weak coupling to the waveguide mode.To better resolve the SQUID resonances, a background response obtained from averaging over the transmission spectra S 21 at different dc fluxes is removed from the individual frequency spectrum.The resulting ∆S 21 = S 21 (f, Φ dc ) − S 21 as a function of applied dc flux Φ dc and microwave frequency f is plotted in Fig. 13 for eight different values of the rf flux amplitude Φ rf .
The yellow color represents nearly complete transmission of the microwaves (i.e.∆S 21 ≲ 0 dB), while the darker green features show conditions where the metamaterial interacts strongly with the passing electromagnetic fields and dissipates power.The green features trace out the resonant response as a function of applied dc magnetic flux with a typical 1Φ 0 periodicity.The red dashed curves mark the single-SQUID resonance and roughly corresponds to the expected boundary between the lowfrequency and highly-tunable SQUID modes, and the high-frequency partial loop modes with lower-tunability introduced in Sec.II E 1.It should be noted that the third kind of mode, the extra-SQUID loops modes, with frequencies around 12f geo ≈ 100 GHz, are beyond the measurement range of our apparatus.Results are presented at eight different applied rf flux amplitudes: [0.033, 0.023, 0.013, 0.01, 0.007, 0.006, 0.005, 0.0006] Φ 0 in panels a) through h).The red dashed curves denote the single-SQUID eigenfrequency.
Figure 14 shows the effects of increasing the driving rf flux amplitude beyond linear response on the spectrum of modes in the overlapping 12 × 12 × 2 sample.The measurement was performed at zero current in the magnet but Φ dc = −0.3Φ0 on the SQUIDs, based on their dc flux tunability curves.The result in Fig. 14 is obtained from two separate power sweeps on the VNA due to its limited dynamic range.The first power sweep was performed from -82 to -42 dBm, while the second from -62 to -35 dBm.In the range where both power sweeps overlap, the average response is shown in Fig. 14  One notes strong tuning of the modes below Ω = 2.5 and relatively small tuning for the higher frequency modes.This behavior is in qualitative agreement with that shown in Fig. 6 for the two corner-coupled SQUIDs.The low-frequency SQUID loop modes are strongly tuned, whereas the partial-loop modes only show modest tuning with rf flux.The experiment appears to just reach the high-power linear limit that is clearly seen in the model results in Fig. 6.However, these large rf flux amplitudes bring about the dangers of sample heating and amplifier saturation.
The response from the overlapping SQUID array SNAP161A is in clear contrast to the measurement on a single layer 12 × 12 × 1 SQUID array (SNAP161D) shown in Fig. 15.The single-layer sample shares the same design with the bottom layer of the two-layer sample (SNAP161A).Unlike the overlapping corner-coupled SQUID array, the single layer SQUID array has only one resonance band tuned with the applied dc flux Φ dc , which is consistent with the coherent response seen in earlier measurements of single-layer SQUID metamaterials [11,20,46,48].The resonant frequency is close to that of a single-SQUID resonance shown as the red dashed curves.At low applied rf flux amplitude, the resonance is above that of the single SQUID since the collective antiferromagnetic coupling reduces the induced flux on each SQUID and thus the effective inductance.According to Eq.( 16), a low L eff < L geo corresponds to a resonance higher than that of the single SQUID as observed experimentally in Fig. 15.As expected for the SQUID modes, their tunability in Φ dc also diminishes with increasing Φ rf , again observed in earlier work on single-layer SQUID metamaterials [11,20,46,48].There are several other interesting features of Fig. 15 worth mentioning.First, note that green blobs appear around the turning points where the resonance tuning curves in applied dc flux from the two adjacent periods meet.These can be attributed to the fact that the metaatom SQUID is strongly hysteretic (β rf > 1), resulting in multiple stable solutions.At points where multiple solutions cross in the frequency-dc-flux space, there can be enhanced resonant responses over a range of frequencies, accounting for these blobs.We also note that the tuning curves are noticeably asymmetric as a function of dc flux at low rf flux amplitudes, and become increasingly symmetric as the rf flux amplitude increases.The asymmetry in the dc tuning curve is related to the hysteresis in dc flux sweep of the rf SQUID metamaterial.For a measurement with decreasing dc flux, the asymmetric tilt of the tuning curve points to the opposite direction.This hysteresis is suppressed at higher applied rf flux amplitude because the strong oscillatory drive can overcome the barriers between local potential minima, preventing the system from becoming stuck in metastable states that are responsible for the hysteresis.

IV. DISCUSSION
Two driven side-by-side inductively-coupled rf SQUIDs have two eigenmodes of oscillation in the linearized case.When the two loops are partially overlapping with significant capacitive coupling, a third eigenmode of oscillation appears, which arises from the partial loop created by the partial overlap between the SQUID loop wires.This new closed circuit is completed by means of displacement currents.The same loop-to-oscillation-mode correspondence has been verified in the calculations for larger systems, i.e. four corner-coupled SQUIDs and the two by two by two array of overlapping SQUIDs.
More generally, for a system of many overlapping corner-coupled SQUIDs, the total number of elementary loops can be determined as follows.If we treat the SQUID loops as vertices and partial loops as edges connecting the corresponding vertices, we can represent the geometry of overlapping SQUIDs as a planar graph, as illustrated in Fig. 16.The most general system can be represented as a disconnected graph, as shown in Fig. 16 where a single SQUID in the lower right corner does not overlap with any other SQUID.Without loss of generality, we can examine each connected sub-graph individually and sum up the number of elementary loops.For each connected planar sub-graph with n vertices (SQUID loops) and m edges (partial loops), the simple cycles of the graph cor-respond to the extra-SQUID loops, the number of which is given by m−n+1, bringing the total number of elementary loops to m+n+m−n+1 = 2m+1 in the sub-graph.The number of unknowns in the voltage formalism for the sub-graph can be obtained as in Sec.II E 3, where 2m − n + 1 unknown voltages and the n phase differences from n junctions, give rise to 2m + 1 unknowns, equal to the number of elementary loops.Tallying up all the sub-graphs, we have shown that for any general system of overlapping SQUIDs, the voltage formalism can result in an exactly determined system that describes the dynamics of the SQUIDs.For example, a single SQUID by itself is a planar connected graph with 1 vertex and 0 edges.The dynamics are described by 1 unknown, the phase difference, and 1 equation of motion, the flux quantization condition for a single SQUID (Eq.( 4)).The N × N × 2 system in Sec.II E 3 is another special case of a system consisting of one planar connected graph.
Although, the extra-SQUID loops studied in the examples above are the smallest circuits corresponding to the very high frequency modes 12f geo ≈ 100 GHz, they can take on larger sizes in a more general geometry as shown in Fig. 16.In fact, the extra-SQUID loop highlighted in blue in Fig. 7 is the smallest realization permitted in our overlapping SQUID design.
The frequency and the tunability in applied dc magnetic flux of the resonance modes can be obtained from Eq.( 16), where the effective inductance depends on the current distribution in the corresponding resonance modes.The modes with very small |L eff | ≪ L geo will have a high resonant frequency Ω res ≈ L geo /L a,eff , and be largely independent of dc flux applied to the metamaterial.On the other hand, the modes with large effective inductance |L eff | ≫ L geo will have a low resonant frequency Ω res ≈ √ β rf cos δ, even lower than the single SQUID resonance, and will be strongly tunable with dc flux.The experimental data on dc-flux tuning of the 12 × 12 × 2 metamaterial is consistent with these expectations.
It should be noted that the only source of loss treated in this work arises from quasiparticle tunneling, which appears as the resistance R in the RCSJ model, and is incorporated into the fully nonlinear equations of motion through the parameter γ.An extension of this treatment would be inclusion of other sources of loss.One candidate source of loss arises from dielectric in the Josephson junction tunnel barrier, the coupling capacitors between SQUIDs, as well as dielectrics surrounding the superconducting wiring.These dielectrics are known hosts for electric-dipole two-level systems [89,90].
The treatment presented here assumes that all of the rf SQUIDs making up the metamaterial are nominally identical, having the same values of geometrical inductance, shunt capacitance, overlap capacitance, critical current, and junction resistance.We also assume that every SQUID experiences the same values of the externallyapplied dc and rf flux, and that the driving rf flux is at a single frequency.It would be interesting to see how the results of this work depend on variations in these quantities due to either statistical or systematic variation in space.
Although our model has mainly focused on the low rf flux amplitude linear limit, we can explore the nonlinear properties of the system at higher rf flux amplitudes from the full numerical solutions to the equations of motion.Thanks to the plethora of resonance bands in the overlapping SQUID metamaterial, one can employ them for broadband parametric amplification or intermodulation generation by harnessing this nonlinearity from the Josephson junctions.Now that capacitive coupling between flux-based superconducting meta-atoms has been established, one can ask whether the capacitive coupling can be varied?For example, there exist many nonlinear dielectric materials whose dielectric properties can be tuned with dc electric field, rf electric field, or with temperature at cryogenic conditions [91,92].In addition, charge qubits can enjoy variable capacitive coupling through a small Josephson junction [64].Thus a certain degree of tunability of capacitive coupling should be possible.

V. CONCLUSIONS
We consider three-dimensional rf SQUID metamaterials with strong capacitive coupling between rf SQUID loops for the first time.Strong displacement currents can flow through the capacitors created by such overlap of the SQUID loops, creating new closed paths for the rf current through the rf SQUID network.The RCSJ model is extended to incorporate the capacitive coupling and the new current paths, leading to the prediction of multiple resonances of comparable frequency to those of the galvanically-closed SQUID loops.The number of resonating loops in our three-dimensional N × N × 2 rf SQUID metamaterial design scales as 8N 2 − 8N + 3. A large N = 12 three-dimensional rf SQUID metamaterial is measured, and is found to behave in a qualitatively different manner from the corresponding singlelayer 12×12×1 metamaterial.The observed multiplicity of resonances are in good agreement with our theory of capacitively-coupled overlapping SQUIDs.Let's now examine the linear limit approximation for this problem.Consider the solutions in the following form ⃗ δ = ⃗ δ rf (t) + ⃗ δ dc , ⃗ δ rf (t) = ⃗ δ rf exp(iΩτ ), and u 1 (t) = û1 exp(iΩτ ), where ⃗ δ rf = ( δa , δb ).Substituting these expressions in the equation of motion Eq. ( 11) and Faraday's law, Eq. ( 6), and rearranging the excitation to the left hand side of the equations, one can obtain a linear system in ( δa , δb , û1 ) : The resonance condition for the system is det( ← → χ = 0), which is a sixth-order equation and can be solved for Ω to obtain the eigen solutions to the two corner-coupled SQUIDs.Both the eigenfrequencies and the nonlinear numerical solutions have been obtained and agree with the results in Secs.II B 5 and II C in the main text.As discussed in Sec.II E, the voltage formalism can be applied to any size of corner-coupled overlapping SQUIDs array.Here, we outline the procedure for setting up the equations of motion for a system of corner-coupled SQUIDs in any general geometry.The first step in studying the dynamics is to identify the independent variables, which are ⃗ δ from the junctions, ⃗ V from the overlapping capacitors, and their time derivatives.We can then express the induced fluxes in the system as where ← → L is the inductance matrix whose i, j the element describes the flux on the loop i induced by the segment j, and ← → I con is the matrix that expresses the currents in each branch in terms of junction currents ⃗ I JJ and the displacement currents through the capacitor nodes c ⃗ V .Next, the same flux relations in the conventional RCSJ model (Sec.II A) are invoked here for each SQUID which can be solved for junction currents ⃗ I JJ ( ⃗ δ, ⃗ V ).The non-SQUID loops, on the other hand, are described by Faraday's law in the following form

FIG. 1 :
FIG. 1: (a) Schematic diagram of a pair of corner-coupled rf SQUIDs with overlapping wiring layers, creating capacitors labeled 1 and 2, along with the two nominally identical Josephson junctions in loop a and loop b represented with "×"s.Junction currents I a0 and I b0 differ from the corresponding non-junction currents I a1 and I b1 , in general.(b) Micrograph of a small section of a 7×7×2 metamaterial made up of corner-coupled SQUIDs (fabricated by the SNEP process), where one representative pair of the capacitively-coupled SQUIDs is highlighted.The different colors of the wiring correspond to different lithographic layers of the device.

1 .
Loops in the two corner-coupled SQUIDs 2 in a closed loop d.The total flux in the closed loop c is simply Φ = ⃗ B • d ⃗ S = c ⃗ A • d ⃗ l where the magnetic field and the vector potential are due to the current in segment 2. We can assign the flux contribution from individual segments as follows: Φ = c ⃗ A•d ⃗ l = i si∈c ⃗ A•d ⃗ l.The partial inductance between segments 1 and 2 is then defined as the ratio, M 1,2 = 1 ⃗ A • d ⃗ l/I 2 .Consequently, it follows that for a closed loop c, M c,2 = si∈c M si,2 , and between two closed loops c and d, M c,d = sj ∈d M c,sj .

FIG. 2 :
FIG. 2: Real part of eigenfrequency solutions Re(Ω) from the characteristic equation det( ← → χ ) = 0 for the two corner-coupled SQUIDs, as a function of dc magnetic flux Φ dc .The parameters for this calculation are given in Appendix B and Table.II.The solid curves with three different colors correspond to the three positive solutions to the characteristic equation, while the black dotted curve is the eigenfrequency for a single SQUID with the same parameters.Due to the hysteretic response of the SQUIDs (β rf > 1), multiple dc flux tuning curves overlap each other in the same range of applied dc flux.

FIG. 3
FIG. 3: a), b) Real and imaginary parts of the solved dimensionless currents ι = 2πL geo I/Φ 0 between Ω = 2.1 and 2.7, for the linearized case of two corner-coupled SQUIDs at zero dc flux.c), d) Real and imaginary parts of the solved currents between Ω = 4.8 and 5.The solid curves are the solutions to the junction currents ι a0,b0 , while the dashed curves are the non-junction currents, ι a1,b1 .Blue curves are for the currents in loop a, and red for loop b.

FIG. 5 :
FIG. 5: Nonlinear solutions to the two corner-coupled SQUIDs at Φ rf ∼ 10 −3 Φ 0 .The parameters for this calculation are given in Appendix B and Table.II.The quantity plotted is |S 21 |(dB) on a logarithmic color scale, as a function of dimensionless frequency Ω = ω/ω geo and applied dc magnetic flux in units of the flux quantum Φ 0 .The white dashed curve corresponds to the eigenfrequency for a single SQUID with the same parameters, and the red curves are the eigenfrequencies from the linear limit solutions for the two corner-coupled SQUIDs as in Sec.II B 5.

FIG. 7 :
FIG. 7: Schematics for larger overlapping SQUID systems.a) The four corner-coupled SQUIDs constitute the minimal system that contains all three different types of loops involved in the dynamics of the gauge-invariant phases.The three types of circuits are galvanicallyconnected SQUID loops (black lines), the partial loops (highlighted in red ), and the extra-SQUID loops (highlighted in blue).b) represents a typical square array of N ×N ×2 SQUIDs (here with N = 2)).The SQUIDs from the two different layers are color coded to illustrate the geometry of this design.Overlap capacitors exist wherever lines of different color cross.An array with N = 12 has been characterized experimentally in this work.

FIG. 8 :
FIG.8: Real part of eigenfrequency solutions Re(Ω) from the characteristic equation det( ← → χ ) = 0 in the linear limit for the four corner-coupled SQUIDs shown in Fig.7 (a), as a function of dc magnetic flux Φ dc in units of Φ 0 .In this case 9 distinct resonance modes can be resolved.The two black horizontal dashed lines delineate the three types of resonant modes.

FIG. 10 :
FIG. 10: Real part of eigenfrequency solutions Re(Ω) from the characteristic equation det( ← → χ ) = 0 in the linear limit for the 2 × 2 × 2 corner-coupled SQUIDs shown in Fig. 7(b), as a function of dc magnetic flux Φ dc in units of Φ 0 .A total of 19 distinct resonance modes can be resolved.The two black horizontal dashed lines delineate the three types of resonant modes.
FIG. 14: Measured change in transmission ∆S 21 (green to yellow color) of the 12 × 12 × 2 overlapping cornercoupled SQUID array as a function of applied rf magnetic flux amplitude from 2 × 10 −4 Φ 0 to 5 × 10 −2 Φ 0 (and the corresponding microwave power incident on the sample in dBm), and rf driving frequency f from 15 to 26.5 GHz (Ω = 1.8 to 3.2), at a temperature of 4.6 K.The measurement was taken at Φ dc = −0.3Φ0 with two separate power sweeps due to the limited 40 dB dynamic range in the VNA power sweep function.The two vertical lines around -62 and -42 dBm are the artifacts from stitching the two sweeps together.

FIG. 16 :
FIG. 16: The correspondence between a general system of overlapping SQUIDs and its graph representation.The SQUID loops are represented by the vertices, and the partial loops connecting the neighboring SQUIDs by the edges.The simple cycles in the graph correspond to extra-SQUID loops in the SQUID array, as highlighted by green and blue arrows in the graph, and corresponding closed loops in the SQUID array.
Appendix D: A general model for a system of overlapping SQUIDs

TABLE II :
Design parameters for the overlapping SQUID metamaterial sample SNAP 161A.See Fig. 1(b) for definitions of d, w, and a SQUID .