# Superconducting Metamaterials

###### Abstract

Metamaterials, i.e. artificial, man-made media designed to achieve properties not available in natural materials, have been the focus of intense research during the last two decades. Many properties have been discovered and multiple designs have been devised that lead to multiple conceptual and practical applications. Superconducting metamaterials made of superconducting metals have the advantage of ultra low losses, a highly desirable feature. The additional use of the celebrated Josephson effect and SQUID (superconducting quantum interference device) configurations enrich the domain of superconducting metamaterials and produce further specificity and functionality. SQUID-based metamaterials are both theoretically investigated but also fabricated and analyzed experimentally in many laboratories and exciting new phenomena have been found both in the classical and quantum realms. The enticing feature of a SQUID is that it is a unique nonlinear oscillator that can be actually manipulated through multiple external means. This domain flexibility is inherited to SQUID-based metamaterials and metasurfaces, i.e. extended units that contain a large arrangement of SQUIDs in various interaction configurations. Such a unit can be viewed theoretically as an assembly of weakly coupled nonlinear oscillators and as such presents a nonlinear dynamics laboratory where numerous, classical as well as quantum complex, spatio-temporal phenomena may be explored. In this review we focus primarily on SQUID-based superconducting metamaterials and present basic properties related to their individual and collective responses to external drives; the work summarized here is primarily theoretical and computational with nevertheless explicit presentation of recent experimental works. We start by showing how a SQUID-based system acts as a genuine metamaterial with right as well as left handed properties, demonstrate that the intrinsic Josephson nonlinearity leads to wide-band tunability, intrinsic nonlinear as well as flat band localization. We explore further exciting properties such as multistability and self-organization and the emergence of counter-intuitive chimera states of selective, partial organization. We then dwell into the truly quantum regime and explore the interaction of electromagnetic pulses with superconducting qubit units where the coupling between the two yields phenomena such as self-induced transparency and superradiance. We thus attempt to present the rich behavior of coupled superconducting units and point to their basic properties and practical utility.

###### keywords:

Superconducting metamaterials, nonlinear metamaterials, superconducting quantum metamaterials, dissipative breathers, chimera states, flat-band localization, self-induced transparency, superradiance, superconducting qubits63.20.Pw, 11.30.Er, 41.20.-q, 78.67.Pt, 05.65.+b, 05.45.Xt, 78.67.Pt, 89.75.-k, 89.75.Kd, 74.25.Ha, 82.25.Dq, 63.20.Pw, 75.30.Kz, 78.20.Ci

^{†}

^{†}journal: Physics Reports\biboptions

super

## Contents

Abstract | 1 | ||

Contents | 2 | ||

1. | Introduction | 3 | |

1.1 | Metamaterials Synthetic Media: Concepts and Perspectives. | 3 | |

1.2 | Nonlinear, Superconducting, and Active Metamaterials. | 4 | |

1.3 | Superconducting Metamaterials from Zero to Terahertz Frequencies. | 5 | |

1.4 | Summary of Prior Work in Superconducting Metamaterials. | 6 | |

1.5 | SQUID Metamaterials. | 9 | |

2. | SQUID-Based Metamaterials I: Models and Collective Properties | 11 | |

2.1 | The rf-SQUID as an artificial magnetic ”atom”. | 11 | |

2.2 | SQUID Metamaterial Models and Flux Wave Dispersion. | 18 | |

2.3 | Wide-Band SQUID Metamaterial Tunability with dc Flux. | 23 | |

2.4 | Energy Transmission in SQUID Metamaterials. | 26 | |

2.5 | Multistability and Self-Organization in Disordered SQUID Metamaterials. | 30 | |

3. | SQUID-Based Metamaterials II: Localization and Novel Dynamic States | 36 | |

3.1 | Intrinsic Localization in Hamiltonian and Dissipative Systems. | 36 | |

3.2 | Dissipative Breathers in SQUID Metamaterials. | 36 | |

3.3 | Collective Counter-Intuitive Dynamic States. | 39 | |

3.4 | Chimera States in SQUID Metamaterials. | 40 | |

3.4.1 SQUID Metamaterials with Non-Local Coupling. | 40 | ||

3.4.2 SQUID Metamaterials with Local Coupling. | 46 | ||

4. | SQUID Metamaterials on Lieb Lattices. | 51 | |

4.1 | Nearest-Neighbor Model and Frequency Spectrum. | 51 | |

4.2 | From flat-Band to Nonlinear Localization. | 52 | |

5. | Quantum Superconducting Metamaterials. | 56 | |

5.1 | Introduction. | 56 | |

5.2 | Superconducting Qubits. | 57 | |

5.3 | Self-Induced Transparency, Superradiance, and Induced Quantum Coherence. | 59 | |

5.3.1 Description of the Model System. | 59 | ||

5.3.2 Second Quantization and Reduction to Maxwell-Bloch Equations. | 60 | ||

5.3.3 Approximations and Analytical Solutions. | 62 | ||

5.3.4 Numerical Simulations. | 64 | ||

6. | Summary | 69 | |

Acknowledgements | 70 | ||

Appendix: Derivation of the Maxwell-Bloch-sine-Gordon equations. | 71 | ||

References | 76 |

## 1 Introduction

### 1.1 Metamaterials Synthetic Media: Concepts and Perspectives

Metamaterials represent a new class of materials generated by the arrangement of artificial structural elements, designed to achieve advantageous and/or unusual properties that do not occur in natural materials. In particular, naturally occurring materials show a limited range of electrical and magnetic properties, thus restricting our ability to manipulate light and other forms of electromagnetic waves. The functionality of metamaterials, on the other hand, relies on the fact that their constitutive elements can be engineered so that they may achieve access to a widely expanded range of electromagnetic properties. Although metamaterials are often associated with negative refraction, this is only one manifestation of their possible fascinating behaviors; they also demonstrate negative permittivity or permeability, cloaking capabilities Schurig et al. (2006), perfect lensing Pendry (2000), high frequency magnetism Linden et al. (2006), classical electromagnetically induced transparency Papasimakis et al. (2008); Kurter et al. (2011a); Jin et al. (2013); Zhang et al. (2017a), as well as dynamic modulation of Terahertz (THz) radiation Li et al. (2017), among other properties. High-frequency magnetism, in particular, exhibited by magnetic metamaterials, is considered one of the ”forbidden fruits” in the Tree of Knowledge that has been brought forth by metamaterial research Zheludev (2010). Their unique properties allow them to form a material base for other functional devices with tuning and switching capabilities Zheludev (2010, 2011); Tong (2018). The scientific activity on metamaterials which has exploded since their first experimental demonstration Smith et al. (2000); Shelby et al. (2001), has led to the emergence of a new, rapidly growing interdisciplinary field of science. This field has currently progressed to the point where physicist, material scientists and engineers are now pursuing applications, in a frequency region that spans several orders of magnitude, from zero Wood and Pendry (2007); Magnus et al. (2008); Navau et al. (2009); Mawatari et al. (2012); Mach-Batlle et al. (2017) to THz Yen et al. (2004); Linden et al. (2004); Withayachumnankul and Abbott (2009); Gu et al. (2010); Jin et al. (2010); Chen et al. (2010); C. H. Zhang et al. (2013) and optical Linden et al. (2006); Shalaev (2007); Soukoulis et al. (2007); Litchinitser and Shalaev (2008). Historically, the metamaterial concept goes back to 1967 Veselago (1967), when V. Veselago investigated hypothetical materials with simultaneously negative permeability and permittivity with respect to their electromagnetic properties. He showed that simultaneously negative permeability and permittivity result in a negative refractive index for such a medium, which would bend the light the ”wrong” way. The realization of materials with simultaneously negative permeability and permittivity, required for negative refractive index, had however to wait until the turn of the century, when D. Smith and his collaborators demonstrated for the first time a structure exhibiting negative refraction in the microwaves Smith et al. (2000). The first metamaterial was fabricated by two interpenetrating subsystems, one them providing negative permittivity while the other negative permeability within the same narrow frequency band. Specifically, an array of thin metallic wires and an array of metallic rings with a slit (split-ring resonators), which were fabricated following the ”recipies” in the seminal works of J. B. Pendry, provided the negative permeability Pendry et al. (1999) and the negative permittivity Pendry et al. (1996), respectively. The wires and the split-rings act as electrically small resonant ”particles”, undertaking the role of atoms in natural materials; however, they are themselves made of conventional materials (highly conducting metals). Accordingly, a metamaterial represents a higher level of structural organization of matter, which moreover is man-made.

