Introduction

The recent advances in the field of quantum information processing has led to a rising demand to explore new systems beyond cavity quantum electrodynamics (QED). One promising candidate is waveguide QED, where quantum systems interact coherently with the mode continuum of a waveguide instead of a cavity. After the pioneering works with single qubits, including the demonstration of resonance flourescence1, quantum correlations of light and single photon routers2, attention shifted to the realization of multiple qubits coupled to a common waveguide. It was derived3 and experimentally verified4 that multiple qubits obtain an infinite range photon mediated effective interaction, which can be tuned with the inter-qubit distance d. Furthermore, the shared collective excitations are of polaritonic nature with lifetimes ranging from extremely sub- to superradiant relative to the radiative lifetime of the individual qubits5,6. The strong intrinsic nonlinearity of the qubits was recently shown to give rise to partially localized polaritons7, topological edge states8,9, and quantum correlations in the scattered light of the array10. The collective quantum properties are exploited in the field of quantum metamaterials11,12. Here, the quantum coherence of the constituting qubits is used to engineer a global optical response which depends on their quantum state13,14,15,16. With respect to quantum information processing, multi-qubit waveguide QED systems could be harnessed in numerous applications such as on demand, highly efficient creation of multi-photon and entangled states17,18,19, storage devices for microwave pulses20, atomic mirrors21, number-resolved photon detection22, slow and even stopped light23. Experimentally, waveguide QED systems have been realized on several platforms including atoms24, quantum dots coupled to nanophotonic waveguides25 and defect centers in diamonds26. Even though superconducting qubits feature advantages such as frequency control, high coherence, near perfect extinction and absence of position and number disorder, superconducting multi-qubit waveguide QED systems have been studied only recently to some extend27,28.

Here, we investigate the mode spectrum of a metamaterial formed by a densely spaced array of eight superconducting transmon qubits coupled to a coplanar waveguide. By employing dedicated flux-bias lines for each qubit, we establish control over their transition frequencies. Thus we are able to alter the number of resonant qubits N at will, allowing us to observe super- and subradiant modes as well as the gradual formation of a bandgap. Our spectroscopic measurements show, that through this control the global optical susceptibility of the metamaterial can be tuned. A demonstration for the collective Autler-Townes splitting (ATS) of eight qubits is presented, which marks an important step toward the implementation of quantum memories.

Results and discussion

Circuit design and properties

The sample under investigation is depicted in Fig. 1(a). The spacing d between adjacent qubits is 400 μm, which is smaller than the corresponding wavelength λ at all accessible frequencies (\(\varphi =\frac{2\pi }{\lambda }d=0.05\!-\!0.16\)). The dense spacing is chosen in order to increase the width of the expected bandgap of Δω = Γ10/φ Γ10 and to fulfil the metamaterial limit of subwavelength dimensions29. Here Γ10 is the radiative decay rate of the individual qubits into the waveguide. The qubits are overcoupled, ensuring a multi-mode Purcell-limited rate Γ10/2π ≈ 6.4 MHz for high extinction and, simultaneously, better subradiant state visibility. All qubits are individually frequency controllable between 3 and 8 GHz by changing the critical currents of the qubit SQUIDs with local flux-bias lines. We compensate unwanted magnetic crosstalk between the flux-bias lines and neighbouring qubits by extracting and diagonalising the full mutual inductance matrix M (see Methods). This allows us to counteract the parasitic crosstalk flux by sending appropriate currents to all qubits, which are not actively tuned. With this calibration scheme we achieve precise control over the individual qubit frequencies. We estimate the residual crosstalk to be smaller than 10−3. The effective Hamiltonian of this system, after formally tracing out photonic degrees of freedom and applying the Markov approximation, is described by6,

