Two-tone spectroscopy of a SQUID metamaterial in the nonlinear regime

E. I. Kiselev, A. S. Averkin, M. V. Fistul, V. P. Koshelets, A. V. Ustinov Institut für Theorie der Kondensierten Materie, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany Physikalisches Institut, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany National University of Science and Technology ”MiSIS”, 119049 Moscow, Russia Center for Theoretical Physics of Complex Systems, Institute of Basic Science (IBS), Daejeon 34126, Republic of Korea Russian Quantum Center, Skolkovo, Moscow 143025, Russia Kotel’nikov Institute of Radio Engineering and Electronics, Moscow 125009, Russia (Dated: June 3, 2019)


I. INTRODUCTION
Superconducting metamaterials have a number of unique properties that are difficult or impossible to achieve in any other media [1]. A generic element of such a metamaterial, its elementary meta-atom, is a superconducting ring split by a Josephson junction. Each meta-atom thus forms a superconducting quantum interference device (SQUID), the most relevant circuit element of superconducting electronics. A Josephson junction acts as a nonlinear inductor, the inductance depending on the superconducting current flowing through it. SQUIDs are used in a variety of applications including sensitive measurements of magnetic fields and the detection of very small currents. In the microwave range, the unique properties of SQUIDs include the tunability of their intrinsic resonance frequency [2][3][4][5], extremely low losses, and bi-and multistable response [6][7][8][9][10], allowing for a wireless switching between opaque and transparent states of the metamaterial [9][10][11].
The nonlinear dynamics of Josephson junctions and SQUIDs plays an essential role in the development of parametric amplifiers [12][13][14][15][16] and even evokes interest in fields as remote as biology [17]. Parametric amplifiers are capable of amplifying one of the two quadratures of a signal without adding extra noise. Their performance is fundamentally limited only by quantum noise-an extraordinary property making them relevant for applications in various fields ranging from astrophysics to the readout of qubits [18]. Driven beyond the parametric regime by a strong microwave driving, a resonant Josephson junction circuit exhibits a bifurcation and enters the bistable regime with two stable dynamical states, which differ by both amplitude and phase. This behavior has been used to build the so-called Josephson bifurcation amplifier [19].
Recently, a lot of attention has been drawn to traveling wave parametric amplifiers (TWPAs), which can be realized with chains of Josephson junctions [20,21] or SQUIDs [14,15,22]. The SQUIDs of a TWPA chain have to be strongly coupled, so that collective oscillations can be induced. In our case, there is no galvanic connection between the SQUIDs. The coupling due to mutual inductance is very small, since the SQUIDs are spaced at large intervals of 92 μm-twice the width of a single SQUID. The distance to the central conductor of the coplanar waveguide (CPW) is approximately 10 μm. Thus, the SQUIDs behave as independent nonlinear oscillators and are excited solely by the signal transmitted through the CPW. This distinguishes our system from the Josephson junction and SQUID chains of Refs. [24][25][26], where a crucial role is played by collective modes. We study the resonant response of the SQUID metamaterial to a weak probe signal, while an additional strong driving (pump) signal pushes the metamaterial into a strongly nonlinear regime. In contrast to studies of intermodulation in SQUID metamaterials [23], where two signals of close frequencies induce sum and difference tones, here we consider signals with remote frequencies and are concerned with the strongly nonlinear behavior induced by the pump. The frequency of the pump FIG. 1. The main idea: the nonlinear response of the metamaterial to a large-amplitude driving tone (blue) is investigated by a weak probe signal (orange). The response to the weak signal is linear; however, it contains information about the large-amplitude response to the driving tone.
is fixed, while the weak signal response is measured over a broad band of frequencies. We note that previously the authors of Refs. [7,8] presented pioneering studies of the response of a single SQUID to a strong driving focusing on the amplitudefrequency and amplitude-phase characteristics. The analyses given in Refs. [7,8] resembles that in Ref. [9]. Here, by extending the experiment and the theoretical framework to two applied signals and performing a mathematical stability analysis of the solutions, we go beyond these works. We find previously unobserved regions of instability in the response of the SQUIDs. Some of these instabilities are associated with SQUID states that generate strong harmonics of the driving signal. Because of the large signal amplitudes, these effects go beyond the scope of a description through Kerr coefficients [25,26]. The key to their understanding lies in the dynamics of strongly driven nonlinear oscillators.
The paper is organized as follows: In Sec. II, we illustrate the main physical ideas behind our study. Section III describes the experimental setup, details of the measurement procedure, and presents our observations. In Sec. IV, we provide the theoretical analysis of the nonlinear dynamics of an rf-SQUID in the presence of both the weak probe signal and the large-amplitude driving tone. We start from revisiting a previously developed analytical approach to the nonlinear SQUID response [9], and then apply its formalism to the two-tone scenario. Section IV concludes the paper.

