Entangling interactions between artificial atoms mediated by a multimode left-handed superconducting ring resonator

Superconducting metamaterial transmission lines implemented with lumped circuit elements can exhibit left-handed dispersion, where the group and phase velocity have opposite sign, in a frequency range relevant for superconducting artificial atoms. Forming such a metamaterial transmission line into a ring and coupling it to qubits at different points around the ring results in a multimode bus resonator with a compact footprint. Using flux-tunable qubits, we characterize and theoretically model the variation in the coupling strength between the two qubits and each of the ring resonator modes. Although the qubits have negligible direct coupling between them, their interactions with the multimode ring resonator result in both a transverse exchange coupling and a higher order $ZZ$ interaction between the qubits. As we vary the detuning between the qubits and their frequency relative to the ring resonator modes, we observe significant variations in both of these inter-qubit interactions, including zero crossings and changes of sign. The ability to modulate interaction terms such as the $ZZ$ scale between zero and large values for small changes in qubit frequency provides a promising pathway for implementing entangling gates in a system capable of hosting many qubits.


I. Introduction
In the field of circuit quantum electrodynamics (cQED), linear resonators are commonly used for readout of artificial atom qubits, and in some cases, for coupling between qubits [1,2].Variants on this architecture are the dominant paradigm for current superconductor-based quantum processors.The resonators are often formed from planar transmission lines with a single mode near the frequency range of the qubits, resulting in physically large footprints of several mm.Utilizing an even larger physical footprint, devices with dense mode spectra near the qubit frequency range have been realized using ultralong linear resonators [3].Similar multimode cQED systems have been studied for implementing quantum memories [4] and quantum simulations [5,6].
The use of metamaterials formed from lumped-element inductors and capacitors allows for the implementation of transmission lines with unconventional wave dispersion.In the case of left-handed dispersion [7], the wave frequency is a falling function of wavenumber above an infrared cutoff frequency f IR , below which waves are unable to propagate [8,9].In the context of cQED, lefthanded metamaterials produce a dense spectrum of orthogonal microwave modes above f IR , which can be engineered to fall in the frequency range of conventional superconducting qubits [10,11].
Ring resonators have been used in integrated photonics systems to form compact optical resonances, or whispering gallery modes, for a broad range of applications * bplourde@syr.eduincluding microwave-optical transducers, microwave frequency combs, and multimode nonlinear optics [12][13][14][15][16]. Superconducting ring resonators with right-handed dispersion have been used in cQED applications resulting in novel properties [17,18], but require a large footprint to ensure a minimum of one wavelength matches the circumference.
In this manuscript, we describe a superconducting ring resonator formed from a left-handed metamaterial transmission line with two transmon qubits coupled at different points around the ring.The left-handed dispersion results in a dense spectrum of modes in the vicinity of the transition frequencies for the qubits.By wrapping the left-handed transmission line into a ring structure, we produce a device with a relatively compact footprint with unique properties that arise due to the symmetry of the ring.We perform a detailed modeling of the standingwave structure and degeneracy breaking in the ring resonator, allowing us to predict the coupling energy scales between the qubits and each ring resonator mode [19].The multimode coupling between the qubits with the ring resonator serving as a bus results in significant variations in both the transverse exchange coupling between the qubits as well as higher order ZZ interactions as the qubits are tuned between different frequency regimes.
Our theoretical modeling of these interactions is in close agreement with our experimental measurements.The ability to modulate the ZZ interaction strength between zero and a large value for a small change in qubit frequency allows for the possibility of implementing fast, high-fidelity entangling gates between pairs of qubits located around the ring resonator.N is depicted with solid (dashed) green lines, labeled m The coupling strength of QA and QB to the mode is determined by the amplitude of the wave at the cell where each qubit is capacitively coupled to the resonator.(b) Theoretical g-coupling values for QA and QB, coupled to a hypothetical 240-cell ring resonator with qubit separation nAB = 6.Simulated data points at which the coupling for QA, g A i /2π, and for QB, g B i /2π, share the same sign for a given mode are shown in purple and follow Eq.(3).Points at which the g A i and g B i coupling to QA and QB, respectively, have opposite signs are shown in gray and are associated with Eq. (4).

