Energy-participation quantization of Josephson circuits

Superconducting microwave circuits incorporating nonlinear devices, such as Josephson junctions, are one of the leading platforms for emerging quantum technologies. Increasing circuit complexity further requires efficient methods for the calculation and optimization of the spectrum, nonlinear interactions, and dissipation in multi-mode distributed quantum circuits. Here, we present a method based on the energy-participation ratio (EPR) of a dissipative or nonlinear element in an electromagnetic mode. The EPR, a number between zero and one, quantifies how much of the energy of a mode is stored in each element. It obeys universal constraints--valid regardless of the circuit topology and nature of the nonlinear elements. The EPR of the elements are calculated from a unique, efficient electromagnetic eigenmode simulation of the linearized circuit, including lossy elements. Their set is the key input to the determination of the quantum Hamiltonian of the system. The method provides an intuitive and simple-to-use tool to quantize multi-junction circuits. It is especially well-suited for finding the Hamiltonian and dissipative parameters of weakly anharmonic systems, such as transmon qubits coupled to resonators, or Josephson transmission lines. We experimentally tested this method on a variety of Josephson circuits, and demonstrated agreement within several percents for nonlinear couplings and modal Hamiltonian parameters, spanning five-orders of magnitude in energy, across a dozen samples.

Superconducting microwave circuits incorporating nonlinear devices, such as Josephson junctions, are a leading platform for emerging quantum technologies. Increasing circuit complexity further requires efficient methods for the calculation and optimization of the spectrum, nonlinear interactions, and dissipation in multi-mode distributed quantum circuits. Here, we present a method based on the energy-participation ratio (EPR) of a dissipative or nonlinear element in an electromagnetic mode. The EPR, a number between zero and one, quantifies how much of the mode energy is stored in each element. The EPRs obey universal constraints and are calculated from one electromagnetic-eigenmode simulation. They lead directly to the system quantum Hamiltonian and dissipative parameters. The method provides an intuitive and simple-to-use tool to quantize multi-junction circuits. We experimentally tested this method on a variety of Josephson circuits, and demonstrated agreement within several percents for nonlinear couplings and modal Hamiltonian parameters, spanning five-orders of magnitude in energy, across a dozen samples.
In this paper, we introduce a circuit quantization method based on the concept of the energy-participation ratio (EPR). We reduce the quantization problem to answering the simple question: what fraction of the energy of mode is stored in element ? This leads to a constrained number between zero and one, the EPR, denoted 45 . This ratio is the key quantity that bridges classical and quantum circuit analysis; we show it plays the primary role in the construction of the system many-body Hamiltonian. Furthermore, dissipation in the system is treated on equal footing by calculating the EPR of lossy element in mode . The EPR method deviates from previous blackbox quantization work 4,6,7 , which uses the impedanceresponse matrix, denoted ( ), where and index ports associated with nonlinear elements. For all pairs of ports, the complex function ( ) is calculated from a finite-element (FE) driven simulation in the vicinity of the eigenfrequency of every mode. Our method re-places these steps with a more economical FE eigenmode simulation, from which one extracts the energy participations and , needed to fully characterize both the dissipative and Hamiltonian properties of the circuit.
To test the method, we compared EPR calculations of circuit parameters to experimentally measured ones for 8 superconducting devices designed with the EPR method, comprising a total of 15 qubits, 8 readout and storage resonator modes, and one waveguide system. The results demonstrate agreement for Hamiltonian parameters spanning over five orders of magnitude in energy. Resonance frequencies were calculated to one percent accuracy, large nonlinear interactions, such as anharmonicities and cross-Kerr frequencies, to five percent, and small, nonlinear interactions to ten percent. This level of accuracy is sufficient for most current quantum information experiments.

II. RESULTS AND DISCUSSION
A. To quantize a simple circuit: qubit coupled to a cavity In this section, we introduce the EPR method of quantum circuit design on a modest, yet informative, example: a transmon qubit coupled to a cavity mode (see Fig. 2). The transmon 46 consists of a Josephson junction shunted by a capacitance. It is embedded in the cavity, which we will consider as a black-box distributed structure. The Hamiltonian of this systemˆf ull can be conceptually separated in two contributions (see Supple- , ì (ì), and ì (ì), respectively, where denotes spatial position. Center inset: ì profile (red: high; blue: low) for the fundamental mode of one of the 3D cavities. Additional FE driven simulations (FE d ) are unnecessary; i.e., the impedance matrix ( ) is not calculated. c The Hamiltonianˆf ull , which includes nonlinear interactions to arbitrary order (see Results ), is computed directly from the eigenanalysis via the EPRs and EPR signs = ±1 of the junctions, . Dissipative contributions due to a lossy element are similarly computed from the loss EPRs ; for linear dissipation, the EPR signs are unnecessary. Direct extraction of Hamiltonian and dissipative parameters from eigensolutions is unique to the EPR method. The geometry of the classical model is modified in an iterative search for the desired dissipative and Hamiltonian parameters (left-pointing arrow). mentary Section A2), whereˆl in consists of all terms associated with the linear response of the junction and the resonator structure, andˆn l consists of terms associated with the nonlinear response of the junction. Restricting our attention to the cavity and qubit modes of the otherwise black-box structure, the analytical form of the Hamiltonian follows from standard circuit quantization 47,48 (see Supplementary Section A4): where and are the angular frequencies of the cavity and qubit eigenmodes defined associated withˆl in , respectively, and whereˆandˆare their annihilation operators, respectively. The Josephson energy can be computed from the Ambegaokar-Baratoff formula adapted to the measured room-temperature resistance of the junction 49 . The junction reduced generalized fluxˆcorresponds to the classical variable ( ) ∫ −∞ ( ) d / 0 , where ( ) is the instantaneous voltage across the junction 47,48 , and 0 ℏ/2 is the reduced flux quantum. The junction flux operator [Eq. (4)] is a linear, real-valued, and non-negative combination of the mode operators (see Supplementary Section A4), and in its expression, and are the quantum zero-point fluctuations of junction flux in the cavity and qubit mode, respectively. It is worth stating that the linear coupling between the cavity and qubit, commonly denoted 46 , is fully factored in our analysis, and is implicitly handled in the extraction of the operators from the electromagnetic simulation.
Our principal aim is to determine the unknown quantities: , , and . As we will show, we extract and compute these quantities from an eigenanalysis of the classical distributed circuit corresponding toˆl in . This includes the qubit-cavity layout, materials, electromagnetic boundary conditions, and a model of the junction as a lumped-element, linear inductor. The eigensolver returns the requested set of eigenmodes and their frequencies, quality factors, and field solutions. By running the eigensolver in the frequency range of interest, we obtain the hybridized cavity and qubit modes, whose eigenfre-quencies and fully determineˆl in (see Supplementary Section C for FE methodology).
To determineˆn l , we need the quantum zero-point fluctuations and , which are calculated from the participation of the junction in the eigenfield solutions. The participation of the junction in mode ∈ { , } is defined to be the fraction of inductive energy stored in the junction relative to the total inductive energy stored in the entire circuit, Inductive energy stored in the junction Total inductive energy stored in mode , (5) evaluated when only mode is excited. Thus, can be computed from the electric ì (ì) and magnetic ì (ì) eigenfields as detailed in Supplementary Section C2; ì denotes spatial position. In the quantum setting, Eq. (5) links ,ˆ, and the state of the circuit, where | denotes a coherent state or a Fock excitation of mode . Note that normal-ordering must be used in Eq. (6); this correct treatment of vacuum fluctuations is detailed in Supplementary Section A6. Simplifying Eq. (6), one expresses the variance of the quantum zeropoint fluctuations and as a function of the classical energy participations , which completely determinesˆn l , and thus completes the description of the system Hamiltonianˆf ull . Here, and can be taken as positive numbers. As we will see in the next section, in the presence of multiple junctions, this is not always true.
Designing experiments with the EPR requires one to further extract fromˆf ull the transition frequencies and nonlinear couplings between modes. Depending on the case, this can be done approximately or exactly using numerical or analytical techniques 20 . This task is easily achieved ifˆn l is a perturbation toˆl in 46 . In this limit, full for our qubit-cavity example can be approximated by the effective, excitation-number-conserving Hamiltonian, see Supplementary Eq. (B.2), whereˆ=ˆ †ˆandˆ=ˆ †ˆdenote the qubit and cavity excitation-number operators, respectively, Δ denotes the 'Lamb shift' of the qubit frequency due to the dressing of this nonlinear mode by quantum fluctuations of the fields, ( ) is the qubit (cavity) anharmonicity, and is the qubit-cavity dispersive shift (cross-Kerr coupling). The Hamiltonian parameters can be calculated directly from the EPR, see Supplementary Section B, Experimentally, the qubit Lamb shift can be obtained as Δ = − /2. Since a single EPR determines the nonlinear interaction for each mode, the parameters and are interdependent, As shown in Supplementary Section A7, the EPRs and obey the constraints 0 ≤ , ≤ 1 and + = 1 .
These relations together with Eqs. (9) and (12) are useful to budget the dilution of the nonlinearity of the junction (see Supplementary Section B2) and to provide insight on the limits of accessible parameters (see Methods). Further, Eq. (13) is used to validate the convergence of the FE simulation.

B. Quantizing the general Josephson system
The simple results obtained in the preceding section will now be generalized to arbitrary nonlinear devices enclosed in a black-box, distributed, electromagnetic structure. While such structures are frequently classified as planar 27,50-52 (2D), quasi-planar 28,29,32,53 (2.5D), or three-dimensional 30,54-56 (3D), we will treat all classes on equal footing. The electromagnetic structure is assumed to be linear in the absence of the enclosed nonlinear devices. For simplicity of discussion, we can consider these devices to be inductive and lumped; distributed nonlinear devices, such as kinetic-inductance transmission lines [57][58][59][60] , can be thought of as a series of lumped ones.
The general nonlinear device that we now consider, referred to as a Josephson dipole, is any lumped, purelyinductive, nonlinear subcircuit with two terminals. The   Figure 3. Schematic representation of the Josephson circuit and its nonlinear elements. a A simple example of a Josephson dipole-a Josephson tunnel junction. b An example of a composite junction, comprising four Josephson junctions in a ring, frustrated by an external magnetic flux Φ ext threading the loop. c Conceptual decomposition of a general Josephson dipole, denoted . For convenience, its potential energy function E (Φ ; Φ ext ) can be Taylor expanded in a sum of nonlinear inductive contributions of increasing order Φ , with relative amplitude , where denotes the index in the series. The energy function can be subjected to external bias parameters Φ ext , such as flux or voltages. d Schematic diagram of a general Josephson circuit conceptually resolved into a purely dissipative (left, ), linear (middle,ˆl in ), and nonlinear (right,ˆn l ) constitutions. key characteristic of the Josephson dipole is that it possesses a characteristic energy function, which encapsulates all details of its constitution. For example, the two-terminal nonlinear device known as the symmetric SQUID 70 is described by the energy func- the generalized flux across the device terminals 47,48 , is the effective Josephson energy, Φ ext is the external flux bias, and the subscript denotes the -th Josephson dipole in the circuit. The flux Φ is defined as the deviation away from the value in equilibrium, as discussed below. To ease the notation, parameters such as Φ ext will be implicit hereafter. Similarly to the example of the single-junction transmon , the energy of a Josephson dipole can be separated in two parts. One part E lin accounts for the linear response of the dipole, while the other E nl accounts for the nonlinear response, where and where the constant sets the scale of the junction energy. This energy scale can be represented by the linear inductance 0 / presented by the Josephson dipole when submitted to a small excitation about its equilibrium.
Frustrated equilibrium. External biases can set up persistent currents in the circuit. These can alter the static [direct-current (dc)] equilibrium of the Josephson system. For example, frustrating a superconducting ring with a magnetic flux sets up a persistent circulating current in the ring. For a Josephson dipole in such a loop, the definition of the flux Φ will differ in Eqs. (14) and (15) as a function of the equilibrium. Due to this adjustment, terms linear in Φ are absent from Eq. (15) by construction. Supplementary Sections A8 and A9 discuss these equilibrium considerations in detail.