The key element for the construction of metamaterials has customarily been the split-ring resonator (SRR), which is a subwavelength ”particle”; in its simplest version it is just a highly conducting metallic ring with a slit. The SRR and all its subsequent versions, i.e., U particles, H particles, or like particles, double and/or multislit SRR molecules, are resonant particles which effectively act as artificial ”magnetic atoms” Caputo et al. (2012). The SRRs can be regarded as inductive-resistive-capacitive () oscillators, featuring a self-inductance , a capacitance , and a resistance , in an electromagnetic field whose wavelength much larger than their characteristic dimension. As long as a metamaterial comprising SRRs is concerned, the wavelength of the electromagnetic field has to be much larger than its unit cell size; then the field really ”sees” the structure as a homogeneous medium at a macroscopic scale and the macroscopic concepts of permittivity and permeability become meaningful. The (effective) homogeneity is fundamental to the metamaterial concept, as it is the ability to structure a material on a scale less than the wavelength of the electromagnetic field of interest. Although in microwaves this is not a problem, downsizing the scale of metamaterial elements to access the optical frequency range may be a non-trivial issue. The advent of metamaterials has led to structures with many different designs of elemental units and geometries, that may extend to one Shamonina and Solymar (2004); Butz et al. (2013a), two Shelby et al. (2001); Mawatari et al. (2012), or three dimensions Zagoskin (2012). One of the most investigated metamaterial designs which does not contain SRRs is the fishnet structure and its versions in two Kafesaki et al. (2007), quasi-two Wuestner et al. (2010), and three dimensions Liu et al. (2008); Valentine et al. (2008). However, all these metamaterials have in common that they owe their extraordinary electromagnetic properties more to their carefully designed and constructed internal structure rather than, e.g., chemical composition of their elements. Metamaterials comprising of split-rings or some other variant of resonant elements, are inherently discrete; discreteness effects do not however manifest themselves as long as the metamaterial responds linearly (low-field intensities) and the homogeneous medium approximation holds. The coupling effects, however, in relatively dense SRR metamaterials are of paramount importance for a thorough understanding of certain aspects of their behavior, since they introduce spectral splitting and/or resonant frequency shifts Gay-Balmaz and Martin (2002); Hesmer et al. (2007); Penciu et al. (2008); Liu et al. (2009); Sersić et al. (2009); Liu and Giessen (2010); Feth et al. (2010). The SRRs are coupled to each other through non-local magnetic and/or electric dipole-dipole interaction, with relative strength depending on the relative orientation of the SRRs in an array. However, due to the nature of the interaction, the coupling energy between neighboring SRRs is already much less than the characteristic energy of the metamaterial; thus in most cases next-nearest and more distant neighbor interactions can be safely neglected. SRR-based metamaterials support a new kind of propagating waves, referred to as magnetoinductive waves, for metamaterials where the magnetic interaction between its units is dominant. They exhibit phonon-like dispersion curves and they can transfer energy Shamonina and Solymar (2004); Syms et al. (2010), and they have been experimentally investigated both in linear and nonlinear SRR-based metamaterials Wiltshire et al. (2003); Dolling et al. (2006); Shadrivov et al. (2007). It is thus possible to fabricate contact-free data and power transfer devices which make use of the unique properties of the metamaterial structure, and may function as a frequency-selective communication channel for devices via their magneto-inductive wave modes Stevens et al. (2010).

Unfortunately, metamaterials structures comprising of resonant metallic elements revealed unfavorable characteristics that render them unsuitable for most practical applications. The SRRs, in particular, suffer from high Ohmic losses at frequencies close to their resonance, where metamaterials acquire their extraordinary properties. Moreover, those properties may only appear within a very narrow band, that is related to the weak coupling between elements. High losses thus hamper any substantial progress towards the practical use of these metamaterials in novel devices. Many applications are also hampered by the lack of tuning capabilities and relatively bulky size. However, another breakthrough came with the discovery of non-resonant, transmission line negative refractive index metamaterials Eleftheriades et al. (2002); Caloz and Itoh (2002), which very quickly led to several applications, at least in the microwaves Caloz (2009). Transmission line metamaterials rely on the appropriate combination of inductive-capacitive () lumped elements into large networks. The tremendous amount of activity in the field of metamaterials since has been summarized in various reviews Smith et al. (2004); Caloz and Itoh (2005); Linden et al. (2006); Padilla et al. (2006); Shalaev (2007); Litchinitser and Shalaev (2008); Anlage (2011); Soukoulis and Wegener (2011); Liu and Zhang (2011); Simovski et al. (2012) and books Eleftheriades and Balmain (2005); Caloz and Itoh (2006); Engheta and Ziokowski (2006); Pendry (2007); Marques et al. (2007); Krowne and Zhang (2007); Ramakrishna and Grzegorczyk (2009); Cui et al. (2010); Cai and Shalaev (2010); Solymar and Shamonina (2009); Noginov and Podolskiy (2012); Tong (2018).

### 1.2 Nonlinear, Superconducting, and Active Metamaterials

Dynamic tunability is a property that is required for applications Fan and Padilla (2015); in principle, one should be able to vary the effective (macroscopic) parameters of a metamaterial in real time, simply by varying an applied field. Tunability provides the means for fabricating meta-devices with switching capabilities Zheludev (2010, 2011), among others, and it can be achieved by the introduction of nonlinearity. Nonlinearity adds new degrees of freedom for metamaterial design that allows for both tunability and multistability - another desired property, that may offer altogether new functionalities and electromagnetic characteristics Lapine (2017), as well as wide-band permeability Lapine et al. (2012). It was very soon after the first demonstration of metamaterials, named at that time as negative refractive index materials, when it became clear that the SRR structure has considerable potential to enhance nonlinear effects due to the intense electric fields which can be generated in their slits OâBrien et al. (2004). Following these ideas, several research groups have demonstrated nonlinear metamaterial units, by filling the SRR slits with appropriate materials, e.g., with a strongly nonlinear dielectric Hand and Cummer (2008), or with a photo-sensitive semiconductor. Other approaches have made use of semiconducting materials, e.g., as substrates, on which the actual metamaterial is fabricated, that enables modulation of THz transmission by Chen et al. (2006). However, the most convenient method for introducing nonlinearity in SRR-based metamaterials was proved to be the insertion of nonlinear electronic components into the SRR slits, e.g., a variable capacitance diode (varactor) Powell et al. (2007); Wang et al. (2008). The dynamic tunability of a two-dimensional metamaterial comprising varactor-loaded SRRs by the power of an applied field has been demonstrated experimentally Shadrivov et al. (2008). Both ways of introducing nonlinearity affect the capacitance of the SRRs which becomes field-dependent; in the equivalent electrical circuit picture, in which the SRRs can be regarded as lumped element electrical oscillators, the capacitance acquires a voltage dependence and in turn a field-dependent magnetic permeability. Nonlinear transmission line metamaterials are reviewed in Ref. Kozyrev and van der Weide (2008).

Nonlinearity does not however help in the reduction of losses; in nonlinear metamaterials the losses continue to be a serious problem. The quest for loss compensation in metamaterials is currently following two different pathways: a ”passive” one, where the metallic elements are replaced by superconducting ones Anlage (2011); Jung et al. (2014a), and an ”active” one, where appropriate constituents are added to metallic metamaterials that provide gain through external energy sources. In order to fabricate both nonlinear and active metamaterials, gain-providing electronic components such as tunnel (Esaki) diodes Esaki (1958) or particular combinations of other gain-providing devices have to be utilized. The Esaki diode, in particular, features a negative resistance part in its current-voltage characteristics, and therefore can provide both gain and nonlinearity in a conventional (i.e., metallic) metamaterial. Tunnel diodes which are biased so that they operate at the negative resistance region of their characteristics may also be employed for the construction of symmetric metamaterials, that rely on balanced gain and loss Lazarides and Tsironis (2013a). symmetric systems correspond to a new paradigm in the realm of artificial or ”synthetic” materials that do not obey separately the parity () and time () symmetries; instead, they do exhibit a combined symmetry Rüter et al. (2010); Regensburger et al. (2012). The notions of symmetric systems originate for non-Hermitian quantum mechanics Bender and Boettcher (1998); Bender et al. (2002), but they have been recently extended to optical lattices El-Ganainy et al. (2007); Makris et al. (2008). The use of active components which are incorporated in metamaterial unit elements has been actually proposed several years ago Boardman et al. (2007), and it is currently recognized as a very promising technique of compensating losses Boardman et al. (2010). Low-loss and active negative index metamaterials by incorporating gain material in areas with high local field have been demonstrated in the optical Xiao (2010). Recently, transmission lines with periodically loaded tunnel diodes which have the negative differential resistance property have been realized and tested as low-loss metamaterials, in which intrinsic losses are compensated by gain Jiang et al. (2011). Moreover, a combination of transistors and a split-ring has been shown to act as a loss-compensated metamaterial element Xu et al. (2012). In the latter experiment, the quality factor for the combined system exhibits huge enhancement compared with that measured for the split-ring alone.

The ”passive” approach to loss reduction employes superconducting materials, i.e, materials exhibiting absence of dc resistance below a particular temperature, known as the critical temperature, . A rough classification of the superconducting materials is made on the basis of their critical temperature; according to that, there are low- and high- superconducting materials. The former include primarily elemental and binary compounds, like Niobium (Nb), Niobium di-Selenide (NbSe) and more recently Niobium Nitride (NbN), while the most known representatives of the latter are the superconducting perovskites such as Yttrium-Barium-Copper-Oxide (YBCO). The latter is the most commonly used perovskite superconductor which typically has a critical temperature , well above the boiling point of liquid Nitrogen. The last few years, there has been an increasing interest in superconducting metamaterials that exploit the zero resistance property of superconductivity, targeting at severe reduction of losses and the emergence of intrinsic nonlinearities due to the extreme sensitivity of the superconducting state to external stimuli Zheludev (2011); Anlage (2011). The direct approach towards fabrication of superconducting metamaterials relies on the replacement of the metallic split-rings of the conventional SRR-based metamaterials by superconducting ones. More shopisticated realizations of superconducting metamaterials result from the replacement of the metallic SRRs by rf SQUIDs (Superconducting QUantum Interference Devices) Likharev. (1986); those SQUID metamaterials are discussed below.