II. Ring resonator modes and qubit coupling: theory
A resonator formed from a conventional transmission line (TL) can be modeled as an array of LC unit cells, each with inductors shunted by capacitors to ground [20].For this discrete model, a linear TL of length d l is partitioned into N unit cells, each with size ∆x = d l /N .The cell capacitance (inductance) are C L = c L ∆x (L L = l L ∆x) with c L (l L ) being the respective quantities per unit length.In the continuum limit N → ∞, the lumpedelement model reproduces the linear dispersion of transverse electromagnetic modes ω(k) = k/ √ c L l L , with k being wave number.This configuration results in a righthanded resonator, as its wave, electric, and magnetic vec-tors make a right-handed set in three dimensions.However, switching the positions of the capacitors and inductors results in a left-handed dispersion relation [8], with ω decreasing for increasing k: In an ideal left-handed metamaterial, varying wave number in the domain of |k∆x| ≤ π shows no upper bound on frequency in the k → 0 limit.In actual devices, a finite upper bound results because these modes hybridize with parasitic right-handed modes.Additionally, left-handed metamaterials exhibit a lower bound at |k∆x| → π, referred to as the infrared cutoff frequency ω IR = 1/2 √ C L L L , which is only slightly perturbed with parasitic modes [21].In between the upper and lower frequency bounds there are a finite number of modes, with those close to the infrared cutoff being densely grouped, and mode separation increasing monotonically with higher frequencies [10].
Next, we consider a left-handed ring resonator design made by connecting the two ends of a left-handed TL.This enforces periodic boundary conditions, similar to Born-von Karman boundary conditions in solid state physics [22].This condition highlights two types of traveling waves in the ring resonator: clockwise and counterclockwise moving waves.Given that the dispersion relation in Eq. ( 1 √ 2, respectively.The parity is indistinguishable at k = 0. Similarly, for the mode at the IR cutoff, the difference between the clockwise and counterclockwise argument k∆x is 2π, thus the two parities are indistinguishable.This is akin to states at the edge of the Brillouin zone in solid state systems.A total of N/2 − 1 degenerate parity pairs can be found in an ideal ring resonator.Figure 1(a) shows two pairs of such degenerate modes in the form of standing waves with a relative phase shift of π/2.In the subsequent sections, we will show that any source of asymmetry between cells or local external coupling to the ring lifts the degeneracy.
Let us consider a ring resonator that is capacitively coupled to two weakly-anharmonic qubits, Q A and Q B , separated by n AB cells around the ring [Fig.1(a)].Using the standard circuit quantization within the rotating wave approximation, the total Hamiltonian reads [23] Approved for Public Release; Distribution Unlimited.PA #: AFRL-2024-0930.
with â † q (â q ) being the creation(annihilation) operator of qubit q = A, B and δ q being the qubit anharmonicity.Here r † k,P (r k,P ) are the creation(annihilation) operator of ring resonator mode k and parity P = E, O, and the qubit-resonator mode coupling strengths g q P (k) are each proportional to the amplitude of the standing wave for the particular mode at the qubit location, as shown in Fig. 1(a), and therefore depend on k and parity P .The even-and odd-parity coupling strengths can be expressed as: with x AB = n AB ∆x being the distance between the qubits (Appendix A).Here, the justification for the terms 'even' and 'odd' parity becomes clear, as even modes couple with the same sign to qubits Q A and Q B , while odd modes couple with opposite signs to the two qubits.The coupling strengths depend periodically on the product of k and x AB .In Fig. 1(a), the solid line is a cosine and the dashed line is a sine mode.The origin is located in the middle between the qubits, thus clarifying the signs in Eqs.(3,4), as Q A is connected to cell number n AB /2, and shows the oscillatory nature of the gcouplings at different mode frequencies.For better visualization of the oscillations, we consider a ring with N = 240 unit cells, but keep the qubit-qubit cell separation n AB = 6.Since the cosine and sine functions do not depend on N , this does not affect the oscillations in coupling strength.We compare the analytical expressions of Eqs.(3,4) with the numerical simulations that predict the same odd (gray) and even (purple) parity behavior.Note that the parity of the g-couplings is crucial for calculation of the effective J-coupling between qubits presented later in the manuscript.

III. Device Layout
The device consists of a 24-cell, left-handed superconducting metamaterial ring-resonator comprised of interdigitated capacitors with double-sided meander-line inductors to ground (Fig. 2).The ring resonator is capacitively coupled to a feedline (R in /R out ) [Fig.2(a)] for probing ring-resonator modes.Two flux-tunable, floating-style asymmetric transmon qubits [24], Q A and Q B , are capacitively coupled to the ring resonator 90 and 180 degrees from the feedline coupling point, respectively.Both qubits have readout resonators inductively coupled to a separate feedline (F in /F out ), as well as separate onchip flux-bias lines for tuning the transition frequency (Appendix B and Appendix C).
Figure 2(e) shows the location of aluminum wirebonds used to connect the central superconducting disk inside Theoretical mode frequencies are shown with green dasheddotted lines when stray inductance due to wirebonds is not included and any degeneracy lifting is due to the qubits or feedline.Solid blue lines show large degeneracy lifting effect of wirebonds on the ring resonator mode frequencies.
the ring resonator to the ground plane of the chip.These ground bonds have non-negligible inductance, which we estimate to be 1 nH/mm [25].The bonds are not spaced symmetrically due to the adjacent locations of the qubits and feedline.These wirebonds break the symmetry of the clockwise and counterclockwise propagating waves and lift the degeneracy of the even and odd ring-resonator modes.The infrared cutoff frequency for this ring resonator device is 4.287 GHz, as shown in Figure 2(e), below which we observe no modes.Between 4-6.5 GHz, we measure thirteen modes to which the qubits may couple.As expected, the modes above the infrared cutoff come in pairs that have broken degeneracy due to the presence Dashed lines in red (blue) show fits to splitting data for QA (QB) obtained from reduced Hamiltonian including one qubit and the three modes (see Appendix F).Horizontal dotted lines show dressed mode frequencies.(c-f) Splittings and reduced Hamiltonian fits for two single modes for each qubit.For all measurements, the spectator qubit is fixed at its upper flux insensitive sweet spot (Φq/Φ0 = 0).As shown in (c), QA has effectively zero coupling to the ring resonator mode at 4.553 GHz, while (d) shows QB has a coupling strength of 26.6 MHz.(g) Extracted and theoretical |g q i |/2π values for each ring-resonator mode i. Error bars on experimental data computed from 95% confidence intervals in g-coupling fits (see also Appendix F).Vertical dashed lines show bare ring resonator mode frequencies.
of the qubits, the feedline, and the wirebonds.The size of the degeneracy lifting ranges from ∼30 MHz near the IR mode to ∼250 MHz at higher frequencies and is in good agreement with theoretical simulations including a model for the wirebonds, as shown in Fig. 2(e) in blue.The same simulation without wirebonds fails to reproduce the observed lifting of degeneracies (Appendix D).

IV. Mode structure and g-coupling
To characterize the coupling strength between each qubit and the different ring resonator modes, we probe the modes by measuring the microwave transmission through the ring resonator feedline (R in /R out ) while scanning the flux for one of the qubits with the other qubit flux-biased at its upper flux-insensitive sweetspot (Appendix E includes details on cabling and readout).We perform these measurements for each of the thirteen modes between 4-6.5 GHz, and observe vacuum Rabi splittings when each qubit hybridizes with the ring-resonator mode photonic state.A selection of these measurements for Q A and Q B is shown in Fig. 3[a-f].To extract the magnitude of g A i (g B i ) for Q A (Q B ), we perform a least-squares minimization to fit each set of splitting data to a reduced Hamiltonian for the ring-resonator mode-qubit system, shown as dashed lines.In regions where the mode spacing and coupling are comparable, we simultaneously fit multiple modes to extract each g A i and g B i , as shown in Fig. 3(a,b) (Appendix F).Near f IR , the mode spacing is 23 MHz between the IR-cutoff mode and the next mode, while g A i /2π = 13 MHz and g B i /2π = 36.5MHz for the IR-cutoff mode.For Q A , one ring-resonator mode near 4.55 GHz [Fig.3(c)] has effectively zero coupling, while for Q B , the extracted magnitude of g B i /2π is 26.6 MHz to the same mode [Fig.3(d)].While theory predicts zerocoupling to the modes at the roots of Eqs.(3,4), this particular mode is not a root.In the simple model without the feedline and wirebonds, this mode has a nonzero coupling strength to Q A ; only in the more complex model does this mode result in a negligible coupling strength (Appendix D).
Experimental and theoretical values for g A i (g B i ) for Q A (Q B ) and ring-resonator modes from 4 GHz to 6.5 GHz are shown in Fig. 3(g).Theoretical calculations include a model of the feedline and the wirebonds, the introduction of which requires the use of numerical calculations instead of the analytical expressions in Eqs.(3,4).Despite the added complexity, we can still observe features of the sine and cosine behavior of g A/B i .For example, we find a pair of modes at 4.344 GHz and 4.475 GHz; the former having a large coupling strength, since it is an odd mode and Eq. ( 4) has a maximum here, while the latter is even and has negligible coupling, following Eq.( 3).
We observe a broad array of coupling magnitudes ranging from 0-80 MHz.While the theoretical model captures the general behavior of g A/B i , a more complete model of the symmetry-breaking perturbations of the ring-resonator circuit is required for better quantitative agreement.Although we cannot directly measure the sign of g, the magnitude and parity of the pairs of g A i and g B i for Q A and Q B to each mode i impacts the qubit-qubit interactions that are mediated by the ring resonator.We utilize these extracted g-couplings for each of the modes in the subsequent sections for analyzing our measured qubit-qubit interactions.

V. J-coupling between qubits
We next explore the transverse exchange coupling between the qubits mediated by virtual photon exchange [26] with the various ring-resonator modes.By using the qubit readout resonators to perform conventional qubit spectroscopy and the local flux-biasing capability to independently tune the qubits, we fix the frequency of one qubit, say Q B , while flux-tuning the other qubit (Q A ).
A nonzero J-coupling between the qubits results in an anti-crossing in the spectrum when the bare qubit frequencies cross.We can adjust the frequency of the crossing point for the bare qubit frequencies relative to the ring-resonator modes and study the variation of the exchange coupling between the qubits (Fig. 4).
To extract the magnitude of J, we fit the spectroscopy data with a reduced Hamiltonian model in which the qubits are restricted to two states and only the nearest ring-resonator modes are accounted for (Appendix G).In the spectroscopy measurements in Fig. 4, we observe a small exchange coupling of 2.1 MHz when the qubits are biased below f IR [Fig.4(a)], while the J-coupling is more than an order of magnitude larger when the qubits are biased in the dense portion of the ring-resonator mode spectrum [Fig.4(b)].
Figure 4(c) shows theoretical predictions for J compared to extracted values for a range of qubit frequencies.In the perturbative Schrieffer-Wolff approximation, J is given by a sum over all modes where ∆ A i (∆ B i ) give the detuning for Q A (Q B ) to each ring resonator mode i.Since J depends on the product g A i g B i for each mode, even-parity modes will have a positive contribution to the sum and odd-parity modes will have a negative one, however this can be changed by the sign of the detuning.When the qubits are biased below f IR and all ∆ q i have the same sign, J will be small, as shown in Fig. 4(c).Although our spectroscopy measurements only provide the magnitude of J, using the perturbative theoretical calculations with different possible parities of g A i and g B i for each mode i, we determine the sign of each qubit-mode coupling that provides the best agreement with the exchange coupling measurements (Appendix G).From the experimentally measured J values and the agreement with the theoretical modeling, it is clear that the exchange coupling can vary in magnitude between zero and at least 41 MHz and can change sign depending on the detuning between the qubits and the various ring-resonator modes.Since the mode spectrum is dense and interactions to modes may be strong, perturbation theory is not guaranteed to be a good choice.Therefore we also compare our results to the non-perturbative Least Action method [27][28][29] at selected frequencies.This method agrees well with both experimental and perturbative data, except for cases where the qubit and resonator mode frequency are quite close, in which case the perturbative method diverges (Appendix H).

VI. ZZ interactions
In general, interactions between pairs of qubits and resonator modes coupled to the qubits lead to ZZ interac- tions between the qubits, where the state of one qubit shifts the transition frequency of the other qubit.Such interactions can be problematic in multi-qubit systems and can generate unwanted and uncontrolled entanglement [28,30].Various approaches have been explored for reducing [31,32] or nulling [27,33] these interactions.On the other hand, the ZZ coupling can also be used for implementing two-qubit entangling gates, provided the interaction strength can be modulated [34][35][36][37].
In the dense multi-mode spectrum of our metamaterial ring resonator system, we explore the ZZ interaction between the two qubits for multiple qubit detunings within different regions of the mode spectrum.We use a standard Ramsey interferometry technique to extract the ZZ interaction strength [38].First, we bias Q B to a particular transition frequency, then we adjust the bias of Q A to various frequencies relative to Q B and the mode spectrum.At each bias point for Q A , we perform a standard Ramsey fringe sequence while stepping through the drive frequency for the X/2 pulses for Q A , thus allowing us to identify the transition frequency for Q A .We then per- In Fig. 5(a,b) we show two series of ZZ measurements as a function of the dressed transition frequency for Q A for two different bias points of Q B .We observe that ζ can vary over a wide frequency range, covering both positive and negative values and crossing through zero, both smoothly in some regions and discontinuously in others.The theoretical curves are calculated by diagonalizing the Hamiltonian and using the definition of the ZZ interaction: , where E ij corresponds to the energy eigenvalue of state i(j) of Q A (Q B ).If two states come into resonance, this definition can be nontrivial, as the identification of which state is, for example, the |11⟩ state and which is a non-computational state, e.g., |20⟩, becomes more difficult (Appendix J).
We observe excellent agreement between our experimental results and theoretical calculations of the variation in the ZZ interaction, capturing both the regions where ζ changes sign smoothly and where it jumps discontinuously.For our multimode ring resonator, a quite small change in the frequency of one qubit, of the order of a few MHz for Q A in this case, results in a change in ζ/2π from zero to tens of MHz.This significant change is due to the large effective coupling between the two qubits, which is mediated by the high mode density.This large coupling causes a strong repulsion of the |11⟩ and |20⟩ states, leading to large changes in ζ (Appendix I).

VII. Discussion
We have demonstrated that the metamaterial ring resonator exhibits a broad range of bipolar interactions between transmon qubits and the ring resonator modes, as well as interqubit interactions mediated by the ring resonator.The symmetry of the ring resonator geometry leads to even or odd parity of the g-coupling for each mode to the two qubits, a property that is unique to the ring geometry, compared to a linear left-handed transmission line configuration.If the degeneracy between the opposite parity g-coupled modes is broken by an asymmetry in the device layout, tuning the qubits between such a pair of modes can lead to the cancellation of ZZ interactions between the qubits.In addition, due to the dense spectrum of the ring resonator modes and the different parities of the g-couplings, there is also a large variation of the J-coupling between qubits as their frequency is varied.Therefore, we can make a large change in both the J-coupling and the ZZ interaction through relatively small changes in qubit frequency.
The large range of tunability for interqubit interactions suggests a potential path for future work based on this device design investigating ZZ-based two-qubit gates.Prior work on two-qubit gates based on in situ modulation of the ZZ interaction using tunable couplers has required tuning of at least one circuit element over a frequency range of the order of GHz [39,40].Fast tuning over such a large range introduces the risk of fidelity degradation through leakage to other modes, both intentional and spurious, within this frequency window.By contrast, the orders of magnitude smaller tuning range required for moving between the on and off regimes of the ZZ interaction for our multimode ring resonator provides a promising pathway for implementing a highfidelity two-qubit entangling gate.Gate errors due to Purcell losses caused by coupling between the qubits and the multiple nearby modes could be reduced by eliminating the separate feedline coupled to the ring resonator, thus making the modes limited by internal rather than external losses.
This platform could also be used for implementing a non-ZZ-based two-qubit gate, for example, based on cross resonance [41,42].By operating the qubits at frequencies where the ZZ-interaction is minimized but the J-coupling is large, one could potentially perform a fast cross resonance entangling operation while minimizing gate errors caused by the parasitic ZZ interaction that occurs in general for systems of coupled transmons [28,43].As an example, for this device, as shown in Fig. 5(a), for a dressed frequency for Q A of 5.195 GHz and dressed frequency for Q B of 5.04 GHz, the ZZ interaction has a zero-crossing; at these same qubit frequencies where ζ/2π = 0, the J-coupling computed from Eq. ( 5) is ∼11 MHz.
In addition, based on the mode structure presented above, the compact ring geometry allows for coupling to more than two qubits, with the possibility of tuning ZZ interactions or J-couplings between different pairs of qubits with selective flux control of each qubit.Coupling a larger number of qubits around the ring, or even coupling together several ring resonators, each coupled to multiple qubits, could be a path for creating a multi-qubit system with long-range and tunable coupling between physically distant qubits.Prior work using a ring resonator architecture for long-range coupling between multiple qubits utilized a 3D ring resonator with righthanded dispersion and required more than an order of magnitude larger physical footprint [18].Our left-handed ring resonator geometry with its multimode spectrum can thus serve as a compact multiqubit bus for an architecture with highly tunable interactions between distant qubits.

Acknowledgments
This work is supported by the U.S. Government under AFOSR grant FA9550-21-1-0020.Portions of this work are supported by AFRL grant FA8750-20-1-1001.Fabrication was performed in part at the Cornell NanoScale Facility, a member of the National Nanotechnology Coordinated Infrastructure (NNCI), which is supported by the National Science Foundation (Grant NNCI-2025233).Any opinions, findings, and conclusions or recommendations expressed in this article are those of the authors and do not necessarily reflect the views of the Air Force Research Laboratory (AFRL).

Appendix A. Theoretical Hamiltonian and derivation of g-coupling
We can derive an analytical expression for the g-coupling values of both qubits from the Lagrangian of the circuit where we ignore the effect of parasitics and therefore drop the index from C L and L L .This expression is given by Here, Φ 0 is the magnetic flux quantum, ϕ j is the flux at node j, j A/B denotes the cell number that is connected to Q A/B , and N is the total number of cells in the ring resonator.In the following, we will assume that the number of cells between the two qubits n AB := j A − j B is even.This matches the setup of the experimental device, where n AB = 6.Since n AB is even, we can place cell number 0 in between the qubits, which now have the indices j A = n AB 2 and j B = − n AB 2 .This choice of reference frame is without loss of generality and can be seen in Fig. 6.The flux-dependent Josephson energy is given by E q J (Φ q ext ) = E q J0 cos , where and E q J1 and E q J2 is the maximum Josephson energy of the two junctions that make up the SQUID loop.The external flux, Φ q ext , gives the flux coupled to the SQUID loop for each qubit q.C S is the shunt capacitance of each qubit and C QM is the qubit-resonator coupling capacitance.
We now introduce the canonical variables of Cooper pair numbers be broken by the qubits, the feedline, imperfections in fabrication and wirebonds (discussion in Appendix D).In our analysis, we include only the symmetry breaking due to the qubits for now, as the complexity increases with more terms.
The values of k∆x/π = 1 and k∆x/π = 0 do not have a corresponding −k value, as they are at the edge and the origin of the first Brillouin zone, respectively.Note that the hypothetical state k∆x/π = −1 differs by 2π from k∆x/π = 1 and is therefore not a unique state.The k∆x/π = 1 state will later be identified as the IR mode, which is not paired with any other −k state and will show no degeneracy with another mode.A similar argument can be made for the k∆x/π = 0 state, which corresponds to the highest frequency mode.In the following we will focus on the paired states, since we are interested in an equation for g(ω).An equation for the g-coupling of the IR mode needs to be calculated separately, which can be done in a similar manner, but will not be shown here.
For the Fourier-transformed canonical variables, we get with the distance between the qubits given by x AB = n AB ∆x, and the left-handed dispersion relation given by . The scaling with ω 2 is unique to left-handed materials and causes a different frequency dependence for the g values than for a system with qubits coupled to a right-handed set of modes.Here φk is expressed in terms of n −k , further highlighting the symmetry between k and −k, which is only broken by the qubit terms.Next, we renormalize the flux to φ = ϕ/Φ 0 , ignore the coupling to the IR mode and the highest frequency mode, and only focus on the paired modes.We transform into even and odd modes using the change of variables and introduce the notation

The total Hamiltonian consists of the parts
All of these terms will be defined in the remainder of this section.Without proof we here assume the interaction between resonator modes H int mode is small compared to the other terms, which is supported by good agreement with numerical results.This approximation breaks down for modes that are very close to the IR mode.For the 24cell ring resonator, no modes are close enough to show a noticeable deviation from theory.In the main text, we also introduce a 240-cell ring resonator.Here, some of the modes are close enough to the IR mode (less than 10 MHz away) and therefore they will show a difference in coupling strength for numerical and analytical results.
For the Hamiltonian of the modes and qubits we get The effective capacitive interaction between the two qubits mediated by the ring resonator capacitances, evaluates to For the interaction between a qubit and an even mode we get Similarly, we can write down the interaction between a qubit and an odd mode as Here, the + branch is for qubit A and − is for B. If we assume both qubits have the same shunt capacitance and the same coupling capacitance, the prefactors are independent of the qubit and given by Approved for Public Release; Distribution Unlimited.PA #: AFRL-2023-2427 We can express the flux and charge operators in terms of creation and annihilation operators a † and a: We .By only considering terms that conserve excitation number based on the rotating wave approximation (RWA), the coupling Hamiltonian reduces to With this Hamiltonian we arrive at an expression for the g-coupling values as a function of k and ω(k) Using the dispersion relation of left-handed materials, we can also express this in terms of the frequency.In the case of a simple dispersion relation without parasitics, this gives us

Note that the proportionality constants α
′ can be either positive or negative, but always have the same sign.Therefore we do not determine the actual sign of the g-couplings to the modes here, only the relative sign between the interactions of each mode to Q A and Q B is given.We refer to this as the parity of the modes or the parity of the g-couplings.Modes with even parity share the same sign between Q A and Q B and modes with odd parity have opposite signs for Q A and Q B .

Appendix B. Device Fabrication
The base layer of the ring resonator device is comprised of 60-nm thick sputtered niobium on a high resistivity (> 10 kΩ-cm) 100-mm silicon wafer.Before sputtering, the silicon wafer undergoes a standard RCA clean process and an oxide etch in a buffered 2% HF solution to remove native oxides.The photolithography patterning is performed by a DUV wafer stepper, and the niobium is etched in an ICP etcher using BCl 3 , Cl 2 , and Ar.All resist is stripped in a TMAH hot strip bath and an oxygen plasma cleaning is performed to remove residual resist residue.We then perform a second 20% buffered HF etch to reduce surface oxides.The device has ground straps along the feedlines labeled R in /R out and F in /F out , as well as the flux-bias lines.Patches of evaporated SiO 2 isolate the straps from the signal traces, and sputtered aluminum provides the electrical ground connection, both patterned via a standard lift-off process.The junctions are pat-100 keV electron-beam lithography of a PMMA/MMA bilayer resist stack, then deposited via a conventional double-angle shadow-evaporation process.
6. Full physical circuit schematic.Right-handed transmission lines labeled Rin/Rout and Fin/Fout in green and yellow, respectively, are used for microwave readout.The 24cell ring-resonator circuit is shown in orange.Capacitively coupled qubits, QA and QB, are shown in red and blue, respectively, and are separated by a distance of nAB = 6.Qubit readout cavities are comprised of quarter-wave coplanar waveguide resonators.

Appendix C. Device parameters
A full circuit diagram of the device is shown in Fig. 6.Parameters for the left-handed ring-resonator are given in Table I.The cell capacitances and inductances were determined via finite element simulation.For all theory calculations and modeling we included the unit-cell stray inductance L R and C R .The circuit diagram for a unit cell including L R and C R is shown in Fig. 7(c).
The resultant band gap at low frequencies due to the left-handed dispersion of the resonator is bounded by the infrared (IR) cutoff frequency, which can be calculated from the cell parameters of the device.Accounting for only the unit cell capacitances and inductances, the IR cutoff for this device is simulated to be 5.6 GHz.Due to the parasitic stray reactances of the ring resonator lumped elements and the grounding wirebonds, the measured IR cutoff is substantially lower at 4.28 GHz.The theoretical infrared cutoff, including the wirebond model and stray reactances, is 4.28 GHz, as shown in Fig. 7.   C L and total center disk capacitance, C C , obtained using Ansys Q3D software.Unit cell inductance and wirebond inductance, L L determined using InductEx software.L W is calculated as 1 nH/mm [25].The stray inductance and capacitance, L R and C R , were found using Sonnet software to identify the self resonances of the lumped-element unit-cell structures.The values in the theory column come from modeling of the device and adjusting parameters to match the measured ring resonator spectrum.

Qubit details are given in
The discrepancy between the theoretical and simulated values for L L are likely due to imperfect grounding and approximations used in estimating the inductive contributions of the wirebonds based on extracted lengths from microscope images.
pacitances are determined via finite element simulation and converting the floating style qubit circuit to the effective transmon circuit using the method described in Ref. [44].The values for E A C and E B C were calculated using Hamiltonian modeling of E C as a function of anharmonicity extracted from independent measurement of the qubit's f 01 and f 02 /2 transition frequencies.We determined E A J0 and E B J0 by fitting spectroscopy data to a characteristic asymmetric transmon curve.This was done because we did not directly measure the qubit maximum transition frequency due to challenges in the chip design.The bare readout cavity frequency for Q A is 6.358 GHz and 6.166 GHz for Q B .From fits, we calculate the maximum transition frequencies for Q A and Q B to be 6.58 GHz and 6.59 GHz, respectively.The proximity of the qubit upper sweetspot to the readout cavity and nearby mode frequencies results in Purcell loss and broad qubit linewidths, making a direct determination of the maximum transition frequencies of the qubits challenging.4 GHz, above which the dispersion becomes righthanded [10].The wirebonds highlighted in blue account for the capacitance of the center disk and an inductance of roughly 1 nH/mm resulting in 1.5 nH per wirebond [25].Ideally wirebonds are short, evenly spaced and densely packed.Due to the structure of our ring resonator, the wirebonds need to avoid the space that the qubits and the feedline occupy, therefore it is impossible to space the bonds evenly, and some areas will have a higher wirebond density than others.This irregular density and positioning of the wirebond connections on the center disk breaks the rotational symmetry of the ring resonator and therefore also breaks the symmetry between clockwise and counterclockwise moving waves and even and odd modes.Because we use numerical methods to simulate the circuit, the position of the wirebond model elements does not exactly follow the experimental placement, as we are focused primarily on achieving an effective model that breaks the symmetry of the circuit with comparable strength.

Appendix E. Cryostat and wiring setup
Measurements are performed on a dilution refrigerator below 15 mK. Figure 8 shows details on the cabling and shielding setup on the dilution refrigerator, as well as the room-temperature electronic configuration for input signals and readout.An infrared-absorbent layer is applied to the inner surfaces of the Cryoperm magnetic shield and the MXC shield [45].A Radiall relay switch on the output line allows switching between measurements via the two feedlines, R in /R out and F in /F out .

Appendix F. g-coupling measurements and fits
To extract g A i and g B i for each mode i, we modeled each mode as an independent harmonic oscillator coupled capacitively to Q A and Q B , respectively.Using standard circuit quantization, we developed a reduced Hamiltonian for the circuit in the basis of transmon charge and resonator excitation number, given by where q is either A or B for Q A and Q B respectively.
The variable E q C is the transmon charging energy and n g is the transmon polarization charge.The operator Îm is the product of the mode identity operators, Îmi is the ith metamaterial mode identity operator, and În is the qubit charge basis identity operator.
A numerical minimization is implemented to fit the relevant eigenvalues of the Hamiltonian to the vacuum Rabi splitting data (main text, Figure 3).We truncate the Hilbert space to include 8 charge states, n, and 4 number states, m, for the resonant modes in the Hamiltonian.For groups of modes with large g-coupling with respect to mode spacing, we include as many as three modes simultaneously in the fit.We do not explicitly include the energy contributions of the spectator qubit in the fits, but we do account for the dispersive χ shift of the ring-resonator modes due to the spectator qubit, which we have fixed at Φ q ext = 0.The free parameters in the fits are the g q i values and E q J0 .We vary E q J0 because the qubit maximum transition frequency is near the qubit readout cavity frequencies and cannot be observed directly, as explained in Appendix C. All other parameters are determined via independent measurements and modeling.
We construct the variance-covariance matrix by taking numerical derivatives of the Hamiltonian with respect to the fit parameters.We compute 95% confidence intervals to obtain the uncertainties in our fit parameters.
Appendix G. Details on J-coupling fits The J-coupling fits were achieved by performing a least squares minimization of the spectroscopy data using a reduced Hamiltonian which reduces the qubits to two-state systems.The summation over modes, i, is reduced to include only the two nearest modes to the crossing of the bare frequencies of Q A and Q B .The dressed qubit frequency for Q A is given by ωA = ω A + i χ A i .To measure J, the flux bias applied to the junction loop of Q B is fixed while we tune the transition frequency of Q A .We obtain ωB from independent measurements in which we detuned Q A .The transition frequency from the ground to the first excited state of Q A is given by ω The amplitude of the exchange term for the two qubits, J, and the maximum Josephson energy for Q A , E A J0 , are the free parameters in the fit.The method for obtaining the magnitudes of the g A i and g B i couplings is outlined in Appendix F. All other parameters are determined via independent measurements.We use the same method described in Appendix F in which we obtain numerical derivatives to compute 95% confidence intervals in our fit parameters.
Table III shows the magnitude of the g-couplings measured for this device.In the last column of Table III the parity of the couplings to the modes is listed.These parities describe the signs of the coupling strength determined by the J-coupling, as explained in the main text.In the basic theoretical model, the mode at 4.927 GHz couples with g = 0 to both qubits and is supposed to have an odd parity.However the experimental J analysis determined that the parity here is even.This can happen for values where theory predicts zero coupling, but symmetry-breaking effects push the g value away from zero.These asymmetric contributions do not need to follow the parity assigned by theory, and in the case of the 4.927 GHz mode, both modes can be pushed to positive g values, even though theory predicts opposite signs.This effect is only relevant for the smallest g values, since the asymmetric contributions are not strong enough for larger values.

Appendix H. Perturbative and non-perturbative J calculations
The theoretical J-coupling strengths are computed using both perturbative and non-perturbative methods.
Here, the parity of the g-coupling values between the two qubits is determined by comparing measured versus perturbative theoretical calculations of the exchange coupling for Q A and Q B , as described in the main text.
The perturbative Schrieffer Wolff method is also highlighted in the main text.Here, it is notable that g A i g B i scales with ω 3 , while ∆ i scales with ω if the resonator frequency is much larger than the qubit frequency.This means that high-frequency modes still contribute some effective coupling strength, even if they are far away from the qubits.Because of this, we sum over all the modes for which we have experimental g values, as well as those that are above the upper sweet spot of the qubits.For those we tune the qubits to the upper sweet spot and observe a change in the dressed resonator frequency.Based on this we can extract a g value for each mode, which is used to model the interaction with these high-frequency modes.
We want to confirm our theoretical predictions of the effective J-coupling also with a non-perturbative method, since this type of 'long-range interaction' is not ideal for perturbative means.Therefore, we use the Least Action method [27][28][29] to block-diagonalize our full-system Hamiltonian.This Hamiltonian includes single and double excitations for both qubits and resonator modes, as well as all possible combinations of one qubit excitation paired with one resonator excitation.We also include two excitations of ring resonator modes which are not the same mode.Since there are many possible combinations for this, we restrict ourselves to those that couple strongly and are relatively close to the frequency domain we are interested in.The resulting Hamiltonian is blockdiagonalized using Least Action and the resulting J 00 coupling strength can be read off.This gives two J values for each anticrossing between the two qubits.One is slightly before the discontinuity of the dressed frequencies and the other one is slightly after it.These two values are usually very close to each other and we average them to get one value for the anticrossing.This is then compared to Schrieffer Wolff and experimental results.

Appendix I. ZZ measurements and fits
When performing the ZZ measurements, we fix the bare frequency of Q B and vary the frequency of Q A , but the dressed frequency of Q B varies depending on the frequency of Q A .To measure the frequency of Q B with high precision, a Ramsey fringe measurement and fit is performed for Q B for each frequency tuning of Q A .The range in fB is 0.7 MHz (1 MHz) for the measurements in which we quote fB to be 5.04 GHz (5.81 GHz).These small variations in fB are due to the effective inter-qubit coupling via the ring resonator modes.
The pulse sequences for extracting the ZZ interaction strength are shown in Fig. 9(a,b).The Ramsey fringe fits were performed by first fitting spectroscopy data of the one-state probability for Q A as a function of idle time for each detuning frequency measured.This corresponds to taking a single horizontal slice of the data shown in Fig. 9(a,b), then applying a damped, oscillatory fit function.We remove fits with error greater than 0.4 MHz in the estimation of f Ramsey .We then fit a line through the f Ramsey data with a slope of one to find the zero point for the Ramsey oscillation frequencies, as shown in Fig. 9(cf).We perform these fits for both pulse sequences, then calculate the difference in transition frequency for Q A to extract ζ/2π.The error for extracting the f Ramsey values is computed from 95% confidence intervals for the Ramsey fits.The error bars for the ζ/2π values as a function of fA in Fig. 9 (main text) are obtained by finding the 95% confidence intervals for a linear fit with a slope of one to the f Ramsey data.FIG. 9. Measurement sequence for obtaining ZZ interaction strength.(a) A Ramsey pulse sequence is performed on QA at variable drive frequencies while QB is idle.(b) A second Ramsey fringe measurement is performed on QA after a π-pulse is applied to QB. (c) Horizontal slices of Ramsey fringe measurements are fit to obtain the Ramsey oscillation frequency, fRamsey, at each drive frequency.The ZZ interaction, denoted ζ/2π, is the change in the dressed transition frequency for QA, fA, depending on whether QB is in the ground or excited state.(d-f) Three example plots of Ramsey frequency as a function of QA drive frequency ( fA drive) for different QA bias points with linear fit lines for ZZ measurements taken when QB is at 5.04 GHz.

Appendix J. Theoretical ZZ calculations
The ZZ interaction is obtained by fully diagonalizing the Hamiltonian that was also used for calculating the J-coupling strengths.Alternatively, the already blockdiagonalized Hamiltonian can be used to only diagonalize the qubit subspace, thus saving computational cost with similar results.The ZZ interaction is then calculated using ζ = E 00 + E 11 − E 10 − E 01 .The biggest impact on the ZZ strength comes from the E 11 level interacting with non-computational states such as E 20 .The high resonator mode density around 4−5 GHz causes the 10 GHz area to be highly populated with double excitations of qubits and resonator modes.For the frequency domain shown in Fig. 10, the most important states are |20; 0000⟩ and |01; 0001⟩, which both couple strongly to |11; 0000⟩, since only one exchange of excitation is required for a transition.The anticrossing with |20; 0000⟩ is especially Since there are multiple other modes close to this region, the exact point of the frequency discontinuity strongly depends on various parameters, making accurate predictions challenging.
interesting here, since it causes a large shift and discontinuity of E 11 , which is also observable in the experimental data in Fig. 10(a).However, the experimental drop-off is positioned at 5.26 GHz, while the theoretical one is at 5.253 GHz.This slight disagreement is probably caused by small inaccuracies in multiple parameters.Most noticeably the anharmonicity of Q A determines the energy of the |20; 0000⟩ state; an inaccuracy here will result in a frequency shift of the drop-off position.There are also multiple other modes nearby, and a change of their frequencies or coupling strengths will impact the position of the discontinuity.

FIG. 1 .
FIG. 1.(a) Mode structure for two pairs of degenerate ring resonator modes.The ring resonator consists of 24 unit cells, shown with gray squares, with each cell containing a capacitor shunted to ground by two inductors.One pair of even (odd) modes with k∆x = 10 2π N is shown in solid (dashed) orange lines, labeled m

O
); a second pair with k∆x = 5 2π ) depends on |k|, waves moving in opposite directions have the same frequency.In a perfect ring with identical cells and no external coupling, equal combinations of clockwise |k⟩ and counterclockwise | − k⟩ waves define the standing wave of frequency ω(k), with even and odd parity |E⟩ = (|k⟩ + | − k⟩)/ √ 2 and |O⟩ = (|k⟩ − | − k⟩)/