Quantum Hamiltonian.
Having conceptually carved out nonlinear contributions from the system Hamiltonianˆf ull , and collected them in the set of E nl functions, we define the linearized Josephson circuit to correspond to everything left over in the system. This linear circuit consists of the electromagnetic circuit external to the Josephson dipoles, combined with their linear inductances . We will use the eigenmodes of the linearized circuit to explicitly constructˆf ull . The eigenmode frequencies and field distributions are readily obtained using a conventional finite-element solver (see Supplementary Section C). The Hamiltonian of the linearized Josephson circuit can thus be expressed as (see Supplementary Section A5) where is the number of modes addressed by the numerical simulation, is the solution eigenfrequency of mode , andˆthe corresponding mode amplitude (annihilation operator), defined by the mode eigenvector. We emphasize that the frequencies will be significantly perturbed by the Lamb shifts Δ , and should be seen as an intermediate parameter entering in the calculation of the rest of the nonlinear Hamiltonian, where is the total number of junctions andˆ Φ / 0 .
In Eq. 17, we have introduced a Taylor expansion of E nl , where the energy and expansion coefficients are known from the fabrication of the Josephson circuit, see Fig. 3(c). For example, for a Josephson junction, the constant is just the Josephson energy, while are the coefficients of the cosine expansion; i.e., is 0 for odd and (−1) /2+1 / ! for even . The expansion is helpful for analytics but does not need to be used in the numerical analysis ofˆn l , see Supplementary Section A.
The Hamiltonianˆf ull is specified since the opera-torsˆare expressed in terms of the mode amplitudes as a linear combination (see Supplementary Section A5). Here, are the dimensionless, real -valued, quantum zero-point fluctuations of the reduced flux of junction in mode . Determination ofˆf ull is now reduced to computing . We achieve this by employing a generalization of the energy-participation ratio.
The energy-participation ratio of junction in eigenmode is defined to be the fraction of inductive energy stored in the junction when only that mode is excited,

Inductive energy stored in junction
Inductive energy stored in mode which is a straightforward extension of Eq. (5), and is similarly computed using normal ordering (see Supplementary Section A6). The EPR is computed from the eigenfield solutions ì (ì) and ì (ì) as explained in Supplementary Section C3. It is a bounded, nonnegative, real number, 0 ≤ ≤ 1. A zero EPR = 0 means that junction is not excited in mode . A unity EPR = 1 means that junction is the only inductive element excited in the mode.
From the EPR , one directly computes the variance of the quantum zero-point fluctuations (see Appendix A6), Equation (21) constitutes the bridge between the classical solution of the linearized Josephson circuit and the quantum Hamiltonianˆf ull of the full Josephson system, up to the sign of .

Universal EPR properties
The quantum fluctuations are not independent of each other, since the EPRs are submitted to three types of universal constraints-valid regardless of the circuit topology and nature of the Josephson dipoles. These are of practical importance, as they are useful guides in evaluating the performance of possible designs and assessing their limitations. As shown in Supplementary Section A7, the EPRs obey one sum rule per junction and one set of inequalities per mode , The total EPR of a Josephson dipole is a quantity that is independent of the number of modes-it is precisely unity for all circuits in which the dipole is embedded. It can only be diluted among the modes. On the other hand, a given mode can accept at most a total EPR of unity from all the dipoles. In practice, this sum rule can be fully exploited only if the bound reaches the total number of relevant modes of the system. The next fundamental property concerns the orthogonality of the EPRs. Rewriting Eq. (21) in terms of the amplitude of the zero-point fluctuation we have where the EPR sign of junction in mode is either +1 or −1. The EPR sign encodes the relative direction of current flowing across the junction. Only the relative value between and for ≠ has physical significance (see Supplementary Figure S8). The EPR sign is calculated in parallel with the process of calculating , from the field solution ì (ì), see Supplementary Section (C.8). We now obtain the EPR orthogonality relationship valid when the sum from 1 to covers all the relevant modes, see Supplementary Eq. (A.60).

Excitation-number-conserving interactions
Thus, as announced, knowledge of the energyparticipation ratios completely specifiesˆn l , through Eqs. (17), (19), and (21). The Hamiltonian can now be analytically or numerically diagonalized using various computational techniques 20 . In this section, our focus will be now to explicitly handle the effect of the nonlinear interactionsˆn l on the eigenmodes. Before treating the case of a general nonlinear interaction, we focus on the leading-order effect ofˆn l in the case of the 'canonical ' which is a generalization of the one found in Eq. (8). In Eq. (25), we have introduced the Lamb shift Δ of mode , the anharmonicity of the mode, and its total dispersive shift (so-called cross-Kerr term) with a different mode, labeled . Each of these parameters is directly calculated from the EPRs. As shown in Supplementary Section (B.3a), for arbitrary and , which we have found useful in handling large circuits, especially for those in excess of 100 modes. We also introduce the diagonal matrices of eigenfrequencies Ω diag ( 1 , . . . , ) and junction energies E J diag ( 1 , . . . , ), which lead to the matrix form of Eq. (26), We have defined the symmetric matrix of dispersive shifts χ, with elements [χ] = . Further discussion of the matrix approach and applications to th-order corrections is deferred to Supplementary Section B2, and the amplitude of an arbitrary multi-photon interaction stemming from the fullˆn l is calculated in Supplementary Section B3.

EPR for dissipation in the circuit
The EPR method treats the calculation of Hamiltonian and dissipation parameters on equal footing. Unlike in the impedance method 4 , one can completely characterize bothˆf ull and the effect of dissipative elements in the  . Comparison between theory and experiment over five-orders of magnitude in energy scale of the system Hamiltonianˆf ull for eight distinct, multi-mode device samples, described in detail in the Methods, including 3D, flip-chip (2.5D), and 3D waveguide architectures incorporating readout and storage resonators and qubit modes. For each device, the dominant parameters inˆf ull , dressed frequencies , bare anharmonicities , and cross-Kerr interactions , were measured and calculated using the EPR method with our open-source pyEPR package 95 . Gray line is of slope one, representing ideal agreement between theory and experiment.
circuit from the eigenfield solutions, ì (ì) and ì (ì). The list of dissipative elements include bulk and surface dielectrics 81-83 , thin-film metals 53,84 , surface interfaces [85][86][87][88][89] , and metal seams 90 . The energy-participation ratio of a dissipative element in mode will be denoted . It is calculated similarly to , as summarized in Supplementary Section D. The participation and the quality factor of the material of this element are used to estimate the total quality factor of mode in the standard way when the fields are not greatly altered by the dissipation ( 1) 30,91-93 , Experimental values of are found in the literature, and some are provided in Supplementary Section D. Equation (29) and the dissipative EPR provide a dissipation budget for the individual influence of each dissipation mechanism in the system, providing a useful tool to optimize design layout for quantum coherence 94 .

C. Comparison between theory and experiment
Applying the EPR method, we designed 8 superconducting samples to test the agreement between the EPR theory and experimental results. We tested several sample configurations, comprising 15 qubits, 8 cavity modes, and one waveguide in three different circuit-QED architectures. The samples were measured in a standard cQED setup, see Methods, at the 15 mK stage of a dilution unit, over multiple cool downs.
Six of the samples were each composed of 2 qubits and one 3D cavity, one sample was composed of 2 qubits and a waveguide, and one sample was a flip-chip, 2.5D system 28 consisting of a flip-chip qubit embedded in a twomode whispering gallery mode resonator 53 (WGMR). The specifics of each sample are discussed in the Methods.
For each sample, we measured the circuit parameters of interest: dressed mode frequencies − Δ , anharmonicities of qubits and high-Q cavities , cross-Kerr frequencies , and input-output quality factors for any readout modes. Our measurement methodology is detailed in the Methods.
The measured parameters were compared to those calculated using the energy-participation method. The linearized Josephson circuit of each sample was modeled in Ansys High-Frequency Electromagnetic-Field Simulator (HFSS). Junctions were modeled as lumped inductors, whose nominal energy was inferred from roomtemperature resistance measurements 49 . To account for the error bars of the measurement and the drift in resistance over time, was adjusted by no more than 10% to fit the measured qubit frequency. To minimize the number of free parameters, we neglect the small junction intrinsic capacitance in our modeling. The tradeoff is a small and estimable systematic offset of the bare simulated mode anharmonicities. We estimate this correction to be on the order of 4% for a = 4 fF. From the eigenfield solutions, we calculated the EPRs and the sign to constructˆand extract its parameters. Detailed steps of the procedure can be found in Supplementary Section C. The results are presented in Tables I-III. Figure 4 summarizes the agreement of the measured and calculated sample parameters, which span five orders of magnitude in frequency. Accounting for , we find that mode frequencies are calculated to one percent accuracy, large nonlinear interaction energies (namely, anharmonicity and cross-Kerr frequencies greater than 10 MHz) are calculated at the 5% level, and small nonlinear interaction energies agree at the 10% level. We highlight that we have used minimal, coarse adjustment to account for shifts in , and otherwise, by neglecting , the calculation is free from adjustable parameters.
The results of Fig. 4 demonstrate the accuracy and applicability of the EPR method. For each device, the EPR results are obtained from a single eigenmode simulation, using full automation of the analysis, provided by our open-source package pyEPR 95 . For current stan-dard applications, we find the agreement sufficient. Further improvements in accuracy would require improved ability to estimate the Josephson dipole energy and its intrinsic capacitance . At the same level of accuracy, improvements in the precision and reproducibility of the implementation and assembly of the Josephson circuit design are needed, such as in chip-clamping techniques, precision machining of the device sample holder and input-output couplers.
Conclusion. An intuitive, easy-to-use and efficient method is needed to design and analyze Josephson microwave quantum circuits. We have described in this article such a method, based on the distribution of the electromagnetic energy in the circuit and its participation in nonlinear and dissipative elements. This so-called EPR method offers physical insight helping the design process, and provides a simple link between the classical circuit and its quantum properties. By comparing our theory to 8 experimental devices incorporating Josephson junctions, we have shown that our method is accurate and applicable to a large range of quantum circuit architectures. It is directly applicable to a broader class of nonlinear inductive elements, such as weak-link nanobridges 38,63 , nanowires 40,41,64,66,67 , and kinetic-inductance thin-films 43,58,59 . While best suited for weakly nonlinear systems, the EPR method is derived within circuit theory without approximations. It can be seen as arising from a change of basis adapted to nonlinear elements, as detailed in Supplementary Section A. In practice, the useful reach of the method is set by the numerical ability to include all relevant electromagnetic modes and to compute the spectrum of the extracted Hamiltonian 20 . We contribute an open-source package pyEPR 95 , which automates the EPR method, and was tested in the design of several further experiments 28,74,96-105 .

Device fabrication.
Unless otherwise noted, samples were fabricated according to the following methodology. Sample patterns, both large and fine features, were defined by a 100 kV electron-beam pattern generator (Raith EBPG 5000+) in a single step on a PMAA/MAA (Microchem A-4/Microchem EL-13) resist bilayer coated on a 430 m thick, double-side-polished, c-plane sapphire wafer, grown with the edge-defined film-fed growth (EFG) technique. Using the bridge-free fabrication technique 106-108 the Al/AlO x /Al Josephson tunnel junctions were formed by a double-angle aluminum evaporation under ultra-high vacuum in a multi-chamber Plassys UMS300 UHV. The two depositions were interrupted by a thermal oxidation step, static 100 Torr environment of 85% argon and 15% oxygen, to form the thin AlO barrier of the tunnel junction. Prior to the first deposition, to reduce junction aging 108 , the exposed wafer surfaces were exposed to 1 minute oxygen-argon plasma clean-  Table I. Two-qubit, one-cavity devices. Summary of measured and calculated Hamiltonian and input-output (I-O) coupling parameters for the six devices described in Methods. Indices D, B, C denote the dark, bright, and cavity modes respectively. The input-output quality factor to the readout cavity is denoted C . For each device, the first (second) row quantifies the measured, , (bare calculated, ) values. The third row quantifies the bare agreement, i.e., ( − ) / . In the anharmonicity column, the bare agreement should be corrected by the systematic shift due to our choice to neglect the junction intrinsic capacitance in our modeling (see Methods). We evaluate the correction to be of order 4%, estimated by taking a nominal junction = 4 fF; hence, an overall corrected agreement of 4.3% for this column.  Table II. Flip-chip (2.5D), one-qubit, one-storage-cavity, one-readout-cavity devices. Summary of measured and calculated Hamiltonian and input-output (I-O) coupling parameters for the device described in Methods. Indices Q, S, C denote the qubit, storage, and readout cavity modes respectively. The input-output quality factor to the readout cavity is denoted C . For each device, the first (second) row quantifies the measured, , (bare calculated, ) values. The third row quantifies the bare agreement, i.e., ( − ) / .
ing, under a pressure of 3 × 10 −3 mbar. After wafer dicing (ADT ProVecturs 7100 ) and chip cleaning, the normal-state resistance of the Josephson junctions was measured to provide an estimate of the Josephson energy, , of the device junctions. The junction energy was to first order estimated by an extrapolation of from room temperature to the operating sample temperature, at approximately 15 mK, using the Ambegaokar-Baratoff relation 109 , where Δ is the superconducting gap of aluminum, is the elementary charge, and ℎ is Planck's constant.
Sample holder. Sample holders were machined in aluminum alloy 6061, seams were formed using thin indium gaskets placed in machined grooves in one of the mating surfaces. Only non-magnetic components were used in proximity to the samples, molybdenum washers, aircraft-alloy 7075 screws (McMaster /Fastener Express) with less than 1% iron impurities, and non-magnetic SubMiniature version A (SMA) connectors.
Cryogenic setup. Samples were thermally anchored to the 15 mK stage of a cryogen-free dilution refrigerator (Oxford Triton 200 ) and were measured using a standard cQED measurement setup 45,51,96 . High-magnetic-permeability, -metal (Amumetal A4K ) shields together with aluminum supercon- Quantum amplifier. The output signal of a sample was processed by a Josephson parametric converter (JPC) anchored at the 15 mK stage and operated in amplification mode 110,111 , before routing to the HEMT. The JPC provide a typical gain of 21 dB with a typical noise-visibility ratio of 6 dB. See Ref. 112 for a review of the parametric amplification.
Frequency and input-output (I-O) coupling measurements. Spectroscopic measurements were used to determine the frequencies of the resonator modes.
Anharmonicities were determined in two-tone spectroscopy 92,113 . Cross-Kerr energies were determined from dressed dephasing measurements 114,115 . In particular, the dressed-dephasing measurement sequence consisted of first preparing the qubit in the ground state, then exciting it to the equator by a /2 pulse. Subsequently, a weak readout tone excited the readout cavity of the qubit for a fixed duration, 10 times the readout cavity lifetime , after which we measure the qubit X and Y Bloch vectors, after waiting for a time 5/ for any photons in the cavity to leak out. By varying the amplitude and frequency of the applied weak-readout tone, we could calibrate both the strength of our readout, in steady-state photon number in the readout cavity, and the value of the cross-Kerr frequency shift between the qubit and readout resonator. The values could be obtained from fits of the X and Y quadratures. For each sample, the coupling quality factor of the readout-cavity mode, denoted , was extracted from the spectroscopic response of the readout cavity at low photon numbers 92,113 , by measuring the scattering parameters, 21 or 11 . To test EPR's robustness to experimental variability and its applicability over wide range of experimental conditions, the presented samples were fabricated in multiple runs and measured in different cooldowns. Some devices were subjected to as many as 6 thermal cycles.
The Hamiltonian parameters and coupling energies for each sample were also calculated, following the EPR method presented in section on the general approach. In particular, we modeled the sample geometry and materials in a FE electromagnetic simulation, as explicated in Supplementary Section C. Our aim in writing this supplementary section has been to provide an easy access point to the practical use of the EPR method, which we hope will benefit the reader, and allow them to adopt it easily. Our choice of simulation software was the Ansys High Frequency Electromagnetic Field Simulation (HFSS), although we emphasize that the EPR ideas translate to any standard EM eigenmode simulation package. Further, we modeled the loss due to the input-output couplers in the simulation as 50 Ω resistive sheets, see Supplementary Section D4. The eigenmode analysis provided the calculated I-O quality factors and Purcell limits. All electromagnetic and quantum analyses, including the extraction of participations form the eigenfields and the numerical diagonalization of the Hamiltonian to extract its quantum spectrum, were performed in a fully automated manner using the freely available pyEPR package 95 .
The mode quality due to the input-output coupling, , was set by the length of the I-O SMA-coupler pin. Its length inside the sample-holder box was measured at room temperature using calipers. This nominal length was used then used in the HFSS model to create a 3D model of the pin inside the sample holder. The quality factor was then obtained from the eigenmode eigenvalue. We remark that the measurement of the pin-length is accurate to no more than 20%; further, it can be affected by various idiosyncrasies, such as bending of the thin SMA center pin. Nonetheless, the predictions of the quality factors for low-Q modes were observed to be very reasonable estimates, and, similarly, the predicted Purcell limits for qubit and high-Q cavity modes were consistent with estimates from measurements.