II. THE MAIN IDEA
The main idea behind our work is that the large-amplitude, steady-state response of a nonlinear system to a strong driving can be investigated by a weak probe signal, which is superimposed on the driving (see Fig. 1). The oscillatory response to the the weak probe signal is expected to be small and hence governed by a linear equation. The parameters of this equation (most relevantly the resonance frequency) will sensitively depend on the amplitude and frequency of the strong driving. Thus, by examining the weak signal response, we gain information about the nonlinear steady-state oscillations.
It is known that the response of (nonhysteretic) rf-SQUIDs-our meta-atoms-in a certain range of driving amplitudes and frequencies, has regions with multiple stable steady states, a feature known as multistability [9]. We find that a distinct resonance dip in the response to the probe signal corresponds to each of the stable states, allowing us to map the multistabilities. Furthermore, unstable regions, in which no stable steady-state solutions exist, are seen as gaps within which no resonant response to the weak probe signal is observed.

A. Experimental setup
The experiments were carried out at a temperature of 4.2 K using a one-dimensional array of 54 single-junction SQUIDs (rf-SQUIDs). In this setting, each SQUID constitutes a metamaterial atom embedded in a CPW. This is shown in Fig. 2(a). The SQUIDs [see Fig. 2(b)] were fabricated using the well-established Nb/AlO x /Nb trilayer process. The SQUID parameters are the same as reported in previous publications and experiments [4]. The value of the Josephson critical current is I c = 1.8 μA, leading to a zero-field Josephson inductance of L j,0 = 82.5 pH. The geometric loop inductance amounts L geo = 183 pH and the shunt capacitance is C = 2.0 pF. These characteristics result in a resonance frequency of approximately 14.8 GHz at zero magnetic field. The ratio of zero-field Josephson and geometric inductances is β L = 2π I c L geo / 0 ≈ 0.45, where 0 = h/(2e) is the magnetic flux quantum. The value β L < 0 indicates that the SQUIDs are in the nonhysteretic regime, thus having each a single stable static state at any magnetic flux.
A strong microwave pump signal generated by a microwave source (the driving tone) was inserted into the CPW to carry out the two-tone spectroscopy of the above-described sample. In the presence of this driving tone, we measured the transmission coefficient S 21 of a weak probe signal FIG. 3. The transmission spectrum of the metamaterial at a driving tone frequency f dr = 15 GHz. Driven into a strongly nonlinear regime by a large-amplitude driving tone, the SQUID metamaterial shows a resonant response to the weak probe signal of the VNA. The resonance frequency is seen as a minimum in the transmission coefficient S 21 . With increasing driving tone power, the resonance frequency drops first and then shows an oscillatory behavior. The blue curve shows the theoretically calculated resonance frequency as given by the harmonic approximation of Eq. (11), while the white curve shows the result obtained by solving Eq. (21).
propagating through the sample. A vector network analyzer (VNA) [see, Fig. 2(c)] was used for these measurements. Typical frequencies for the two signals were 1-10 GHz and 14-20 GHz for the driving tone and 10-16 GHz for the probe signal. No external magnetic field was applied.