Superconducting metamaterials are not however limited to the above mentioned realizations, but they also include other types of artificial metamaterials; thin superconducting plates have been used in a particular geometrical arrangement to ”beat the static” Narimanov (2008) and make possible a zero frequency metamaterial (dc metamaterial) Wood and Pendry (2007); Magnus et al. (2008); Navau et al. (2009); Gömöry et al. (2012); Mawatari et al. (2012). Other types of superconducting metamaterials in the form of heterostructures, where superconducting layers alternate with ferromagnetic or non-magnetic metallic layers have been shown to exhibit electromagnetically induced transparency Kurter et al. (2011a); Wu et al. (2011a); Jin et al. (2013), switching capabilities Kurter et al. (2012), magnetic cloaking, and concentration Prat-Camps et al. (2013). Recently, tunable electromagnetically induced transparency has been demonstrated in a Niobium Nitride (NbN) terahertz (THz) superconducting metamaterial. An intense THz pulse is used to induce nonlinearities in the NbN thin film and thereby tune the electromagnetically induced transparency-like behavior Zhang et al. (2017a). Furthermore, the dynamic process of parity-time () symmetry breaking was experimentally demonstrated in a hybridized metasurface which consists of two coupled resonators made from metal and NbN Wang et al. (2017). Negative refraction index metamaterials in the visible spectrum, based on MgB/SiC composites, have been also realized Limberopoulos et al. (2009), following prior theoretical investigations Kussow et al. (2007). Moreover, there is substantial evidence for negative refraction index in layered superconductors above the plasma frequency of the Josephson plasma waves Golick et al. (2010), that was theoretically investigated by several authors Pimenov et al. (2005); Rakhmanov et al. (2010). Other types of superconducting metamaterials include those made of magnetically active planar spirals Kurter et al. (2011b), as well as those with rather special (”woodcut”) geometries Savinov et al. (2012a), two-dimensional arrays of Josephson junctions Adams (2013), as well as superconducting ”left-handed” transmission lines Salehi et al. (2005); Wang and Lancaster (2006). Recently, in a novel one-dimensional Josephson metamaterial composed of a chain of asymmetric SQUIDs, strong Kerr nonlinearity was demonstrated Zhang et al. (2017b). Moreover, the Kerr constant was tunable over a wide range, from positive to negative values, by a magnetic flux threading the SQUIDs.

### 1.3 Superconducting Metamaterials from Zero to Terahertz Frequencies

There are several demonstrations of superconducting metamaterial elements which exhibit tunability of their properties by varying the temperature or the applied magnetic field Ricci et al. (2005, 2007); Gu et al. (2010); Fedotov et al. (2010); Chen et al. (2010); Jung et al. (2013); Trang et al. (2013). We should also mention the theoretical investigations (nonlinear circuit modeling) on a multi-resonant superconducting split-ring resonator Mazdouri and Javadzadeh (2017), and on a ”meta-atom” composed of a direct current (dc) SQUID and a superconducting rod attached to it, which exhibits both electric and magnetic resonant response Shramkova et al. (2017). Superconducting split-rings combined into two-dimensional planar arrays form superconducting metamaterials exhibiting tunability and switching capabilities at microwave and THz frequencies Gu et al. (2010); Ricci et al. (2007); Kurter et al. (2010); Jin et al. (2010); Chen et al. (2010); Zhang et al. (2011); Wu et al. (2011b); Zhang et al. (2012, 2013); C. H. Zhang et al. (2013); Grady et al. (2013). Up to the time of writing, metamaterials comprising superconducting SRRs employ one of the following geometries:

(i) square SRRs with rectangular cross-section in the double, narrow-side coupled SRR geometry Ricci et al. (2005); Ricci and Anlage (2006); Ricci et al. (2007);

(iii) square SRRs with rectangular cross-section in the single SRR geometry Gu et al. (2010);

(iv) electric inductive-capacitive SRRs of two different types Singh et al. (2012).

Also, novel metamaterial designs including a ”woodcut” type superconducting metamaterial, and niobium-connected asymmetrically split-ring metamaterials were demonstrated Savinov et al. (2013). All these metamaterials were fabricated in the planar geometry, using either conventional, low superconductors such as niobium (Nb) and niobium nitride films, or the most widely used member of the high superconductor family, i.e., the yttrium-barium-copper-oxide (). The experiments were performed in microwaves and in the (sub-)THz range ().

All these superconducting metamaterials share a common feature: they all comprise resonant sub-wavelength superconducting elements, that exhibit a strong response at one particular frequency, i.e., the resonance frequency, . That resonance frequency is tunable under external fields, such as temperature, constant (dc) and time-periodic (ac) magnetic fields, and applied current, due to the extreme sensitivity of the superconducting state to external stimuli. (Note however that for some geometries there can be more than one strong resonances.) The experimental investigation of the resonances and their ability for being shifted either to higher or lower frequencies relies on measurements of the complex transmission spectrums, with dips signifying the existence of resonances. However, not only the frequency of a resonance but also its quality is of great interest in prospective applications. That quality is indicated by the depth of the dip of the complex transmission magnitude in the corresponding transmission spectrum, as well as its width, and quantified by the corresponding quality factor . In general, the quality factor increases considerably as the temperature decreases below the critical one at . Other factors, related to the geometry and material issues of the superconducting SRRs that comprise the metamaterial, also affect the resonance frequency . Thus, the resonance properties of a metamaterial can be engineered to achieve the desired values, just like in conventional metamaterials. However, for superconducting metamaterials, the thickness of the superconducting film seems to be an important parameter, because of the peculiar magnetic properties of superconductors. Using proper design, it is possible to switch on and off the resonance in superconducting metamaterials in very short time-scales, providing thus the means of manufacturing devices with fast switching capabilities.

### 1.4 Summary of earlier work in superconducting metamaterials

In this Subsection, a brief account is given on the progress in the development and applications of superconducting (both classical and quantum) metamaterials, i.e., metamaterials utilizing either superconducting materials or devices, is given. A more detailed and extended account is given in two review articles on the subject Anlage (2011); Jung et al. (2014a), as well as in Chapter of a recently published book Tong (2018). The status of the current research on SQUID metamaterials is discussed separately in the next Subsection (Subsection ). In the older of the two review articles Anlage (2011), the properties of superconductors which are relevant to superconducting metamaterials, and the advantages of superconducting metamaterials over their normal metal counterparts are discussed. The author reviews the status of superconductor-ferromagnet composites, dc superconducting metamaterials, radio frequency (rf) superconducting metamaterials, and superconducting photonic crystals (although the latter fall outside the domain of what are usually called metamaterials). There is also a brief discussion on SQUID metamaterials, with reference to the theoretical works in which it was proposed to use an array of rf SQUIDs as a metamaterial Du et al. (2006); Lazarides and Tsironis (2007). In the second review article Jung et al. (2014a), a more detailed account on the advantages of superconducting metamaterials over their normal metal counterparts was given, along with an update on the status of superconducting metamaterials. Moreover, analogue electromagnetically-induced transparency superconducting metamaterials and superconducting SRR-based metamaterials are also reviewed. In this review article, there is also a discussion on the first experiments on SQUID metamaterials Butz et al. (2013a); Trepanier et al. (2013); Jung et al. (2014b) which have confirmed earlier theoretical predictions. However, a lot of experimental and theoretical work on SQUID metamaterials has been performed after the time of writing of the second review article in this field. The present review article aims to give an up-to-date and extended account of the theoretical and experimental work on SQUID metamaterials and reveal their extraordinary nonlinear dynamic properties. SQUID metamaterials provide a unique testbed for exploring complex spatio-temporal dynamics. In the quantum regime, a prototype model for a ”basic” SCQMM which has been investigated by several authors is reviewed, which exhibits novel physical properties. Some of these properties are discussed in detail.

Superconductivity is a macroscopic quantum state of matter which arises from the interaction between electrons and lattice vibrations; as a results, the electrons form pairs (Cooper pairs) which condense into a single macroscopic ground state. The latter is the equilibrium thermodynamic state below a transition (critical) temperature . The ground state is separated by a temperature-dependent energy gap from the excited states with quasi-particles (quasi-electrons). As mentioned earlier, the superconductors are roughly classified into high and low critical temperature ones (high and low, respectively). In some circumstances, the Cooper pairs can be described in terms of a single macroscopic quantum wavefunction , whose squared magnitude is interpreted as the local density of superconducting electrons (), and whose phase is coherent over macroscopic dimensions. Superconductivity exhibits several extraordinary properties, such as zero dc resistance and the Meissner effect. Importantly, it also exhibits macroscopic quantum phenomena including fluxoid quantization and the Josephson effects at tunnel (insulating) barriers and weak links. When two superconductors and are brought close together and separated by a thin insulating barrier, there can be tunneling of Cooper pairs from one superconductor to the other. This tunneling produces a supercurrent (Josephson current) between and , , where is the critical current of the Josephson junction and is the gauge-invariant Josephson phase, with and being the phases of the macroscopic quantum wavefunctions of and , respectively, the electromagnetic vector potential in the region between and , and the flux quantum ( is Planck’s constant and the electron’s charge). Depending on whether is time-dependent or not, the appearance of the supercurrent is referred to as the ac or the dc Josephson effect.