Two-qubit, one-cavity devices
Device description. We measured 6 samples that were each comprised of two qubits and one cavity. The cavity was a standard, machined aluminum cavity 54 . It housed either one or two sapphire chips, which were either patterned with transmon qubits or simply blank. Each transmon consisted of two thin-film aluminum pads connected by a Josephson junction. We tested two configurations of chips and patterns. Configuration A consisted of one chip with two orthogonal qubits, as depicted in Fig. 5(a). Similarly, configuration B consisted of one chip with two parallel qubits, depicted in Fig. 5(b). The two qubits were aligned parallel to each other; however, unlike configuration B, there was no galvanic connection between them. The results of the measurements are presented in Table I. Samples R1C9, R2C1, R7C1, and R2C1, R3C2 were fabricated in configuration A, sample DT3 was fabricated in configuration B. Three of the sample (R2C1, R7C1, R3C1) were fabricated simultaneously on the same sapphire wafer, all with nominally identical dimensions. Additionally, R2C1 and R7C1 were designed to have nominally the same Josephson junctions energy, . The rest of the samples (R1C9, DT3, and R3C2) were fabricated at different times and on different wafers. The dimensions of their transmons and the inductance of the junctions were designed to be different. Only sample R3C2 was designed to be very similar to the nominally identical sample R2C1 and R7C1, but with adjusted . For samples R2C1, R7C1, R3C1, and R3C2 a second, un-processed, un-patterned, blank sapphire chip was placed in parallel with the qubit carrying chip [see Fig. 5(c)] to purposefully lower the readout cavity frequency, thus bringing it within the JPC amplification band.
Configurations A and B were designed to test the ability of the EPR method to calculate the mixing between strongly coupled modes. The strong coupling was achieved in two distinct ways. First, configuration A used the spatial proximity of the two qubits to yielded a strong capacitive coupling between them, which resulted in large qubit-qubit mixing. Second, instead of spatial proximity, configuration B used a galvanic connection between the qubits to yield strong hybridization. Our two-qubit designs share some similar-inspirit characteristics with the promising recent developments reported in Refs. 116-121, but our implementation is distinct and is designed to provide several unique advantages. Mode structure and interesting physical insights. Configuration A is characterized by strong capacitive coupling between the two transmons, which have different pad sizes, see Fig. 5(c), and hence different normal-mode frequencies. Due to the strong hybridization, each qubit normal mode consists of some excitation in the vertical and some in the horizontal transmon. With some foresight, we will label the vertical mode bright (B), and the horizontal dark (D). The brightmode resonance is higher in frequency, and thus is closer to the resonance of the readout cavity mode (C). This smaller detuning made it a natural choice for designing stronger cou-pling between it, (B), and the readout mode (C). This was implemented by orienting the transmon design that participates in mode (B) vertical.
To understand this design choice, let us first consider the popular analogy 46,122 between circuit-QED and cavity-QED, often used to discuss mode couplings. In this atomic analogy, the transmon qubit is analogous to a real atom inside the cavity. Thus, it can described by an electric dipole moment ì d . Meanwhile, its coupling, cross-Kerr, etc. to the cavity mode are derived from the electric-dipole coupling interaction. In particular, the coupling amplitude is proportional to ì d · ì , where ì is the cavity electric field at the transmon junction. From this analogy, one can infer that the coupling is maximized when the two are parallel, ì d ì , and one could hope to measure a strong cross-Kerr between the bright qubit and the cavity. This successful conclusion is true, but a coincidence. We will shortly discuss how this popular analogy fails spectacularly for the dark mode in configuration B. Instead, we will argue that a correct way to understand the nonlinear coupling between the two modes is through the participation ratio, which will provide the correct coupling for both configuration A and B.
Before proceeding to configuration B, we note one further useful features that configuration A exhibits. In particular, while the bright qubit mode can be Purcell limited 116,123 , the dark mode is simultaneously Purcell protected. Thus, one can potentially achieve a high ratio in the I-O bath coupling of the two qubits.
Configuration B has two qubit modes, which we will also label dark (D) and bright (B). Since both transmons are designed with the exact same transmon pad geometry and junction energy , see Fig. 5, we can expect that no single junction is preferred, due to the symmetry of the sample. This is in sharp contrast to the asymmetric energy distribution in configuration A. Returning to configuration B, we can estimate that in each qubit mode, both junctions participate equally and with near maximal allowed participation, If the two transmons were well-separated spatially and not connected, they would be uncoupled. However, the galvanic connection between the two lower pads, see Fig. 5(a), results in a very strong hybridization and splitting between the nominally identical transmons. The result of the strong hybridization is a symmetric and antisymmetric combination of the two bare transmons. In other words, the hybridization results in a common mode, namely (B), where both junctions oscillate in-phase, and a differential mode namely (D), where both junctions oscillate out-of-phase. These phase relationships are captured by the signs: In an attempt to understand how these hybridized qubit modes will couple to the cavity mode (C), let us first consider the atomic analogy again. When the two junctions oscillate in phase, in the (B) mode, the net dipole moment of the bright mode, ì B , must be large, since it is the sum of the two junction dipole contributions. Secondly, ì B must be oriented in the vertical direction, parallel to the cavity electric field ì . Hence, we would conclude that the bright mode coupling is large, ì B · ì 0, and there should be a strong cross-Kerr interaction between the cavity and bright qubit. Continuing the analogy in the case of the dark mode, we would deduce that the net dipole moment of the dark mode is zero, since the two junctions oscillate out of phase, and cancel each other's contribution, ì D = 0. Thus, we should not expect any coupling between the dark qubit and the cavity mode, ì D · ì = 0. To the contrary of this conclusion, as can be seen in the measured results in Table I, the nonlinear coupling of the dark and bright qubit to the cavity is nearly equal. The atomic analogy and the dipole argument have failed completely. We can understand the origin of this failure and how to arrive at the correct conclusion by using the energy-participation ratio. As embodied in Eq. (26), in the dispersive regime, the nonlinear coupling between two modes, in this case a qubit and cavity, is given by the overlap of the EPR distribution. In particular, the cross-Kerr amplitude between the dark qubit and the readout cavity mode is given by where both junctions have the same junction energy , and (resp: ) denotes the dark qubit (resp: cavity) mode frequency. The signs, used in the atomic dipole logic do not factor into the coupling, because the Josephson mechanics is fundamentally different. To obtain BC , one can replace the label 'D' with 'B' in Eq. (34). Then, it is easy to use Eq. (31) to show that the ratio of two Kerr couplings is not zero, but rather of order unity, Failure of the conventional dipole approach. We showed that although the heuristic atomic analogy seems seductively accurate, it fails completely in some cases to predict the nonlinear couplings. Instead, one can use the intuition and calculation method provided by the EPRs.
As an added note, we observe that Eqs. (32) and (32) embodies the orthogonality of the participations, see Eq. (24). We also remark that although the atomic analogy fails in the case of the nonlinear couplings, it can yield some guidance when considering the linear mixing of the modes, useful for discussing the Purcell effect. To illustrate, let us briefly extend the atomic analogy. The dipole-like coupling between the bright mode and the cavity suggests that the bright mode will inherit some coupling to the environment, mediated by the cavity. Thus, since ì B · ì 0, we can expect the bright qubit to potentially be Purcell limited. In contrast, since ì D · ì = 0, we could expect the dark qubit to be Purcell protected. Both of these qualitative Purcell predictions are valid, but to quantify them, we will use the EPR method and FE eigenmode simulation of the sample, as will be discussed shortly. Experimental results. Table I summarizes the results of the agreement between the measured and calculated Hamiltonian parameters for all two-qubit, one-cavity samples. The three modes in each sample are labeled dark (D), bright (B), and cavity (C); the reason for this convention is described above. In all samples, the qubits were designed to be in the dispersive regime with respect to the cavity, which was detuned by 2-4 GHz. However, in a large number of the samples, the two qubits were strongly hybridized, often necessitating higher-order nonlinear corrections to be included in the calculation. This strong hybridization was used as a test of the theory in this more challenging and fickle regime. In total, for each sample we measured and calculated 8 frequency parameters and one dimensionless, coupling quality factor, , of the readout cavity mode. In particular, in the low-excitation limit, the nonlinear interactions among the modes were characterized by the effective dispersive Hamiltonian whereˆ,ˆ, andˆdenote the dark, bright, and cavity photon-number operator, respectively. The coupling of the resonator mode to the bath is given by the Lindblad superoperator term D [ˆ] , where = / , and is the density operator.
We remark that all samples in configuration A demonstrated a large asymmetry in the Kerr coupling between the bright-to-cavity and dark-to-cavity coupling, . In contrast, samples in configuration B demonstrated near equal coupling, ≈ . In both configurations A and B, the dark mode was Purcell protected, we calculated a Purcell coupling factor of Purcell 10 7 , using the eigenmode method described in Supplementary Section D4. On the other hand, the bright mode was somewhat Purcell limited, Purcell ≈ 10 6 . From the relative Rabi amplitudes of the dark and bright qubits, we could verify the order of magnitude scaling calculated for the Purcell effect.

Two-qubit, single-waveguide devices
We measured a two qubit sample inside of a waveguide. of the waveguide, and placed /4 away from the termination wall, at the measurement frequency. The rest of the experimental setup was identical to that described in section 'Methods of the experiment'. The two qubit modes were labeled dark and bright, similarly to the samples discussed in the section 'Two-qubit, one-cavity devices.'. Table III presents the agreement between the measured and calculated key Hamiltonian parameters of the sample. These consist of the two mode frequencies, two qubit anharmonicities, and the strong cross-Kerr interaction between the two qubits.
3. Flip-chip (2.5D), one-qubit, one-storage-cavity, one-readout-cavity devices We also designed a multilayer planar 28 (2.5D) circuit-QED sample, depicted in Fig. 7, with the EPR method. It consisted of high-Q storage mode (S), one low-Q readout cavity (C), and one control transmon qubit (Q). The two cavity modes were formed in the footprint of a single whispering gallery mode resonator 53 . The three modes were in the dispersive regime, and the storage mode was used to encode and decode quantum information, as well as to observe parity revivals. Details of the sample design have been reported in Ref. 28.
The agreement between the measured parameters of the sample and those obtained by the EPR calculation methods are presented in Table II.
Data availability. Data are available from the authors on reasonable request.