B. Observations
We observed that the frequency-dependent transmission coefficient S 21 of the weak probe signal shows, even in the presence of a strong pump signal, a narrow resonant drop at the resonance frequency f res . The resonance frequency f res itself strongly varies with the power of the externally applied pump tone. Typical experimentally measured dependencies of S 21 ( f pr , P) on the frequency f pr of the weak probe signal and the power P of the driving tone are presented in Figs. 3-5, for driving tone frequencies of 15, 6, and 1 GHz, respectively.
For all driving tone frequencies, we observed a substantial decrease of the resonance frequency f res with increasing driving tone power. At high powers, the resonance frequency f res showed an oscillatory behavior.
An additional feature seen in Fig. 3 is the splitting of the resonance curves of different SQUIDs. This splitting is due to the slightly different strength of coupling of each single SQUID to the CPW. If the coupling of a SQUID is comparatively weak, the effect of the driving tone on this particular SQUID also will be smaller. Thus, its resonance curve will be shifted to higher values of the driving tone power.
At lower frequencies of the driving tone ( f dr 7 GHz), we observe gaps in the probe signal's resonant response. Such a gap is seen in Fig. 4 at 12 GHz and a driving power of FIG. 4. The transmission spectrum of the metamaterial at a driving tone frequency f dr = 6 GHz. Apart from the typical oscillatory behavior of the resonance frequency, we observe a gap in the resonance curve around 12 GHz, at approximately 0 dBm driving power. This gap indicates an unstable region in the response of the SQUID. In the vicinity of the gap, the resonance curve splits into two branches corresponding to two stable states that can be occupied by the SQUIDs-the response is bistable. The blue curve shows the resonance frequency as calculated using Eq. (11) and is not plotted in the region of instability predicted by Eq. (16). The white curve shows the resonance frequency obtained from Eq. (21). In the bistable region, the green curve, which was calculated using Eq. (18), depicts the resonance associated with a state that is characterized by a strong response at twice the driving frequency (see Sec. IV B 3). around 0 dBm. It indicates an instability in the response of the strongly driven SQUIDs, as will be explained below. In the vicinity of the gap, the resonance curve splits. This is due to the existence of two stable SQUID states in this range of parameters. Both states are occupied by some of the 54 SQUIDs of our metamaterial, allowing us to image the bistability directly.
Lowering the driving frequency even further, we observe resonant sidebands in the probe signal's transmission coefficient S 21 (see Fig. 5). These sidebands are located above and below the main resonance curve, at frequencies f res ± 2 f dr , where f dr is the frequency of the driving tone and f res is the resonance frequency at a given power of the driving tone.
All above observations can be explained by an analysis of the nonlinear rf-SQUID dynamics in the presence of a large amplitude microwave signal. The oscillations of the resonance frequency with the power of a driving tone are obtained by treating the response to the probe signal as that of a harmonic oscillator with parameters depending on the power of the driving tone. To obtain the resonance gaps and sidebands, we have to go beyond this harmonic model and turn to a more detailed treatment of the generic SQUID nonlinearity.

A. The nonlinear response of rf-SQUIDs to a harmonic driving
In order to explain the effects reported in the previous section, we assume that the interaction between adjacent 033096-3 FIG. 5. Transmission spectrum of the metamaterial at a driving tone frequency f dr = 1 GHz. Sideband resonances appear at f ms,0 ± 2 GHz, where f ms,0 is the resonance frequency at a given power of the driving tone. The blue curve shows the main resonance as given by Eq. (11). The white curve shows the resonance frequencies obtained from Eq. (21). Equation (21) includes the effect of sum and difference tones of the driving and probe signals. Including these sum and difference tones gives a natural explanation for the resonant sidebands.
rf-SQUIDs is negligible and that therefore the electrodynamic response of the superconducting metamaterial is determined by the response of individual SQUIDs. The dynamics of an rf-SQUID is governed by a simple cirquit model: the resistively and capacitively shunted Josephson junction model [27] in parallel with the geometric inductance of the SQUID ring L geo [see Fig. 2(b)]. In this model, the time-dependent Josephson phase ϕ(t ) satisfies the nonlinear equation where ϕ ext consists of both an externally applied driving tone and a probe signal, β L = L geo /L j , L j = 0 /2π I c with I c being the critical current and R being the normal resistance of the Josephson junction. We start by considering the nonlinear dynamics of an rf-SQUID in the presence of a strong driving tone ϕ ext = ϕ dr sin(ω dr t ). Approximate solutions to Eq. (1) are obtained by exploiting the fact that the response of an rf-SQUID to a driving tone of frequency ω dr will mostly be dominated by harmonic oscillations at the same frequency ω dr . This approach has already been successfully applied to a study of dynamic metastable states that can be excited in an rf-SQUID metamaterial [9] and uses the following ansatz: The validity of this ansatz has been analyzed in the methods section of Ref. [9]. It was found that the ratio between the first and third harmonics is numerically small, which justifies the monochromatic approach of Eq. (2). Notice that ϕ dr ∝ √ P for the driving tone amplitude holds, where P is the power of the driving tone. When we insert the ansatz of Eq. (2) into Eq. (1), the nonlinear term of Eq. (1) is expanded into a Fourier series with the help of a Jacobi-Anger identity [28]: sin (ϕ a sin (ω dr t + δ)) where J n (x) are Bessel functions of nth order. All terms of this expansion except for n = 0 are neglected according to our assumptions and Eq. (1) reads where ω c = R/L geo and the geometrical frequency is given by ω geo = (L geo C) −1/2 . By rewriting Eq. (4) in terms of sin (ω dr t ) and cos (ω dr t ) and comparing the coefficients of these functions, we obtain an equation for the amplitude of the driving tone ϕ dr in terms of the amplitude of the Josephson phase response ϕ a : For the SQUIDs used here, the resistance R in is of the order of 10 k and hence ω c ≈ 50×10 12 /s. Thus, to a good approximation, the last term in Eq. (5) can be neglected. From Eq. (5), one can see that if the SQUID is driven at a frequency near the geometric resonance, i.e., ω dr ≈ ω geo , we have ϕ dr ≈ 2β L |J 1 (ϕ a )|. Since J 1 (ϕ a ) is an oscillating function, the inverse function ϕ a (ϕ dr ) is multivalued. This phenomenon is referred to as multistability and was experimentally observed by Jung et al. [9], as well as by Refs. [7,8], however, without putting emphasis on it. The response amplitude of a SQUID as a function of the driving amplitude ϕ a (ϕ dr ) for ω dr ≈ 0.998ω geo is shown in Fig. 6. Some parts of the curve ϕ a (ϕ dr ) are unstable in the sense that any small perturbation to the amplitude ϕ a will cause a change of the state of the SQUID. These parts are marked with red dotted lines in Fig. 6.
To investigate the stability of the solutions of Eq. (5), we followed a standard procedure [29], where the time evolution of a small perturbation to the phase and amplitude of the SQUID response as given by Eq. (2) is studied. The response is written ϕ a sin (ω dr t + δ) = x sin ω dr t + y cos ω dr t with ϕ a = x 2 + y 2 , δ = arctan(y/x) and the perturbations λ i (t ) are added: x → x + λ 1 (t ), y → y + λ 2 (t ), λ i (t ) 1. Linear differential equations are then derived for λ i (t ) from Eq. (1) assuming the validity of Eq. (5). In unstable regions, the perturbations λ i (t ) show an exponential growth, while stability is assured if only oscillatory solutions exist for λ i (t ). The results agree with the general expectations [9,23,29] that those sections of ϕ a (ϕ dr ) are unstable, where the curve's slope is negative.