Superconductors bring three unique advantages to the development of metamaterials in the microwave and sub-THz frequencies which have been analyzed in Refs. Anlage (2011); Jung et al. (2014a). Namely, (i) low losses (one of the key limitations of conventional metamaterials), (ii) the possibility for higher compactification of superconducting metamaterials compared to other realizations (superconducting SRRs can be substantially miniaturized while still maintaining their low-loss properties), and (iii) strong nonlinearities inherent to the superconducting state, which allow for tunability and provide switching capabilities. The limitations of superconducting metamaterials arise from the need to create and maintain a cryogenic environment, and to bring signals to and from the surrounding room temperature environment. Quite fortunately, closed-cycle cryocoolers have become remarkably small, efficient and inexpensive since the discovery of high superconductors, so that they are now able to operate for several years unattended, and moreover they can accommodate the heat load associated with microwave input and output transmission lines to room temperature. Superconductors can also be very sensitive to variations in temperature, stray magnetic field, or strong rf power which can alter their properties and change the behavior of the metamaterial. Thus, careful temperature control and high quality magnetic shielding are often required for reliable performance of superconducting metamaterials.

Superconducting metamaterials exhibit intrinsic nonlinearity because they are typically made up of very compact elementary units, resulting in strong currents and fields within them. Nonlinearity provides tunability through the variation of external fields. For example, a connected array of asymmetrically-split Nb resonators shows transmission tunable by current at sub-THz frequencies due to localized heating and the entrance of magnetic vortices Savinov et al. (2012b). The change in superfluid density by a change in temperature was demonstrated for a superconducting thin film SRR Ricci and Anlage (2006). Later, it was demonstrated that the resonant frequency of a Nb SRR changes significantly with the entry of magnetic vortices Ricci et al. (2007). Similar results showing complex tuning with magnetic field were later demonstrated at microwaves using high superconducting SRRs Trang et al. (2013) and at sub-THz frequencies with similarly designed Nb SRRs Savinov et al. (2013). The nonlinearity associated with the resistive transition of the superconductor was exploited to demonstrate a bolometric detector at sub-THz frequencies using the collective properties of an asymmetrically split-ring array made of Nb Savinov et al. (2013). Thermal tuning has been accomplished at THz frequencies by varying the temperature in high (YBCO) metamaterial Chen et al. (2010) and NbN electric inductive-capacitive thin film resonators Wu et al. (2011b). Enchancement of thermal tunability was accomplished by decreasing the thickness of the high superconducting films which make up square SRRs Singh et al. (2013). Fast nonlinear response can be obtained in superconducting films in THz time-domain experiments. In such an experiment it was found that intense THz pulses on a NbN metamaterial could break significant number of Cooper pairs to produce a large quasi-particle density which increases the effective surface resistance of the film and modulates the depth of the SRR resonance C. H. Zhang et al. (2013); Zhang et al. (2013). The tuning of high (YBCO) SRR metamaterial with variable THz beam intensity has demonstrated that the resonance strength decreases and the resonance frequency shifts as the intensity is increased Grady et al. (2013).

A natural opportunity to create a negative real part of the effective magnetic permeability is offered by a gyromagnetic material for frequencies above the ferromagnetic resonance Chui and Hu (2002). However, the imaginary part of is quite large near the resonance and limits the utility of such a metamaterial. A hybrid metamaterial, resulting from the combination of the gyromagnetic material mentioned earlier with a superconductor can help to reduce losses Nurgaliev (2008). A superlattice consisting of high superconducting and manganite ferromagnetic layers (YBCO/LSMO) was created and it was shown to produce a negative index band in the vicinity of the ferromagnetic resonance () at magnetic fields between and Pimenov et al. (2005). More recently, a metamaterial composed of permanent magnetic ferrite rods and metallic wires was fabricated. This metamaterial exhibits not only negative refraction but also near-zero refraction, without external magnetic field Bi et al. (2014).

The concept of a dc metamaterial operating at very low frequencies that could make up a dc magnetic cloak has been proposed and investigated in Ref. Wood and Pendry (2007). The first realization of such a metamaterial, which is based on non-resonant elements, consists of an array of superconducting plates Magnus et al. (2008). The superconducting elements exclude a static magnetic field, and provide the foundation for the diamagnetic effect, when that field is applied normal to the plates. The strength of the response depends on the ratio between the dimension of the plates and the lattice spacing. An experimental demonstration of a dc metamaterial cloak was made using an arrangement of Pb thin film plates Magnus et al. (2008). Subsequent theoretical work has refined the dc magnetic cloak design and suggested that it can be implemented with high superconducting thin films Navau et al. (2009). It was later demonstrated experimentally that a specially designed cylindrical superconductor-ferromagnet bilayer can exactly cloak uniform static magnetic fields Gömöry et al. (2012).

Superconducting rf metamaterials consisting of two-dimensional Nb spirals developed on quartz substrates show strong tunability as the transition temperature is approached Kurter et al. (2010, 2011b). Rf metamaterials have great potential in applications such as magnetic resonance imaging devices for non-invasive and high resolution medical imaging Wiltshire et al. (2001). The superconducting rf metamaterials have many advantages over their normal metal counterparts, such as reducing considerably the Ohmic losses, compact structure, and sensitive tuning of resonances with temperature or rf magnetic field, which makes them promising for rf applications. Similar, high rf metamaterials (in which the spirals are made by YBCO) were also fabricated, which enable higher operating temperatures and greater tunability Ghamsari et al. (2013).

The elementary units (i.e., the meta-atoms) in a metamaterial can be combined into meta-molecules so that the interactions between the meta-atoms can give rise to qualitatively new effects, such as the classical analogue of the electromagnetically induced transparency (EIT). This effect has been observed in asymmetrically-split ring metamaterials in which Fano resonances have been measured as peaks in the transmission spectrum, corresponding to metamaterial induced transparency Fedotov et al. (2010). Metamaterials consisting of normal metal - superconductor hybrid meta-molecules can create strong classical EIT effects Kurter et al. (2011a). The meta-molecule consists of a gold (Au) strip with end caps and two superconducting (Nb) SRRs (the ”bright” and the ”dark” element, respectively). A tunable transparency window which could even be switched off completely by increasing the intensity of the signal propagating through the meta-molecule was demonstrated Kurter et al. (2011a, 2012). EIT effects were also observed in the THz domain utilizing NbN bright and dark resonators to create a transparency window Wu et al. (2011a). Further experiments on all-superconducting (NbN) metamaterials utilizing strongly coupled SRR-superconducting ring elements showed enhanced slow-light features Jin et al. (2013).

A quantum metamaterial is meant to be an artificial optical medium that (a) comprise quantum coherent unit elements whose parameters can be tailored, (b) the quantum states of (at least some of) these elements can be directly controlled, and (c) maintain global coherence for sufficiently long time. These properties make a quantum metamaterial a qualitatively different system Zagoskin (2011); Zagoskin et al. (2016). Superconducting quantum metamaterials offer nowadays a wide range of prospects from detecting single microwave photons to quantum birefringence and superradiant phase transitions Jung et al. (2014a). They may also play a role in quantum computing and quantum memories. The last few years, novel superconducting devices, which can be coupled strongly to external electromagnetic field, can serve as the quantum coherent unit elements of superconducting quantum metamaterials (SCQMMs). For example, at ultra-low temperatures, superconducting loops containing Josephson junctions exhibit a discrete energy level spectrum and thus behave in many aspects as quantum meta-atoms. It is very common to approximate such devices as two-level quantum systems, referred to as superconducting qubits, whose energy level splitting corresponds to a frequency of the order of a few GHz. The interaction between light and a SCQMM is described by photons coupling to the artificial two-level systems, i.e., the superconducting qubits. The condition of keeping the energy of thermal fluctuations , where is Boltzmann’s constant and the temperature, below the energy level splitting of the qubit, where is Planck’s constant and the transition frequency, requires temperatures well below . In the past few years, research on superconducting qubits has made enormous progress that paves the way towards superconducting qubit-based quantum metamaterials.