$$\frac{{H}_{{\rm{eff}}}}{\hslash }=\mathop{\sum }\limits_{j}^{8}({\omega }_{j}{b}_{j}^{\dagger }{b}_{j}+\frac{{\chi }_{j}}{2}{b}_{j}^{\dagger }{b}_{j}^{\dagger }{b}_{j}{b}_{j})+{\rm{i}}\frac{{{{\Gamma }}}_{10}}{2}\mathop{\sum }\limits_{k\ne j}^{8}{b}_{k}^{\dagger }{b}_{j}{{\rm{e}}}^{-{\rm{i}}\frac{\omega d}{c}| k-j| },$$
(1)

with the bosonic creation operator \({b}_{j}^{\dagger }\left|0\right\rangle =\left|{{\rm{e}}}_{j}\right\rangle\), exciting the j-th qubit at frequency ωj; where we assume a e+iωt time dependence of the excitations. Here, χj/2π is the qubit anharmonicity for which we find spectroscopically weakly varying values around −275 MHz. The last term of Heff describes the effective qubit–qubit coupling. Due to the specific choice of small d its imaginary part dominates over the real part, leading to a suppressed exchange-type interaction between the qubits. The expected eigenfrequencies ωξ of Heff in the single excitation limit and χ → 0 are shown in Fig. 2(a).

Fig. 1: Sample and transmission.
figure 1

a Optical micrograph of the metamaterial. It is composed of eight superconducting transmon qubits capacitively coupled to a coplanar waveguide. Local flux-bias lines provide individual qubit frequency control in the range 3–8 GHz. b Transmission S21 for different numbers N of resonant qubits at ωr/2π = 7.898 GHz and low drive powers. With increasing N, the emergence of subradiant states (visible as peaks in transmission) can be observed. Black dotted lines are fits to the expected transmission using a transfer matrix calculation. Black arrows mark calculated frequencies of the two brightest subradiant states for N = 8.

Fig. 2: Bandgap and polariton relaxation rates.
figure 2

a Dependence of absolute transmission S212 on the number of resonant qubits N. Crosses mark calculated eigenfrequencies ωξ of Heff and their corresponding radiative decay rates Γξ10 (color coded). With increasing N a bandgap of strongly supressed transmission opens up. b Measured radiative decay rates Γξ of the brightest (blue crosses) and second brightest (red crosses) subradiant states. Orange and blue solid triangles are the corresponding calculated rates, the green line is the analytical result Γξ = 8N3φ2/π4. For the brightest subradiant state, a fit to a power-law with exponent b = 2.93 ± 0.13 (black dashed line) confirms the scaling of ΓξN3.

Bandstructure and collective metamaterial excitations

