Realization and modeling of rf superconducting quantum interference device metamaterials

We have prepared meta-atoms based on radio frequency superconducting quantum interference devices (RF SQUIDs) and examined their tunability with dc magnetic field, rf current, and temperature. RF SQUIDs are superconducting split ring resonators in which the usual capacitance is supplemented with a Josephson junction, which introduces strong nonlinearity in the rf properties. We find excellent agreement between the data and a model which regards the Josephson junction as the resistively and capacitively-shunted junction. A magnetic field tunability of 80 THz/Gauss at 12 GHz is observed, a total tunability of 56$%$ is achieved, and a unique electromagnetically-induced transparency feature at intermediate excitation powers is demonstrated for the first time. An RF SQUID metamaterial is shown to have qualitatively the same behavior as a single RF SQUID with regards to DC flux and temperature tuning.


I. INTRODUCTION
Metamaterials are artificially structured media designed to have electromagnetic properties not found in nature. These properties arise from both the structure of individual meta-atoms and the interactions between them, resulting in interesting collective behavior. A wide variety of properties and applications have been pursued using metamaterials, including negative index of refraction [1] [2] [3], super-resolution imaging [4] [5], cloaking [6] [7], transformation optics [8] [9], and perfect absorption [10].
Most of the applications invoking novel optical properties impose three stringent constraints. First, the metamaterial must have low attenuation. Features such as evanescent wave amplification [11] and negative refraction are strongly suppressed by even small amounts of loss [12] [13] [14] [15] [11] [16]. For example, the enhanced loss in metamaterials approaching the plasmonic limit has imposed a severe limitation on visible wavelength metamaterials composed of noble metal nano-structures [17] [18] [19] [20]. Secondly, the meta-atoms must be of deep sub-wavelength dimensions to achieve the metamaterial limit, as opposed to the photonic crystal limit. This has been an issue both in the visible and microwave regimes, where meta-atom sizes often approach the scale of the wavelength to minimize losses [21] [22]. Thirdly, it is desirable to make metamaterials that have textured properties in space (e.g. for cloaking or transformation optics), or that can be tuned and re-configured after the metamaterial has been fabricated [23]. The ability to tune the electromagnetic response over a wide range, and on short time scales, is desirable for applications such as * mctrep@umd.edu; these authors contributed equally † dmchang@umd.edu; these authors contributed equally software-defined radio [24] and filters for digital rf receivers [25] [26]. Superconducting metamaterials have been proposed to address all three of these constraints [27]. In addition, the quantum coherent nature of the superconducting state leads to qualitatively new phenomena such as fluxoid quantization and Josephson effects.
Here we focus on the tunability and reconfigurable nature of superconducting Josephson metamaterials.

A. Tunable Superconducting Metamaterials
Superconductors are fundamentally nonlinear due to the nonlinear Meissner effect [28] [29]. In general this intrinsic nonlinearity will be encountered near the limits of parameter space spanned by temperature, applied magnetic field, and applied current.
Superconducting metamaterials are tunable with temperature due to the temperature dependence of the superfluid density and kinetic inductance [30] [31]. This type of tunability has been exploited in various superconducting metamaterials with some success [32] [33] [34] [35] [36]. However, temperature tuning is slow since the thermal inertia of the meta-atoms can be large, even at low temperatures. Typical estimates for temperature tuning response times are on the order of 10 µs [37].
Applied currents can also be used to tune superconducting metamaterials; they can cause superfluid depairing, which increases the kinetic inductance [38]. For example, applied currents can tune sub-THz transmission of a metamaterial composed of a network of resonators connected by a superconducting wire loop [37]. However, applied currents often induce magnetic vortices in the superconductor before the de-pairing critical current is reached [39] [37]. These vortices move under the influence of the high-frequency currents, creating enhanced inductance and dissipation. The additional dissipation renders the superconductors less attractive for applications.
Magnetic field tuning of superconductors is attractive for applications because it can have a large effect on metamaterial properties. For example, superconducting split-ring resonators (SRRs) have been tuned by both DC and RF magnetic fields [40] [41]. A DC magnetic field was shown to add magnetic flux to a thin film Nb SRR, increasing its inductance and loss [40]. The RF magnetic field created enhanced RF screening currents at discrete locations in the SRR, resulting in enhanced inductance and dissipation as magnetic flux moved into and out of the superconducting film at high frequency [40] [42] [43] [44] [45] [46]. It was also found that a superconducting resonator exhibiting an analog of electromagneticallyinduced transparency showed a strong switching behavior at high excitation power [36], for similar reasons. However, the insertion of magnetic flux into superconducting materials is often too slow and too dissipative for tuning applications; a low-dissipation quantum effect would be better. For example, flux quantization could be used to discretely tune a superconducting meta-atom [47].The addition of the Josephson effect to superconducting metamaterials adds a mechanism of tunability that offers a large degree of high speed tunability and low dissipation.