References 36
A. Theoretical foundation of the energy-participation method Supplementary Table I: Table of key symbols and relationships used in the derivation of the energy-participation-ratio method.

Basic Josephson circuit variables and parameters
Symbol Value Description Superconducting magnetic-flux quantum-the ratio of Planck's quantum of electromagnetic action ℎ and the charge of a Cooper pair 2 .
Generalized magnetic flux of circuit branch , at time . The instantaneous voltage ( ) across the terminals of branch at time is equivalently denoted Φ ( ). In general, , but we take the initial flux Φ (−∞) to be zero, corresponding to the circuit in an equilibrium state (see Sec. A8). With this convention, Φ is a deviation away from equilibrium. This variable is analogous to the elongation of a mechanical spring.
Φ ( ) Generalized magnetic-flux deviation of a non-linear device, the -th Josephson dipole.
The nonlinear component of the energy function of junction . For certain situations, it may be favorable to select the partition of E such that E nl contains linear interactions. If possible, it is often convenient to expand E nl in a Taylor series around the circuit operating where are the dimensionless coefficients of the expansion.
Column vector consisting of the flux deviations of all circuit branches enclosed in the minimum spanning tree. The roman subscript t denotes the tree. The flux of individual tree-branches are denoted Φ t 1 , Φ t 2 , . . .

Linearized Josephson eigenmodes
, The label indexes sets associated with the eigenmodes of the linearized Hamiltonian H lin . We consider ∈ {1, . . . , }, where is the total number of eigenmodes of relevance determined by the context. That is, is either the total number of eigenmodes of H lin or of the relevant eigenmodes included in the finite-element simulation.
Hamiltonian function of the Josephson system. It can be partitioned into H lin , a part comprised of purely linear terms (quadratic in Φ t , or equivalently Φ m ), and H nl , a part generally comprised of non-linear terms (higher-than quadratic in Φ t , or equivalently Φ m ).
Φ m , Q Column vectors (of length ) whose elements are the eigenmode flux Φ (considered as the generalized position) and charge (considered as the generalized momentum) canonical variables, associated with the Hamiltonian H lin .
Eigenfrequency of the -th mode of H lin , carries a dimension of circular frequency.  The quantum zero-point fluctuation (ZPF) of the reduced magnetic fluxˆof Josephson dipole due to mode ; i.e., for a given mode, the ZPF magnitude gives the nonzero standard deviation of the magnetic flux in the ground state. The amplitude of the ZPF is determined by the physical structure of mode , and can be understood in terms of the energy-participation ratio and its sign .
In this section, we derive the EPR method from first principles. We first review quantum electromagnetic circuit theory (Sec. A1), then use it to find the quantum eigenmodes of the circuit (Secs. A2-A5). In Sec. A6, we define the EPR of a Josephson dipole in a quantum mode, and use it to find the quantum zero-point fluctuations (ZPF). The universal properties and sum rules of the EPRs are detailed in Sec. A7.
In Section A8, we detail use of the EPR method for biased systems-those that incorporate active elements (such as current and voltage sources) or external bias conditions that result in persistent currents (such as a frustration by an external magnetic flux).
The EPR derivation consists of a series of exact transformation. No approximation are made in deriving the EPR or in using it to find the quantum Hamiltonian of a general Josephson system. In this sense, the results are universal. Approximation are made in practice when using numerical methods.