FIG. 2 .
FIG. 2. (a) Chip layout for metamaterial ring-resonator device.(b, c) Optical micrographs of device with falsecolor highlighting.(d) SEM image of Josephson junction in qubit.Device has two flux-tunable transmon qubits, QA and QB, coupled to a ring resonator at the positions indicated.Rin/Rout connections to upper feedline used to probe ring resonator modes.Fin/Fout connections used for measuring readout resonators coupled to each qubit.(e) Measured ring resonator mode frequencies are shown with gray dashed lines.Theoretical mode frequencies are shown with green dasheddotted lines when stray inductance due to wirebonds is not included and any degeneracy lifting is due to the qubits or feedline.Solid blue lines show large degeneracy lifting effect of wirebonds on the ring resonator mode frequencies.

FIG. 3 .
FIG. 3. Vacuum Rabi splittings measured via Rin/Rout, between three ring resonator modes and (a) QA, (b) QB.Dashed lines in red (blue) show fits to splitting data for QA (QB) obtained from reduced Hamiltonian including one qubit and the three modes (see Appendix F).Horizontal dotted lines show dressed mode frequencies.(c-f) Splittings and reduced Hamiltonian fits for two single modes for each qubit.For all measurements, the spectator qubit is fixed at its upper flux insensitive sweet spot (Φq/Φ0 = 0).As shown in (c), QA has effectively zero coupling to the ring resonator mode at 4.553 GHz, while (d) shows QB has a coupling strength of 26.6 MHz.(g) Extracted and theoretical |g q i |/2π values for each ring-resonator mode i. Error bars on experimental data computed from 95% confidence intervals in g-coupling fits (see also Appendix F).Vertical dashed lines show bare ring resonator mode frequencies.

