Observation of a dissipative phase transition in a one-dimensional circuit QED lattice

Condensed matter physics has been driven forward by significant experimental and theoretical progress in the study and understanding of equilibrium phase transitions based on symmetry and topology. However, nonequilibrium phase transitions have remained a challenge, in part due to their complexity in theoretical descriptions and the additional experimental difficulties in systematically controlling systems out of equilibrium. Here, we study a one-dimensional chain of 72 microwave cavities, each coupled to a superconducting qubit, and coherently drive the system into a nonequilibrium steady state. We find experimental evidence for a dissipative phase transition in the system in which the steady state changes dramatically as the mean photon number is increased. Near the boundary between the two observed phases, the system demonstrates bistability, with characteristic switching times as long as 60 ms -- far longer than any of the intrinsic rates known for the system. This experiment demonstrates the power of circuit QED systems for studying nonequilibrium condensed matter physics and paves the way for future experiments exploring nonequilbrium physics with many-body quantum optics.

Condensed matter physics has been driven forward by significant experimental and theoretical progress in the study and understanding of equilibrium phase transitions based on symmetry and topology. However, nonequilibrium phase transitions have remained a challenge, in part due to their complexity in theoretical descriptions and the additional experimental difficulties in systematically controlling systems out of equilibrium. Here, we study a one-dimensional chain of 72 microwave cavities, each coupled to a superconducting qubit, and coherently drive the system into a nonequilibrium steady state. We find experimental evidence for a dissipative phase transition in the system in which the steady state changes dramatically as the mean photon number is increased. Near the boundary between the two observed phases, the system demonstrates bistability, with characteristic switching times as long as 60 ms -far longer than any of the intrinsic rates known for the system. This experiment demonstrates the power of circuit QED systems for studying nonequilibrium condensed matter physics and paves the way for future experiments exploring nonequilbrium physics with many-body quantum optics.
Over the past decades, there has been remarkable progress in studying both real and synthetic quantum materials. Advances in nanoscale fabrication and cryogenics have allowed for exquisite control of electronic systems -unlocking strongly correlated electronic states and topological materials 1 . Simultaneously, the ability to model desired Hamiltonians with ultra-cold Fermi and Bose gases has allowed unprecedented access to synthetic material properties 2 . As a whole, much of the development of condensed matter physics has focused on the study of (quasi-)equilibrium physics, which is more accessible both experimentally and theoretically. However, the constant presence of dissipation, noise, and decoherence belie the fact that, ultimately, the world is nonequilibrium.
A phase transition indicates a sometimes sudden change in the physical properties of a system as a function of some external system parameter. Thermal phase transitions are well-understood in the context of statistical mechanics and occur when the free energy becomes nonanalytic. At zero temperature, the role of quantum fluctuations gives rise to a new set of quantum phase transitions which involve a sudden change in the ground state of a Hamiltonian H; a phase transition occurs when the gap between the first excited state and the ground state closes. These concepts need to be extended to consider nonequilibrium steady states, as the system is no longer in its ground state but rather in a state that balances drive and dissipation. In a dissipative phase transition, the steady state abruptly changes as a system parameter is varied 26 . Whenever the system is describable in terms of a Lindblad master equation 3 ,ρ = Lρ, then such a transition is signalled by the closing of the lowest excitation gap in the spectrum of the Liouvillian superoperator L.
In recent years, interacting photons have emerged as an excellent candidate for studying nonequilibrium condensed matter physics due to the lack of particle number conservation 4 . In cavity quantum electrodynamics, strong coupling between atoms and a cavity can mediate effective photon-photon interactions [5][6][7] . Arrays of coupled microwave 8 or optical 9 cavities can be fabricated by conventional lithographic techniques, and the competition between on-site interactions and hopping between neighboring cavities can give rise to quantum phase transitions of light [10][11][12] . A wide range of many-body effects have been predicted in these systems, including a Mott insulator-superfluid phase transition 10-12 and fractional quantum Hall-like states of light [13][14][15] . Experiments on small systems have demonstrated low-disorder lattices 8 , a dynamical quantum phase transition in a cavity dimer 16 , and chiral ground state currents in a cavity trimer 17 . Circuit QED lattices are inherently open systems, with dissipation an ever-present force that leads both to qubit relaxation as well as the inevitable loss of photons from microwave cavities. While dissipation presents an obstacle for quantum information processing, it is of fundamental interest in the study of nonequilibrium phase transitions. Just as excitations inevitably leak from the system, it is easy to add photons back and to drive into a steady state, making these systems particularly amenable to the study of dissipative phase transitions.
In this paper we present experimental evidence for a dissipative phase transition in a circuit QED lattice. We observe that at drive frequencies between the low-power resonance frequencies of the system, there exists a region of hysteresis and bistability where the steady state of the system switches stochastically between two states ρ 1 and ρ 2 . By determining the corresponding switching rates, we can obtain the so-called asymptotic decay rate which characterizes the closing of the spectral gap of L. At the transition between the two states, the characteristic switching times become exceptionally long, a key characteristic of a dissipative phase transition. A similar observation has recently been made in a single cavity system with multiple qubits 18 .
FIG. 1. 72-site circuit QED lattice. a, Coplanar waveguide resonators, each with a bare cavity frequency of 7.5 GHz, are capacitively coupled to form a linear chain on a 2.5 × 2.5 cm 2 chip. Each resonator is coupled to its neighboring resonators to yield a hopping rate t/2π ≈ 144 MHz and has an average photon loss rate of κ/2π ≈ 1.6 MHz. At three intermediate chain sites, three-way coupling capacitors provide ports for secondary input and output lines (arrows on sides). b-c, A transmon qubit is capacitively coupled to the center pin near the edge of each resonator in the lattice, ensuring coupling to the fundamental mode of each resonator with strength g. The coupled resonator-qubit system forms the fundamental unit cell of the lattice. d, The circuit can be modeled as a linear chain of coupled oscillators, each dispersively coupled to a weakly anharmonic multi-level system.
Our device, shown in Fig. 1a, consists of a linear chain of 72 lattice sites. Each site comprises a coplanarwaveguide resonator with fundamental-mode frequency ω/2π = 7.5 GHz, coupled to a transmon qubit 19 (Fig. 1b) placed at one of the resonator's voltage antinodes. Resonators are capacitively coupled to neighboring resonators (Fig. 1c), so that photons can hop between nearest-neighbor sites. Variations in transmon qubit frequencies in fabrication are a likely source of uncontrolled disorder that is difficult to compensate for in our lattice. We therefore use an asymmetric SQUID-loop geometry allowing each qubit to be tuned over a finite frequency range via an applied magnetic flux. Because individually tuning 72 qubits is currently infeasible in our system, we instead employ a global magnetic field to simultaneously tune all qubit frequencies. Because each qubit is intentionally fabricated with a SQUID-loop of random area, this randomizes the frequency of all qubits within a band of frequencies near 8.5 GHz. In this way, we can ensure that features of interest are universal to the system rather than artifacts of a particular instance of disorder (see Supplementary Information I).
To experimentally study the nonequilibrium behav-ior of the device, we monitor the homodyne transmission across the lattice while varying the drive frequency and scanning the drive power over more than five orders of magnitude ( Fig. 2a). At low drive powers, we find the expected discrete transmission peaks associated with the interaction-shifted eigenmode frequencies of the resonator lattice. As we vary the mean photon number in the system by increasing the strength of the drive, we observe that a sudden change in system behavior occurs: transmission peaks split and then, at around −10 dB of drive power, abruptly give way to a region of strongly suppressed transmission. In this high-power region, peak-like features are completely absent.
The transition between the low-and high-power phases can be more thoroughly explored by measuring the transmission at a single drive frequency while sweeping the drive power either from low to high (2c) or from high to low (Fig. 2d). Doing so reveals a significant region exhibiting hysteresis, which is located at the top of the low-power lobes where the transition to the high-power phase occurs. Subtracting the transmission signals for the two different sweep directions clearly marks the hysteretic regime, as shown in Fig. 2e.
To gain insight into this behavior, we model the system as a one-dimensional chain of identical circuit QED elements, as illustrated in Fig. 1d. The corresponding Hamiltonian includes terms for the resonator, qubit, and the resonator-qubit coupling on each site j, hopping of photons between nearest-neighbor resonators, and a coherent drive (acting on site 1 only). Each resonator contributes a single harmonic mode, H r j = ω a † j a j , where a † j and a j are the creation and annihilation operators for photons on site j. The low-lying transmon levels are described as an anharmonic oscillator H q and projectors P N that truncate the Hilbert space to levels N E J /2E C within the transmon's cosine well. (Interestingly, the sign of U only affects the system dynamics, but not the steady state given that the qubit and drive frequencies are tuned accordingly -see Supplementary Information II.A.) The operators b † j and b create and annihilate qubit excitations, E C is the single electron charging energy and E J the effective Josephson energy of the transmon qubit. Note that we neglect disorder effects within this model. Qubit-resonator coupling and photon hopping take the simple forms H rq j = g(a j b † j + h.c.) and H hop j,j = t(a j a † j + h.c.). Within rotating-wave approximation, the microwave drive acting on site 1 is given by H d = (t)a 1 e iω d t +h.c. In our model, we account for qubit relaxation and intrinsic photon loss (at rates Γ and κ, respectively) by employing the standard Lindblad master ρ} is the usual action of the Lindblad damping operator.
The experimentally observed transition is remarkably well captured by simple, quasi-classical mean-field theory that decouples the sites, but allows for mean-field parameters to differ from site to site. Allowing for sitedependent parameters is particularly relevant for our case in which the drive only acts on one end of the resonator chain, rather than on every site.) Within the quasi-classical treatment 20 , the quadrature amplitudes α j = a j and β j = b j play the role of mean-field parameters and obey the equations From these equations, we obtain the steady-state transmission signal S 21 ∼ a j and the second-order coherence function g (2) (0) = a † j a † j a j a j /| a † j a j | 2 , choosing j as the label of the output port resonator (see Methods).
The steady-state transmission (Fig. 2b) reproduces all of the qualitative features of the experimental data. We find that the transition occurs beyond the point where the dispersive approximation holds, and further observe that the mean-field solution predicts an increasing accu-mulation of transmon excitations. Higher transmon levels are a crucial model ingredient, as calculations based on the simpler Jaynes-Cummings lattice do not yield results consistent with experiment. Quasi-classically, the drop in transmission in the high-power phase is associated with chaotic dynamics with parallels to results previously obtained for a driven, dissipative Bose-Hubbard chain 20 .
While bistability obtained within nonequilibrium mean-field theory generally has to be considered with care, it is interesting to note that the mean-field solution reveals a region of bistability and hysteresis consistent with that detected experimentally (see Supplementary  Information, Fig. 2). Finding multiple steady states appears at odds with Spohn's theorem 21 : the steady-state solution of the Lindblad master equation is unique as long as Hilbert space is finite (or can be safely truncated), and minimal conditions for nature and number of relaxation channels are satisfied. Thus, dissipative phase transitions and stationary bistability can only occur in the thermodynamic limit -when the number of lattice sites tends to infinity and/or when truncation fails due to accumulation of excitations in the strong drive limit, such as in the breakdown of photon blockade on a single Jaynes-Cummings site 18,22 .
Bistability and hysteresis can, however, be produced dynamically 23 as recently discussed by Casteels and coworkers in the context of a Bose-Hubbard dimer 24 . As shown in their paper, hysteresis arises from parameter sweeps across a point where the spectral gap of the Liouvillian superoperator L (nearly) closes, in a manner 7  GHz is plotted as a function of power. κ and γ are included for reference to indicate that ADR can be as large as five orders of magnitude slower than relevant timescales of the device. When either γ1→2 or γ2→1 are slower than the duration of the measurement pulse, τm, we cannot reliable extract a characteristic switch rate. In these cases we select the smallest extracted switching time which is larger than 1/τm. Upward (downward) pointing triangles indicate when γ1→2 (γ2→1) are less than 1/τm, circles indicate when both rates are used.
analogous to the Kibble-Zurek mechanism 25 . Similar to the case studied by Casteels et al., we find that meanfield theory can capture certain qualitative aspects of the bistability and hysteresis, but necessarily fails in aspects related to sweep times and quantum fluctuations. The dramatic suppression of transmission and loss of all resonance peaks beyond a certain drive power threshold are indicative of a dissipative phase transition, arising from the intricate interplay of dissipation, driving, and nonlinearity of the system. The crucial quantity for such a transition is the gap in the spectrum of the Lindblad superoperator L. If the real part of one of its eigenval-ues approaches zero, then deviations of the steady state along the "direction" of the corresponding L-eigenstate become increasingly long-lived and ultimately allow for a new steady state to emerge. The negative real part of the eigenvalue λ closest to zero, − Re λ, is known as the asymptotic decay rate (ADR) 26 . An approximation for the ADR can be extracted by single-shot measurements of the dynamics in the bistable region as follows.
We apply a drive with constant frequency and amplitude, and record single-shot time traces of the homodyne amplitude and phase. Our measurements show that the system undergoes switching between two metastable states on timescales large compared to system-intrinsic timescales (Fig. 3a). The state of the system at each point along a single-shot trajectory is classified as either ρ 1 or ρ 2 , and characteristic dwell times are extracted. The statistics acquired from many single-shot trajectories allow us to extract average rates γ 1→2 , γ 2→1 for the switching between the two metastable states ρ 1 and ρ 2 observed at low power and high power, respectively (labeled as 1 and 2 in the figure).
The extracted switching rates allow us to estimate the asymptotic decay rate by adopting a simplified rateequation model 27 describing the probabilities p 1 and p 2 for the system to be in metastable state ρ 1 or ρ 2 (see Methods section for details): Diagonalization of this system yields the stationary and purely decaying eigenmodes ρ s = (γ 2→1 ρ 1 + γ 1→2 ρ 2 )/γ Σ and ρ ADR = γ 2→1 ρ 1 − γ 1→2 ρ 2 with corresponding eigenvalues zero and λ ADR = −γ Σ = −(γ 1→2 + γ 2→1 ). Hence, this simplified model predicts an asymptotic decay rate of − Re λ ADR = γ Σ . Remarkably, the asymptotic decay rate, shown in Fig.  3b-c, reaches a minimum value as low as ∼ 10 Hz, which is five orders of magnitudes lower than the rates set by photon decay and transmon relaxation in our system. This vast timescale discrepancy delivers strong evidence for the onset of a dissipative phase transition. Similar to the situation of equilibrium phase transitions, it is only in the thermodynamic limit that the the spectral gap can fully close and turn the crossover between two steadystate phases into a phase transition in the strict sense 28 .
We gather additional evidence for the approach to a dissipative phase transition by measuring fluorescence power spectra and second-order coherence functions in our system. To this end, two different driving pulse shapes, Fig. 4a-b, are used to access the distinct states of the system. Within the region of bistability, we can perform state initialization either in the low-power phase 1 or the high-power phase 2 by approaching the final drive amplitude ξ either from a lower or a higher drive amplitude. After this initialization period, the two pulses maintain the constant drive amplitude ξ, during which time, the transmitted signal is detected using heterodyne detection with a 32 MHz intermediate frequency. The power spectrum is then obtained by performing a Fourier transform on the heterodyne signal. The second-order correlation function g (2) (0) is measured using techniques outlined in Ref. 29 . Figure 4 indicates that the low-power state can be characterized by a single, coherent drive tone (g (2) = 1) and that the high-power state can be characterized by broadband and multimode (see supplement) emission and bunching (g (2) ≈ 2). In addition, the onset of the high-power state has a stark linewidth broadening of the drive tone and has a region of strong bunching g (2) ≈ 5 for the down pulse as the system transitions to the high-power state. Experimental measurements of g (2) in Fig. 4e agree well with theory results shown Fig.  4f, barring the strong bunching observed with the 'down pulse' (Fig. 4b) at the high ξ side of the bistability region. Based on our modeling, the experiment involves both large numbers of photons and excitations of higher trans-mon levels, and hence may indeed approach the thermodynamic limit necessary for the observation of a dissipative phase transition. This work demonstrates the potential for circuit QED lattices as a controllable platform that can guide a deeper theoretical and experimental understanding of nonequilibrium condensed matter physics.