B. Superconducting Meta-Atoms Employing the Josephson Effect
A superconductor can be described by a macroscopic phase-coherent complex quantum wavefunction Ψ = √ n s e iθ [48]. This wavefunction inherits its phase coherence from the underlying microscopic BCS (Bardeen, Cooper, Schrieffer) wavefunction describing the Cooper pairing of electrons in the metal. The absolute square of the macroscopic quantum wavefunction is the local superfluid density n s ∼ |Ψ| 2 [49]. When two superconductors are brought close together and separated by a thin insulating barrier, there can be tunneling of Cooper pairs between the two materials [50] [48]. This tunneling results in two types of Josephson effect.
The DC Josephson effect produces a DC current between the two superconductors which depends on the phase difference between their macroscopic quantum wavefunctions as A( r, t) · d l is the gauge-invariant phase difference between superconductors 1 and 2, A( r, t) is the magnetic vector potential in the region between the superconductors, I c is the critical (maximum) current of the junction, and Φ 0 = h 2e ∼ = 2.07 × 10 −15 Tm 2 is the flux quantum (h is Planck's constant and e is the electronic charge).
The AC Josephson effect relates a DC voltage drop across the junction to a time-varying gauge-invariant phase difference, and therefore to an AC current in the junction: dδ dt = 2π Φ0 V , where V is the DC voltage across the junction. In general the AC impedance of a Josephson junction contains both resistive and reactive components [51].
In the presence of both DC and AC currents, the inductance of the junction can be described approximately as [49] The Josephson inductance is a strong function of applied currents and fields. Tuning of Josephson inductance in a transmission line metamaterial geometry was considered theoretically by Salehi, Majedi, and Mansour [52] [53].

C. SQUID Metamaterials
In this work we focus on meta-atoms comprised of a superconducting loop interrupted by a single Josephson junction, commonly known as a radio frequency superconducting quantum interference device (RF SQUID). Since the RF SQUID meta-atom (shown in Fig. 1(b)) is a natural quantum analog of the SRR, we examine the electromagnetic properties of the SQUID near its selfresonant frequency.
The original purpose of the RF SQUID was to measure small magnetic fields and operate as a flux to frequency transducer [54] [55]. The first proposal to use an array of RF SQUIDs as a metamaterial was made by Du, Chen, and Li [56] [57]. Their calculation assumes the SQUID has quantized energy levels and considers the interaction of microwave photons with the lowest states of the SQUID potential. For small detuning of the photon frequency above the transition from the ground state to the first excited state, the medium presents a negative effective permeability. The frequency region of negative permeability is diminished by a non-zero de-phasing rate and negative permeability will disappear for de-phasing rates larger than a critical value.
RF SQUIDs interacting with classical electromagnetic fields were modelled by Lazarides and Tsironis [58]. They considered a two-dimensional array of RF SQUIDs in which the Josephson junction was treated as a parallel combination of resistance, capacitance, and Josephson inductance. Near resonance, a single RF SQUID can have a large diamagnetic response. An RF SQUID array displays a negative real part of effective permeability for a range of a frequencies and magnetic fields. The permeability is oscillatory as a function of applied magnetic flux, and is suppressed with applied fields that induce currents in excess of the critical current of the Josephson junction. Similar calculations, which included interactions between the RF SQUIDs, were carried out by Maimistov and Gabitov [59]. The response of a twodimensional RF SQUID metamaterial under an applied oscilating magnetic field was also numerically investigated by Lazarides and Tsironis [60]. The weakly cou- pled elements of the metamaterial are predicted to show bistability and synchronization.
Related work on a one-dimensional array of superconducting islands that can act as quantum bits (qubits) was considered by Rakhmanov, et al. [61]. When interacting with classical electromagnetic radiation, the array can create a quantum photonic crystal, which can support a variety of nonlinear wave excitations. A similar idea based on a SQUID transmission line was implemented to perform parametric amplification of microwave signals [62] [63]. A one-dimensional array of dc SQUIDs has been utilized as a tunable composite right/left-handed transmission line [64]. A dc-field tunable one-dimensional SQUID metamaterial embedded in a co-planar waveguide structure has been recently demonstrated [65] [66] [67].
It is expected that SQUID metamaterials will be able to perform functions similar to galvanically connected SQUID arrays which have been developed in recent years. One possible application takes advantage of the SQUID's extreme sensitivity to magnetic flux to create compact, wideband antennas, sensitive to high frequency magnetic fields, as opposed to electric fields [68] [69] [70] [71]. Other possible applications include low noise amplifiers for RF sensing and qubit readout [72] [73] [63], and highly sensitive magnetometers and filters [74] [75] [76].
Our objective in this paper is to demonstrate the first SQUID metamaterial in a geometry that can be easily extended to a fully three-dimensional structure which can interact with free-space electromagnetic waves. We demonstrate that RF SQUID meta-atoms are remarkably tunable with very small DC magnetic fields, as well as temperature and RF currents.
The remainder of the paper is organized as follows. We first present models of RF SQUID linear and nonlinear behavior. Then we discuss the design and fabrication of a single RF SQUID meta-atom and a dense array of such SQUIDs. Section V presents our experimental results of the RF SQUID meta-atom and array in waveguide geometries, demonstrating tunability of the meta-atom properties with DC flux, temperature, and RF flux. Section VI is a comparison of data and model predictions, and discussion of the results, and Section VII serves as a conclusion.

II. MODELING
An ideal Josephson junction can be treated as an inductor with a variable inductance given by Eq. 1. A more realistic model (valid for low frequencies and currents) shunts the junction with a sub-gap resistance R (representing the tunneling of normal state electrons across the junction) and capacitance C (the capacitance of two overlapping conductors separated by an insulator) [49]. We model the RF SQUID as a resistively and capacitively shunted junction (RCSJ) in parallel with an inductor representing the inductance, L, of the superconducting loop. The result is an RLC circuit (see circuit diagram in Fig.  1(d)) with a resonant frequency given by where T is the temperature and Φ app is the magnetic flux applied to the SQUID loop. The quality factor for this parallel RLC circuit is proportional to the resonant frequency and given by, To understand the behavior of the resonance as a function of temperature and applied flux it is necessary to solve for the time-dependent gauge-invariant phase difference δ(T, Φ app ).
We use the RCSJ model to solve for δ subject to the condition that the total flux through the loop must be an integer number (n) of flux quanta, i.e. Φ = nΦ 0 . The total flux, Φ, is related to the applied flux, Φ app , and the current-induced flux by where I is the current through the loop, which is the sum of the currents through the junction, the resistor, and the capacitor. We can rewrite Eq. 4 in the terms of δ as where Φ DC is the DC flux bias of the applied magnetic field and Φ rf and ω are the amplitude and angular frequency of the applied rf field respectively. By making the substitution Φ = Φ0δ 2π , we have taken n = 1 because the choice of integer shifts δ (which has 2π periodicity) by 2π. Equation 5 can be recast in non-dimensional form as where Φ0 , and f DC = ΦDC Φ0 . A well-studied mechanical analogue of the Josephson junction is the driven and damped pendulum [77]. Equation 6 resembles the nonlinear differential equation that governs that system. However, incorporating the junction into a superconducting loop subject to flux quantization introduces an additional term linear in δ on the right-hand side of Eq. 6. Unlike the driven, damped pendulum, this equation is linear in the high power limit in addition to the low power limit. For our chosen parameters, δ is always periodic in time and its variation is dominated by the same frequency as the driving flux. We do not observe chaos for parameter values relevant to the experiments discussed below.
At intermediate powers the full non-linear Eq. 6 must be solved for δ, and therefore f 0 (T, f DC , f rf ), but the equation can be linearized and simplified in the low and high power limits. In the low rf power limit (where it is assumed that the time-varying component of the gauge invariant phase difference is very small i.e. δ rf (τ ) << 1) Eq. 6 can be linearized by separating the phase difference into a DC and rf components, i.e.
Equation 6 simplifies to the following time-independent and time-dependent equations: where α = 1 + β rf cos δ DC . In this case the DC flux dictates δ DC and Eq. 9 has an analytic solution for δ rf given by From this result one finds that δ and consequently f 0 (T, f DC ) are periodic functions of DC flux, f DC , with maximum and minimum values at f DC = n, and f DC = n + 1/2 respectively, where n is any integer. The resonant frequency of the SQUID, f 0 , is also temperature dependent due to the temperature dependence of the critical current, I c . As temperature increases the maximum frequency decreases and the minimum frequency increases resulting in a reduction of total flux tunability of the SQUID. When the temperature reaches the critical temperature of niobium (T c = 9.2 K), |L JJ | is expected to diverge (since I c → 0). In this case the Josephson junction loses the ability to modify the resonance and the resonant frequency reduces to This resonance depends solely on the non-junction properties of the SQUID and is insensitive to applied flux, Φ app , and temperature, T . Equation 6 is also linear in the high power limit (where it is assumed δ rf >> 1 since β rf sin δ+δ ≈ δ for β rf < 1), where it separates into the following time dependent and time independent equations: 2πf rf sin (Ωτ ) = δ rf + 1 Q geo Equation 13 has the same analytic solution as Eq. 10 with α replaced by 1, so unlike Eq. 10 it has no dependence on temperature or DC flux. In the high power limit the resonant frequency reduces to the same value as in the high temperature limit, Eq. 11.

III. DESIGN AND FABRICATION
We have designed and measured both a single RF SQUID meta-atom and a dense 27x27 array of metaatoms, which were manufactured using the Hypres 0.3 µA/µm 2 Nb/AlO x /Nb junction process on silicon substrates. However, measurements of other junctions from this run suggest the critical current density is closer to 0.2 µA/µm 2 at 4.2 K. The superconducting loop is composed of two Nb films (135 nm and 300 nm thick) which are connected by a via and a Nb/AlO x /Nb Josephson junction (see Fig. 1(b)). There is additional capacitance where these layers overlap (with SiO 2 dielectric) which is necessary to bring the resonant frequency within the measurable range. When designing the SQUIDs we have control over the inner and outer radius of the loop and thus the loop inductance L, the critical current of the junction I c , and the overlap capacitance C. Values for L, I c , and C were chosen to maximize tunability within the measurable frequency range 6.5-22 GHz (dictated by the available waveguides) while remaining in the low noise (Γ = 2πkB T Φ0Ic < 1 and L F = 1 kB T Φ0 2π 2 >> L [78]) and non-hysteretic (β rf = 2πLIc Φ0 < 1) limits. For the single meta-atom shown in Fig. 1(b), the inner and outer radii of the SQUID loop are 100 µm and 400 µm respectively with 3 µm diameter holes in the film every 10 µm to pin vortices [79]. We estimate L to be 0.33 nH based on FastHenry calculations [80]. The area of the Josephson junction and the area of the capacitor are 5.3 µm 2 and 600 µm 2 respectively, yielding nominal design values of I c (4.2 K) = 0.97 µA and C = 0.32 pF. Figure 1(c) shows a portion of a 27x27 array of nominally identical RF SQUID meta-atoms with a center to center separation of 83 µm. The inner and outer radii of the SQUID loops are 30 µm and 40 µm respectively. We estimate L to be 0.12 nH. The area of the Josephson junction and the area of the capacitor are 13 µm 2 and 1800 µm 2 respectively, yielding nominal design values of I c (4.2 K) = 3.7 µA and C = 0.84 pF.

IV. EXPERIMENT
The single RF SQUID (and later the 27x27 RF SQUID array) is oriented in a 7.6 cm long Ku (K) rectangular waveguide which has a single propagating mode operating frequency range from 9.5-19 GHz (15-26 GHz), as shown in Fig. 1(a). An attenuated microwave signal from an Agilent E8364C network analyzer excites a TE 10 mode in the waveguide producing an rf magnetic field, B rf , perpendicular to the plane of the SQUID. The transmission, S 21 , is measured after a cryogenic low noise amplifier (LNA) and a room-temperature (RT) amplifier. The experiment is conducted in a three-stage pulsed-tube cryostat with a sample base temperature of 6.5 K. The sample is magnetically shielded by a mu-metal cylinder and a superconducting niobium open cylinder inside the cryostat. Superconducting coils surrounding the waveguide bias the SQUID by generating a perpendicular DC magnetic field, B DC in the range from 0 to ±1 µT.
The resonance is detected as a dip in the frequency dependent transmission magnitude through the waveguide |S 21 (ω)|. We have considered two methods for calculating scattering parameters from the solution for δ(t). One method is to consider the SQUID an effective medium with an effective relative permeability, µ r [81] [65], where Φ ac is the AC flux response of the loop, the angle brackets represent time averaging, and F is the filling fraction of the SQUID in the medium [58] [57]. S 21 is proportional to the ratio of the transmitted electric field, E T , to the incident field E 0 , and can be calculated using the continuity of E and H fields at the boundaries of the effective medium and the empty waveguide as where l is the length of the medium, γ = k µr k0 , k = µ r ω c 2 − π a 2 is the wave number in the medium, and k 0 = k(µ r = 1) is the wave number in the empty waveguide. For the single SQUID (as opposed to the array) the choice of l and F is not straightforward. They are used as fitting parameters (along with resistance R) for the width and depth of the measured resonances. The important parameter for the fit is the relationship between F and l, which we find to be F = A loop 0.79l , where A loop is the area of the SQUID loop (calculated from the average radius).
An alternate method for estimating S 21 focuses on the power dissipated in the resistor, R in Fig. 1(d).
where P T is the transmitted power and P 0 is the incident power. This assumes that the only power not transmitted through the waveguide is the power dissipated in the resistor. It does not account for reflection or other loss mechanisms.

V. RESULTS
The resonance of an RF SQUID responds to three tuning parameters: DC magnetic field, temperature, and rf power. Modifying the resonance with DC magnetic field is the most straightforward. Figure 2 shows the experimental results for |S 21 | of an RF SQUID meta-atom as a function of DC flux and frequency. The incident rf power and the temperature were fixed at -80 dBm and 6.5 K, respectively. Resonance dips in |S 21 (ω)| appear as the red features against a yellow background of unaffected signals (|S 21 | = 0dB). The signal was extracted by subtracting |S 21 | at 16 K (which is well above the critical temperature) from |S 21 | at 6.5 K, removing a background variation, and applying a threshold to identify the resonance. The resonance shows good periodicity with DC flux, with a maximum value of 16.9±0.3 GHz and a minimum at or below 10 GHz. The cutoff frequency of the Ku waveguide is 9.5 GHz which imposes a lower frequency limit on this measurement. The same sample measured in an X band waveguide (single propagating mode operating frequency range from 6.6-13 GHz) has a minimum resonant frequency of 9.5 ± 0.5 GHz. The small magnitude of the resonance dips can be attributed to the small size of the SQUID relative to the waveguide. This becomes less of a factor when considering arrays of SQUIDs which occupy a more appreciable fraction of the waveguide cross section. The quality factor, Q, is extracted by fitting the |S 21 (ω)| data to a model of an RLC-resonator inductively coupled to a transmission line [82]. Q also shows periodicity with DC flux, following the same trend as the resonant frequency (Eq. 3), with a maximum value of 54 and minimum at or below 20.
At an rf power of -80 dBm, (Φ rf ≈ 0.003Φ 0 ) the low power linearized equation (Eq. 9) is valid (since δ rf << 1). Treating the RF SQUID as an effective medium in the waveguide, |S 21 (ω, Φ DC )| is calculated from Eqs. 8, 10, 14, and 15 and plotted as the contour lines in Fig.  2. The capacitance C and critical current I c (6.5 K) are adjusted from their nominal values and taken to be 0.42 pF and 0.55 µA, respectively. These values were chosen so the resonant frequency at T = 6.5 K and Φ DC = 0 agree with the measurements in both the low and high power limits i.e. f 0 (6.5 K, Φ DC = 0) = 16.8 GHz for low power (Eq. 2) and f 0,geo = 13.5 GHz (Eq. 11) for high power. The model and the data agree both in terms of the flux dependence and magnitude of the tunability (∼ 7 GHz), and the magnitude of the resonance dip.
The DC flux tunability is modified by changes in temperature. The red features in Fig. 3 show the DC flux tuned resonance at 6.5 K, 7.6 K, and 8.3 K for an rf power of -80 dBm. The flux tunability is reduced from ∼ 7 GHz at 6.5 K, to ∼ 3 GHz at 7.6 K, and ∼ 1 GHz at 8.3 K. The quality factor, Q, is 50 at Φ DC = 0 and does not show significant variation with temperature. The grey lines follow the resonance by using the solution for δ DC from Eq. 8 and calculating the resonant frequency using Eqs. 1 and 2 (assuming δ ≈ δ DC ). The reduction in tunability with increased temperature in consistent with the model. The only parameter varied to fit this data is the temperature dependent critical current.
We also measured |S 21 | as a function of frequency and dc flux at low power for various fixed temperatures ranging from 6.5 to 9.2 K (not shown). For f DC = 0 the resonant frequency decreased from 17 to 13.5 GHz with increasing temperature and for f DC = 1/2 it increased from 9.5 to 13.5 GHz over the same temperature range. Consistent with the model, the resonant frequency saturates in the high temperature limit at f 0,geo = 13.5 ± 0.2 GHz (Eq. 11).
The temperature dependence of DC flux tuning for the 27x27 RF-SQUID array is shown in Fig. 4. The res- onance (red features) shows good periodicity with DC flux with a maximum of 21.5 ± 0.2 GHz and a minimum below 15 GHz at 6.5 K and -60 dBm rf power. As temperature increases, the flux tunability is reduced from > 6.5 GHz for 6.5 K to 2.5 GHz at 7.9 K. A flux offset of about 0.1Φ 0 shifts the curve along the flux axis without affecting the periodicity. The periodic horizontal features result from standing waves in the system due to an impedance mismatch in the external circuit.
The meta-atoms in the array respond to the tuning parameters coherently, creating a stronger signal and a slightly wider resonance dip compared to a single RF SQUID. The flux dependence of the array is qualitatively the same as that of a single SQUID. This is demonstrated by the grey lines in Fig. 4 which are solutions for the resonant frequency of a single SQUID with the following parameters: L = 0.12 nH (nominal value), C = 0.65 pF, I c (6.5K)= 1.2 µA, and I c (7.9K)= 0.47 µA.
The resonance of an RF SQUID can also be modified by current induced by Φ rf . Figure 5 shows |S 21 | as a function of rf power and frequency for a single RF SQUID meta-atom in the representative case of f DC = 1/6 and T = 6.5 K. The resonant response is represented by the red features. When the rf power is low, the resonant frequency remains at a fixed value, 16 GHz, and the resonance dip has a constant depth. As the rf power increases, the resonant frequency decreases and the resonance dip becomes shallower. The resonance disappears at intermediate powers before reappearing as a strong resonance at high power (-55 dBm).
Calculating the rf power dependence of the resonance requires a full numerical solution of the nonlinear Eq. 6. Instead of regarding the RF SQUID as an effective medium, we employ Eq. 16 to estimate |S 21 |. The results of this calculation are plotted as contour lines in Fig.  5 and the calculated minimum of |S 21 | is plotted as a dashed line. Both show excellent agreement with the data, including the depth of the |S 21 | dip.
There are two transition points in the rf power dependent resonance at −75dBm (Φ rf ≈ 0.006Φ 0 ) and -57 dBm (Φ rf ≈ 0.05Φ 0 ). We can understand the significance of these powers and the behavior in the low, intermediate, and high power regions using the model discussed in Sec. II. In the low power limit, the resonance does not change with increasing rf power because S 21 from Eq. 16 does not depend on input rf power, P 0 , since P 0 ∝ V 2 (V ∝ f rf from Eq. 10 and f rf ∝ √ P 0 because of the relationship between input power and magnetic fields in a waveguide).
The first transition occurs at an rf power of −75 dBm, where the resonant frequency begins to decrease and the resonance dip becomes shallower. At this rf power the oscillation amplitude of δ during an rf period is π/2 on resonance. The resonant frequency is modified because above −75 dBm rf power, the time averaged phase difference δ deviates from δ DC . In the case of f DC = 1/6, δ > δ DC resulting in a decrease in cos δ from its low power value, a corresponding increase in L JJ (Eq. 1), and an ultimate decrease in resonant frequency f 0 (Eq. 2). However in the case of f DC = 1/2, the modification in δ causes an increase in cos δ from its low power value and an ultimate increase in resonant frequency, f 0 . This is consistent with the experimental results (not shown). The resonance dip becomes shallower because V is no longer proportion to f rf ; instead the input power, P 0 , increases more than the dissipated power, V 2 /R.
The second transition occurs at an rf power of −57 dBm where δ oscillates with an amplitude exceeding 2π on resonance and the Josephson junction undergoes multiple phase slips in each rf period. This is the high rf power limit (Eq. 13), where the tunable resonant frequency reduces to a fixed value (Eq. 11), i.e. f 0,geo = 13.5 GHz. We expect the losses to be greater in this regime because the phase slips are dissipative, resulting in a deeper resonance dip, which is clearly evident in the data.
A similar evolution of a transition frequency with RF power has been observed in single qubits coupled to microwave cavities containing a small number of photons [83] [84]. The Janes-Cummings Hamiltonian shows such behavior at high photon excitation number, reproducing the observed dispersion and amplitude variation of the transition with increasing RF power [85].

VI. DISCUSSION
The design parameters of the SQUID meta-atom are the critical current I c (T ), the capacitance C, the sub-gap resistance R, and the loop inductance L. In our comparison of data and model for the single meta-atom we use FIG. 6. Critical current temperature dependence Ic(T ) calculated from measured resonant frequencies for a single SQUID at zero DC flux and corresponding fixed temperatures the nominal design value for loop inductance because the lithographic process dictates the loop geometry to high precision. We adjusted the capacitance C and critical current I c (6.5 K) from their nominal design values to fit the experimentally observed resonant frequencies at P = −80 dBm using Eq. 2 and P = −50 dBm using Eq. 11 for T = 6.5 K and Φ DC = 0. The fit value for capacitance is 30% higher than the design value, which exceeds the 20% tolerance quoted for the Nb/AlO x /Nb process [86]. (The same method was used to determine I c (6.5 K), I c (7.9 K), and C for the fit to the array data in Fig. 4.) The sub-gap resistance R = 820 Ω, and the filling fraction F and medium length l in Eqs. 14 and 15 were adjusted to match the width and depth of the observed |S 21 (ω)| minimum at Φ DC = 0, T = 6.5 K, and -80 dBm. With these parameters we quantitatively explain the DC flux and rf power dependence of the RF SQUID meta-atom.
We extracted the temperature dependence of the critical current I c (T ) from the flux dependence of the resonant frequency at different fixed temperatures. The critical current of the junction at each temperature was calculated by substituting the maximum resonant frequency as a function of flux into Eqs. 1 and 2. The results for I c (T ), shown in Fig. 6, are consistent with previous results on Nb/AlO x /Al/Nb tunnel junctions [87].
We have demonstrated tunability of the RF SQUID meta-atom with temperature, RF current, and DC magnetic flux. The SQUID meta-atom has an unusual response to changes in temperature. The resonant frequency can either increase or decrease, depending on the DC magnetic flux applied to the meta-atom; the slope of temperature tuning of f 0 can be large or small, and either positive or negative.
Tuning with rf power can also either increase or decrease the resonant frequency, depending on the applied DC magnetic flux. However, the strength of the resonance response varies with rf power while it remains constant with temperature tuning. At low pow-ers (< −75 dBm), and high powers (> −57 dBm), the meta-atom has a resonant interaction with electromagnetic radiation, but at intermediate powers it becomes essentially transparent. This is similar to demonstrations of metamaterial-induced transparency in which a tunable "spectral hole" is created by interfering resonant processes in two or more meta-atoms making up a meta-molecule [88] [89] [90] [91]. There, the transparency condition in the superconducting meta-molecule could be suppressed in a switching transition at high power [36]. In our case, the transparency process is due to the nonlinear Josephson effect and is self-induced, making transparency in the RF SQUID meta-atom simpler than previous implementations. Also unlike other realizations, the transparency condition arises from a decrease in dissipation without enhancement in loss at nearby frequencies [92]. As such it is potentially useful as a power limiter for sensitive front-end receivers [93].
The DC flux tuning of the SQUID meta-atom is remarkably sensitive. At its maximum value, the flux tunability (defined as the frequency change divided by the change in magnetic field) is approximately 80 THz/Gauss, at 12 GHz and 6.5 K. Hence only extremely small magnetic fields are required to make very substantial changes in the properties of the meta-atom.
The next question is how quickly can tuning by each parameter be achieved. The shortest time scale for tuning a superconductor is ∼ /∆ ∼ 1 ps, where ∆ is the energy gap. In this case, the RC time constant of 0.3 ns (L/R time constant 0.4 ps) imposes an upper limit on the intrisic switching speed. Tuning with temperature depends on changes to the critical current of the junction and is relatively slow, with typical estimates on the order of 10 µs [37]. A pulsed rf power measurement of the SQUID meta-atom is consistent with a response time less than 500 ns. Hence the rf power tuning time is in the sub-µs range perhaps only limited by the RC time. Flux tuning of SQUID-like superconducting qubits has been accomplished on nano-second time scales limited only by the rise-time of the applied current pulse [94] [95].
The tunability of the meta-atom in normalized units (defined as tuning frequency range divided by the center frequency of tuning) is 56%. The meta-atom loss varies by about 64% (defined as the loss range divided by maximum loss) over this entire range. We estimate the figure of merit (defined as the ratio of the real part of µ to the imaginary part) to be 1700 on resonance in the case of low power and zero DC flux. The goal of achieving wide meta-atom tunability on short time scales with low loss is well within reach.
The coherent response of the meta-atoms in the 27x27 array demonstrates the feasiblity of a widely tunable RF SQUID metamaterial. The SQUIDs are close enough that their interactions play a role in their dynamics. We define a coupling parameter which is the ratio of mutual inductance to self-inductance; κ ∼ = µ0πr 4 4d 3 L = 0.02, where r is the radius of the loop and d is the center-to-center spacing of the SQUIDs. The ratio of the wavelength to the lattice parameter for the measured array is ∼ 350, which is well-within the meta-material limit. The array displays an overall tunability of 46% and a flux tunability of 4 THz/Gauss at 19 GHz and 6.5K. Future work will further explore arrays of RF SQUID meta-atoms in both the interacting and non-interacting limits.

VII. CONCLUSION
We have designed, fabricated, and presented an rf superconducting quantum interference device meta-atom in a geometry that can be scaled for free-space interactions with electromagnetic fields. The meta-atom proves to be highly tunable with dc magnetic field, rf currents, and temperature. The RF SQUID meta-atom proves to have low losses over its wide (56%) range of tunability. A novel form of electromagnetically-induced transparency has been observed in this meta-atom with increasing in-duced rf current. We have obtained a detailed quantitative understanding of this nonlinear meta-atom with first principles modeling. The data and model are in excellent agreement. Finally, a two-dimensional array of RF SQUIDs has shown coherent tuning with DC flux and temperature in a manner similar to a single SQUID.