B. Two-tone resonant spectroscopy of rf-SQUIDs
In this subsection, we analyze the resonant response of the superconducting metamaterial to both the weak probe signal and the strong driving tone. First, we will treat the response within a harmonic approximation; i.e., we will approximate the probe signal response as that of a harmonic oscillator with parameters depending on the driving tone amplitude. This approximation is sufficient to explain the oscillations of the SQUID's resonance frequency f res with the increasing power of the driving tone. Next, the sidebands and unstable regions seen in the experiments will be described within a more detailed analysis.

Harmonic approximation
The two signals applied to the SQUID are the driving tone ϕ dr sin (ω dr t ) and the probe signal of the VNA ϕ pr,ext . The response is assumed to be of the form ϕ = ϕ a sin (ω dr t + δ) + ϕ pr . Equation (1) then takes the following form: For a weak probe signal, the response ϕ pr will be small. Expanding up to the first order in ϕ pr , the nonlinear term of (6) becomes sin (ϕ a sin(ω dr t + δ)) + cos (ϕ a sin(ω dr t + δ))ϕ pr . (7) At this point, a key assumption enters the calculation: Since the VNA signal is very weak and |ϕ pr | ϕ a , its presence does not significantly alter the state of the SQUID. It is hence assumed that the Eq. (4) is still valid. By substituting Eqs. (4) and (7) into Eq. (6), we obtain an effective equation for the response ϕ pr to the VNA signal: ϕ pr,ext = ϕ pr + ω −2 geoφ pr + ω −1 cφ pr + β L cos (ϕ a sin(ω dr t + δ))ϕ pr .
This equation with a time-dependent, periodic coefficient is known as Hill's equation [30]. Solutions of the homogeneous part of this equation can be unstable in the sense that their amplitude exponentially grows in time [30,31]. These instabilities, however, only occur for specific values of ω dr , ϕ a , and ω c and will be studied in next subsection. Here, solely the special solution to (8) is of interest, and only the time average of cos (ϕ a sin(ωt + δ)) over one period, i.e., the zeroth Fourier mode, will be taken into account. We use the Jacobi-Anger identity [28] cos (ϕ a sin (ω dr t + δ)) and approximate Eq. (8) by The frequency of the resonant drop of the probe signal's transmission coefficient S 21 then is This result gives the shift of the resonance frequency as a function of the amplitude of the driving response ϕ a [see Eq. (5)]. In the experimental data, however, the transmission coefficient S 21 is plotted against the power of the driving tone, P. We thus introduce a coupling constant C by means of When fitting Eq. (11) to experimental data, C is the only free parameter. It has to be obtained independently for every driving tone frequency, since the transmission through microwave cables and devices is frequency dependent. In Figs. 3-5, the SQUID resonance frequency predicted by Eq. (11) is plotted blue. Unsurprisingly, the analysis presented so far describes the experiment best at high driving tone frequencies (Fig. 3), since in this regime the time-dependent terms in the expansion of Eq. (9) are expected to average out. However, even at low driving tone frequencies Eq. (11) is in a good accord with experimental results (Figs. 4 and 5).