METHODS
Experimental methods. The cavities of the circuit QED lattice were etched using standard optical lithography and plasma etching techniques from a 200 nm thick Nb film on a 25 × 25 cm 2 sapphire substrate. Transmon qubits were designed to have Josephson junctions with dimensions 200 × 180 nm 2 and 450 × 450 nm 2 and were fabricated according to the "Manhattan" technique outlined by Potts et al. 30 , using electron beam lithography and aluminum evaporation. Similar transmon qubits have coherence times T 1 = 1 µs and coupling constant of g/2π = 265 MHz. Measurements were performed at a temperature of 7.5 mK in a dilution refrigerator, and inside a superconducting solenoid magnet controlled by a room-temperature DC voltage source. Transmission measurements are performed using a network analyzer, switching-rate measurements using standard homodyne detection techniques. All power-spectrum measurements were done by taking the Fourier transform of a heterodyne signal, and g (2) measurements were implemented using the homodyne techniques described by Eichler et al. 29 (see Supplementary Information for further details).

Numerical solution of the mean-field equations.
We solve for the stationary state of the mean-field equations (3) by time evolution and extracting the long-time limit, since root-finding methods are difficult to handle for the large system of nonlinear equations 20 . In the high-power phase the dynamics is chaotic, so that additional time averaging in the long-time limit is required. For instance, the second-order coherence function is obtained by evaluating g (2) (0) = |α(t)| 2 t /| |α(t)| t | 2 , where the time average · t is carried out over a time interval that excludes any initial transient behavior.
Model underlying the ADR estimate. First, consider stochastic switching between two pure states |1 and |2 . The simplest description is based on a two-level Hamiltonian H = E 21 |2 2| where E 21 is the energy difference between the two states, and the master equatioṅ (5) with D[L] denoting the usual Lindblad damping superoperator for jump operator L. The resulting 4 × 4 Liouvillian L is block-diagonal, where one of the two blocks fully captures the dynamics of density matrices of the form ρ(t) = p 1 (t)|1 1| + p 2 (t)|2 2|, where the probabilities p 1,2 obey the rate equation (4). This model can be extended and made more realistic by considering subsets of pure states that make up the two metastable states ρ 1 and ρ 2 , which are likely to be mixed states rather than pure states.