There are several theoretical investigations on the physics of one-dimensional arrays of superconducting qubits coupled to transmission-line resonators Rakhmanov et al. (2008); Zagoskin et al. (2009); Ian et al. (2012); Viehmann O and F (2013); Shvetsov et al. (2013); Volkov and Fistul (2014); Asai et al. (2015); Ivić et al. (2016); Asai et al. (2018). Moreover, two-dimensional Zueco et al. (2012) and three-dimensional Zagoskin (2012) SCQMMs based on Josephson junction networks were proposed. A more extended discussion of the theoretical works on SCQMMs is given in Subsection . Still there is little progress in the experimental realization of such systems. The first SCQMM which was implemented in 2014 Macha et al. (2014), and comprises flux qubits arranged in a double chain geometry. In that prototype system, the dispersive shift of the resonator frequency imposed by the SCQMM was observed. Moreover, the collective resonant coupling of groups of qubits with the quantized mode of a photon field was identified, despite of the relatively large spread of the qubit parameters. Recently, an experiment on an SCQMM comprising an array of twin flux qubits, was demonstrated Shulga et al. (2018). The qubit array is embedded directly into the central electrode of an Al coplanar waveguide; each qubit contains Josephson junctions, and it is strongly coupled to the electromagnetic waves propagating through the system. It was observed that in a broad frequency range, the transmission coefficient through that SCQMM depends periodically on the external magnetic field. Moreover, the excitation of the qubits in the array leads to a large resonant enhancement of the transmission. We undoubtedly expect to see more experiments with arrays of superconducting qubits placed in transmission lines or waveguides in the near future.

### 1.5 SQUID Metamaterials

The rf SQUIDs, mentioned above, are highly nonlinear superconducting devices which are long known in the Josephson community and encompass the Josephson effect Josephson (1962). The simplest version of a SQUID is made by a superconducting ring which is interrupted by a Josephson junction (JJ); the latter is typically formed by two superconductors separated by a thin insulating (dielectric) layer. The current through the insulating layer and the voltage across the junction are then determined by the celebrated Josephson relations and crucially affect the electromagnetic behavior of the rf SQUID. SQUIDs have found numerous technological applications in modern science Kleiner et al. (2004); Clarke and Braginski (2004a, b); Fagaly (2006); they are most commonly used as magnetic field sensors, since they can detect even tiny magnetic fields and measure their intensity with unprecedented precision. SQUID metamaterials constitute the direct superconducting analogue of conventional (metallic) nonlinear (i.e., varactor loaded) SRR-based magnetic metamaterials, which result from the replacement of the nonlinear SRRs by rf SQUIDs. The latter possess inherent nonlinearity due to the Josephson element. Similarly to the conventional (metallic), SRR-based magnetic metamaterials, the SQUIDs are coupled magnetically to each other through magnetic dipole-dipole interactions. Several years ago, theoretical investigations have suggested that rf SQUID arrays in one and two dimensions can operate as magnetic metamaterials both in the classical Lazarides and Tsironis (2007) and in the quantum regime Du et al. (2006), and they may exhibit negative and/or oscillating effective magnetic permeability in a particular frequency band which encloses the resonance frequency of individual SQUIDs. Recent experiments on single rf SQUIDs in a waveguide demonstrated directly the feasibility of constructing SQUID-based thin-film metasurfaces Jung et al. (2013). Subsequent experiments on one-dimensional, quasi-two-dimensional, and truly two-dimensional SQUID metamaterials have revealed a number of several extraordinary properties such as negative diamagnetic permeability Jung et al. (2013); Butz et al. (2013a), broad-band tunability Butz et al. (2013a); Trepanier et al. (2013), self-induced broad-band transparency Zhang et al. (2015), dynamic multistability and switching Jung et al. (2014b), as well as coherent oscillations Trepanier et al. (2017). Moreover, nonlinear localization Lazarides et al. (2008a) and nonlinear band-opening (nonlinear transmission) Tsironis et al. (2014), as well as the emergence of dynamic states referred to as chimera states in current literature Lazarides et al. (2015); Hizanidis et al. (2016a), have been demonstrated numerically in SQUID metamaterial models. Those counter-intuitive dynamic states, which have been discovered numerically in rings of identical phase oscillators Kuramoto and Battogtokh (2002), are reviewed in Refs. Panaggio and Abrams (2015); Yao and Zheng (2016). Moreover, numerical investigations on SQUID metamaterials on Lieb lattices which posses a flat band in their frequency spectrum, reveal the existence of flat-band localized states in the linear regime and the more well-known nonlinearly localized states in the nonlinear regime Lazarides and Tsironis (2017). The interaction of an electromagnetic wave with a diluted concentration of a chain of SQUIDs in a thin film suggests a mechanism for the excitation of magnetization waves along the chain by a normally incident field Maimistov and Gabitov (2010). In the linear limit, a two-dimensional array of rf SQUIDs acts as a metasurface that controls the polarization of an electromagnetic wave Caputo et al. (2015).

SQUID arrays have been also integrated in larger devices in order to take advantage of their extraordinary properties; notably, amplification and squeezing of quantum noise has been recently achieved with a tunable SQUID-based metamaterial Castellanos-Beltran et al. (2008). Other important developments demonstrate clearly that SQUID-based metamaterials enable feedback control of superconducting cubits Risté et al. (2012), observation of Casimir effects Lähteenmäki et al. (2013), measurements of nanomechanical motion below the standard quantum limit Teufel et al. (2009), and three-wave mixing Roch et al. (2012). At sufficiently low (sub-Kelvin) temperatures, SQUID metamaterials provide access to the quantum regime, where rf SQUIDs can be manipulated as flux and phase qubits Poletto et al. (2009); Castellano et al. (2010). Recently, the technological advances that led to nano-SQUIDs make possible the fabrication of SQUID metamaterials at the nanoscale Wernsdorfer (2009).

From the above discussion it should be clear that the field of superconducting metamaterials, in which superconductivity plays a substantial role in determining their properties, has expanded substantially. In this review, we focus on the SQUID metamaterials, that represent an area of the field of superconducting metamaterials, which however has already reached a level of maturity. We also focus on SCQMMs, and in particular on a prototype model for a chain of charge qubits in a transmission-line resonator Rakhmanov et al. (2008). The SCQMMs are related to the (classical, i.e., not truly quantum) SQUID metamaterials in that they also encompass the Josephson effect. In Section 2, we describe the SQUID metamaterial models used for simulating real systems in current research, we provide the corresponding dispersion of flux waves which can propagate in SQUID metamaterials, and we present numerical results (along with selected experimental ones), which reveal novel properties such as wide-band tunability, energy transmission, and multistability. In Section 3, we present and discuss results on nonlinear localization in SQUID metamaterials, which leads to the generation of states referred to as discrete breathers. In that Section, we also emphasize the possibility for the emergence of chimera states in SQUID metamaterial models with either nonlocal or local (nearest-neighbor) couplng between their elements (i.e., the SQUIDs). In Section 4, the dynamical model for SQUID metamaterials on a Lieb lattice is presented, along with its full frequency spectrum. The latter contains a flat band, which allows for the formation of flat-band localized states in the linear regime. The case of nonlinearly localized states, which can be formed in the nonlinear regime, as well as the transition between the two regimes, is investigated. In Section 5, we describe a model SCQMMs (a chain of charge qubits in a superconducting transmission-line resonator) and discuss the possibility for having propagating self-induced transparent or superradiant pulses in that medium. Most importantly, those pulses induce quantum coherence effects in the medium itself, by exciting population inversion pulses in the qubit subsystem. Moreover, the speed of the propagating pulses can be controlled by proper engineering of the parameters of the qubits. The most important points made in this review are summarized in Section 6.

## 2 SQUID-Based Metamaterials I: Models and Collective Properties

### 2.1 The rf-SQUID as an artificial magnetic ”atom”

The Superconducting QUantum Interference Device (SQUID) is currently one of the most important solid-state circuit elements for superconducting electronics Anders et al. (2010); among many other technological applications Kleiner et al. (2004); Clarke and Braginski (2004b); Fagaly (2006), SQUIDs are used in devices that provide the most sensitive sensors of magnetic fields. Recent advances that led to nano-SQUIDs Wernsdorfer (2009) makes the fabrication of SQUID metamaterials at the nanoscale an interesting possibility. The radio-frequency (rf) SQUID, in particular, shown schematically in Fig. 1(a), consists of a superconducting ring of self-inductance interrupted by a Josephson junction (JJ) Josephson (1962). A JJ is made by two superconductors connected through a ”weak link”, i.e., through a region of weakened superconductivity. A common type of a JJ is usually fabricated by two superconductors separated by a thin dielectric oxide layer (insulating barrier) as shown in Fig. 2(a); such a JJ is referred to as a superconductor-insulator-superconductor (SIS) junction. The fundamental properties of JJs have been established long ago Barone and Patternó (1982); Likharev. (1986), and their usage in applications involving superconducting circuits has been thoroughly explored. The observed Josephson effect in such a junction, has been predicted by Brian D. Josephson in 1962 and it is of great importance in the field of superconductivity as well as in physics. That effect has been exploited in numerous applications in superconducting electronics, sensors, and high frequency devices. In an ideal JJ, whose electrical circuit symbol is shown in Fig. 2(b), the (super)current (Josephson current) through the JJ and the voltage across the JJ are related through the celebrated Josephson relations Josephson (1962)

(1) |

where is the critical current of the JJ, is the flux quantum, with with and being the Planck’s constant and the electron’s charge, respectively, and is the difference of the phases of the order parameters of the superconductors at left () and right () of the barrier and , respectively (Fig. 2(a)), i.e., (the Josephson phase). In the presence of an electromagnetic potential , the corresponding gauge-invariant Josephson phase is