FIG. 4 .
FIG. 4. (a) and (b) show qubit spectroscopy for QA as a function of flux when QB is at a fixed frequency, denoted by a blue dotted line.Red dashed lines show fits to data obtained via minimization of reduced Hamiltonian from which we determine J/2π (see Appendix G).(c) Experimental and theoretical J-coupling values as a function of frequency.Diamonds show experimental values for J/2π, with error bars computed from 95% confidence intervals in fits (see also Appendix G).Gray lines and purple crosses show theoretical J/2π values calculated using Schrieffer-Wolff and Least Action, respectively.Dotted vertical lines show bare ring-resonator mode frequencies.

FIG. 5 .
FIG. 5. ζ/2π as a function of the dressed frequency of QA, fA, when the dressed frequency of QB, fB, is (a) 5.04 GHz, (b) 5.81 GHz (see Appendix I).Measured and theoretical values are shown as diamonds and gray lines, respectively.Vertical dashed lines in blue and gray show the location of QB and the ring resonator modes in frequency space, respectively.Insets show Ramsey oscillation frequency fit data as a function of QA drive frequency taken from Ramsey fringe measurements.Red data points and fit lines generated by performing a simple Ramsey pulse sequence (shown in red box) on QA at multiple detunings and extracting oscillation frequency.Data resulting from pulse sequence in which a π-pulse is applied to QB, followed by a Ramsey on QA (shown in purple box).Error bars computed from 95% confidence intervals for both Ramsey oscillation fits (see also Appendix I).Intersection of fit lines where Ramsey oscillation frequency vanishes indicate QA frequency, from which we compute ζ/2π.