I. DEVICE PARAMETERS
The device consists of 72-coplanar waveguide cavities coupled by capacitors which are formed using gaps in the center pin of the resonator. Each tunneling capacitor is designed to have a capacitance of 20.7 pF, resulting in a tunneling matrix element t/2π = 144 MHz. The threeway couplers are integrated in the lattice to maintain the cavity-cavity hopping rate while introducing a weak coupling to an input/output transmission line to weakly probe the internal behavior of the lattice.
One qubit is coupled to each cavity. From finiteelement simulations, we predict a charging energy E c /h = 180 MHz. Similar qubits measured in other devices have measured coupling rates g/2π = 265 MHz. Each transmon SQUID has a fixed width of 9 µm but the heights are chosen to be a random number between 8 µm and 20 µm, ensuring that a global magnetic field can continuously change the qubit frequencies but will never return to the exact same qubit frequency distribution. Together with previous measurements of E J and simulated values for E c , we expect qubit frequencies to be between 8 and 8.8 GHz. This is confirmed in Fig. 5 which shows a measurement of transmission as a function of external magnetic field. Each mode exhibits frequency shifts of different flux periodicity, indicating the existence of numerous qubits. Based on optical microscope inspection of the sample, it is expected that after the experiment was complete, a minimum 60% of our qubits are fully functional, with another 20% for which the smaller Josephson junction of the SQUID did not make a physical connection; it is unknown if this damage occured due to handling the device after the experiment was complete.