(2) |

In an ideal JJ, in which the (super)current is carried solely by Cooper pairs, there is no voltage drop across the barrier of the JJ. In practice, however, this can only be true at zero temperature (), while at finite temperatures there will always be a quasi-partice current. The latter is carried single-electron excitations (quasi-electrons) resulting from thermal breaking of Cooper pairs due to the non-zero temperature, and it is subjected to losses. In superconducting circuits, a real JJ is often considered as the parallel combination of an ideal JJ (described by Eqs. (1)), a shunting capacitance , due to the thin insulating layer, and a shunting resistance , due to quasi-electron tunneling through the insulating barrier, i.e., the quasi-particle current (Fig. 2(c)). This equivalent electrical circuit for a real JJ is known as the Resistively and Capacitively Shunted Junction (RCSJ) model, and is widely used for modeling SIS Josephson junctions. The ideal JJ can be also described as a variable inductor, in a superconducting circuit. From the Josephson relations Eqs. (1) and the current-voltage relation for an ordinary inductor , it is easily deduced that the Josephson inductance is

(3) |

Note that Eq. (3) describes a nonlinear inductance, since depends both on the current and the voltage through the Josephson phase .

Due to the Josephson element (i.e., the JJ), the rf SQUID is a highly nonlinear oscillator that responds in a manner analogous to a magnetic ”atom”, exhibiting strong resonance in a time-varying magnetic field with appropriate polarization. Moreover, it exhibits very rich dynamic behavior, including chaotic effects Fesser et al. (1983); Ritala and Salomaa (1984); Sørensen et al. (1985) and tunability of its resonance frequency with external fields Lazarides and Tsironis (2012). The equivalent electrical circuit for an rf SQUID in a time-dependent magnetic field threading perpendicularly its loop, as shown schematically in Fig. 1(b), comprises a flux source in series with an inductance and a real JJ described by the RCSJ model.

The dynamic equation for the flux threading the loop of the rf SQUID is obtained by direct application of Kirchhoff’s laws, as

(4) |

where is the magnetic flux quantum, is the critical current of the JJ, and is the temporal variable. Eq. (4) is derived from the combination of the single-SQUID flux-balance relation

(5) |

and the expression for the current in the SQUID provided by the RCSJ model

(6) |

Eq. (4) has been studied extensively for more than three decades, usually under an external flux field of the form

(7) |

i.e., in the presence of a time-independent (constant, dc) and/or a time-dependent (usually sinusoidal) magnetic field of amplitude and frequency . The orientation of both fields is such that their flux threads the SQUID loop. In the absence of dc flux, and very low amplitude of the ac field (, linear regime), the SQUID exhibits resonant magnetic response at

(8) |

where

(9) |

is the inductive-capacitive () SQUID frequency and SQUID parameter, respectively. Eq. (4) is formally equivalent to that of a massive particle in a tilted washboard potential

(10) |

with being the Josephson energy. That potential has a number of minimums which depends on the value of , while the location of those minimums varies with the applied dc (bias) flux . For (non-hysteretic regime) the potential is a corrugated parabola with a single minimum which moves to the right with increasing , as it is shown in Fig. 3(a). For (hysteretic regime) there are more than one minimums, while their number increases with further increasing . A dc flux can move all these minimums as well (Fig. 3(b)) The emergence of more and more minimums with increasing at is illustrated in Fig. 3(c). It should be stressed here that in general the external flux is a sum of a dc bias and an ac (time-periodic) term. In that case, the potential as a whole rocks back and forth at the frequency of the external driving ac flux, . Then, in order to determine the possible stationary states of the SQUID, the principles of minimum energy conditions do not apply; instead, the complete nonlinear dynamic problem has to be considered.

Normalization.- For an appropriate normalization of the single SQUID equation (4) and the corresponding dynamic equations for the one- and two-dimensional SQUID metamaterials discussed below, the following relations are used

(11) |

i.e., frequency and time are normalized to and its inverse , respectively, while all the fluxes and currents are normalized to and , respectively. Then, Eq. (4) is written in normalized form as

(12) |

where the overdots denote derivation with respect to the normalized temporal variable , is the normalized external flux, and

(13) |

is the rescaled SQUID parameter and dimensionless loss coefficient, respectively. The term which is proportional to in Eq. (12) actually represents all of the dissipation coupled to the rf SQUID.

The properties of the many variants of the SQUID have been investigated for many years, and they can be found in a number of review articles Gallop and Petley (1976); Clarke (1986); Jenks et al. (1997); Koelle et al. (1999); Kleiner et al. (2004); Fagaly (2006); Clarke (2010), textbooks Gallop (1991), and a Handbook Clarke and Braginski (2004a, b). Here we focus on the multistability and the tunability properties of rf SQUIDs, which are important for our later discussions on SQUID metamaterials. As it was mentioned earlier, the rf SQUID is a nonlinear oscillator which exhibits strong resonant response at a particular frequency to a sinusoidal (ac) flux field. For low amplitudes of the ac flux field, i.e., in the linear regime, the single-SQUID resonance frequency is given in Eq. (8); in units of the inductive-capacitive () SQUID frequency, the single-SQUID resonance frequency is

(14) |

It was theoretically demonstrated that the SQUID resonance can be tuned within a broad band of frequencies either by a dc flux bias or by the amplitude of an ac flux field . Soon after these predictions, these tunability properties were confirmed experimentally Butz et al. (2013a); Trepanier et al. (2013). Further experiments showed that the single-SQUID resonance frequency is also tunable by the power of the ac flux field, as well as with temperature Zhang et al. (2015). In particular, with increasing the amplitude of the ac flux field from low to high values, the single-SQUID resonance frequency shifts from to (i.e., towards lower frquencies). Note that for , a typical value for and very close to those obtained in the experiments, the single-SQUID resonance frequency may vary from to , that is more than of variation (broad-band tunability). Moreover, the single-SQUID resonance curve, i.e., the oscillation amplitude of the flux through the SQUID loop as a function of the driving frequency , changes dramatically its shape. Such a resonance curve for relatively high ac flux amplitude is shown in Fig. 4(a) Hizanidis et al. (2016a). Resonance curves like that are calculated from Eq. (12). As it can be observed, the curve ”snakes” back and forth within a narrow frequency band around the geometrical resonance frequency . The solid-blue branches of the resonance curve indicate stable solutions, while the dashed-black branches indicate unstable ones (hereafter referred to as stable and unstable branches, respectively). The stable and unstable branches merge at particular points (turning points) denoted by SN; at all these points, in which , saddle-node bifurcations of limit cycles occur. Clearly, by looking at Fig. 4(a) one can identify that there are certain frequency bands in which more than one simultaneously stable solutions exist. For the sake of illustration, a vertical line has been drawned at . For that frequency, there are five (5) simultaneously stable solutions which are marked by the letters , and four (4) unstable ones. This can be seen more clearly in the inset of Fig. 4(a) which shows an enlargement of the main figure around . Actually, at this frequency the number of possible solutions for the chosen set of simulation parameters is maximum. Thus for that set of parameters, is the maximum multistability frequency. The corresponding trajectories of the stable solutions at are shown in the phase portraits in Fig. 4(b). Multistability is enchanced (i.e., more solution branches appear) with increasing the ac flux amplitude or lowering the loss coefficient .

An approximation to the resonance curve for is given by Hizanidis et al. (2016a)

(15) |

where , , , and , which implicitly provides as a function of . The approximate flux amplitude - driving frequency curves from Eq. (15) are shown in Fig. 4(a) in turquoise color; clearly, they show excellent agreement with the numerical snaking resonance curve for . From Eq. (15), the frequency of the first saddle-node bifurcation (the one with the lowest ) can be calculated accurately. For simplicity, set and into Eq. (15) and then use the condition to obtain , where is the flux amplitude at which the first saddle-node bifurcation occurs. By substitution of into the simplified Eq. (15), we get , where is the frequency at which the first saddle-node bifurcation occurs. For the parameters used in Fig. 4(a), i.e., and , we get and which agree very well with the numerics.

The dynamic complexity for frequencies around the single-SQUID resonance increases immensely in a SQUID array with a relatively large number of SQUIDS. Although the coupling between SQUIDs is discussed in detail in the next Subsection (Subsection ), we believe it is appropriate to show here the solutions around the single-SQUID multistability frequency for a system of two coupled SQUIDs (Fig. 4(c)). The two SQUIDs, and , are identical, and they are coupled magnetically with strength through their mutual inductances. In Fig. 4(c), the flux through the loop of SQUID , , is plotted as a function of frequency ; since the presence of chaotic solutions was expected, the value of was plotted at the end of fifty (50) consecutive driving periods for each (after the transients have died out). The frequency interval of Fig. 4(c) is the same as that in the inset of Fig. 4(a). Different colors have been used to help distinguishing between different solution branches; also, unstable solutions have been omitted for clarity. It can be observed that the number of stable states for the two-SQUID system is more than two times larger than the stable states of the single-SQUID. Moreover, apart from the periodic solutions, a number of (coexisting) chaotic solutions has emerged. The dynamic complexity increases with increasing the number of SQUIDs which are coupled together. This effect, which has been described in the past for certain arrays of nonlinear oscillators is named as attractor crowding Wiesenfeld and Hadley (1989); Tsang and Wiesenfeld (1990). It has been argued that the number of stable limit cycles (i.e., periodic solutions) in such systems scale with the number of oscillators as . As a result, their basins of attraction crowd more and more tightly in phase space with increasing . The importance of this effect for the emergence of counterintuitive collective states in SQUID metamaterials is discussed in Subsection .