Parametric instabilities of the SQUID response
At certain values of the driving frequency ω dr and the driving amplitude ϕ dr , the solution to the homogeneous part of Eq. (8) becomes unstable [30]. To describe this phenomenon, called parametric instability or parametric resonance [31], we take into account the n = 1 term in the expansion of Eq. (9), which oscillates with the frequency of 2ω dr . The effective probe signal equation (8) then becomes ϕ pr,ext = [1 + β L J 0 (ϕ a ) + J 2 (ϕ a ) cos (2ω dr t )]ϕ pr 033096-5 It is well known [30,31] that the parametrically unstable regions of such a Mathieu-type equation occur at driving frequencies of where m is an integer number and ω ± is a small frequency interval whose value depends on the driving tone amplitude ϕ a . By applying the stability criterion of Ref. [31] to Eq. (13), estimations on ω ± are obtained. For m = 1, we have .
The m = 1 instability occurs when ω ≈ ω res ≈ ω geo holds for the driving frequency. It is therefore not surprising that the intervals given by Eq. (15) coincide with the instability regions already discussed in the Sec. IV A. For the m = 2 instability, however, we observe a different behavior. In this case, the interval of instability is not symmetric: (16) Figure 4 shows the experimental data obtained at a driving tone frequency of 6 GHz. The m = 2 instability is clearly seen as a gap in the resonance curve around 12 GHz at a driving power of approximately 0 dBm. The gap in the blue curve which shows the resonance frequency as given by Eq. (11) indicates the unstable region predicted by Eq. (16). Furthermore, higher order instabilities were observed for m = 3 and m = 4 at driving frequencies of 4 and 3 GHz respectively. While reaching an instability in the m = 1 case just indicates a jump to another stable state predicted by Eq. (5), in the higher order cases (5) does not predict any multistabilities. Thus, the m = 2, 3, . . . instabilities suggest that the SQUID response in the region of these instabilities is more complicated than assumed in deriving Eq. (5). In fact, we demonstrate below that the resonant drop seen in Fig. 4 above the gap and at powers above 0 dBm indicates that a fraction of the SQUIDs occupies a state where harmonics with frequencies ω dr and 2ω dr are approximately equally strong. Such a state indeed goes beyond the approximations, leading to (5).

A new, symmetry-breaking bistability in the SQUID response
To explain the bistability observed at a driving frequency of 6 GHz (see Fig. 4), it is necessary to go beyond the approximate ansatz of Eq. (2). The bistability appears when, while increasing the driving amplitude, the resonance frequency f res is tuned to twice the driving frequency f dr . At the same time, the frequency of the second harmonic of the response with a frequency of 2 f dr , too, becomes approximately equal to the resonance frequency. The amplitude of the second harmonic can therefore be expected to grow in this regime. In fact, numerics show that the resonant drop in Fig. 4, which is not described by Eqs. (11) or (21), results from a SQUID state where the amplitudes ϕ a1 , ϕ a2 , of the first and second harmonics with frequencies ω dr and 2ω dr have the same order of magnitude. We therefore use the ansatz ϕ(t ) = ϕ a1 sin (ω dr t + δ 1 ) + ϕ a2 sin (2ω dr t + δ 2 ) (17) to describe this state. The frequency of the resonant drop is obtained analogously to (11) and given by The amplitudes ϕ a1 and ϕ a2 are obtained numerically as Fourier coefficients of the first and second harmonics of the solution to Eq. (1). For the numerical analysis, we used the standard fourth-order Runge-Kutta method and generated a solution ϕ(t ) of Eq. (1). To find the amplitudes ϕ a1 and ϕ a2 , we projected the solution ϕ(t ) onto the Fourier modes, e.g., where T = 2π/ω and n is the number of periods over which the solution was averaged. The green solid line in Fig. 4 shows the resonance as described by Eq. (18) in comparison with experimental results. The fact that a strong second-order harmonic can be generated by the SQUIDs without a dc magnetic field being applied is rather remarkable. Such a state dynamically breaks the ϕ → −ϕ centrosymmetry of Eq. (1). Dynamical symmetry breaking is a phenomenon known from the theory of nonlinear oscillators [32,33]. In fact, the observation of this bistability demonstrates the advantages of the two-tone spectroscopy over frequency-amplitude and frequency-phase measurements, in which it cannot be imaged with such ease. We expect that other irregularities in the weak signal response, e.g., the ones seen in Fig. 5, can similarly be explained in terms of states with a strong response at harmonics of the driving tone. Figure 5 shows the response to the probe signal at a low driving tone frequency ω dr ω pr . We observe sideband resonances approximately 2 GHz above and below the main resonance. These sidebands appear due to the term J 2 (ϕ a ) cos (2ω dr t ) in the probe signal equation (13). The term mixes the oscillation at ω pr with modes at frequencies separated from it by a multiple of 2ω dr . To calculate the positions of the sidebands, we therefore use the ansatz ϕ(t ) = a sin (ω pr t ) + b 1 sin (ω pr t + 2ω dr t )