n k = 1 √ 2 N , −1 + 4 N
N j e ik∆xj n j , where k∆x/π ∈ [−1 + , ..., 1 − 2 N , 1] form a Brillouin zone.We can immediately see that for most k values there is a corresponding −k mode, where positive k values correspond to modes traveling around the ring resonator clockwise and negative k values correspond to counterclockwise moving modes.This symmetry is essential to the device, but can Approved for Public Release; Distribution Unlimited.PA #: AFRL-2023-2427 can identify the zero point fluctuation of the charge as n zpf = ℏ 2eLω , which scales with 1/ √ ω.Similarly, the zero point fluctuations of the charge for each qubit scales with √ ω A/B .We redefine the proportionality factors α A/B E/O to include the zero point fluctuations, except for their frequency dependence, and label them as α

FIG. 7 .
FIG. 7. (a) Image of wirebonds used for grounding the device, as well as attachments to signal traces.(b) GDS of device layout including approximate wirebond positions shown with white lines.(c) Theoretical wirebond model used for modeling mode resonances and coupling strengths.

FIG. 8 .
FIG. 8. Wiring schematic for room temperature and cryogenic cabling down to sample box.Standard heterodyne setup is used for measurement.All filters are low-pass unless otherwise indicated.

FIG. 10 .
FIG. 10.(a) Experimental and theoretical ZZ interaction ζ/2π as a function of QA frequency fA shown with diamonds and grey lines, respectively.(b) Multi-photon energies as a function of fA.The light grey regions in (a) and (b) denote a region of uncertainty in the theoretical ZZ calculation.The discontinuity of the ZZ interaction around 5.25 GHz is caused by an anticrossing of the |11⟩ and |20⟩ states sible in (b).Since there are multiple other modes close to this region, the exact point of the frequency discontinuity strongly depends on various parameters, making accurate predictions challenging.

Table II .
The effective qubit ca-

TABLE I :
Ring resonator cell parameters obtained via finite element simulation.Unit cell capacitance,

TABLE II :
Parameters for two transmon qubits, Q A and Q B .

TABLE III :
Magnitude of qubit-ring resonator g-coupling parameters and parity for Q A and Q