II. THEORY
A. Mapping of positive -negative U The sign of the Hubbard interaction U has a significant impact on the system's energy spectrum. However, if we are merely interested in the steady-state behavior, the sign of U is relatively less important as long as the system remains in the transmon's cosine well. We begin our discussion first with the quasi-classical treatment. For the steady state, the time derivative of any expectation value vanishes. Hence, the quasi-classical dynamical equations reduce to with the qubit-resonator detuning ∆ = Ω − ω. Now, we consider a mapping such that This mapping physically corresponds to tuning both the drive and qubit frequencies in a specific way.
(For our open chain, the sign flip of g and t is allowed by a local gauge transformation.) By taking the complex conjugate on both sides of the dynamical equations together with the mapping, we find: The mapping thus effectively flips the sign of U . In other words, if both the drive and qubit frequencies are tunable, the sign of U is not a concern and the truncation to low-lying transmon levels is justified. The above discussion can be easily generalized to the Lindblad master equation. In that case, instead of mapping α j and β j to their complex conjugates, we map the steady-state density matrix ρ s to its complex conjugate ρ * s . Note that any steady-state expectation value can be calculated using ρ * s through the equation

B. Simulation of hysteresis
Hysteresis appears in the simulation if the quasiclassical dynamical equations have more than one attractor. Hence, we can numerically investigate the region of hysteresis by picking different initial states and checking whether the system evolves to the same or different attractors in the long-time limit. We illustrate this by simulating the transmission for a 20-site chain in Figure  6. (The region of hysteresis is qualitative the same for Transmission as a function of external magnetic field. By changing the strength of the external magnetic field, individual qubit frequencies can be continuously changed. The non-periodic behaviour of the transmission peaks is due to the random areas of the transmons in the cQED lattice. both 20-site and 72-site chains.) Subtraction of the simulated transmission with two different initial states shows a region with more than one attractor. The shape of this hysteresis region is consistent with the experimental data. Readers may notice that the result becomes 'noisy' in the high power regime. This is the result of the chaotic dynamics which makes the time average strongly sensitive to the initial state and the averaging window. However, the hysteresis is reassuring by direct inspection of the system time evolution. In the region of hysteresis, one initial state gives stationary behavior (stable fixed point) in the long time limit while the other gives chaotic dynamics (strange attractor).