Sideband resonances
and ϕ pr,ext = ϕ pr sin (ω pr t ). To describe the sidebands nearest to the main resonance, modes with frequencies further away from ω pr , i.e., ω pr ± 4ω dr , ω pr ± 6ω dr , and so on can be  neglected. Since the dissipation is rather small, one can neglect phase shifts δ in the above ansatz. We substitute (20) into (13) and solve for a(ϕ pr ). With dissipation being neglected, resonances occur at those frequencies ω pr at which a(ϕ pr ) ∝ ϕ pr diverges, or, equivalently, ϕ pr (a) = 0.
Solving Eq. (21) involves finding the roots of a sixth-order polynomial in ω pr , which is done numerically. For most values of the driving tone frequency ω dr and power P, Eq. (21) has three positive roots, corresponding to the main resonant response and two sidebands. The resonance and sideband curves predicted by Eq. (21) for the driving tone frequency of 1 GHz are plotted in Fig. 5 (the white solid line) and are in good agreement with the experiment. For higher driving tone frequencies (Fig. 4), the sidebands are outside of the experimentally accessible frequency range. The main resonance curves given by Eq. (21) (15) and (16).

V. CONCLUSIONS
We presented the results of a two-tone spectroscopy of an rf-SQUID metamaterial consisting of 54 single SQUIDs placed in a transmission line. The power of the pump (driving) tone was typically much higher than the power of the probe signal used to measure the transmission spectrum.
We observed that the resonance frequency of the metamaterial, seen as a drop in the transmission spectrum of the probe signal at frequencies between 10 and 15 GHz, shows a characteristic oscillatory dependence on the power of a pump tone. The frequency of the pump tone was varied between 1 and 20 GHz, changing the shape of the resonance curve.
For pump tone frequencies below or in the region of the resonance frequencies, i.e., for our parameters f dr < 14.5 GHz, we observed the gaps in the resonance curves of the transmission spectrum. A bistability in the response was directly observed at an intermediate range of the pump frequency (6 GHz). At low pump frequencies f dr sidebands located approximately 2 f dr above and below the main resonance curve were found.
Most of our observations are well described theoretically by an approach based on an approximate analytical model for the response of a SQUID to strong external driving [9]. Whereas the shape of the resonance curves naturally followed from an extension of the model to two tones of different amplitudes, the gaps could be explained as parametric instabilities of the SQUID response. The observed bistability is shown to appear as a result of the existence of a symmetry-breaking state with strong second-harmonic generation.
We believe that the described effects have multiple applications. The instability-induced gaps can be used to create metamaterials with tunable transparancy. The symmetry breaking second-harmonic generation can be employed in frequency doublers. Last but not least, the remarkably clear small signal response at even the highest pump powers shows that a study of parametric amplification in the highly nonlinear regime could be a worthwhile endeavor.
To summarize, at high driving power levels SQUID metamaterials exhibit a rich spectrum of features tunable by the power and the frequency of the pump tone. The nonlinear twotone spectroscopy is a powerful method for studying these features and can be applied to other systems. The high degree of nonlinearity and the unique multistable behavior make superconducting circuits an exciting playground for studies of tunable metamaterials.