A1. Review of electrical circuit theory, and the Josephson junction
We briefly recount the basic formulation of classical electrical theory. This formulation will be used as a stepping stone in the following derivation. We focus on the lumped-element regime. An electrical element is said to be in the lumped-element regime when its physical dimensions are negligible with respect to the electromagnetic wavelengths considered in the analysis. In other words, self-resonances and parasitic internal dynamical (A.1) By convention, the reference orientation of voltage is opposite to that of current. This is not the convention typically used in Lenz's law; in electromagnetism, the current-density vector ì and the electric-field vector ì are usually projected on the same orientation.
From charge conservation, it follows that the timeinstantaneous current ( ) and the charge ( ) having passed through the element, obey a relation similar to that of the voltage and flux, Reference state and initial conditions. In Eqs. (A.1) and (A.2), we now assume zero-valued initial condi-tions, Φ (−∞) = (−∞) = 0. In the case of a circuit frustrated by sources such that the equilibrium system state is non-zero, we can define the variables , , , and Φ to denote deviations away from the global equilibrium of the circuit, as discussed in more detail in Sec. A8. By way of analogy, imagine a mechanical spring. The spring is stretched from its rest position to a new equilibrium by a second stretched spring. The deviations of the spring are measured from the equilibrium position of the combined system (this is what we mean by global equilibrium), not the spring in isolation.
Power and energy. The instantaneous power ( ) delivered to an element and the total energy absorbed by the element E ( ) are Capacitive and inductive elements. A capacitive (resp., inductive) element is described by an algebraic relationship between charge and voltage (resp., flux and current). Let us introduce the simplest linear, passive, time-invariant elements. The simplest capacitor is defined by the constitutive relationship ( ) = ( ), where is its capacitance-a positive, real constant. The dual relationship Φ ( ) = ( ) defines the simplest inductor, where is its inductance, also a positive, real constant. In terms of the flux across the element, the energies of the simple capacitor and inductor are respectively (and assuming Φ (−∞) = (−∞) = 0). Since , ≥ 0 and Φ ( ) ∈ R for all , we can verify that the total energy gained by these elements is always positive, as required for passive elements.
Josephson tunnel junction. The chief non-linear element used in circuit quantum electrodynamics is the Josephson tunnel junction 37,48,126 , characterized by the flux-controlled inductive relationship ( ) = 0 sin (Φ ( ) / 0 ), where 0 ℏ/2 is the reduced magnetic flux quantum and 0 is the junction critical current. The energy function of the junction in terms of flux is where 0 0 denotes the Josephson energy. For small deviations about Φ = 0, ignoring the constant energy offset, To lowest order, the junction responds as a linear inductor with inductance / 2 0 [compare this to Eq. (A.4)] It is useful to introduce the reduced magnetic flux Φ/ 0 of the junction. Single-vs. multi-valued energy functions. The Josephson junction exhibits a fundamental asymmetry. Inverting the current-flux relationship leads to a multivalued function with an infinite number of branches; i.e., = sin −1 ( ( ) / 0 ) + 2 or = − sin −1 ( ( ) / 0 ) + 2 , where is some integer. For flux-controlled inductors, the multi-valued situation can be avoided by favoring a description in terms of flux, rather than charge.
Junction in a frustrated circuit. Embedding the junction in a frustrated circuit can lead to a non-zero equilibrium value for . For deviations − eq away from the equilibrium value eq , see also Eq. (A.64), We can identify the differential inductance of the junction at = eq to be eq = /cos eq . In Sec. A8, we discuss how can be taken as a deviation away from the equilibrium value eq , in effect mapping − eq ↦ → , see Eq. (A.63). Also in Sec. A8, we discuss sources terms such as the one presented by sin eq − eq . Such energy terms linear in turn the junction into an active component, capable of supplying current; i.e., we observe that the current relationship of the junction to lowest order ( ) = E Φ = 0 sin eq + −1 eq Φ ( ) + O Φ 2 contains the constant term 0 sin eq , which acts like a current source.
Relationship between the gauge-invariant superconducting phase difference of a tunnel junction and the reduced magnetic flux (compact vs. non-compact variables). In superconductivity, the Josephson energy coupling two small superconducting islands has a cosine dependence on the gauge-invariant phase difference of the superconducting phases of the two islands 70,127 . This macroscopic variable is a phase angle-a compact variable in the half-open interval ∈ [0, 2 [; in contrast with the non-compact variable ∈ [−∞, +∞], which must be used in circuits where a superconducting wire connects the two sides of the junction.
In our treatment of the Josephson circuits so far, we have completely ignored this subtlety. Rather, we have based our discussion on the non-compact variable . The relationship between these two collective, macroscopic variables is = Φ/ 0 mod 2 . Although the variable is non-compact, the associated wavefunction ( ) is submitted in practice to constraints (like confinement in one or a few potential wells) such that it is decomposed onto a basis set of wavefunctions that are indexed by a single discrete (rather than continuous) index; for example, the Fock basis ( ), with ∈ N. Quantum-mechanically, the representation of ( ) and ( ) is therefore not very different since ( ) is represented by the discrete rotor basis, ( ), with ∈ Z.
A broken symmetry. The compact support of corresponds to a symmetry that is usually broken in most circuits-that of the impossibility to distinguish between different values of differing by 2 . Losses associated with the junction, or coupling to other elements, such as in the RF-SQUID or fluxonium qubit 128 , render 2 turns of macroscopically distinguishable-hence demanding a description in terms of non-compact support corresponding to a point on an open-ended line rather than a circle. For certain special cases, however, it may be advantageous to retain the compact support version. But even in such cases, one can start with the noncompact version of and recover the compact version by a limit procedure 46 Similar results can be obtained for current controlled inductors and for generalized capacitors.
A network of circuit elements. An electrical circuit is an interconnected collection of circuit elements. The connectivity of the elements can be described by an oriented graph. Each branch in the graph corresponds to one element. The -th element is associated with the instantaneous voltage Kirchhoff 's two universal circuit laws. The following two laws are universal and topological in nature. They describe relationship among the branch variables, independent of the constitution of the branch elements. In other words, they apply to nonlinear, time-dependent, and even hysteric elements.
Kirchhoff where the sum runs over all branches that form the -th loop. For a given branch , the positive (resp., negative) sign in Eq. (A.9) is selected if its flux reference direction aligns (resp., is opposite to) the loop orientation. Algebraically, KVL leads to a set of constraints among the network variables. Thus, in the context of a Lagrangian description of the circuit in which the generalized position variables are taken to be fluxes Φ , the KVL conditions express a set of holonomic constraints that need to be eliminated in order to obtain a Lagrangian of the second kind 129 .
Kirchhoff 's current law (KCL) is a statement of the conservation of charge: at every node in the circuit, the algebraic sum of all the current leaving or entering the node is equal to zero. Recast another way, for all branches connected to node , The negative sign is chosen for branches whose current reference direction points toward the -th node. In the flux-based Lagrangian description of the circuit, the KCL algebraic conditions become the Lagrangian equations of motion.
Eliminating the KVL constraints using the method of the minimum spanning-tree. The set of KVL algebraic equations defined in Eq. (A.9) reduce the number of independent branch fluxes Φ . We can systematically choose a minimal set of independent branch fluxes using the minimum spanning-tree graph method 48,126,130,131 . For our derivation, it is not necessary to explicitly construct the tree. A set of branches from the graph can be selected to form a complete and minimum spanning tree. In general, there are many satisfactory tree sets of branches. Different trees are related by a simple algebraic transformation; similar to a basis change. The branches that belong to the spanning tree can be labeled t 1 , t 2 , . . . The flux of the -th spanning-tree branch is denoted Φ t ( ). In subscripts, roman (resp., italic) symbols denote labels (resp., variables). The spanning-tree can be organized in a column vector which serves the purposes of a basis for the description of the circuit. The branch-fluxes not in the spanning tree (links or chords of the graph) are obtained by a linear transformation of Φ t . We define the energy functions and Lagrangian of the system in terms of Φ t in the Sec. A3.

A2. The Josephson system and its non-linear Josephson dipoles
In the main text, under Quantizing the general Josephson system, we introduced the notion of a Josephson system-a general electromagnetic environment that incorporates nonlinear devices, referred to as Josephson dipoles. The Josephson system is treated as a distributed black-box structure.
Discretization. We aim to model the Josephson system as realistically as possible. To account in detail for its physical layout, materials, boundary conditions, and dipole structures, we aim to leverage conventional electromagnetic analysis techniques, such as the finiteelement (FE) method. The FE method subdivides the physical layout of the system. Using a set of basis functions, the electromagnetic circuit is discretized 132,133 . The discretized circuit can be represented by a lumpedelement model. In principle, we can take the limit of infinite subdivision. The Josephson dipoles are assumed to be the only non-linear elements in the circuit. All other elements are linear, as representative of the linear nature of Maxwell's equations.
Dissipation. As the object of interest is the control of quantum information in the system, in this section, we focus on systems with low dissipation. This condition requires that the quality factor of all modes of relevance is high, 1, where is the mode index. In Sec. D, we treat dissipation as a perturbation to the lossless solutions.
Josephson dipole. The Josephson dipole was introduced in the main text, under Quantizing the general Josephson system. For simplicity of discussion, here, we treat the dipole as a lumped, two-terminal, flux-controlled element.
The -th Josepson dipole in the system is fully specified by its energy function E Φ ; Φ ,ext , see Eq. (A.8). The energy depends only on the magnetic flux Φ ( ) across the terminals of the dipole and on any external parameters Φ ,ext that control the energy landscape; these can include a voltage bias, a current bias, and an external magnetic field bias. The index runs from unit to the total number of Josephson dipoles in the circuit. This formulation is rather general. It encapsulate the wide span of Josephson dipoles discussed in the main text. These include simple devices, such as the Josephson tunnel junction or a nanowire, and also composite devices, such as SQUIDs, SNAILs, or more general subcircuits. The underlying physical phenomenon giving rise to the low-loss, non-linearity of the dipole is immaterial.
Partition of the Josephson dipole energy-function. It is always possible to partition E Φ ; Φ ,ext into a linear and non-linear part, (A.12) This division is purely a conceptual one-the Josephson dipole cannot be physically divided into a linear and nonlinear part. By selecting an equilibrium point of the circuit and defining the branch fluxes Φ as deviations away from the equilibrium [discussed in more detail in Sec. (A8)], the partitions take the concrete form where is an overall scaling factor of the energy function, defined in Eq. (A.13a); we refer to it as the Josephson dipole energy scale. In Eq. (A.13c), we have introduced the Taylor series expansion of E nl around the equilibrium state of the circuit and have introduced the dimensionless expansion coefficients . We stress that the expansion is not needed for the EPR method. It is merely a convenient tool for working analytically with weakly non-linear circuits, such as the transmon qubit. For notation simplicity, in Eq. (A.13) the dependance of E lin , E nl , , and on the external bias parameter Φ ,ext is made implicit, and we will continue to do so henceforth; in other words, keep in mind that Φ ,ext and Φ ,ext . Example of a Josephson dipole: the Josephson junction. We illustrate the partitioning construction defined in Eq. (A.13) using the example of the Josephson tunnel junction. For an un-frustrated junction, it follows from Eq. (A.5) that is the Josephson energy. The energy function E lin is associated with the linear response of the junction. It presents the inductance = 2 0 / . The energy E nl is associated with the response of non-linear inductor. The expansion coefficient of E nl as defined in Eq. (A.13c) are In partitioning E , we included all of the linear response of the junction in E lin , leaving none for E nl ; i.e., E nl lacks quadratic terms in Φ . However, this is not required. There are certain cases for which retaining some part of the linear response in E nl is advantageous. The equilibrium flux of the junction shifts to a nonzero equilibrium value eq, , as determined by the circuit equilibrium considerations (see Sec. A8). Applying the partition defined in Eq. (A.14) to Eq. (A.7), we find that in terms of the out-of-equilibrium flux deviation ,

A3. Energy of the Josephson circuit and its Lagrangian
Capacitive energy. The total capacitive energy E cap of the Josephson system is simply the algebraic sum of the total energy of all its capacitive elements. Using Eq. (A.4), and summing over all capacitive branches, where is the capacitance of branch . Each of these fluxes can be expressed in terms of the linearly-independent spanning-tree fluxes Φ t . The energy function is quadratic in Φ t , where C is the capacitance matrix of the circuit 47,48,126 . It follows from KVL and the constitutive relationships of the capacitors that C is a positive-definite, real, symmetric (PDRS) matrix. In the continuous limit of space, the total capacitive energy E cap can be found using Eq. (C.3). Inductive energy. The total inductive energy in the circuit E ind is similarly the algebraic sum of the total energy of all circuit inductive branches; i.e., E ind ∈ind. E (Φ ). In the Josephson system, inductive branches come in two distinct flavors: linear and nonlinear. Physically, linear inductive branches are associated with the geometry and magnetic fields. We denote their total energy E mag . The non-linear inductive branches (Josephson dipoles) are generally associated with the kinetic inductance of electrons. We denote their total energy E kin . Hence, The energy of the linear branches E mag is the dual of Eq. (A.17). It is also a quadratic form, where the inductance matrix L −1 mag completely describes all linear, magnetic-in-origin inductances in the circuit 47,48,126 . Due to its nature, L −1 mag is PDRS. In the continuous limit of space, the total magnetic inductive energy E mag can be found using Eq. (C.4). The total inductive kinetic energy of the circuit, associated with the non-linear dipoles, is This energy is not stored in the magnetic fields. However, we can group the magnetic energy and the linear part of E kin together to express the inductive energy as a partition of linear and non-linear contributions where we have introduced the total inductance matrix of the circuit L −1 and the total non-linear energy function of the circuit E nl . For later use, we can obtain the series expansion of E nl in terms of that of E nl , defined in Eq. (A.13b), Generalized coordinates for the Lagrangian. We have expressed E cap and E ind in terms of the independent set of spanning-tree fluxes Φ t . We could equivalently have expressed E cap and E ind in terms of the charge variables . However, as discussed in Sec. A1, the flux-controlled Josephson dipoles present an asymmetry which favors treatment in the flux basis. We hence employ Φ t as the generalized position coordinates in the Lagrangian description of the circuit, and Φ t as the generalized velocity.
Lagrangian of the Josephson circuit. The Lagrangian function of the Josephson circuit follows from KCL. Energy functions with Φ t as their argument (resp., Φ t ) play the role of potential (resp., kinetic) energies. The system Lagrangian is the difference of the total kinetic and potential energy functions, where we have partitioned the Lagrangian into a linear L lin and nonlinear L nl part.
Equation (A.24) explicitly constructs the Josephson system Lagrangian. It partitions it into a linear and nonlinear part. In Sec. A4, we diagonalize L lin to find the eigenmodes of the linearized system. In Sec. B, we treat the effect of L nl .

A4. Eigenmodes of the linearized Josephson circuit
In this sub-sectino, we diagonalize L lin and find its eigenmodes, eigenfrequencies and eigenvectors (i.e., spatial-mode profiles). These intermediate result provide a key stepping stone on our path to quantizing the Josephson system and treating L nl . The process conceptually parallels that taken by the finite-element (FE) electromagnetic (EM) solver in an eigenanalysis of the linearized Josephson system.
The Lagrangian L lin is the sum of two quadratic forms, see Eq. (A.24). We use the standard method for their simultaneous diagonalization 129 , based on a series of principle-axis transforms. We then transform the Lagrangian into a diagonalized Hamiltonian.
Diagonalizing the inductance matrix. Since the inverse inductance matrix L −1 is a PDRS matrix, we can diagonalize it with a real orthogonal matrix O L , obey- which motivates the principle-axis transformation of the magnetic flux defined by whereΦ is a rotated and then scaled version of Φ t . Under this transformation, the transformed Lagrangian function becomes Since the capacitance matrix C is PDRS and it is transformed by a rotation and then a dilation, it follows thatC is also PDRS. More generally, the eigenvalues of a matrix are invariant under a similarity transform, such as the one employed in Eq. (A.29).
Diagonalizing the capacitance matrix. SinceC is PDRS, we diagonalize it with a real, orthogonal transformation OC, such that where ΛC is the dimensionless, diagonal matrix constructed from the eigenvalues ofC and I C is the identity matrix with physical dimensions of capacitance.
Employing the orthogonality transformation O T C in a manner similar to the one used with O L , we rotate (but do not scale) the coordinates for a second time. Under this second principle-axis transform, we define the eigenmode magnetic flux variable in terms of which the Lagrangian L lin is diagonal, where the nonlinear part under the transformation is

Equations of motion. The Lagrangian equations of motion,L
where the linear and nonlinear parts of the Hamiltonian expressed in the eigenmode coordinates are where the diagonal eigenfrequency matrix Ω of the Hamiltonian H lin is and 1 is unity carrying physical dimensions of inductance. The entries of Ω are the eigenmode frequencies of the linearized circuit. These correspond to the eigenfrequencies solved for by the FE eigenanalysis.
Eigenvectors of the Josephson system. The eigenvector matrix E relates the spanning-tree fluxes Φ t to the eigenmode ones Φ m , It is found by concatenating the principle-axis transformations defined in Eqs. (A.27) and (A.31), The eigenvector matrix E is real and positive-definite, since it is the produce real, positive-definite transforms. It is dimensionless and, in general, non-symmetric. It is related to the square root of the inductance matrix, EE = LI −1 H and E −1 E −1 = L −1 I H . The eigenvector matrix E represents the eigenfield solutions found in the FE analysis. It is key in determining the quantum zero-point fluctuations of the mode and dipole fluxes, as shown in the following section.

A5. Quantizing the Josephson circuit
We quantize H full using Dirac's canonical approach 134 . Before taking passage from classical to quantum, we introduce the complex mode amplitude operator -the classical analog of the bosonic amplitude operatorˆ. This provides a direct path to second quantization in the eigenmode basis of H lin .
Complex action-angle variables. We define the vector of action-angle variables α = ( 1 , . . . , ) T by the noncanonical, complex transformation where 1 H is unity with dimensions of inductance. The normalization 1/ √ 2ℏΩ is chosen so that the Poisson bracket of the action-angles is , * P = 1/( ℏ) . In terms of the action-angles, the Hamiltonian remain diagonal, (A. 43) whereˆ Φ / 0 is the reduced magnetic-flux operator for Josephson dipole . This is Eq. (14) of the main text. In the following, we decomposeˆin terms ofˆ; i.e., in second quantization with respect to the eigenmodes of H lin . In writing Eq. (A.42), we have omitted the 1 2 ℏ ground-state energy of every mode. In the limit of infinite discretization, the sum of these ground energies tends to infinity, a standard conceptual difficulty in quantum field theory. Since physical experiments observe changes in the energy of the field, the vacuum energy can be neglected (except for special cases; e.g., Casimir effect). To proceed, we introduce helpful notation.
Notation for vectors of operators. A bold symbol typeset in roman, such as Φ t , denotes a vector or matrix, whose elements are constants or variables. Since variable, such as Φ t 1 are promoted to quantum operators,Φ t 1 , we can accommodate vectors of operators in our notation with a hat symbol; e.g., The spanning-tree flux operatorΦ t corresponding to the -th spanning-tree-branch flux variable Φ t . Similarly, the eigenmode flux operatorΦ m corresponding to the -th eigenflux variable Φ m . Zero-point fluctuations of the eigenoperators. Inverting Eq. (A.39), one finds the vectors of eigenflux and eigencharge operators, respectively, where the diagonal matrices of the quantum ZPF of the operators are We recall that the operatorΦ = Φ ZPF ˆ † +ˆ has a zero-mean, gaussian-distributed distribution in the ground state; i.e., its mean is 0 Φ 0 = 0 and its vari- Interpretation of the eigenmode ZPFs and effective mode inductances, and impedances. The mode impedance = Φ ZPF / ZPF = 1 H and ZPF Φ ZPF = √︁ ℏ /2 × 1 H depend only on . These quantities are in general removed from direct physical meaning. Physical meaning is extracted by using them as computational tools to calculate the quantum ZPF of the branch fluxes and charges. The EPR method negates the need to explicitly compute them; rather, we only directly work with physically measure quantities, such as the eigenfields of the modes and the Josephson dipole EPRs, as discussed in the next sub-section. The ( , ) element of Φ ZPF t is the ZPF of the -th spanning-tree branch due to mode and is either positive or negative. The overall sign is arbitrary. Hover, the relative sign between two branches and in the same mode determines if they are excited by the mode in-phase or out-of phase with each other. The eigenvalue matrix E is constructed canonically so as to correctly relate the fluctuations. The eigenvector matrix of the standard product C −1 L −1 , obtained in the Lagrangian equations of motion, is not canonical and cannot be used in Eq. (A.47).
Second quantization of the dipole flux operator. Since E is real, the reduced-magnetic-flux ZPF amplitudes can be chosen to be real-valued numbers; Eq. 14(c) of the main text follows,

A6. Energy-participation ratio (EPR)
We define energy-participation ratio in the context of quantum circuits. We motivate the omission of vacuum energy contributions and use the EPR to find the quantum zero-point fluctuations . Definition of EPR in plain english. The EPR of Josephson dipole in eigenmode is the fraction of inductive energy allocated to the Josephson dipole when mode is excited.
Interpretation of the energy-participation ratio. The EPR quantifies how much of the inductive energy of a mode is allocated to a Josephson dipole. The lowest possible participation is zero-the Josephson dipole inductor does not participate in the mode. When the mode is excited, none of the excitation energy flows to the dipole. The largest possible participation is unity. When the mode is excited, all of its excitation flows to the Josephson dipole; none of it is distributed to any other inductor. As we show in the following, a larger participation leads to larger quantum vacuum fluctuations.
Transmon-qubit example. Consider the transmon qubit coupled to a readout cavity mode (see Methods). The Josephson junction has an EPR of near unity in the qubit eigenmode, since nearly all of the inductance of the transmon is due to the kinetic inductance of the tunnel junction. On the other hand, the EPR of the junction in the cavity mode is on the order of 10 −2 . The junction kinetic inductance contributes only on the order of one percent to the total inductance of the cavity, associated with the large cavity box. This leads to the large ZPF of the junction flux in the transmon mode, ZPF ∼ 1, and the therefore much smaller fluctuations of the junction in the cavity mode ZPF ∼ 10 −1 (see also Sec. A7).
Energy offset. The EPR is defined in terms of the energy excitation of a mode, rather than in terms of its absolute energy.
Vacuum energy. Mathematically, the vacuum energy is a consequence of the non-commutativity ofˆandˆ †, which in effect adds a constant offset to the energy of each mode. In calculating the EPR, we take the reference energy of an element relative to its vacuum energy.
Equipartition theorem. Classically, for a linear circuit, the energy of an eigenmode oscillates in time between inductive and capacitive energy. At periodic intervals given by / , the inductive energy of the mode is exactly zero. For this reason, we use the time-average energy of the Josephson dipole and the eigenmode. We recall that the time-averaged energy is half of the peak energy. This can be exploited in expressing the total inductive energy as half the total mode energy, We will discuss what states are used in calculating the expectation values below. The ground state energy is taken to be zero.
Calculating the EPR. These ideas lead us to the following definition of the EPR in the quantum setting, where the over-line denotes a time average and the expectation value is taken over a state with an excitation only in mode . As discussed above, all energies are referenced to their ground-state expectation values, to disregard vacuum-energy contributions (for special notation to accommodate this, see remark on Wick ordering at the end of this section). Using Eqs.  The EPR is independent of the excitation amplitude . EPR for a coherent state excitation. A coherent state excitation of mode is denoted | ˆ( ,ˆ) |vac , whereˆis the displacement operator of the th mode, ( ,ˆ) exp ˆ † − * ˆ , and is a non-zero complex number. The coherent state time-evolution is a rotation, | ( ) = (0) − . Using Eq. (A.51), we find the EPR to again be excitation independent and exactly equal to that given in Eq. (A. 53).
Quantum fluctuations in terms of EPR. The EPR for a Fock or coherent state excitation is given by Eq. (A.53). Inverting this expression, we find the ZPF in terms of the EPR, where is the EPR sign, ∈ {−1, +1}. In practice, the sign is calculated using Eq. (C.8). Algebraically, the sign can be found from the eigenvector matrix E. If the th row of the eigenvector matrix corresponds to the th spanning-tree branch flux and the th column corresponds to the th eigenmode, then = sign [E] , where [E] is the ( , )th entry of the matrix.
EPR sign: sign freedom. The value of an individual sign is completely arbitrary. It does not have measurable consequences in the same way that the amplitude of a standing mode can be taken to be positive or negative, indicating a -phase shift. The sign derives its physical meaning in relationship to other signs of the same mode but different Josephson dipoles (see Supplementary Figure S8). If the signs of two dipoles are the same (resp., different), then the two dipoles oscillate in-phase (resp., out-of-phase). The sign of the EPR can flipped by either flipping the reference direction of the junction or the overall phase of the mode.
EPR for circuit design. The EPR in Eq. (A.54) is essentially the only free parameter in engineering the quantum ZPF and the amplitudes of the non-linear couplings. It is readily calculated from the classical eigenmode FE simulation of the distributed physical layout of the circuit, as detailed in Sec. C. In this way, the EPR serves as a bridge between the linearized classical circuits treated with FE solvers and the Josephson circuit in the quantum domain.
Remark: Equivalent formulation using Wick normal-ordered expectation value for energies. Mathematically enforcing the energy to be referenced to that of the vacuum state can be equivalently accomplished by disregarding the commutation relationships in expectation values of the energy. The compact notation that accomplishes this was introduced in quantum optics and is also used in quantum field theory: for an operatorˆ, one takes the (Wick) normal-ordered form of the operator 136 , denoted by a colon on either side of the operator, : :. Thus, :ˆˆ † : =ˆ †ˆand : ˆ †ˆ 2 : =ˆ † 2ˆ2 . The normal-ordered form of the operators should not be confused with the normal-ordering transformation applied to the operator, which takes the commutation relations into account; e.g., In the presence of non-linear interactions, the vacuum energy leads to observable effects, such as the Lamb shift, introduced in Eq. (7) of the main text. The definition given in Eq. (A.50) can be equivalently stated using Wick notation: (A.56)

A7. Universal EPR properties
The energy-participation ratios obey four universal properties. These properties are valid regardless of the circuit topology and nature of the Josephson dipoles. They follow directly from the normal-mode structure of the eigenmodes of lin , obtained in Sec. A4 and are linked to the ZPF, as discussed in Sec. A6. In the following, denotes the total number of modes.
The EPR is bounded. The EPR is an energy fraction. It follows from its definition in Eq. (A.50) that it is a real number comprised between zero and one-since the Josephson dipole energy is always positive and equal to or smaller than the total inductive energy of the mode, The total EPR of a dipole follows a sum rule-it is conserved while being diluted among the modes. The total EPR of a dipole is conserved-it is exactly unity across all modes, which follows from Eqs. (A.38) and (A.51). Increasing the EPR of a Josephson dipole in one mode will proportionally reduce its participation across other modes. The EPR is neither created, nor destroyed. It is distributed among the modes. Adding additional modes to the circuit or removing modes from the circuit does not increase or decrease the total EPR of a dipole.
The total EPR of a mode is at most unity. A single mode can at most have a total EPR of unity, and no where is the EPR sign. Remark: Use of these four properties in quantum circuit design. These universal EPR constraints are useful in the design of quantum circuits, especially for designs that require weak and strong nonlinear interactions simultaneously. For example, see Methods, where we discuss the impossibility of a design we targeted due to the EPR constraints. We subsequently use the constraints to obtain a best approximation of the design, and to gain insight into the range of possible non-linear couplings. In our experience, the EPR has allowed us to thus circumvent the need to run a prohibitively expensive finite-element sweep to explore the many possible design-parameter regimes.
A8. The biased Josephson system: equilibrium state in the presence of persistent currents The magnetic flux across a Josephson dipole Φ , introduced in Sec. A1, is defined as a deviation away from the equilibrium configuration of the Josephson system. A system that incorporates an active element-one that can act as a voltage or current source-or a system subjected to a frustrated constraint-such as found in a conducting loop frustrated by an external magnetic flux-can have a non-zero equilibrium state. Referencing the Josephson dipole and spanning-tree fluxes Φ t from equilibrium guarantees that the total inductive energy E ind (Φ t ) is at a minimum for Φ t = 0.
Static-equilibrium conditions. In static equilibrium, the net generalized forces (currents) at each node in the system vanish and, consequently, the generalized acceleration vanishes, d 2 Φ t d 2 = 0. The corresponding magnetic flux of the spanning-tree branches in static equilibrium Φ t,eq can be found by extremizing the Lagrangian with respect to the generalized position variables 129 , Solving these conditions for the equilibrium flux Φ t,eq amounts to solving for the direct-current (dc) operating point of the circuit-a standard, classical problem; although one that is non-local in nature and, in the presence of Josephson dipoles, one that involves solutions to non-linear equations. In general, classical numerical methods can be used to find the equilibrium of the circuit, especially when the equilibrium equations are transcendental. However, for many practical situations, Eq. (A.61) simplifies or need not be evaluated at all, as discussed in Sec. A9. Equilibrium flux of a Josephson Dipole in isolation vs. in a circuit. In general, a Josephson dipole has a different equilibrium flux across its terminals when in isolation vs. when embedded in a system. When in isolation, the native equilibrium flux of the Josephson dipole is found from the simple local condition E Φ = 0, where E is the energy function of the Josephson dipole and Φ is the magnetic flux across the dipole. When embedded in a system, the equilibrium flux of the Josephson dipole, found from Eq. (A.61), is in general not a local property of the dipole anymore-but one of the system, as illustrated by the following example.
Example: Josephson junction in a frustrated ring. Consider a Josephson tunnel junction in isolation-not connected to any other elements. It's energy function, see Eq. (A.5), is E (Φ ) = − cos (Φ / 0 ), where Φ is the total magnetic flux across the junction in isolation. (Here, the subscript denotes the junction in isolation.) The energy is minimized at the native equilibrium-flux value Φ nat eq = 0. Now, consider embedding this junction in a superconducting ring frustrated by an external magnetic flux Φ ext , as depicted in Supplementary Figure S1(a). Suppose the geometric ring inductance is , see Supplementary Figure S1(b). The equilibrium condition of the circuit, given by Eq. (A.61), reduces to Kepler's transcendental equation where the effective dc current sourced to the loop is eff = −1 Φ ext . In the additional presence of a voltage source and current source as depicted in see Supplementary Figure S1(b), the net effective current is eff = −1 Φ ext + . Without applied forces on the spring, the native equilibrium length of the spring is 0 eq . Stretching or compressing the spring from equilibrium is measured by the deviation away from the native equilibrium. The spring constant is . (b) Depiction of two connected springs, with spring constant and , constrained between two walls separated by a fixed distance ext . The equilibrium lengths of the two springs in the system are eq and eq . These differ from their equilibrium lengths in isolation; i.e., eq ≠ eq . Deviations from the system equilibrium are denoted by .
A solution Φ eq (Φ ext , Φ , ) of Eq. (A.62) can be obtained numerically or using Lagrange inversion. In the EPR method, the canonical flux deviation away from the selected equilibrium Φ eq is then The energy function of the junction now in terms of the deviation Φ is Using the series expansion of Eq. (A.64), given in Eq. (A.7), the Josephson dipole energy of the junction, as defined in Eq. (A.13), is (Φ ext , Φ , ) cos Φ eq / 0 and the leading-order non-linear coefficient is 3 (Φ ext , Φ , ) = − 1 6 tan Φ eq / 0 . Mechanical analogy: spring in isolation. By way of analogy, a Josephson dipole is like a mechanical spring. Imagine a spring in isolation, see Supplementary Figure S2(a). No external forces act on it. Its native equilibrium length (when the spring is in isolation) is 0 eq . Stretching the spring to length results in the force = − ( ) on its ends, where − 0 eq denotes a deviation away from equilibrium. In general, ( ) is a nonlinear function; however, in equilibrium, (0) = 0. For a linear spring, the energy, E ( ) = 1 2 2 , is minimized for = 0.
Mechanical analogy: spring in a frustrated system. Imagine embedding the spring in a simple system as depicted in Supplementary Figure S2(b). The system imposes the KVL-like constraint condition + = ext , where is the length of the second spring, and ext is the total distance between the two walls. For linear springs, the length of the spring in the equilibrium of the system is eq = eq + ext − eq /( + ), where and eq are the spring constant and native equilibrium length of the second spring, respectively. The energy of the spring is E ( ) = 1 2 2 , where we omit constant energy terms and use − eq to denote deviations away from the equilibrium length of the spring in the system (not the spring in isolation).
Remark: voltage source and equilibrium. The voltage bias does not appear in the effective current bias of the junction. The capacitor isolates it in dc from the rest of the circuit and prevents it from establishing a dc current in any loop. In Sec. A9, we discuss further situation in which a bias has no effect on the equilibrium.
Remark: voltage source and dc loops. A dc voltage source is absent from the dc-conducting superconducting loop formed by the junction and inductor, since its magnetic flux grows linearly in time, Φ = , and it would hence supply a current to the loop that grows infinitely large. At some time, this would break the abstraction of the ideal voltage source or quench the superconductivity.

A9. The biased Josephson system: simplifications
In Sec. A8, we reviewed the equilibrium conditions of the Josephson system. Here, we develop several common situations in which the analysis of these conditions simplifies from a global to a local one or one that need not be done at all. In the region of the gap, we define the segment of the linecontour for the magnetic flux Φ ext to be the a minimal-length one across the gap-i.e., in terms of electrical schematics, we take the flux across the effective capacitance associated with the gap in the loop. (e) Example similar to (d) but with a SNAIL element instead of tunnel junction.
A dc-conducting loop incorporating a Josephson dipole with no internal loops. If a simple Josephson dipole with no loops internal to it (such as the Josephson tunnel junction; not a SQUID or SNAIL), is embedded in a dc-conducting distributed and frustrated loop, see Supplementary Figure S3(b), then the flux of the dipole in the equilibrium of the total circuit cannot be determined locally. Rather, calculating the equilibrium requires incorporating information about the geometric-inductance response of the distributed loop at dc. Classical simulation methods can be used to obtain an equilibrium point.
A dc-conducting loop incorporating a Josephson dipole with internal loops. In the most general case, a Josephson dipole is part of a dc-conducting loop; see Supplementary Figure S3(c). If the loop is not frustrated, then the equilibrium of flux of the dipole can be taken as the native equilibrium flux of the dipole in isolation. However, in practice, even if no frustration is intended, there may be spurious magnetic flux that will frustrate the loop. If the loop is frustrated by sources (which can include other Josephson dipoles able to source current, such as a biased SNAIL), the equilibrium condition of the global loop should be calculated accounting for the geometric inductances.
Voltage offset. A voltage offset will have also have a frustrating effect on the eigenspectrum of the Hamiltonian. However, for the purposes of calculating the potential energy minima for use in the EPR linearization, it is sufficient to consider only the frustration due to currents.
Multiple equilibrium states. The equilibrium conditions can yield multiple solutions that locally minimize the energy function E ind . It is preferable to choose a lowest-energy, stable equilibrium state for the operating point of the circuit. In principle, any of the minima can be used in the EPR method as the operational point; the physics of all other minima can be reconstructed from the operating point. However, we note that the presence of multiple minima can lead to phase slips 72,77,[137][138][139][140] .

B. Nonlinear interactions, effective Hamiltonians, and the EPR
In this section, we find the amplitude of a general nonlinear interaction due toˆn l explicitly. In Sec. B1, we find the effective interaction Hamiltonian of weakly nonlinear systems in the dispersive regime to leading-order. In Sec. B2, to more systematically handle large systems, we introduce the EPR matrix P and use it find the Kerr matrix and the vectors of anharmonicities and Lamb shifts. In Sec. B3, we find the general, normal-ordered form ofˆn l and an expression to calculate the amplitude of any many-body interaction contained within. In Sec. B4, we use these results to describe the parametric activation of nonlinearities in a pumped Josephson circuit.

B1. Excitation-number-conserving interactions of weakly-nonlinear systems
Josephson circuits that are weakly non-linear and in the dispersive regime, such as the transmon-cavity and transmon-transmon circuits, have, in the absence of drives, dominant non-linear interactions that conserve excitation numbers.
Perturbative and dispersive. For such weakly nonlinear systems with small zero-point fluctuations 1, for all modes and junctions j, the Hamiltonianˆn l exerts only a perturbative effect on the eigenspectrum ofˆl in . Hence, we treat its effect order-by-order using the expansionˆn l = =1 ∞ =3ˆ, obtained in Eq. (A. 43). For this approach, we focus on the regime in which the energy difference between two modes is much larger than the phase excitation of the Josephson dipoles, ℏ ( − ) ˆ for all , , , and . Example: transmon coupled to a readout cavity mode. Recall the simple example circuit quantized at the beginning of the main text. A qubit ( ) is coupled to a cavity ( ). To leading order in , the circuit Hamiltonianˆn l is, using Eq. (A.43), where andˆare the Josephson energy and flux operator of the tunnel junction, and and are the qubit and cavity mode ZPF, respectively.
Expanding the multinomial. Expanding the = 4 multinomial term in Eq. (B.1), we find a weighted sum of all possible four-body interactions between the qubit and cavity; example terms includeˆ ˆ † Higher-order nonlinearity yields lower-order coupling. Due to the non-commutativity of the operators, the normal-ordering procedure results in terms that have lower polynomial order than those in original unordered expression; these new terms includeˆ †ˆand 3ˆ. Such quadratic terms (e.g.,ˆ †ˆ,ˆ †ˆ, andˆ †ˆ) dress the modes ofˆl in in a linear manner-both renormalizing their frequencies and hybridizing them. These new linear couplings cannot be a-priori straightforwardly included inˆl in , since they occur only as a non-classical consequence of the ZPF in the non-linearity; they do not appear in a classical treatment of the Josephson system, which lacks the non-commuting operator aspect. Hence, their amplitudes are determined by the quantum ZPF of the eigenmodes. Remark: A self-consistent approach to include these linear couplings in lin is possible under certain assumptions 17 .
Rotating-wave approximation in the dispersive regime. To leading-order, under the rotating-wave approximation (RWA), only excitation-number conserving interactions inˆn l (those with an equal number of raising and lowering operators in each mode) contribute. Thus, Eq. (B.1) reduces to the effective, dispersive qubit-cavity Hamilto-nianˆ= where ∈ { , } is the mode label, Δ is the effective Lamb shift, is the anharmonicity, and is the total cross-Kerr frequency shift induced between modes and . The bar over denotes the RWA; the subscript = 4 denotes the highest the power of the Taylor expansion included in the effective Hamiltonian. The excitation spectrum of this Hamiltonian is illustrated in Supplementary Figure S4. Remark: Using first-order time-independent perturbation theory in place of the RWA yields the same result as the RWA.
Hamiltonian parameters. The Hamiltonian parameters are obtained from the normal ordering and using Eq. (A.54), Constraints in the Hamiltonian. First, the structure of the non-linearity imposes certain constraints. For example, ≤ 2 √ , where the equality holds only for the one-junction case, = 1. Second, the universal EPR properties impose a strict limit on physically realizable designs, subject to Eqs. (A.57)-(A.60).
Degeneracies. So far, we implicitly assumed the system does not have accidental degeneracies, such as those that occur when the frequency of a mode is an integer multiple of that of another, ≈ × , where ∈ Z. In the case of such a degeneracy, additional terms survive the RWA; for example, if = 3 , then the nonexcitation conserving termˆ3ˆ † survives the RWA.
EPR sign & example. Due to the even-power nature of the nonlinear interactions, the EPR sign drops out altogether from Eq. (B.3). The interaction strength between the qubit and cavity modes only depends on the overlap int the EPRs of the two modes alone. The significance of this is curiously showcased by Device DT3, presented in Methods. The qubit eigenmodes of DT3 are the equally-hybridized symmetric (bright, B) and antisymmetric (dark, D) combinations of the two bare transmon qubit modes. Naively, using the analogy of atoms in a cavity, one reasons that the symmetric mode couples to the cavity (C), | BC | 0, and the antisymmetric modes does not, | DC | = 0, due to the out-of-phase oscillation of the junction current dipole-a cancelation effect. However, from the EPR expression, it is seen that, to the contrary, the signs are irrelevant-no such cancelation is possible. Rather, according to the EPR method, since both modes have equal participation in the two junctions ( 1 = 2 = 1 = 2 ), the couplings to the cavity are essentially equal, DC ≈ BC ; indeed, this is observed in the experiment (see the corresponding data table in the Methods section).
Generalizing to many modes. Equations (B.2) and (B.3) were derived for the example of a qubit-cavity circuit; however, they generalize straightforwardly to modes under the simple extension ∈ {1, . . . , }.

B2. EPR matrices and many-body interactions
To more easily and systematically handle large-scale circuits and higher-order non-linearities, we introduce the EPR matrix P, comprising the × EPRs as its entries, where Ω is the diagonal eigenfrequency matrix and where 1 denotes a column vector of length with all entries equal to unity.
Dilution of the nonlinearity. The dilution of the nonlinearity of the Josephson dipole elements among the eigenmodes is neatly expressed in Eq. (B.6). The Josephson dipole energies E J −1 are diluted by ΩP through a congruence transform to the Kerr coefficients. While the eigenfrequencies Ω weighs the contribution to each mode, the participation matrix P dictates the dilution of the junction energies and their nonlinearity among the modes.
Higher-order nonlinear corrections and dilution. Using the results of Sec. B3, the -th other correction to the Kerr matrix is In general, the Lamb shift correction depends on the total ZPF frustration of the junctions.

B3. General many-body interactions
So far, we explicitly calculated the leading-order correction on the spectrum ofˆl in due toˆn l . We expressed the eigenmode interactions using normal-ordered many-body terms, such asˆ †ˆ †ˆˆ. Here, we extend the analysis and compute the amplitude of any general normal-ordered term inˆn l .
The form a general many-body interaction. A general normal-ordered many-body interaction has the form α,βˆ † 1 1 · · ·ˆ †ˆ1 1 · · ·ˆ Importantly, the commutator ˆ,ˆ † = =1 2 is scalar-valued, which allows us to use the normal-ordering non-commutative binomial theorem to expand the -th power term ofˆn l 141-143 . Using the non-commuting expansion, where floor 2 gives the greatest integer that is less than or equal to 2 and the total variance of the ZPF of the -th dipole is Sinceˆis the sum of operators that commute, see Eq. (B.15), we can now expand the powers ofˆusing the classical multinomial theorem, where the multi-index shorthand α =1 , the multinomial coefficient is 20) and the sum condition |α| = means that the sum includes terms with all possible tuples α such that their 1-norm has the value .
The general normal-ordered form. Combining Eqs. (B.16)-(B.19), we find the general normal-ordered many-body form ofˆn l to all orders in and without approximations, The energy-dimensioned amplitude of an interaction due to the -th power ofˆn l for a general many body modeinteraction term is the sum of the individual junction contributions, The Josephson system can be subjected to strong external drives used to parametrically activate or enhance nonlinear couplings. We illustrate the use of Equation (B.22) in this context by calculating the amplitude of a excitation-swapping beam-splitter interaction.
Example: parametrically-activated beam-splitter interaction. Motivated by the setup of device WG1, we aim to parametrically activate a beam-splitter (BS) interaction between two non-resonant modes of a Josephson system; for example, such an interaction can be used as a Q-switch 96,145 . The system has a high-quality, storagecavity mode ( ) and another low-quality, readout-cavity mode ( ). The two modes are far detuned and obey the conditions outlined at the start of Sec. B1. From the sys-temˆn l , we aim to obtain the effective BS Hamiltonian Parametric activation and interaction rate. A third mode of the system, indexed by , is off-resonantly driven at frequencies and with amplitudes 1 and 2 , respectively; within the RWA, the drive Hamiltonian is General parametrically-activated interaction. The procedure illustrated with the BS example generalizes straightforwardly to the parametric activation of most other interactions and to finding their rates using Eq. (B.22). The procedure is useful for more complex and even cascaded processes 97,144 used in dissipation engineering 147 . It is reported that the limit of procedure is typically reached when approaches unity, and other activated nonlinear processes become non-negligible 148-151 .

C. Finite-element electromagnetic-analysis methodology
We detail the EPR methodology for the finite-element (FE) analysis of the Josephson circuit. In Sec. C1, we model a Josephson dipole in the FE simulation as a rectangular sheet with an inductive lumped-element boundary condition with inductance . In Secs. C2 and C3, we extract the EPR and EPR signs from the result of an eigenanalysis simulation of the linearized Josephson circuit, corresponding toˆl in . This step completes the classical analysis part of the EPR method; from here, the Hamiltonian is fully specified, as described in Secs. A and B. The eigenresult also provides complete information on the dissipation and input-output coupling of the circuit, as described in Sec. D.
These steps, and the calculations detailed in this text, are automated by the open-source software package pyEPR 152 , which we offer to the community.

C1. Modeling the Josephson dipole
In the FE model of the Josephson circuit, we model a Josephson dipole as a simple rectangular sheet with a lumped-element boundary condition 4 , see Supplementary Figure S6. The sheet abstracts away the physical layout of the Josephson dipole and its wiring leads.
Idealization of the deep-sub-wavelength features. This idealization of the Josephson dipole is justified in that its physical size is in the deep-sub-wavelength regime with respect to the eigenmodes of interest. For example, for the devices described in the Methods section, the separation between the Josephson dipole size and the mode wavelengths of interest is approximately five orders of magnitude. We hence treat a Josephson dipole in the FE model as lumped-element inductor with inductance , given by Eq. (A.13a); the linearization is taken with respect to the circuit equilibrium, see Sec. A8.
Electromagnetic model of the lumped inductance. The -th Josephson dipole is modeled as a twodimensional sheet , see Supplementary Figure S6. The sheet is assigned a surface-impedance boundary condition, which imposes ì = (ˆ× ì ) across the sheet, where ì and ì are the tangential electric and magnetic fields of the sheet, respectively,ˆis the sheet normal vector, and is the complex-valued surface impedance corresponding to a total sheet inductance of . The hat symbol over denotes a unit vector in the context of electromagnetic fields; it is not to be confused with the hat notation used for quantum operators.
Reducing the model complexity: ignoring leads. If the geometric inductance of the wire leads connecting a Josephson dipole to larger distributed surfaces (such as the pads of a transmon qubit) is negligible in comparison to and the wire leads are deeply sub-wavelength in size, then the wire leads can be ignored, as depicted in Supplementary Figure S6. If the wire leads have nonnegligible inductance, then they can either be included in the simulation or their linear inductance can be included in the definition of the Josephson dipole. In practice, it is preferable to be able to abstract away the wire leads as much as possible, since their feature sizes are typically orders of magnitude smaller than other all other design features. The inclusion of such fine detail in the model can result in very large geometric aspect ratios, which in turn result in increased computational costs. However, while this finer detail is more representative of the physical design, we have found that it does not lead to noticeable corrections for typical cQED devices.
Performance tip: mesh operations. The sheet is generally one of the smallest features of the FE model. Due to this small size but critical role, one can gently seed a higher level of mesh on to speed up the eigenanalysis. However, caution should be used to avoid seeding too heavy of an initial mesh, which can instead lead to poor convergence of the simulation. Convergence can be verified by plotting as a function of the simulation adaptive pass number and by verifying the conditions de- If a Josephson circuit incorporates exactly one Josephson dipole, then we use the global electric and magnetic eigenmode energies to directly calculate the EPR of the dipole in the mode.
Energy balance. The time-averaged electromagnetic energy in a resonantly excited mode is equally split into an inductive E ind and capacitive E cap contribution 125 . This detailed balance, E ind = E cap , holds even in the presence of dissipation and defines the eigenmode condition. In the presence of a Josephson dipole, the inductive energy is split into a magnetic E mag and a kinetic E kin contribution; E ind = E mag + E kin . The magnetic contribution is associated with magnetic fields and geometric inductance. The kinetic contribution is associated with the Josephson dipole kinetic inductance and the flow of electrons, and their inertia. From the point of view of the FE analysis, the magnetic energy is stored in the magnetic eigenfields ì and the kinetic energy is stored in the lumped-element boundary condition on . If lumpedelement capacitive boundary conditions are absent from the model, then the capacitive eigenmode energy is stored entirely in electric eigenfields ì , E cap = E elec ; hence, Calculating EPR from the energy balance. Using Eq. (C.1) and Eq. (5) of the main text, = E kin /E ind , the EPR is calculated from the ratio of global energy quantities, where the total magnetic-and electric-field energies are computed from the eigenfields phasors, where ì max ( , , ) [resp., ì max ( , , )] is the eigenmode electric (resp., magnetic) phasor, and ← → (resp., ← → ) denotes the electric-permittivity (resp., magnetic-permeability) tensor. The spatial integrals are performed over total volume of the device. The electric eigenfield is related to the phasor by C3. Calculating the EPR in the case of multiple Josephson dipoles

EPR
of the -th Josephson dipole. In the case of multiple Josephson dipoles, the total kinetic energy E kin is itself split among the dipoles, see Eq. (A.20). We calculate the EPR of junction in mode using the EPR definition, Eq. (A.50), where is the peak value of the Josephson dipole current in mode . The current is calculated from the integral of the mode surface-current density ì , ( , , ) over the dipole sheet , where is the length of the sheet, see Supplementary Figure S7.
EPR sign. The EPR sign is calculated from the orientation of the current . The absolute orientation is relative, as described in the main text, and is determined by defining a directed line DL across . If the phasor ì , is aligned with DL then = +1; otherwise, = −1. Hence, we extract the sign using The direction of the line merely establishes a convention for a positive EPR sign.
Enforcing energy balance. The convergence of , a quantity extracted from the local eigenfield solutions, can be enhanced by re-normalizing the set of mode EPR to ensure energy balance, see Eq. (C.1). If lumped-element capacitive boundary conditions are absent, then the EPR can be renormalized to ensure that their total sum is equal to the ratio =1 = E kin /E ind , which is a globally calculated quantity, and is therefore expected to converge quicker.
The above calculations are automated by pyEPR 95 .

C4. Remarks on the finite-element eigenmode approach
Finding the eigenfrequencies. Rather than searching for the location of an unknown pole in an impedance response of the circuit and then honing in on it to perform finer sweeps, the eigenmode analysis returns the lowest modes above a minimum frequency of interest. This typically lifts the requirement for a prior knowledge of the mode frequencies.
Single simulation. The FE eigenmode performs a single simulation from which complete information of the Josephson circuit is extracted. This is in contrast to the typical process flow used in an impedance analysis, which requires that mode frequencies be first identified so that a series of individual narrow-frequency-range impedanceresponse sweeps (one for each mode and junction) can be performed.
Closed-circuit optimization. We have found that the above two feature (see remarks) can speed up the iterative refinement of a quantum design and can circumvents difficulties associated with finding and fitting unknown, narrow-line poles in an impedance analysis.

D. Dissipation budget and input-output coupling
In this section, we summarize the methodology used to fully characterize dissipation and input-output coupling in the Josephson system. Loss of energy results from material losses and radiative boundaries which guide energy away from the system. Additionally, control of the system is achieved by means of the radiative boundaries, such as input-output (I-O) coupling. The dissipation budget, comprising the individual losscontribution bound of each lossy and radiative element, is extracted from the same eigensolution used to calculate the EPR , in essentially the same way-by calculating the fraction of the -th mode energy stored in the -th lossy element-the lossy EPR . While lossy EPR signs can also be calculated for each element, these are not needed for linear dissipation; however, their role in I-O coupling is detailed in Sec. D4.

D1. Dissipation budget
The dissipation budget comprises the estimated energy-loss-rate contributions due to each loss mechanism and each lossy object in the Josephson circuit 85,88,91,92,113,125,[153][154][155][156] . We classily losses as capacitive, inductive, or radiative, and summarize their calculation here. Each loss mechanism is described by a corresponding lossy EPR and an intrinsic quality factor . Assuming the dissipation is linear and the mode of interest is underdamped, the loss rates of each mechanism simply add; so, the total quality factor of the -th mode is 92,113,154 where cap and ind are the total mode quality factors due to capacitive and inductive losses (i.e., those proportional to the intensity of the electric ì 2 and magnetic field ì 2 , respectively, see Secs. D2 and D3), and rad is the total radiative mode quality factor, see Sec. D4. In this section, the mode index is implicit.
Remark. Beyond these intrinsic mechanisms, extrinsic factors, such as ionizing radiation, can play a significant role in understanding loss mechanisms, such as those due to quasiparticle in superconducting quantum circuits 157,158 .

D2. Capacitive loss
The total capacitive mode quality cap is the weighted sum of intrinsic quality factors cap (i.e., the inverse of the dielectric loss tangent) of all lossy dielectrics in the Josephson circuit 92,154 , where cap is the lossy energy-participation ratio of theth dielectric in the mode-i.e., cap is the fraction of capacitive energy stored in the dielectric element . We classify lossy capacitive elements as either bulk or surface. For example, the volume of three-dimension dielectric object, such as a chip substrate, is associated with bulk capacitive loss 83,94,159 . On the other hand, the surface of the substrate, which may be a surface-dielectric layer is classified as a surface-loss element 88,94 . The lossy EPR for bulk capacitive loss is calculated from the eigenfield solutions, where the integral is carried over the volume of theth bulk dielectric object; the total electric energy E elec is defined in Eq. (C.3). The lossy EPR for a surface dielectric is approximated by where the surface-dielectric layer thickness is and its permittivity is .

D3. Inductive loss
Physically, inductive losses originate from the dissipative flow of electrical current in metals or through metalmetal seams. Supercurrent loss due to quasiparticles and vortices can be accounted for in an effective quality factor of the conducting surface. We denote the intrinsic inductive quality factor of a lossy object ind . The bound on the total inductive-loss quality factor ind of mode is a weighted sum of ind , where ind is the lossy EPR for inductive element . We classify inductive losses as either surface, bulk, or seam. Surface conductive loss (in the skin-depth). The fraction of eigenmode energy stored in the skin depth 0 of metal surface , denoted surf , if the lossy EPR of the surface, and is obtained from the eigenfield solutions, where is the magnetic permeability of the surface; typically, = 0 , and ì max, is the magnetic field phasor parallel to the surface. For superconductors, ind,surf is the kinetic inductance fraction 91,154 , commonly denoted . The total magnetic energy E mag is defined in Eq. (C.4).
Surface conductive loss: intrinsic quality. In the normal state, a metal, such as copper or aluminum, has an intrinsic inductive quality factor of order unity 125 , ind,surf ≈ 1. However, in the superconducting state, the lower bound on the quality factor is typically found to be in the range of several thousand 160 , ind,surf 10 3 . The lower bound for thinfilm superconducting aluminum has been measured to exceed ind,surf > 10 553 . Bulk magnetic loss. The mode magnetic field can couple to bulk magnetic impurities present in the volume of lossy object . The lossy EPR for the bulk-inductiveloss mechanism of volume is This coupling is typically negligible in current superconducting quantum circuits. Seam loss. The seam formed by pressing two metals together provides an electrical bridge for the current to flow across but also introduces a key dissipation mechanism 90 . A common example of a seam is the one formed by the two halves of the metal sample holders used to house a cQED chip. In the FE analysis, we model the seam by a line path seam tracing out the location of the seam at the surface of the two mating walls. The inductive lossy EPR participation for the seam is where the seam thickness is denoted , its magnetic permeability , and its the penetration depth 0 . It is convenient to rewrite the total seam loss due to seam in terms of an effective seam admittance seam , defined in Ref. 90, (D.9)

D4. Radiative loss and input-output coupling
The Josephson circuit incorporates radiative boundaries, typically purposefully introduced to serve as inputoutput ports of the circuit 125 , thus providing a means to perform system measurement and control. For cQED devices, the port exposes the circuit to an external transmission line conduit, such as a coaxial cable or a coplanar waveguide. The total radiative quality factor rad â j Φ a in,1âin,2 a out,2 a out,1 (a) Figure S8. Schematic of a transmonqubit circuit (mode operatorˆ), comprising a Josephson tunnel junction (flux operatorΦ ) and a capacitor, coupled to two input-output ports. The input and output fields of the left (resp., right) transmission line areˆi n,1 andˆo ut,1 (resp.,ˆi n,2 andˆo ut,2 ). (a) Ports 1 and 2 both coupled to the top qubit node; hence, both port EPR signs are equal, 1 = 2 . (b) Port 1 couples to the top node, but port 2 couples to the bottom qubit node. The port EPR signs have opposite signs, 1 = − 2 .
of mode is the sum of the individual port contributions −1 rad = =1 −1 , where is the total number of ports and is the quality factor due to port . Below, we describe how port is modeled in the FE simulation to extract from the eigensolutions. Radiative energy loss. Energy stored in mode can leak at a rate through port and be guided away by the transmission line. While for certain modes, such as readout ones, this coupling is desired, for other modes, such as a qubit one, this coupling is often considered spurious. In the case of a qubit mode, the energy loss to a readout port is seen as a manifestation of the Purcell effect 123 ; while, for the readout mode of the same structure, the energy loss to the readout port sets the rate of information gain 161 . The rate is calculated from port EPR , as detailed in the following for both wanted and spurious terms. The port EPR sign is calculated concurrently and is important for the system drive configuration.
FE model of the port. In the presence of a port, the boundary of the Josephson circuit is somewhat ambiguous-we can include more or less of the port and conduit structure in the model. We choose to include a minimal but sufficiently large portion of these to faithfully model the disturbing effect of the boundary condition on the eigenmodes. For example, in the case of the qubit-cavity structure of Figure 1(a) of the main text, we include a short stub of the I-O coaxial cable in the FE model. The length of the stub can be determined from the effect of a sweep of its length on the target parameters. We have found that an alternative heuristic measure is to use the decay of the eigenfields inside the port structure and to make sure that the end of the port structure is at least several exponential decay lengths long (the field in the port structure decays exponentially since the eigenmodes are generally below its cutoff). The port structure termination surface is treated as a resistive sheet with to model the effect of the transmission line guiding waves away from the Josephson circuit. The sheet effective resistance is . In the case of a 50 Ω port line, = 50 Ω. While a resistive boundary condition can be assigned to the sheet, in the case of high-quality modes, it is possible to use a perturbative approach and to perform a lossless simulation of the FE model, from which the loss can be extracted (see also remark at end of this section).
Calculating the input-output (I-O) coupling rate from the port EPR . From the lossless eigensolutions of mode , the energy lost to the effective resistor of port during one mode oscillation period where is the peak current across the port due to the excitation of mode ; see Eq. (C.7). The total mode energy E ( ) at time decays at rate Hence, assuming a high quality mode, the energy loss during one oscillation period, between times = 0 and , is 12) This shows that the loss to port during one period is = E (0), which we equate to the expression of Eq. (D.10) to find the I-O coupling rate in terms of quantities calculated from the eigenfields, . (D.14) Sign of the I-O participation. If there are multiple ports, the port EPR sign , calculated using Eq. (C.8), is important in a manner similar to that explained for the case of the EPR sign in the case of multiple Josephson dipoles; see the main text. To illustrate, consider a simple transmon-qubit circuit coupled to two transmission lines in two different ways as depicted in the two panels of Supplementary Figure S8. While in both configuration the eigenmode frequency and quality factor is identical, the configurations are inequivalent. Consider driving both transmission lines with the same amplitude and in-phase, then the circuit of Supplementary Figure S8(a) is excited but that of Supplementary Figure S8(b) is not.
Remark on modeling the port termination as resistive vs. lossless. The port sheet can be treated as either resistive or lossless. In the former, the sheet is assigned a lumped-element boundary condition with impedance matching the port input impedance as seen from the system; typically designed to be = 50 Ω 10,162 . The eigenresults fully account for the effect of the dissipation on the mode profile. These effects are negligible in the case of high-quality modes, but become significant as approaches unity. For lossless treatment of , suitable when 1, the termination simply assigned a perfectly conducting boundary condition. In both treatment, the eigenmode fields can be used to calculate the loss due to . Reactive ports. In addition to being resistive, ports can have a reactive component, typically occurring in the presence of non-idealities. For example, a port coupled to transmission line suffering from a down-line reflection will have some of the energy it leaks out to the line come back to it 163 . The reactive part of the port structure thus houses energy in its internal modes. There can be treated by modifying the boundary condition on or for a more general treatment can be accounted for by including the port, line, and scatterer structure in the FE model. This larger model will account for hybridization between the line modes and those of the system 100 .