In Fig. 5, a number of flux amplitude - frequency and current amplitude-frequency curves are presented to demonstrate the tunability of the resonance frequency by varying the amplitude of the ac field or by varying the dc flux bias . Since the properties of a SQUID-based metamaterial are primarily determined by the corresponding properties of its elements (i.e., the individual SQUIDs), the tunability of a single SQUID implies the tunability of the metamaterial itself. In Figs. 5(a) and (b) the flux amplitude - frequency curves are shown for a SQUID in the non-hysteretic and the hysteretic regime with and , respectively. The ac field amplitude increases from top to bottom panel. In the top panels, the SQUID is close to the linear regime and the resonance curves are almost symmetric without exhibiting visible hysteresis. A sharp resonance appears at as predicted from linear theory. With increasing the nonlinearity becomes more and more appreciable and the resonance moves towards lower frequencies (middle and lower panels). Hysteretic effects are clearly visible in this regime. The resonance frequency of an rf SQUID can be also tuned by the application of a dc flux bias, , as shown in Figs. 5(c) and (d). While for the resonance appears close to (although slightly shifted to lower frequencies due to small nonlinear effects), it moves towards lower frequencies for increasing . Importantly, the variation of the resonance frequency does not seem to occur continuously but, at least for low , small jumps are clearly observable due to the inherently quantum nature of the SQUID which is incorporated to some extent into the phenomenological flux dynamics Eq. (12). The shift of the resonance with a dc flux bias in a single SQUID has been observed in high critical temperature (high) rf SQUIDs Zeng et al. (2000), as well as in low critical temperature (low) rf SQUIDs Trepanier et al. (2013); Zhang et al. (2015). In Fig. 5, there was no attempt to trace all possible branches of the resonance curves for clarity, and also for keeping the current in the SQUID to values less than the critical one for the JJ, i.e., for .

Another illustration of the multistability in an rf SQUID is shown in Fig. 6, which also reveals an anti-resonance effect. In Figs. 6(a) and (b), the current amplitude - frequency curves are shown in two cases; one close to the weakly nonlinear regime and the other in the strongly nonlinear regime, respectively. In Fig. 6(a), the curve does not exhibit hysteresis but it is slightly skewed; the resonance frequency is , slightly lower than the SQUID resonance frequency in the linear regime, (for ). In Fig. 6(b), the ac field amplitude has been increased by an order of magnitude with respect to that in Fig. 6(a), and thus strongly nonlinear effects become readily apparent. Five (5) stable branches can be identified in a narrow frequency region around , i.e., around the geometrical (inductive-capacitive, ) resonance frequency (unstable brances are not shown). The upper branches, which are extremely sensitive to perturbations, correspond to high values of the current amplitude , which leads the JJ of the SQUID to its normal state. The red arrows in Figs. 6, point at the location of an anti-resonance Hizanidis et al. (2016b) in the current amplitude - frequency curves. Such an anti-resonance makes itself apparent as a well-defined dip in those curves, with a minimum that almost reaches zero. The effect of anti-resonance has been observed in nonlinearly coupled oscillators subjected to a periodic driving force Woafo et al. (1998) as well as in parametrically driven nonlinear oscillators Chakraborty and Sarkar (2013). However, it has never before been observed in a single, periodically driven nonlinear oscillator such as the rf SQUID. In Figs. 6(c) and (d), enlargements of Figs. 6(a) and (b), respectively, are shown around the anti-resonance frequency. Although the ”resonance” region in the strongly nonlinear case has been shifted significantly to the left as compared with the weakly nonlinear case, the location of the anti-resonance has remained unchanged (eventhough in Figs. 6(a) and (b) differ by an order of magnitude). The knowledge of the location of anti-resonance(s) as well as the resonance(s) of an oscillator or a system of oscillators, beyond their theoretical interest, it is of importance in device applications. Certainly these resonances and anti-resonances have significant implications for the SQUID metamaterials whose properties are determined by those of their elements (i.e., the individual SQUIDs). When the SQUIDs are in an anti-resonant state, in which the induced current is zero, they do not absorb energy from the applied field which can thus transpass the SQUID metamaterial almost unaffected. Thus, in such a state, the SQUID metamaterial appears to be transparent to the applied magnetic flux as has been already observed in experiments on two-dimensional SQUID metamaterials Zhang et al. (2015); the observed effect has been named as broadband self-induced transparency. Moreover, since the anti-resonance frequency is not affected by , the transparency can be observed even in the strongly nonlinear regime, for which the anti-resonance frequency lies into the multistability region. In that case, the transparency of the metamaterial may be turned on and off as it has been already discussed in Ref. Zhang et al. (2015). Thus, the concept of the anti-resonance serves for making a connection between an important SQUID metamaterial property and a fundamental dynamical property of nonlinear oscillators.

The tunability of the SQUID resonance with a dc magnetic field and the temperature has been investigated in recent experiments Jung et al. (2013); Trepanier et al. (2013). Those investigations rely on the measurement of the magnitude of the complex transmission as one or more external parameters such as the driving frequency, the dc flux bias, and the temperature vary. Very low values of indicate that the SQUID is at resonance. In Fig. 7, the resonant response is identified by the red features. In the left panel, it is observed that the resonance vary periodically with the applied dc flux, with period . In the middle panel of Fig. 7, the effect of the temperature is revealed; as expected, the tunability bandwidth of the resonance decreases with increasing temperature. In the right panel of Fig. 7, the variation of the resonance frequency with the applied rf power is shown. Clearly, three different regimes are observed; for substantial intervals of low and high rf power, the resonance frequency is approximatelly constant at and , respectively, while for intermediate rf powers the resonance apparently dissapears. The latter effect is related to the broad-band self-induced transparency Zhang et al. (2015).

### 2.2 SQUID Metamaterials Models and Flux Wave Dispersion

Conventional (metallic) metamaterials comprise regular arrays of split-ring resonators (SRRs), which are highly conducting metallic rings with a slit. These structures can become nonlinear with the insertion of electronic devices (e.g., varactors) into their slits Lapine et al. (2003); Shadrivov et al. (2006a); Wang et al. (2008); Boardman et al. (2010); Lapine et al. (2014). SQUID metamaterials is the superconducting analogue of those nonlinear conventional metamaterials that result from the replacement of the varactor-loaded metallic rings by rf SQUIDs as it has been suggested both in the quantum Du et al. (2006) and the classical Lazarides and Tsironis (2007) regime. Recently, one- and two-dimensional SQUID metamaterials have been constructed from low critical temperature superconductors which operate close to liquid Helium temperatures Butz et al. (2013a, b); Trepanier et al. (2013); Zhang et al. (2015, 2016); Trepanier et al. (2017). The experimental investigation of these structures has revealed several novel properties such as negative diamagnetic permeability Jung et al. (2013); Butz et al. (2013a), broad-band tunability Butz et al. (2013a); Trepanier et al. (2013), self-induced broad-band transparency Zhang et al. (2015), dynamic multistability and switching Jung et al. (2014b), as well as coherent oscillations Trepanier et al. (2017), among others. Some of these properties, i.e., the dynamic multistability effect and tunability of the resonance frequency of SQUID metamaterials, have been also revealed in numerical simulations Lazarides and Tsironis (2013b); Tsironis et al. (2014). Moreover, nonlinear localization Lazarides et al. (2008a) and the emergence of counter-intuitive dynamic states referred to as chimera states in current literature Lazarides et al. (2015); Hizanidis et al. (2016a, b) have been demonstrated numerically in SQUID metamaterial models. The chimera states have been discovered numerically in rings of identical phase oscillators Kuramoto and Battogtokh (2002) (see Ref. Panaggio and Abrams (2015) for a review).

The applied time-dependent magnetic fields induce (super)currents in the SQUID rings through Faraday’s induction law, which couple the SQUIDs together through dipole-dipole magnetic forces; although weak due to its magnetic nature, that interaction couples the SQUIDs non-locally since it falls-off as the inverse cube of their center-to-center distance. Consider a one-dimensional linear array of identical SQUIDs coupled together magnetically through dipole-dipole forces. The magnetic flux threading the th SQUID loop is Lazarides et al. (2015)

(16) |

where the indices and run from to , is the external flux in each SQUID, is the dimensionless coupling coefficient between the SQUIDs at positions and , with being their corresponding mutual inductance, and

(17) |

is the current in each SQUID given by the RCSJ model Likharev. (1986), with and being the flux quantum and the critical current of the JJs, respectively. Recall that within the RCSJ framework, , , and are the resistance, capacitance, and self-inductance of the SQUIDs’ equivalent circuit, respectively. Combination of Eqs.(16) and (17) gives