III. SWITCHING RATE EXTRACTION
In order to extract switching rates necessary for the calculation of the asymptotic decay rate (ADR), it is necessary to gather statistics on the lifetimes of each state in the region of bistability. To do this, a pulse shown in Fig.  4 (CW) microwave source. After the initialization phase of the pulse, the homodyne amplitude and phase are extracted from the output signal by first passing the raw signal through a series of amplifiers before being mixed down to using a local oscillator (LO) at the drive frequency using an IQ mixer. The I and Q outputs of the IQ mixer are then sent through a 1.9 MHz low-pass filter (LPF) before being amplified, digitized, and sent to the measurement computer (see Fig. 7). The amplitude and phase of the homodyne signal are then extracted by taking The 1.9 MHz filter was chosen to reduce noise in the single-shot signal in an effort to reduce the number of false-counts in the steady state transitions. This filter frequency is well above observed transition rates. The digitizer sampling rate was chosen to be 50 MS/s and the data was down-sampled to include only every tenth data point to avoid memory constraints. After all processing, each trajectory consists of 1.5 × 10 6 data points. Figure 7 outlines the thresholding algorithm used to determine switching rates γ 1→2 and γ 2→1 . For each drive amplitude and frequency, seven 0.3 s trajectories are acquired and the data is placed in histograms for amplitude and phase. It is then determined whether the resultant distributions are gaussian, meaning there is no bistability, or bimodal, meaning that the there is bistability. In the case of a bimodal distribution, data are classified as either being in the high or low-power state with a threshold given by the mean of the peak locations. The measured quantity (amplitude or phase) which has the fewest his- togram counts at the threshold is then used as the measured quantity for the state lifetime determination. New data is then acquired and categorized according to the threshold. Once the new data has been categorized, state lifetimes are extracted. The lifetimes are then binned in a nonuniform histogram bins shown in Fig. 8.
Depending on the population of the histogram bins, the bins can be summed to create histograms with larger time intervals. This scheme ensures that histograms are sufficiently populated for situations involving short or long lifetimes. After the binning procedure has been performed, the resultant distribution is fit to an exponential to extract a characteristic switching time, τ c . From this, the switching rate for the state of interest is γ = 1/τ c .
The physical measurement setup for the jump rate extraction is shown in Fig. 10.