First, we characterize the mode spectrum of the metamaterial in dependence on the number of resonant qubits N by measuring the transmission coefficient S21(ω) while the qubits are consecutively tuned to a common resonance frequency at ωr/2π = 7.898 GHz, compare Fig. 1 (b) and Fig. 2 (a). The incident photon power Pinc is kept below the single photon level (Pinc\({\hbar}\)ωΓ10) to avoid saturation of the qubits. For a single qubit, the well known resonance-fluorescence was observed as a single dip in transmission1. By fitting the complex transmission data to the expected transmission function S21(ω) (see Methods) the individual coherence properties of all eight qubits at ωr can be extracted (see Supplementary Note 1). We find the average radiative rates Γ10/2π = 6.4 MHz and the intrinsic non-radiative rates Γnr/2π = 240 – 560 kHz. For N ≥ 2 resonant qubits the system obtains multiple eigenmodes and the super- and subradiant polariton modes start to emerge. The superradiant mode is manifested as a wide transmission dip above ωr, the subradiant modes can be identified as transmission peaks below ωr. We note that the peak shape is created by Fano-interference of the sub- and the superradiant modes (see Supplementary Note 5), which also causes the calculated eigenfrequencies of Heff to not exactly coincide with the maximum of the peaks in Fig. 1 (b) and Fig. 2(a). For N ≥ 6 qubits a second darker subradiant mode is visible between ωr and the brightest subradiant mode. The limiting factor for the observation of the subradiant states is the intrinsic qubit coherence as given by Γnr, setting an upper threshold for the maximum observable lifetime. Darker subradiant states with Γξ < Γnr decay in the qubits into dielectric channels or dephase due to flux noise in the SQUIDs before they are remitted into the waveguide. The observed transmission coefficient S21 is in good agreement with the calculated transmission based on a transfer matrix approach30. The asymmetric lineshape of the resonances is a parasitic effect caused by interference of the signal with low-Q standing waves in the cryostat31. We account for this effect in the transfer matrix calculation by adding semi-reflective inductances in front and after the qubit array. As shown in Fig. 2(a), a frequency region of strongly suppressed transmission with S212 <−25 dB is opening up above ωr with increasing N. This effect is associated with the emergence of a polaritonic bandgap, where the effective refractive index becomes purely imaginary, as expected for any kind of resonant periodic structures29,32. For N = 8 qubits we extract a bandgap bandwidth of Δω ≈ 1.9Γ10. Compared to the expected bandgap width of Δω = Γ10/φ ≈ 6.3Γ10 of the structure for N →  this places our system size in the transitioning regime between a single atom and a fully extended metamaterial with a continuous mode spectrum. The radiative decay rates Γξ of the observed subradiant states shown in Fig. 2(b) are extracted by fitting Lorentzians to the corresponding modes in the reflection data (not shown). It can be generally shown, that S212 indeed obtains a Lorentzian shape in the vicinity of each ωξ30,33,34. From a fit to a power law Nb with exponent b = 2.93 ± 0.13 we find that the rate of the brightest subradiant states scales with ΓξN3, which we also find analytically for densely packed qubit structures with φ 1 from Heff (see Supplementary Note 4). The found scaling is the complementary asymptotic of the theoretically predicted ΓξN−3 law for the darkest subradiant modes in references5,6,32. Small deviations between calculated and measured values of Γξ are caused by imperfect qubit tuning and distortions of the observed reflection coefficient due to the microwave background.

The control over the mode spectrum can be further elaborated with the off-resonant situation, where one qubit has a finite detuning Δ from the common resonance frequency ωr of the residual qubits by sweeping it through the common resonance as shown in Fig. 3(a). Here, qubit 8 is tuned through the collective resonance of qubits 1–7. For large detunings Δ Γ10 the modes of the ensemble and qubit 8 are not hybridized (not shown). For smaller detunings an additional partially hybridized subradiant mode appears, which becomes for Δ ≈ 0 the brightest subradiant state of the fully hybridized 8-qubit system. The result is in good agreement with the transfer matrix calculation and direct diagonalization of Heff as shown in Fig. 3(b). The level repulsion between the eigenstates is caused by the residual exchange-type interaction between the qubits due to the finite inter-qubit distance d. In Fig. 3 there are several blind spots where the subradiant states turn completely dark, occurring when the frequency of the detuned qubit matches the frequency of a dark mode. This is explained by the Fano-like interferences35 between the detuned qubit resonance and the modes of the resonant qubits, which are analyzed in more detail in the Supplementary Note 5.

Fig. 3: Off-resonant qubits.
figure 3

a Measured absolute reflection S22 for N = 7 resonant qubits at ωr/2π = 7.897 GHz (solid black line) and qubit 8 being tuned through the collective resonance (black dashed line). b Calculated reflection with transfer matrix method and relevant eigenfrequencies of Heff (dash dotted lines) show good agreement with experimental result in (a). The vertical line indicates the resonance with ω8 = ωr, corresponding to the situation observed in Fig. 1.

Collective Autler-Townes splitting