(18) |

where is the inverse of the coupling matrix

(19) |

with being the coupling coefficient betwen nearest-neighboring SQUIDs.

The dimensionless coupling strength between the SQUIDs at site and (normalized center-to-center distance ), can be calculated either analytically or numerically. The self-inductance of the SQUIDs, for example, can be either estimated by an empirically derived equation Brojeny and Clem (2003) or it can be calculated using commercially available software (FastHenry). The mutual inductance can also be obtained numerically from a FastHenry calculation or it can be approximated using basic expressions from electromagnetism. The magnetic field generated by a wire loop, which is the approximate geometry of a SQUID, at a distance greater than its dimensions, is given by the Biot-Savart law as , where is the current in the wire, is the radius of the loop, and is the distance from the center of the loop. The mutual inductance between two such (identical) loops lying on the same plane is given by

(20) |

where it is assumed that the field is constant over the area of each loop, . For square loops of side , the radius should be replaced by . Eq. (20) explains qualitatively the inverse cube distance-dependence of the coupling strength between SQUIDs.

In normalized form Eq. (18) reads ()

(21) |

where the relations given in Eq. (11) have been used. Specifically, frequency and time are normalized to and its inverse , respectively, the fluxes and currents are normalized to and , respectively, the overdots denote derivation with respect to the normalized temporal variable, , , with being the normalized driving frequency, and , are given in Eq. (13). The (magnetoinductive) coupling strength between SQUIDs, which can be estimated from the experimental parameters in Ref. Kirtley et al. (2005), as well as from recent experiments Jung et al. (2013); Butz et al. (2013a), is rather weak due to its magnetic nature (of the order of in normalized units). Since that strength falls-off approximatelly as the inverse-cube of the distance between SQUIDs, a model which takes into account nearest-neighbor coupling only is sufficient for making reliable predictions. In that case, the coupling matrix assumes the simpler, tridiagonal and symmetric form

(22) |

For , the inverse of the coupling matrix is approximated to order by

(23) |

Substituting Eq. (23) into Eq. (21), the corresponding dynamic equations for the fluxes through the loops of the SQUIDs of a locally coupled SQUID metamaterial are obtained as

(24) |

where is the ”effective” external flux, with being the normalized external flux. The effective flux arises due to the nearest-neighbor approximation. For a finite SQUID metamaterial (with SQUIDs), is slightly different for the SQUIDs at the end-points of the array; specifically, for those SQUIDs since they interact with one nearest-neighbor only.

Linearization of Eq. (21) around zero flux with and gives for the infinite system

(25) |

By substitution of the plane-wave trial solution into Eq. (25), with being the wavevector normalized to ( is the side of the unit cell, see Fig. 1) , and using

(26) |

where is the ”distance” from the main diagonal of , we get

(27) |

where . Note that for the infinite system the diagonal elements of the inverse of the coupling matrix have practically the same value which is slightly larger than unity. The frequency is very close to the resonance frequency of individual SQUIDs, . Eq. (27) is the nonlocal frequency dispersion. By substitution of the same trial solution into Eq. (24), we get the nearest-neighbor frequency dispersion for flux waves, as

(28) |

Eqs. (27) and (28) result in slightly different frequency dispersion curves as can be observed in the lower panel of Fig. 8 for two different values of the coupling coefficient .

In order to increase the dimensionality and obtain the dynamic equations for the fluxes through the loops of the SQUIDs arranged in an orthogonal lattice, as shown in the left panel of Fig. 9 in which the unit cell is enclosed inside the green-dotted line, we first write the corresponding flux-balance relations Lazarides et al. (2008a); Tsironis et al. (2009); Lazarides and Tsironis (2013b),

(29) |

where is the flux threading the th SQUID of the metamaterial, is the total current induced in the th SQUID of the metamaterial, and are the magnetic coupling coefficients between neighboring SQUIDs, with and being the mutual inductances in the and directions, (). The subscripts and run from to and to , respectively. The current is given by the RCSJ model as

(30) |

Eq. (29) can be inverted to provide the currents as a function of the fluxes, and then it can be combined with Eq. (30) to give the dynamic equations Lazarides et al. (2008a)

(31) |

In the absence of losses (), the earlier equations can be obtained from the Hamiltonian function

(32) |

where

(33) |

is the canonical variable conjugate to , and represents the charge accumulating across the capacitance of the JJ of each rf SQUID. The above Hamiltonian function is the weak coupling version of that proposed in the context of quantum computation Roscilde et al. (2005). Eqs. (2.2) can be written in normalized form as

(34) |

where the overdots denote differentiation with respect to the normalized time , and

(35) |

The frequency dispersion of linear flux-waves in two-dimensional SQUID metamaterials can be obtained with the standard procedure, by using plane wave trial solutions into the linearized dynamic equations (34). That procudure results in the relation

(36) |

where is the normalized wavevector. The components of are related to those of the wavevector in physical units through with and being the wavevector component and center-to-center distance between neighboring SQUIDs in and direction, respectively. A density plot of the frequency dispersion equations (36) on the plane is shown in the right panel of Fig. 8 for a tetragonal (i.e., ) SQUID metamaterial. Assuming thus that the coupling is isotropic, i.e., , the maximum and minimum values of the linear frequency band are then obtained by substituting and , respectively, into Eq. (36). Thus we get

(37) |

that give an approximate bandwidth .

Note that the dissipation term appearing in Eq. (2.2) may result from the corresponding Hamilton’s equations with a time-dependent Hamiltonian Lazarides et al. (2008a)

(38) |

where is the Josephson energy, , and

(39) |

is the new canonical variable conjugate to which represents the generalized charge across the capacitance of the JJ of each rf SQUID. The Hamiltonian in Eq. (2.2) is a generalization in the two-dimensional lossy case of that employed in the context of quantum computation with rf SQUID qubits Roscilde et al. (2005); Corato et al. (2006).

### 2.3 Wide-Band SQUID Metamaterial Tunability with dc Flux

An rf SQUID metamaterial is shown to have qualitatively the same behavior as a single rf SQUID with regard to dc flux and temperature tuning. Thus, in close similarity with conventional, metallic metamaterials, rf SQUID metamaterials acquire their electromagnetic properties from the resonant characteristics of their constitutive elements, i.e., the individual rf SQUIDs. However, there are also properties of SQUID metamaterials that go beyong those of individual rf SQUIDs; these emerge through collective interaction of a large number of SQUIDs forming a metamaterial. Numerical simulations using the SQUID metamaterial models presented in the previous sub-section confirm the experimentally observed tunability patterns with applied dc flux in both one and two dimensions. Here, numerical results for the two-dimensional model are presented; however, the dc flux tunability patterns for either one- or two-dimensional SQUID metamaterials are very similar. Due to the weak coupling between SQUIDs, for which the coupling coefficient has been estimated to be of the order of in normalized units Zhang et al. (2015), the nearest-neighbor coupling between SQUIDs provides reliable results. The normalized equations for the two-dimensional model equations (34) are Lazarides et al. (2008a); Lazarides and Tsironis (2013b)

(40) |

where and . The total (symmetrized) energy of the SQUID metamaterial, in units of the Josephson energy , is then

(41) |

For , the average of that energy over one period of evolution

(42) |

where with being the normalized driving frequency, is constant when the obtained solution is locked to that driving frequency. In the following, the term ”tunability of the resonance” for the SQUID metamaterial is used; however, the notion of the resonance is rather appropriate for a single SQUID only. The SQUID metamaterial is capable of absorbing substantial amount of energy for frequencies within its linear frequency band, given by Eq. (36); however, the energy absorption in that band is far from being uniform. Thus, the term ”tunability of the resonance” for the SQUID metamaterial refers to that frequency at which the highest absorption of energy occurs (also note that for strong nonlinearities, the SQUID metamaterial can absorb significant amount of energies in other frequency bands as well, see next Section). Typical resonance tunability patterns are shown in Fig. 10 as density plots of as a function of the driving frequency (in nutural units) and the dc flux bias for several combinations of and . Note that the energy background in Fig. 10 has been removed for clarity. The thick lines with darker color depict the regions of the map with high energy absorption. The parameters used in the calculations are similar to the experimental ones Butz et al. (2013a); Trepanier et al. (2013); Zhang et al. (2015), although no attempt was made to exactly fit the observed patterns. These parameter values are consistent with a single-SQUID resonance frequency , which is also used in the calculations and to express the frequency in natural units. In Fig. 10, the ac flux amplitude increases from left to right, while the coupling increases from top bottom. The resonance become stronger as we move from the left to the right panels, as the nonlinear effects become more and more important with increasing . When going from top to bottom panels, with increasing , a smearing of the resonance is observed, along with the appearence of secondary resonances. The latter manifest themselves as thin red (dark) curves that are located close to the main shifting pattern, and they are better seen around half-intefger values of the applied dc flux. In order to obtain accurate tunability patterns, a few hundreds of absorbed energy-frequency curves (one for each ) have been calculated. Eqs. (40) are typically integrated with a fourth-order Runge-Kutta algorithm with constant time-step.