IV. EMISSION PROPERTIES
The transmission shown in Figure 2 of the main text exhibits an abrupt change as the system crosses into the high power regime. As illustrated in Figure 4 of the main text, the high power state can be characterized by broadband emission around the drive tone. When a power spectrum measurement is performed (in this case using a spectrum analyzer) observing the entire range of frequencies with low-power transmission peaks, the system emits at nearly all the low-power transmission peaks of the system. This can be understood qualitatively by con-  9. Multimode Emission. By driving a single mode of the system, in this case at 7.535 GHz, emission is observed both below (a) and above (b) the drive tone. Emission peaks (shown in red) qualitatively match low-power transmission peaks (shown in blue). Emission is plotted in units of dB above the background.
sidering the system Hamiltonian in the eigenmode basis of the cavity chain.
We can rewrite the Hamiltonian in the cavity eigenmode basis by the substitution: a j = N µ=1 W jµãµ , where W jµ is the weight of the µth eigenmode at site j. This gives whereω µ is the frequency of the µth eigenmode. This equation can be thought of as a collection of modes coupled to a bath of transmon qubits. In this basis, each cavity eigenmode interacts with all other cavity eigenmodes through the communal qubit bath. This effective interaction explains the presence of multimode emission shown in Figure 9.  Fig. 2 are performed using the same experimental setup but using a network analyzer (Keysight PNA-X N5241A) to source the RF signal and the returning signal after the two Miteq amplifiers is sent directly back into the network analyzer.