The intrinsic quantum nonlinearity of the system due to the anharmonic nature of the qubits can be probed by increasing the microwave power beyond the single photon regime (Pinc > \({\hbar}\)ωΓ10). At higher power the 0 → 1 transition of a single qubit will start to saturate and transmission at resonance will increase back to unity1, see Fig. 4(b). We could experimentally verify the prediction of the saturation of an ensemble of resonant qubits at higher drive rates of refs. 1,3, which can be observed in Fig. 4(a). We point out that all spectroscopic features of the metamaterial, such as super- and subradiant modes as well as the bandgap, are saturable. The power P50% needed to saturate the transmission at ωr to 50% (S212 = 0.5) grows with an approximate scaling \(\propto \mathrm{ln}\,(N)\). The observed quantum nonlinear behavior establishes a border to previously studied metamaterials consisting of harmonic resonators rather than qubits36. Furthermore the anharmonic level structure of the transmon can be used to electromagnetically open a transparency window around the frequency of the 0 → 1 transition, by employing the ATS37. As indicated in Fig. 4(c), a coherent control tone with Rabi strength Ωc and frequency ωc drives the 1 → 2 transition, while a weak microwave tone with Ωp Ωc is probing the transmission at frequencies ωp around the 0 → 1 transition. The control tone is dressing the first qubit level, creating two hybridized levels, which are separated proportionally to its amplitude, and is therefore creating a transparency window with respect to the probe. To the best of our knowledge this effect was so far only demonstrated for a single qubit in superconducting waveguide QED, leading to applications like single photon routers2. Here, we observe the collective ATS for up to N = 8 resonant qubits (Fig. 4(d)). Analogously to the single qubit case, the observed level splitting is proportional to Ωc. At control tone powers >−110 dBm the bandgap is rendered fully transparent and transmission close to unity is restored. In the experiment we find slightly smaller than unity values of S21 ≈ 0.75 due to the interference of the signal with the microwave background. In the bandstructure picture the collective ATS can be understood as the dressed states of the individual qubits giving rise to two independent Bloch bands38. Therefore, the two branches of the collective ATS with supressed transmission are collectively broadened bandgaps and have a larger linewidth than the dressed states of the single qubit ATS. As pointed out in reference38 the collective ATS demonstrates active control over the bandstructure of the metamaterial via the parameters ωc and Ωc. The observed splittings are in agreement with a transfer matrix calculation and a full master equation simulation (see Supplementary Note 6). Minor deviations are caused by imperfect qubit tuning and differing qubit anharmonicities of about σχ ≈ 10 MHz. We argue here that the single photon router concept with a single qubit can significantly be improved with multiple qubits, which form a much wider stop-band with higher saturation power, thus permitting to route also multiple photons.

Fig. 4: Power saturation and Autler-Townes splitting.
figure 4

a Absolute power transmission S212 for N = 8 resonant qubits shows saturation with increasing power due to the anharmonic level structure of the transmons. b Comparison between the saturation of S212 at ωr for a single qubit and N = 8 qubits. c Schematic illustration of the Autler-Townes effect for a ladder-type three-level system. A control tone is driving the 1 → 2 transition with Rabi strength Ωc and frequency ωc, while a weak tone is probing the transmission of the 0 → 1 transition. d Experimental demonstration of the collective Autler-Townes splitting for N = 8. The red dashed line marks the fitted level separation to Ωc.

In conclusion, we demonstrated a fully controllable quantum metamaterial consisting of eight densely packed transmon qubits coupled to a waveguide. Such intermediate system size combined with individual qubit control allowed us to explore the transition from a single mode regime to a continuous band spectrum. By tuning the qubits consecutively to resonance, we observed the emergence of a polaritonic bandgap, and confirmed the scaling of the brightest dark mode decay rate with the qubit number. Active control over the band structure of the ensemble was demonstraded by inducing a transparency window in the bandgap region, using the Autler-Townes effect. Our work promotes further research with higher qubit numbers to realize a long-living quantum memory.

Methods

Fabrication

The sample is fabricated with two consecutive lithography steps from thermally evaporated aluminum in a Plassys MEB550s shadow evaporator on a 500 μm sapphire substrate. In a first step solely the qubits are patterned with 50 keV electron beam lithography. The Josephson junctions are patterned with a bridge-free fabrication technique39 using a PMMA/PMMA-MAA double resist stack. Before double-angle evaporation, the developed resist stack is cleaned for 6 min with an oxygen-plasma to remove resist residuals in the junction area to reduce the impact of junction aging39. In a second optical lithography step we pattern the CPW-waveguide and the ground-plane in a liftoff-process on S1805 photo resist. The SQUIDs are formed by two Josephson junctions, enclosing an area of 560 μm2. By design, the junction areas are 0.12 and 0.17 μm2 with a designed asymmetry of 17%.

Calibration of magnetic crosstalk

We extract the full 8 × 8 mutual inductance matrix M between the qubits and the bias coils in 28 consecutive measurements. For that, the transmission through the chip is observed at a fixed frequency while tuning the currents of two qubits Ix and Iy such that they get tuned through the observation frequency. Fitting the slope of the observed qubit traces gives access to Mxy/Mxx and Myx/Myy. When M is known, compensation currents which are send to all qubits which are not actively tuned, can be calculated to compensate unwanted crosstalk. We estimate the residual crosstalk to be smaller than 0.1%. Further information is provided in the Supplementary Note 2.

Normalization of spectroscopic data

We normalize the transmission data by dividing the raw data \({S}_{21}^{\,\text{raw}\,}\) by the transmission data at high powers \({S}_{21}^{\text{sat}}{:}{S}_{21}^{\,\text{calib}\,}(\omega )={S}_{21}^{\,\text{raw}\,}(\omega )/(a{S}_{21}^{\,\text{sat}\,}(\omega ))\), where a ≈ 1 is a constant factor accounting for weak fluctuations of the amplifier gain. Reflection data is normalized by dividing the raw data \({S}_{22}^{\,\text{raw}\,}\) with its maximum value at the qubit resonance frequency fr: \({S}_{22}^{\,\text{calib}}={S}_{22}^{\text{raw}}/{S}_{22}^{\text{raw}\,}({f}_{\text{r}})\). This approximation is justified for the sample under investigation since the extinction of single and multiple qubits is very close to 1. A more rigorous normalization of the data based on energy conservation S212 + S112 ≈ 1 is not applicable here, due to differing signal paths for reflection and transmission measurements, as pointed out in the main text. In addition, in any experimental system the insertion losses due to impedance mismatches of the feedline can not be avoided. Since they are in general not symmetric on both sides of the sample, S11 ≠ S22 and therefore S212 + S222 ≠ 1.

Characterization of individual qubits

The amplitude transmission coefficient of a driven two-level system side-coupled to a waveguide S21 is given by2:

$${S}_{21}=1-\frac{{{{\Gamma }}}_{10}}{2{\gamma }_{10}}\frac{1-{\rm{i}}\frac{\omega -{\omega }_{\text{r}}}{{\gamma }_{10}}}{1+{\left(\frac{\omega -{\omega }_{\text{r}}}{{\gamma }_{10}}\right)}^{2}+\frac{{{{\Omega }}}_{{\rm{p}}}^{2}}{({{{\Gamma }}}_{10}+{{{\Gamma }}}_{l}){\gamma }_{10}}}$$
(2)

The decoherence rate γ10 = Γ10/2 + Γnr is the sum of radiative decay Γ10 and non-radiative decay rates Γnr = ΓΦ + Γl/2. Here ΓΦ accounts for pure dephasing of the qubit and Γl for all non radiative relaxation channels. We define the extinction coefficient as \(1-{(1-{{{\Gamma }}}_{10}/2{\gamma }_{10})}^{2}\), measuring the suppression of power-transmission at very low drive powers. A circle fitting procedure40 is used to fit equation (2) to the measured complex transmitted signal S21 in the limit of weak driving \({{{\Omega }}}_{{\rm{p}}}^{2}\ll ({{{\Gamma }}}_{10}+{{{\Gamma }}}_{l}){\gamma }_{10}\). The decoherence rates of the individual qubits and further details on the fitting procedure are provided in the Supplementary Note 1.