Energy-participation quantization of Josephson circuits

Abstract

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.

I INTRODUCTION

Quantum information processing based on the control of microwave electromagnetic fields in Josephson circuits is a promising platform for both fundamental physics experiments and emerging quantum technologies (Devoret2013, ; GoogleSupremacy2019, ; Blais2020, ). Key to the success of this platform is the ability to quantitatively model the distributed quantized electromagnetic modes of the system, their nonlinear interactions, and their dissipation (see Fig. 1). This challenge is the subject of intensifying interest (Nigg2012, ; Bourassa2012, ; Solgun2014, ; Solgun2015, ; Smith2016, ; Gely2017, ; Malekakhlagh2017-Cutoff-Free, ; Pechal2017, ; Parra-Rodriguez2018, ; Parra-Rodriguez2018a, ; Ansari2018, ; Krupko2018, ; Malekakhlagh2018, ; Solgun2017, ; Petrescu2019, ; You2019-Koch, ; DiPaolo2019, ; Menke2021, ; Gely2020, ; Kyaw2020, ; Yan2020, ; Minev2021-lom, ; qiskit-metal-code, ), as experimental architectures (Barends2013, ; Minev2016, ; Brecht2016, ; Reagor2016-cavity, ; Gambetta2017, ; Rosenberg2017, ; Versluis2017, ; Naik2017, ; Krantz2019, ; Kjaergaard2019, ) and nonlinear devices (Josephson1962, ; Vijay2010, ; Kerman2010, ; Larsen2015, ; DeLange2015-nanowire, ; Janvier2015-atomic-point, ; Maleeva2018, ; Wang2018, ) scale in both complexity and diversity.

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 $m$ is stored in element $j$? This leads to a constrained number between zero and one, the EPR, denoted $p_{mj}$ (Minev2019-Thesis, ). 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 $p_{ml}$ of lossy element $l$ in mode $m$.

The EPR method deviates from previous black-box quantization work (Nigg2012, ; Solgun2014, ; Solgun2015, ), which uses the impedance-response matrix, denoted $Z_{jj^{\prime}}\left(\omega\right)$, where $j$ and $j^{\prime}$ index ports associated with nonlinear elements. For all pairs of ports, the complex function $Z_{jj^{\prime}}\left(\omega\right)$ is calculated from a finite-element (FE) driven simulation in the vicinity of the eigenfrequency of every mode. Our method replaces these steps with a more economical FE eigenmode simulation, from which one extracts the energy participations $p_{ml}$ and $p_{mj}$, 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.

Figure 1: Conceptual overview. a Illustration of the physical model of an example quantum device, which comprises two three-dimensional (3D) cavities (grey enclosures), each housing several qubit chips (green boxes). A close-up view of one of the chips is depicted in the inset. The dotted box in the center of the chip schematically outlines a non-linear inductive sub-circuit, referred to as a Josephson dipole. b Results of a finite-element eigenmode analysis ($\mathrm{FE}_{\mathrm{e}}$) of the Josephson circuit linearized about its equilibrium. The $m$-th mode eigenfrequency and electric and magnetic fields are $\omega_{m}$, $\vec{E}_{m}\left(\vec{r}\right)$, and $\vec{H}_{m}\left(\vec{r}\right)$, respectively, where $r$ denotes spatial position. Center inset: $\left|\vec{E}_{m}\right|$ profile (red: high; blue: low) for the fundamental mode of one of the 3D cavities. Additional FE driven simulations ($\mathrm{FE}_{\mathrm{d}}$) are unnecessary; i.e., the impedance matrix $Z_{jj^{\prime}}\left(\omega\right)$ is not calculated. c The Hamiltonian $\hat{H}_{\mathrm{full}}$, which includes nonlinear interactions to arbitrary order (see Results ), is computed directly from the eigenanalysis via the EPRs $p_{mj}$ and EPR signs $s_{mj}=\pm 1$ of the junctions, $j$. Dissipative contributions due to a lossy element $l$ are similarly computed from the loss EPRs $p_{ml}$; for linear dissipation, the EPR signs $s_{ml}$ 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).

II RESULTS AND DISCUSSION

Figure 2: Quantizing a simple circuit. a Illustration of a 3D cavity enclosing a transmon qubit chip. The cross symbol marks the location of a Josephson junction. Vertical blue arrows depict the electric field $\vec{E}_{m}\left(\vec{r}\right)$ of the fundamental cavity mode, TE₁₀₁. b Equivalent two-mode lumped-element- representation of the distributed circuit. Operators $\hat{a}_{q}$ and $\hat{a}_{c}$ denote the qubit and cavity mode operators, respectively.

II.1 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 (Koch2007, ) 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 $\hat{H}_{\mathrm{full}}$ can be conceptually separated in two contributions (see Supplementary Section A2),

$$\hat{H}_{\mathrm{full}}=\hat{H}_{\mathrm{lin}}+\hat{H}_{\mathrm{nl}}\;,$$ (1)

where $\hat{H}_{\mathrm{lin}}$ consists of all terms associated with the linear response of the junction and the resonator structure, and $\hat{H}_{\mathrm{nl}}$ 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 (Yurke1984, ; Devoret1995, ) (see Supplementary Section A4):

$$\hat{H}_{\mathrm{lin}} = \hbar\omega_{c}\hat{a}_{c}^{\dagger}\hat{a}_{c}+\hbar\omega_{q}\hat{a}_{q}^{\dagger}\hat{a}_{q}\;,$$ (2)

$$\hat{H}_{\mathrm{nl}} = -E_{J}\left[\cos\left(\hat{\varphi}_{J}\right)+\hat{\varphi}_{J}^{2}/2\right]\;,$$ (3)

$$\hat{\varphi}_{J} = \varphi_{q}\left(\hat{a}_{q}+\hat{a}_{q}^{\dagger}\right)+\varphi_{c}\left(\hat{a}_{c}+\hat{a}_{c}^{\dagger}\right)\;,$$ (4)

where $\omega_{c}$ and $\omega_{q}$ are the angular frequencies of the cavity and qubit eigenmodes defined associated with $\hat{H}_{\mathrm{lin}}$, respectively, and where $\hat{a}_{c}$ and $\hat{a}_{q}$ are their annihilation operators, respectively. The Josephson energy $E_{J}$ can be computed from the Ambegaokar-Baratoff formula adapted to the measured room-temperature resistance of the junction (Gloos2000, ). The junction reduced generalized flux $\hat{\varphi}_{J}$ corresponds to the classical variable $\varphi_{J}\left(t\right)\coloneqq\int_{-\infty}^{t}v_{J}\left(\tau\right)\,\mathrm{d}\tau/\phi_{0}$, where $v_{J}(\tau)$ is the instantaneous voltage across the junction (Yurke1984, ; Devoret1995, ), and $\phi_{0}\coloneqq\hbar/2e$ 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, $\varphi{}_{c}$ and $\varphi_{q}$ 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 $g$ (Koch2007, ), 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: $\omega_{c},\omega_{q},\varphi_{q}$ and $\varphi_{c}$. As we will show, we extract and compute these quantities from an eigenanalysis of the classical distributed circuit corresponding to $\hat{H}_{\mathrm{lin}}$. 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 eigenfrequencies $\omega_{c}$ and $\omega_{q}$ fully determine $\hat{H}_{\text{lin}}$ (see Supplementary Section C for FE methodology).

To determine $\hat{H}_{\mathrm{nl}}$, we need the quantum zero-point fluctuations $\varphi_{q}$ and $\varphi_{c}$, which are calculated from the participation of the junction in the eigenfield solutions. The participation $p_{m}$ of the junction in mode $m\in\left\{c,q\right\}$ is defined to be the fraction of inductive energy stored in the junction relative to the total inductive energy stored in the entire circuit,

$$p_{m}\coloneqq\frac{\text{Inductive energy stored in the junction}}{\text{Total inductive energy stored in mode }m}\;,$$ (5)

evaluated when only mode $m$ is excited. Thus, $p_{m}$ can be computed from the electric $\vec{E}_{m}(\vec{r})$ and magnetic $\vec{H}_{m}(\vec{r})$ eigenfields as detailed in Supplementary Section C2; $\vec{r}$ denotes spatial position. In the quantum setting, Eq. (5) links $p_{m}$, $\hat{\varphi}_{J}$, and the state of the circuit,

$$p_{m}=\frac{\langle\psi_{m}|\frac{1}{2}E_{J}\hat{\varphi}_{J}^{2}|\psi_{m}\rangle}{\langle\psi_{m}|\frac{1}{2}\hat{H}_{\mathrm{lin}}|\psi_{m}\rangle}\;,$$ (6)

where $\left|\psi_{m}\right\rangle$ denotes a coherent state or a Fock excitation of mode $m$. 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 zero-point fluctuations $\varphi_{c}$ and $\varphi_{q}$ as a function of the classical energy participations $p_{m}$,

$$\varphi_{c}{}^{2}=p_{c}\frac{\hbar\omega_{c}}{2E_{J}}\quad\text{and}\quad\varphi_{q}{}^{2}=p_{q}\frac{\hbar\omega_{q}}{2E_{J}}\;,$$ (7)

which completely determines $\hat{H}_{\mathrm{nl}}$, and thus completes the description of the system Hamiltonian $\hat{H}_{\mathrm{full}}$. Here, $\varphi_{c}$ and $\varphi_{q}$ 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 $\hat{H}_{\mathrm{full}}$ the transition frequencies and nonlinear couplings between modes. Depending on the case, this can be done approximately or exactly using numerical or analytical techniques (DiPaolo2019, ). This task is easily achieved if $\hat{H}_{\mathrm{nl}}$ is a perturbation to $\hat{H}_{\text{lin}}$ (Koch2007, ). In this limit, $\hat{H}_{\mathrm{full}}$ for our qubit-cavity example can be approximated by the effective, excitation-number-conserving Hamiltonian, see Supplementary Eq. (B.2),

$$\hat{H}_{\text{eff}}=\left(\omega_{q}-\Delta_{q}\right)\hat{n}_{q}+\left(\omega_{c}-\Delta_{c}\right)\hat{n}_{c}-\chi_{qc}\hat{n}_{q}\hat{n}_{c}\\ -\frac{1}{2}\alpha_{q}\hat{n}_{q}\left(\hat{n}_{q}-\hat{1}\right)-\frac{1}{2}\alpha_{c}\hat{n}_{c}\left(\hat{n}_{c}-\hat{1}\right)\;,$$ (8)

where $\hat{n}_{q}=\hat{a}_{q}^{\dagger}\hat{a}_{q}$ and $\hat{n}_{c}=\hat{a}_{c}^{\dagger}\hat{a}_{c}$ denote the qubit and cavity excitation-number operators, respectively, $\Delta_{q}$ denotes the ‘Lamb shift’ of the qubit frequency due to the dressing of this nonlinear mode by quantum fluctuations of the fields, $\alpha_{q}$ ($\alpha_{c}$) is the qubit (cavity) anharmonicity, and $\chi_{qc}$ is the qubit-cavity dispersive shift (cross-Kerr coupling). The Hamiltonian parameters can be calculated directly from the EPR, see Supplementary Section B,

$$\alpha_{q} =\frac{1}{2}\chi_{qq}=p_{q}^{2}\frac{\hbar\omega_{q}^{2}}{8E_{J}}\;,$$ (9)

$$\alpha_{c} =\frac{1}{2}\chi_{cc}=p_{c}^{2}\frac{\hbar\omega_{c}^{2}}{8E_{J}}\;,$$ (10)

$$\chi_{qc} =p_{q}p_{c}\frac{\hbar\omega_{q}\omega_{c}}{4E_{J}}\;.$$ (11)

Experimentally, the qubit Lamb shift can be obtained as $\Delta_{q}=\alpha_{q}-\chi_{qc}/2$. Since a single EPR $p_{m}$ determines the nonlinear interaction for each mode, the parameters $\chi_{qc}$ and $\alpha_{q}$ are interdependent,

$$\chi_{qc}=\sqrt{\chi_{qq}\chi_{cc}}=2\sqrt{\alpha_{q}\alpha_{c}}\;.$$ (12)

As shown in Supplementary Section A7, the EPRs $p_{c}$ and $p_{q}$ obey the constraints

$$0\leq p_{q},p_{c}\leq 1\quad\text{and}\quad p_{q}+p_{c}=1\;.$$ (13)

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.

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 $\Phi_{j}^{\mathrm{ext}}$ threading the loop. c Conceptual decomposition of a general Josephson dipole, denoted $j$. For convenience, its potential energy function $\mathcal{E}_{j}(\Phi_{j};\Phi_{j}^{\mathrm{ext}})$ can be Taylor expanded in a sum of nonlinear inductive contributions of increasing order $\Phi_{j}^{p}$, with relative amplitude $c_{jp}$, where $p$ denotes the index in the series. The energy function can be subjected to external bias parameters $\Phi_{j}^{\mathrm{ext}}$, such as flux or voltages. d Schematic diagram of a general Josephson circuit conceptually resolved into a purely dissipative (left, $\kappa_{m}^{l}$), linear (middle, $\hat{H}_{\mathrm{lin}}$), and nonlinear (right, $\hat{H}_{\mathrm{nl}}$) constitutions.

II.2 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 (Blais2004, ; Wallraff2004, ; Barends2013, ; FYan2016, ) (2D), quasi-planar (Minev2013, ; Minev2016, ; Brecht2016, ; Rosenberg2017, ) (2.5D), or three-dimensional (Paik2011, ; Rigetti2012, ; Reagor2016-cavity, ; Axline2016, ) (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 (Yurke2006, ; HoEom2012, ; Vissers2015, ; Mortensen2016, ), can be thought of as a series of lumped ones.

The simplest nonlinear device comprises a single element, such as a Josephson tunnel junction [see Fig. 3(a)], an atomic-point contact (Koops1996, ; Bretheau2013, ), a nanobridge (Vijay2010, ; Peltonen2016, ), a semiconducting nanowire (Mooij2006, ; Ku2010, ; Abay2014, ; Larsen2015, ; DeLange2015-nanowire, ; Casparis2016, ), or another hybrid structure (Shim2014, ). A multi-element device, such as a SQUID (Zimmerman1966, ; Clarke2004, ), a SNAIL (Frattini2017, ) [see Fig. 3(b)], a superinductance (Manucharyan2012, ; Pop2014, ; Muppalla2017, ), or a junction array (Corlevi2006, ; Hutter2011, ; Bourassa2012, ; Masluk2012, ; Bell2012, ; Weibl2015, ; Macklin2015, ) refers to a subcircuit composed of purely inductive lumped elements. This subcircuit can also be subjected to external controls, such as voltage or flux biases.

The general nonlinear device that we now consider, referred to as a Josephson dipole, is any lumped, purely-inductive, nonlinear subcircuit with two terminals. The 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 (Clarke2004, ) is described by the energy function $\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)=-E_{j}\left(\Phi_{j}^{\mathrm{ext}}\right)\cos\left(\Phi_{j}/\phi_{0}\right)$, where $\Phi_{j}$ is the generalized flux across the device terminals (Yurke1984, ; Devoret1995, ), $E_{j}$ is the effective Josephson energy, $\Phi_{j}^{\mathrm{ext}}$ is the external flux bias, and the subscript $j$ denotes the $j$-th Josephson dipole in the circuit. The flux $\Phi_{j}$ is defined as the deviation away from the value in equilibrium, as discussed below. To ease the notation, parameters such as $\Phi_{j}^{\mathrm{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 $\mathcal{E}_{j}^{\mathrm{lin}}$ accounts for the linear response of the dipole, while the other $\mathcal{E}_{j}^{\mathrm{nl}}$ accounts for the nonlinear response,

$$\mathcal{E}_{j}\left(\Phi_{j}\right)=\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)+\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)\;,$$ (14)

where

$$\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right) \coloneqq \frac{1}{2}E_{j}\left(\frac{\Phi_{j}}{\phi_{0}}\right)^{2}\;,$$ (15)

and where the constant $E_{j}$ sets the scale of the junction energy. This energy scale can be represented by the linear inductance $L_{j}\coloneqq\phi_{0}/E_{j}$ 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 $\Phi_{j}$ will differ in Eqs. (14) and (15) as a function of the equilibrium. Due to this adjustment, terms linear in $\Phi_{j}$ are absent from Eq. (15) by construction. Supplementary Sections A8 and A9 discuss these equilibrium considerations in detail.

II.2.1 Quantum Hamiltonian.

Having conceptually carved out nonlinear contributions from the system Hamiltonian $\hat{H}_{\mathrm{full}}$, and collected them in the set of $\mathcal{E}_{j}^{\mathrm{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 $L_{j}$. We will use the eigenmodes of the linearized circuit to explicitly construct $\hat{H}_{\mathrm{full}}$. 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)

$$\hat{H}_{\mathrm{lin}} =\sum_{m=1}^{M}\hbar\omega_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}\;,$$ (16)

where $M$ is the number of modes addressed by the numerical simulation, $\omega_{m}$ is the solution eigenfrequency of mode $m$, and $\hat{a}_{m}$ the corresponding mode amplitude (annihilation operator), defined by the mode eigenvector. We emphasize that the frequencies $\omega_{m}$ will be significantly perturbed by the Lamb shifts $\Delta_{m}$, and should be seen as an intermediate parameter entering in the calculation of the rest of the nonlinear Hamiltonian,

$$\hat{H}_{\mathrm{nl}} =\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}=\sum_{j=1}^{J}E_{j}\left(c_{j3}\hat{\varphi}_{j}^{3}+c_{j4}\hat{\varphi}_{j}^{4}+\cdots\right)$$ (17)

$$=\sum_{j=1}^{J}E_{j}\sum_{p=3}^{\infty}c_{jp}\hat{\varphi}_{j}^{p}\;,$$ (18)

$$\hat{\varphi}_{j} =\sum_{m=1}^{M}\varphi_{mj}\left(\hat{a}_{m}^{\dagger}+\hat{a}_{m}\right)\;,$$ (19)

where $J$ is the total number of junctions and $\hat{\varphi}_{j}\coloneqq\hat{\Phi}_{j}/\phi_{0}$. In Eq. 17, we have introduced a Taylor expansion of $\mathcal{E}_{j}^{\mathrm{nl}}$, where the energy $E_{j}$ and expansion coefficients $c_{jp}$ are known from the fabrication of the Josephson circuit, see Fig. 3(c). For example, for a Josephson junction, the constant $E_{j}$ is just the Josephson energy, while $c_{jp}$ are the coefficients of the cosine expansion; i.e., $c_{jp}$ is $0$ for odd $p$ and $\left(-1\right)^{p/2+1}/p!$ for even $p$. The expansion is helpful for analytics but does not need to be used in the numerical analysis of $\hat{H}_{\mathrm{nl}}$, see Supplementary Section A.

The Hamiltonian $\hat{H}_{\mathrm{full}}$ is specified since the operators $\hat{\varphi}_{j}$ are expressed in terms of the mode amplitudes as a linear combination (see Supplementary Section A5). Here, $\varphi_{mj}$ are the dimensionless, real-valued, quantum zero-point fluctuations of the reduced flux of junction $j$ in mode $m$. Determination of $\hat{H}_{\mathrm{full}}$ is now reduced to computing $\varphi_{mj}$. We achieve this by employing a generalization of the energy-participation ratio.

The energy-participation ratio $p_{mj}$ of junction $j$ in eigenmode $m$ is defined to be the fraction of inductive energy stored in the junction when only that mode is excited,

$$p_{mj} \coloneqq\frac{\text{Inductive energy stored in junction }j}{\text{Inductive energy stored in mode }m}$$

$$=\frac{\langle\psi_{m}|\frac{1}{2}E_{j}\hat{\varphi}_{j}^{2}|\psi_{m}\rangle}{\langle\psi_{m}|\frac{1}{2}\hat{H}_{\mathrm{lin}}|\psi_{m}\rangle}\;,$$ (20)

which is a straightforward extension of Eq. (5), and is similarly computed using normal ordering (see Supplementary Section A6). The EPR $p_{mj}$ is computed from the eigenfield solutions $\vec{E}_{m}(\vec{r})$ and $\vec{H}_{m}(\vec{r})$ as explained in Supplementary Section C3. It is a bounded, non-negative, real number, $0\leq p_{mj}\leq 1$. A zero EPR $p_{mj}=0$ means that junction $j$ is not excited in mode $m$. A unity EPR $p_{mj}=1$ means that junction $j$ is the only inductive element excited in the mode.

From the EPR $p_{mj}$, one directly computes the variance of the quantum zero-point fluctuations (see Appendix A6),

$$\boxed{\varphi_{mj}^{2}=p_{mj}\frac{\hbar\omega_{m}}{2E_{j}}\;.}$$ (21)

Equation (21) constitutes the bridge between the classical solution of the linearized Josephson circuit and the quantum Hamiltonian $\hat{H}_{\mathrm{full}}$ of the full Josephson system, up to the sign of $\varphi_{mj}$.

II.2.2 Universal EPR properties

The quantum fluctuations $\varphi_{mj}$ 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 $j$ and one set of inequalities per mode $m$,

$$\sum_{m=1}^{M}p_{mj}=1\quad\text{and}\quad 0\leq\sum_{j=1}^{J}p_{mj}\leq 1\;.$$ (22)

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 $M$ 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

$$\varphi_{mj}=s_{mj}\sqrt{p_{mj}\hbar\omega_{m}/2E_{j}}\;,$$ (23)

where the EPR sign $s_{mj}$ of junction $j$ in mode $m$ is either $+1$ or $-1$. The EPR sign encodes the relative direction of current flowing across the junction. Only the relative value between $s_{mj}$ and $s_{mj^{\prime}}$ for $j\neq j^{\prime}$ has physical significance (see Supplementary Figure S8). The EPR sign $s_{mj}$ is calculated in parallel with the process of calculating $p_{mj}$, from the field solution $\vec{H}(\vec{r})$, see Supplementary Section (C.8). We now obtain the EPR orthogonality relationship

$$\sum_{m=1}^{M}s_{mj}s_{mj^{\prime}}\sqrt{p_{mj}p_{mj^{\prime}}}\, =0\;,$$ (24)

valid when the sum from $1$ to $M$ covers all the relevant modes, see Supplementary Eq. (A.60).

II.2.3 Excitation-number-conserving interactions

Thus, as announced, knowledge of the energy-participation ratios completely specifies $\hat{H}_{\mathrm{nl}}$, through Eqs. (17), (19), and (21). The Hamiltonian can now be analytically or numerically diagonalized using various computational techniques (DiPaolo2019, ). In this section, our focus will be now to explicitly handle the effect of the nonlinear interactions $\hat{H}_{\mathrm{nl}}$ on the eigenmodes. Before treating the case of a general nonlinear interaction, we focus on the leading-order effect of $\hat{H}_{\mathrm{nl}}$ in the case of the ‘canonical’ Josephson system. In this case, the $J$ Josephson dipoles are all Josephson tunnel junctions, characterized by Eq. (A.14b), and the dispersive regime is satisfied for all pairs of modes $k$ and $m$; i.e., $\omega_{k}-\omega_{m}\gg E_{j}c_{jp}\left<\hat{\varphi}_{j}^{p}\right>$ for $p\geq 3$ and in the absence of strong drives. The leading-order nonlinear terms are the subset of $p=4$ terms that conserve excitation number. After normal ordering, see Supplementary Section B1, one finds the effective Hamiltonian

$$\hat{\overline{H}}_{4}=-\hbar\sum_{m=1}^{M}\Delta_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}+\frac{\alpha_{m}}{2}\hat{a}_{m}^{\dagger 2}\hat{a}_{m}^{2}+\sum_{n<m}\chi_{mn}\hat{a}_{m}^{\dagger}\hat{a}_{n}^{\dagger}\hat{a}_{m}\hat{a}_{n}\;,$$ (25)

which is a generalization of the one found in Eq. (8). In Eq. (25), we have introduced the Lamb shift $\Delta_{m}$ of mode $m$, the anharmonicity $\alpha_{m}$ of the mode, and its total dispersive shift $\chi_{mn}$ (so-called cross-Kerr term) with a different mode, labeled $n$. Each of these parameters is directly calculated from the EPRs. As shown in Supplementary Section (B.3a), for arbitrary $m$ and $n$,

$$\chi_{mn}=\sum_{j=1}^{J}\frac{\hbar\omega_{m}\omega_{n}}{4E_{j}}p_{mj}p_{nj}\;,$$ (26)

while $\alpha_{m}=\chi_{mm}/2$ and $\Delta_{m}=\sum_{n=1}^{M}\chi_{mn}/2$. Equation (26) implements, mathematically, the idea that the amplitude of these nonlinear couplings is the result of a spatial-mode scalar product of the EPRs. Remarkably, from Eq. (26) it is seen that the EPRs are essentially the only free parameters subject to design, when determining the nonlinear couplings, since $\omega_{m}$, $\omega_{n}$, and $E_{j}$ are generally tightly constrained by experimental considerations.

Equation (26) can be cast in matrix form by introducing the EPR matrix

$$\mathrm{\mathbf{P}}\coloneqq\left(\begin{array}[]{ccc}p_{11}&\cdots&p_{1J}\\ \vdots&\ddots&\vdots\\ p_{M1}&\cdots&p_{MJ}\end{array}\right)\;,$$ (27)

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 $\mathrm{\mathbf{\Omega}}\coloneqq\operatorname{diag}\left(\omega_{1},\ldots,\omega_{M}\right)$ and junction energies $\mathrm{\mathbf{E_{J}}}\coloneqq\operatorname{diag}\left(E_{1},\ldots,E_{J}\right)$, which lead to the matrix form of Eq. (26),

$$\begin{array}[]{rcclc}\text{Kerr matrix:}&&\bm{\chi}&=&\frac{\hbar}{4}\mathrm{\mathbf{\Omega PE_{J}^{-1}P^{\intercal}}}\bm{\Omega}\;,\\ \text{Anharmonicity:}&&\alpha_{m}&=&\frac{1}{2}\left[\bm{\chi}\right]_{mm}\;,\\ \text{Lamb shift:}&&\Delta_{m}&=&\frac{1}{2}\sum_{m^{\prime}=1}^{M}\left[\bm{\chi}\right]_{mm^{\prime}}\;.\end{array}$$ (28)

We have defined the symmetric matrix of dispersive shifts $\bm{\chi}$, with elements $\left[\bm{\chi}\right]_{mm^{\prime}}=\chi_{mm^{\prime}}$. Further discussion of the matrix approach and applications to $p$th-order corrections is deferred to Supplementary Section B2, and the amplitude of an arbitrary multi-photon interaction stemming from the full $\hat{H}_{\text{nl}}$ is calculated in Supplementary Section B3.

II.2.4 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 (Nigg2012, ), one can completely characterize both $\hat{H}_{\mathrm{full}}$ and the effect of dissipative elements in the circuit from the eigenfield solutions, $\vec{E}_{m}(\vec{r})$ and $\vec{H}_{m}(\vec{r})$. The list of dissipative elements include bulk and surface dielectrics (Martinis2005, ; Patel2013, ; Dial2016, ), thin-film metals (Vissers2012, ; Minev2013, ), surface interfaces (Wenner2011, ; Geerlings2012, ; Sandberg2013, ; Wang2015, ; Bruno2015-loss, ), and metal seams (Brecht2015, ). The energy-participation ratio of a dissipative element $l$ in mode $m$ will be denoted $p_{ml}$. It is calculated similarly to $p_{mj}$, as summarized in Supplementary Section D. The participation $p_{ml}$ and the quality factor $Q_{l}$ of the material of this element are used to estimate the total quality factor of mode $m$ in the standard way when the fields are not greatly altered by the dissipation ($Q_{m}\gg 1$) (Gao2008, ; Geerlings2013, ; Reagor2016-cavity, ; Brecht2017-micromachined, ),

$$Q_{m}^{-1}=\sum_{l}{p_{ml}Q_{l}^{-1}}\;.$$ (29)

Experimental values of $Q_{l}$ are found in the literature, and some are provided in Supplementary Section D. Equation (29) and the dissipative EPR $p_{ml}$ 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 (Martinis2014, ).

Figure 4: Comparison between theory and experiment over five-orders of magnitude in energy scale of the system Hamiltonian $\hat{H}_{\mathrm{full}}$ 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 $\hat{H}_{\mathrm{full}}$, dressed frequencies $\omega_{m}^{\prime}$, bare anharmonicities $\alpha_{m}$, and cross-Kerr interactions $\chi_{mn}$, were measured and calculated using the EPR method with our open-source pyEPR package (Note1, ). Gray line is of slope one, representing ideal agreement between theory and experiment.

II.3 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 (Minev2016, ) consisting of a flip-chip qubit embedded in a two-mode whispering gallery mode resonator (Minev2013, ) (WGMR). The specifics of each sample are discussed in the Methods.

For each sample, we measured the circuit parameters of interest: dressed mode frequencies $\omega_{m}-\Delta_{m}$, anharmonicities of qubits and high-Q cavities $\alpha_{m}$, cross-Kerr frequencies $\chi_{mn}$, and input-output quality factors $Q_{C}$ 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 $E_{j}$ was inferred from room-temperature resistance measurements (Gloos2000, ). To account for the error bars of the measurement and the drift in resistance over time, $E_{j}$ 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 $C_{J}$ 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 $C_{J}=4\,\mathrm{fF}$. From the eigenfield solutions, we calculated the EPRs $p_{mj}$ and the sign $s_{mj}$ to construct $\hat{H}$ and extract its parameters. Detailed steps of the procedure can be found in Supplementary Section C. The results are presented in Tables 1–3.

Figure 4 summarizes the agreement of the measured and calculated sample parameters, which span five orders of magnitude in frequency. Accounting for $C_{J}$, 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 $E_{j}$, and otherwise, by neglecting $C_{j}$, 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 (Note1, ). For current standard applications, we find the agreement sufficient. Further improvements in accuracy would require improved ability to estimate the Josephson dipole energy $E_{j}$ and its intrinsic capacitance $C_{J}$. 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 (Vijay2010, ; Peltonen2016, ), nanowires (Mooij2006, ; Abay2014, ; Larsen2015, ; DeLange2015-nanowire, ; Casparis2016, ), and kinetic-inductance thin-films (HoEom2012, ; Vissers2015, ; Maleeva2018, ). 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 (DiPaolo2019, ). We contribute an open-source package pyEPR , which automates the EPR method, and was tested in the design of several further experiments (Leghtas2015, ; Minev2016, ; Mundhada2017, ; Touzard2017, ; Campagne-Ibarcq2018-RemoteEnt, ; Muppalla2017, ; Wang2019-cav-atten, ; Grimm2019, ; Minev2019Nature, ; Campagne-Ibarcq2019, ; Winkel2019, ; Winkel2019a, ).

III METHODS

III.1 Methods of the experiment

Device	Frequency (MHz)			Anharmonicity (MHz)		Cross-Kerr (MHz)						I-O coupling
	$\omega_{\mathrm{D}}/2\pi$	$\omega_{\mathrm{B}}/2\pi$	$\omega_{\mathrm{C}}/2\pi$	$\alpha_{\mathrm{D}}/2\pi$	$\alpha_{\mathrm{B}}/2\pi$	$\chi_{\mathrm{DB}}/2\pi$		$\chi_{\mathrm{BC}}/2\pi$		$\chi_{\mathrm{DC}}/2\pi$		$Q_{\mathrm{C}}$
R9C1	4951	5664	9158	138	170	92.		4.	7	0.	4	$5.20\times 10^{3}$
	4866	5691	9154	150	185	99.		4.	2	0.	55	$7.40\times 10^{3}$
	-1.7%	0.5%	-0.04%	8%	8%	7%		-12%		27%		29%
R2C1	4823	5567	8947	150	192	64.	5	4.	8	0.	3	$4.97\times 10^{3}$
	4770	5640	8950	161	211	67.	7	5.	88	0.	46	$5.44\times 10^{3}$
	-1.1%	1.3%	0.03%	6.8%	9%	4.	7%	18%		35%		9 %
R7C1	4726	5475	8999	156	189	67.		4.	8	0.	34	$2.68\times 10^{3}$
	4770	5640	8950	161	211	67.	7	5.	88	0.	46	$3.07\times 10^{3}$
	0.9%	2.9%	-0.55%	3.1%	10%	1%		18%		26%		13%
R3C2	4845	5620	8979	152	195	61		5.	1	0.	3	$2.11\times 10^{3}$
	4770	5640	8950	161	211	67.	7	5.	88	0.	46	$1.78\times 10^{3}$
	-1.5%	0.4%	-0.3%	5.6%	7.6%	9.	9%	13%		35%		-19%
R3C1	4688	5300	9003	148	174	85		5.		0.	33	$2.43\times 10^{3}$
	4745	5265	8922	159	198	73		5.	1	0.	37	$5.65\times 10^{3}$
	1.2%	-0.7%	-0.9%	6.9%	12.1%	-16%		2%		9%		57%
DT3	6160	7110	9170	130	150	278		3.		2.	5	$9.17\times 10^{3}$
	6100	7141	9155	140	177	312		3.	9	3.	1	$7.33\times 10^{3}$
	-1.0%	0.4%	-0.15%	7%	15%	11%		23%		19%		-25%

Table 1: 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 $Q_{\mathrm{C}}$. For each device, the first (second) row quantifies the measured, $m$, (bare calculated, $c$) values. The third row quantifies the bare agreement, i.e., $\left(c-m\right)/c$. 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 $C_{J}=4\,\mathrm{fF}$; hence, an overall corrected agreement of 4.3% for this column.

Device	Frequency (MHz)			Anharmonicity (MHz)	Cross-Kerr (MHz)		I-O coupling
	$\omega_{\mathrm{Q}}/2\pi$	$\omega_{\mathrm{S}}/2\pi$	$\omega_{\mathrm{C}}/2\pi$	$\alpha_{\mathrm{D}}/2\pi$	$\chi_{\mathrm{QS}}/2\pi$	$\chi_{\mathrm{QC}}/2\pi$	$Q_{\mathrm{C}}$
WG1	4890	7070	7267	310	0.25	0.30	$20\times 10^{3}$
	4820	7020	7340	325	0.29	0.33	$16\times 10^{3}$
	-1.4%	-0.7%	1.0%	4.6%	13%	9%	-22%

Table 2: 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 $Q_{\mathrm{C}}$. For each device, the first (second) row quantifies the measured, $m$, (bare calculated, $c$) values. The third row quantifies the bare agreement, i.e., $\left(c-m\right)/c$.

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 $\mu$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 (Rigetti2009, ; Lecocq2011-bridge-free, ; Pop2012-Junction, ) the Al/Al$\text{O}_{\text{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_x barrier of the tunnel junction. Prior to the first deposition, to reduce junction aging (Pop2012-Junction, ), the exposed wafer surfaces were exposed to 1 minute oxygen-argon plasma cleaning, under a pressure of $3\times 10^{-3}$ mbar. After wafer dicing (ADT ProVecturs 7100) and chip cleaning, the normal-state resistance $R_{N}$ of the Josephson junctions was measured to provide an estimate of the Josephson energy, $E_{J}$, of the device junctions. The junction energy was to first order estimated by an extrapolation of $R_{N}$ from room temperature to the operating sample temperature, at approximately $15$ mK, using the Ambegaokar-Baratoff relation (Ambegaokar1963, ),

$$E_{J}=\frac{1}{2}\frac{h\Delta}{\left(2e\right)^{2}}R_{N}^{-1}\;,$$ (30)

where $\Delta$ is the superconducting gap of aluminum, $e$ is the elementary charge, and $h$ is Planck’s constant.

Frequency (MHz)		Anharmonicity (MHz)		Cross-Kerr (MHz)
$\omega_{\mathrm{D}}/2\pi$	$\omega_{\mathrm{B}}/2\pi$	$\alpha_{\mathrm{D}}/2\pi$	$\alpha_{\mathrm{B}}/2\pi$	$\chi_{\mathrm{DB}}/2\pi$
6010	8670	85	180	278
5824	8878	97	206	281
-3.2%	2.3%	12%	13%	1.1%

Table 3: Two-qubit, one-waveguide devices. Summary of measured and calculated Hamiltonian parameters for the device described in Methods . Indices D and B denote the dark and bright modes, respectively. For each device, the first (second) row summarizes the measured, $m$, (bare calculated, $c$) values. The third row quantifies the bare agreement, i.e., $\left(c-m\right)/c$.

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 (Wallraff2004, ; Leghtas2015, ; Minev2019-Thesis, ). High-magnetic-permeability, $\mu$-metal (Amumetal A4K) shields together with aluminum superconducting shields enclosed all samples. Microwave input and output lines were filtered with Eccosorb CR-110 infrared-frequency filters (Rigetti2012, ; Geerlings2013, ), thermally anchored at the 15 mK stage. Output lines were additionally filtered with cryogenic isolators (Quinstar CWJ1019-K414) and 12 GHz K&L multi-section lowpass filters. Output lines leading up to the high-electron-mobility transistor (HEMT) amplifier (Low Noise Factory), anchored at 4 K, were superconducting (CoaxCo Ltd. SC-086/50-NbTi-NbTi PTFE).

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 (Bergeal2010, ; Abdo2013, ), 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. (Roy2017-Review, ) 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 (Geerlings2013, ; Reagor2016, ). Cross-Kerr energies were determined from dressed dephasing measurements (Gambetta2006-dephasing, ; Gambetta2008-qm-traj, ). 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 $\pi/2$ pulse. Subsequently, a weak readout tone excited the readout cavity of the qubit for a fixed duration, 10 times the readout cavity lifetime $\kappa_{r}$, after which we measure the qubit X and Y Bloch vectors, after waiting for a time $5/\kappa_{r}$ 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 $Q_{C}$, was extracted from the spectroscopic response of the readout cavity at low photon numbers (Geerlings2013, ; Reagor2016, ), by measuring the scattering parameters, $S_{21}$ or $S_{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\,\Omega$ 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 (Note1, ).

The mode quality due to the input-output coupling, $Q_{C}$, 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 $Q_{C}$ 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.

Figure 5: Two–qubit, one-cavity devices. a and b Not-to-scale diagram illustrating chip configurations A and B, respectively. Vertical blue arrow indicates cavity electric field orientation. Crosses mark the location of Josephson tunnel junctions. c Optical photograph of sample R2C1. Bottom half of aluminum sample holder is visible; top half is removed. The two-qubit chip (outlined by the dashed green box) is housed in the middle of the readout cavity (highlighted in blue). Cavity fundamental mode electric field profile $\vec{E}$ depicted by arrows. Input-output SMA pin coupler labeled I-O.

III.2 Devices

III.2.1 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 (Paik2011, ). 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 1.

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, $E_{j}$. 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 $E_{j}$. 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-in-spirit characteristics with the promising recent developments reported in Refs. Gambetta2011-Purcell, ; Srinivasan2011, ; Diniz2013, ; Dumur2015, ; Zhang2017, ; Roy2017-3qubits, , 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 bright-mode 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 coupling 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 (Koch2007, ; Devoret2007-HowStrongCanCouplingBe, ) 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 $\vec{\mathrm{d}}_{B}$. 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 $\vec{\mathrm{d}}_{B}\cdot\vec{E}$, where $\vec{E}$ 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, $\vec{\mathrm{d}}_{B}\parallel\vec{E}$, 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 (Houck2008-Purcell, ; Gambetta2011-Purcell, ), 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 $E_{J}$, 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,

$$p_{\mathrm{D}1}=p_{\mathrm{D}2}=p_{\mathrm{B}1}=p_{\mathrm{B}2}\approx\frac{1}{2}\;.$$ (31)

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:

$$s_{\mathrm{D}1} =1\;, s_{\mathrm{D}2} =-1\;,$$ (32)

$$s_{\mathrm{B}1} =1\;, s_{\mathrm{B}2} =+1\;.$$ (33)

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, $\vec{d}_{\mathrm{B}}$, must be large, since it is the sum of the two junction dipole contributions. Secondly, $\vec{d}_{\mathrm{B}}$ must be oriented in the vertical direction, parallel to the cavity electric field $\vec{E}$. Hence, we would conclude that the bright mode coupling is large, $\vec{d}_{\mathrm{B}}\cdot\vec{E}\gg 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, $\vec{d}_{\mathrm{D}}=0$. Thus, we should not expect any coupling between the dark qubit and the cavity mode, $\vec{d}_{\mathrm{D}}\cdot\vec{E}=0$. To the contrary of this conclusion, as can be seen in the measured results in Table 1, 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

$$\chi_{\mathrm{DC}} =\frac{\hbar\omega_{D}\omega_{C}}{4E_{J}}\left(p_{\mathrm{D}1}p_{\mathrm{C}1}+p_{\mathrm{D}2}p_{\mathrm{C}2}\right)\;,$$ (34)

where both junctions have the same junction energy $E_{J}$, and $\omega_{B}$ (resp: $\omega_{C}$) 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 $\chi_{\mathrm{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,

$$\chi_{\mathrm{BC}}/\chi_{\mathrm{DC}}=\omega_{B}/\omega_{D}\;.$$ (35)

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 $\vec{d}_{\mathrm{B}}\cdot\vec{E}\gg 0$, we can expect the bright qubit to potentially be Purcell limited. In contrast, since $\vec{d}_{\mathrm{D}}\cdot\vec{E}=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 1 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, $Q_{C}$, 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

$$\hat{\overline{H}}/\hbar=\omega_{D}\hat{n}_{D}+\omega_{B}\hat{n}_{B}+\omega_{C}\hat{n}_{C}\\ -\frac{1}{2}\alpha_{D}\hat{n}_{D}\left(\hat{n}_{D}-\hat{1}\right)-\frac{1}{2}\alpha_{B}\hat{n}_{B}\left(\hat{n}_{B}-\hat{1}\right)\\ -\chi_{DB}\hat{n}_{D}\hat{n}_{B}-\chi_{DC}\hat{n}_{D}\hat{n}_{C}-\chi_{BC}\hat{n}_{B}\hat{n}_{C}\;,$$ (36)

where $\hat{n}_{D},\hat{n}_{B}$, and $\hat{n}_{C}$ 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 $\kappa_{C}\mathcal{D}[\hat{a}_{C}]\rho$, where $\kappa_{C}=\omega_{C}/Q_{C}$, and $\rho$ 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, $\chi_{BC}\gg\chi_{DC}$. In contrast, samples in configuration B demonstrated near equal coupling, $\chi_{BC}\approx\chi_{DC}$. In both configurations A and B, the dark mode was Purcell protected, we calculated a Purcell coupling factor of $Q_{\text{Purcell}}^{D}\gg 10^{7}$, using the eigenmode method described in Supplementary Section D4. On the other hand, the bright mode was somewhat Purcell limited, $Q_{\text{Purcell}}^{B}\approx 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.

Figure 6: Schematic representation (not to scale) of superconducting waveguide device DTW1. Two-qubit chip is placed $\lambda/4$ from a short termination in the waveguide. Guided input-output waves are launched and monitored through the I-O port connecting to a standard SMA adapter (not shown).

III.2.2 Two-qubit, single-waveguide devices

We measured a two qubit sample inside of a waveguide. Figure 6 presents the setup, and depicts the sample, which was of the configuration B type presented in Fig. 5. The sample chip was positioned inside an aluminum WR90 waveguide. The waveguide was terminated in a short at one side, and attached to an impedance-matched SMA coupler port on the input-launcher side, which was used to drive and measure the waveguide. The chip was centered inside the cross-section of the waveguide, and placed $\lambda/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 3 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.

Figure 7: Illustration of flip-chip (2.5D) device WG1. a Depiction (not to scale) of chip stack consisting of two chips separated by a 100 $\upmu\text{m}$ vacuum gap. The inner face of each chip supports part of the pattern of a multi-layer whispering gallery-mode resonator (WGMR) resonator (Minev2013, ). The lower layer contains a 2.5D, aperture transmon qubit (Minev2016, ) embedded in the WGMR. b Zoomed-in view of the lower layer of the aperture transmon. The cross marks the location of the Josephson tunnel junction device, which connects the lower center trace to an island embedded in the lower aperture. The capacitance to the top layer significantly participates in determining qubit parameters.

III.2.3 Flip-chip (2.5D), one-qubit, one-storage-cavity, one-readout-cavity devices

We also designed a multilayer planar (Minev2016, ) (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 (Minev2013, ). 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. Minev2016, . The agreement between the measured parameters of the sample and those obtained by the EPR calculation methods are presented in Table 2.

Data availability.

Data are available from the authors on reasonable request.

Code availability.

The source code for the EPR method is open-sourced and can be found at http://github.com/zlatko-minev/pyEPR.

Acknowledgments.

We thank S.M. Girvin, R.J. Schoelkopf, S. Nigg, H. Paik, A. Blais, H.E. Türeci, F. Solgun, V. Sivak, S. Touzard, N. Frattini, S. Shankar, C. Axline, V.V. Albert, K. Chou, A. Petrescu, D. Cody, A. Eickbusch, and E. Flurin for valuable discussions, and the I. Siddiqi and B. Huard groups for using the early versions of pyEPR. This research was supported by the US Army Research Office (ARO) Grant No. W911NF-18-1-0212. Z.K.M. acknowledges partial support from the ARO (W911NF-16-1-0349). M.H.D. acknowledges partial support from the ARO (W911NF-18-1-0020) and the Air Force Office of Scientific Research (FA9550-19-1-0399). The view and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office or the US Government. The US Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

Author contributions.

Z.K.M. conceived and developed the method and the theory, performed the experiments, and analyzed the data. Z.L. and M.H.D. contributed to the theory. I.M.P. and S.O.M. contributed to the measurement and fabrication of the devices L.C. assisted with the simulations. Z.K.M. and M.H.D. wrote the manuscript. All authors provided suggestions, discussed the results and contributed to the manuscript.

Competing interests.

The authors declare no competing interests.

Supplementary Information:
Energy-participation quantization of Josephson circuits

Contents

1. I INTRODUCTION
2. II RESULTS AND DISCUSSION
   1. II.1 To quantize a simple circuit: qubit coupled to a cavity
   2. II.2 Quantizing the general Josephson system
   3. II.3 Comparison between theory and experiment
3. III METHODS
   1. III.1 Methods of the experiment
   2. III.2 Devices
4. A Theoretical foundation of the energy-participation method
   1. A1 Review of electrical circuit theory, and the Josephson junction
   2. A2 The Josephson system and its non-linear Josephson dipoles
   3. A3 Energy of the Josephson circuit and its Lagrangian
   4. A4 Eigenmodes of the linearized Josephson circuit
   5. A5 Quantizing the Josephson circuit
   6. A6 Energy-participation ratio (EPR)
   7. A7 Universal EPR properties
   8. A8 The biased Josephson system: equilibrium state in the presence of persistent currents
   9. A9 The biased Josephson system: simplifications
5. B Nonlinear interactions, effective Hamiltonians, and the EPR
   1. B1 Excitation-number-conserving interactions of weakly-nonlinear systems
   2. B2 EPR matrices and many-body interactions
   3. B3 General many-body interactions
   4. B4 The driven Josephson system: parametrically-activated interactions
6. C Finite-element electromagnetic-analysis methodology
   1. C1 Modeling the Josephson dipole
   2. C2 Calculating the EPR $p_{m}$ in the case of a single Josephson dipole
   3. C3 Calculating the EPR $p_{mj}$ in the case of multiple Josephson dipoles
   4. C4 Remarks on the finite-element eigenmode approach
7. D Dissipation budget and input-output coupling
   1. D1 Dissipation budget
   2. D2 Capacitive loss
   3. D3 Inductive loss
   4. D4 Radiative loss and input-output coupling

Supplementary Section A Theoretical foundation of the energy-participation method

Supplementary Table 1: 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
$\Phi_{0}$	$h/(2e)$	Superconducting magnetic-flux quantum—the ratio of Planck’s quantum of electromagnetic action $h$ and the charge of a Cooper pair $2e$.
$\phi_{0}$	$\Phi_{0}/\left(2\pi\right)$	Reduced flux quantum.
$\Phi_{b}\left(t\right)$	$\int_{-\infty}^{t}v_{b}\left(\tau\right)\mathrm{d}\tau$	Generalized magnetic flux of circuit branch $b$, at time $t$. The instantaneous voltage $v_{b}\left(\tau\right)$ across the terminals of branch $b$ at time $\tau$ is equivalently denoted $\dot{\Phi}_{b}\left(\tau\right)$. In general, $\Phi_{b}\left(t\right)=\int_{-\infty}^{t}v_{b}\left(\tau\right)\mathrm{d}\tau+\Phi_{b}\left(-\infty\right)$, but we take the initial flux $\Phi_{b}\left(-\infty\right)$ to be zero, corresponding to the circuit in an equilibrium state (see Sec. A8). With this convention, $\Phi_{b}$ is a deviation away from equilibrium. This variable is analogous to the elongation of a mechanical spring.
$j,\,J$		The subscript or index $j$ denotes the $j$-th Josephson dipole, where $j\in\left\{1,\,\ldots,\,J\right\}$, and $J$ denotes the total number of Josephson dipoles.
$\Phi_{j}\left(t\right)$		Generalized magnetic-flux deviation of a non-linear device, the $j$-th Josephson dipole.
$\varphi_{j}\left(t\right)$	$\Phi_{j}\left(t\right)/\phi_{0}$	Reduced magnetic-flux variable of Josephson dipole $j$.
$\mathcal{E}_{j}\left(\Phi_{j}\right)$	$\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)+\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)$	Energy function of Josephson dipole $j$. Typically split into two component: $\mathcal{E}_{j}^{\mathrm{lin}}$ associated with linear interactions and $\mathcal{E}_{j}^{\mathrm{nl}}$ associated with nonlinear ones.
$\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)$	$\frac{1}{2}E_{j}\left(\Phi_{j}/\phi_{0}\right)^{2}$	Linear component of the energy function of junction $j$, defined by the energy scale $E_{j}$. Terms linear in $\Phi_{j}$ are absent since we selected an equilibrium operating point of the circuit. The energy function $\mathcal{E}_{j}^{\mathrm{lin}}$ may never contain non-linear terms.
$\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)$	$\mathcal{E}_{j}\left(\Phi_{j}\right)-\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)$	The nonlinear component of the energy function of junction $j$. For certain situations, it may be favorable to select the partition of $\mathcal{E}_{j}$ such that $\mathcal{E}_{j}^{\mathrm{nl}}$ contains linear interactions. If possible, it is often convenient to expand $\mathcal{E}_{j}^{\mathrm{nl}}$ in a Taylor series around the circuit operating point: $\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)=E_{j}\sum_{p=3}^{\infty}c_{jp}\left(\frac{\Phi_{j}}{\phi_{0}}\right)^{p}$, where $c_{jp}$ are the dimensionless coefficients of the expansion.
$\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\left(t\right)$	$\left(\Phi_{\mathrm{t}_{1}},\,\Phi_{\mathrm{t}_{2}},\,\ldots\right)^{\intercal}$	Column vector consisting of the flux deviations of all circuit branches enclosed in the minimum spanning tree. The roman subscript $\mathrm{t}$ denotes the tree. The flux of individual tree-branches are denoted $\Phi_{\mathrm{t}_{1}},\,\Phi_{\mathrm{t}_{2}},\,\ldots$
Linearized Josephson eigenmodes
$m,\;M$		The label $m$ indexes sets associated with the eigenmodes of the linearized Hamiltonian $\mathcal{H}_{\mathrm{lin}}$. We consider $m\in\left\{1,\,\ldots,\,M\right\}$, where $M$ is the total number of eigenmodes of relevance determined by the context. That is, $M$ is either the total number of eigenmodes of $\mathcal{H}_{\mathrm{lin}}$ or of the relevant eigenmodes included in the finite-element simulation.
$\mathcal{H}_{\mathrm{full}}$	$\mathcal{H}_{\mathrm{lin}}+\mathcal{H}_{\mathrm{nl}}$	Hamiltonian function of the Josephson system. It can be partitioned into $\mathcal{H}_{\mathrm{lin}}$, a part comprised of purely linear terms (quadratic in $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$, or equivalently $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$), and $\mathcal{H}_{\mathrm{nl}}$, a part generally comprised of non-linear terms (higher-than quadratic in $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$, or equivalently $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$).
$\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\,\mathrm{\mathbf{Q}}_{m}$		Column vectors (of length $M$) whose elements are the eigenmode flux $\Phi_{m}$ (considered as the generalized position) and charge $Q_{m}$ (considered as the generalized momentum) canonical variables, associated with the Hamiltonian $\mathcal{H}_{\mathrm{lin}}$.
$\omega_{m}$		Eigenfrequency of the $m$-th mode of $\mathcal{H}_{\mathrm{lin}}$, carries a dimension of circular frequency.
$\mathrm{\mathbf{\Omega}}$	$\operatorname{Diag}\left(\omega_{1},\omega_{2},\ldots\right)^{\intercal}$	Diagonal $M\times M$ matrix comprising the eigenfrequencies.
$\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$	$\mathrm{\mathbf{E}}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$	The spanning-tree branch fluxes $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$ are a linear combination of the eigenmode $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$ fluxes. The two are related by an affine transformation given by the properly-constructed eigenvector matrix $\mathrm{\mathbf{E}}$, see Eqs. (A.38) and (A.37).
Quantum operators and energy-participation ratios
$\hat{H}_{\mathrm{lin}}$	$\sum_{m=1}^{M}\hbar\omega_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}$	Hamiltonian operator corresponding to the Hamiltonian function $\mathcal{H}_{\mathrm{lin}}$, expressed in second quantization with respect to the eigenmodes of $\mathcal{H}_{\mathrm{lin}}$. The eigenmodes (expressed in terms of $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$ and $\mathrm{\mathbf{Q}}_{m}$) define the bosonic amplitude (lowering) operators $\hat{a}_{m}$.
$\hat{\varphi}_{j}$	$\sum_{m=1}^{J}\varphi_{mj}\hat{a}_{m}+\text{H.c.}$	Operator corresponding to the reduced flux $\varphi_{j}$ of Josephson dipole $j$; $\hat{\varphi}_{j}$ admits a linear decomposition in terms of the mode operators $\hat{a}_{m}$. The coefficients of this expansion are the quantum zero-point fluctuations $\varphi_{mj}$.
$\varphi_{mj}$	$s_{mj}\sqrt{\frac{\hbar\omega_{m}}{E_{j}}p_{mj}}$	The quantum zero-point fluctuation (ZPF) of the reduced magnetic flux $\hat{\phi}_{j}$ of Josephson dipole $j$ due to mode $m$; i.e., for a given mode, the ZPF magnitude gives the non-zero standard deviation of the magnetic flux in the ground state. The amplitude of the ZPF $\varphi_{mj}$ is determined by the physical structure of mode $m$, and can be understood in terms of the energy-participation ratio $p_{mj}$ and its sign $s_{mj}$.
$p_{mj},\;p_{ml}$		Energy-participation ratios (EPR) of Josephson dipole $j$ and lossy element $l$ in mode $m$, comprised between zero and unity, $0\leq p_{mj},p_{ml}\leq 1$.
$s_{mj},\;s_{ml}$		Sign $s_{mj}$ (resp., $s_{ml}$) of the energy-participation ratio $p_{mj}$ (resp., lossy EPR $p_{ml}$). The value of the sign is either $+1$ or $-1$. See Supplementary Figure S8.

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 degrees of freedom of the element are neglectable.

Electrical theory describes the physical laws governing four basic manifestations of electricity. For each element, these are captured in the variables for voltage $v\left(t\right)$, current $i\left(t\right)$, charge $Q\left(t\right)$, and magnetic flux $\Phi\left(t\right)$, as functions of time $t$. Maxwell’s equations and conservation of charge lead to the six universal relationships presented in Eqs. (A.1)–(A.3) (feynmanLec-Vol-II, ; Pozar, ).

The time-instantaneous voltage $v\left(t\right)$ and generalized magnetic flux across an element are related by the lumped-element version of Faraday’s law of induction (Yurke1984, ; Devoret1995, ; feynmanLec-Vol-II, ; Pozar, ),

$$v\left(t\right)=\frac{\mathrm{d}}{\mathrm{d}t}\Phi\left(t\right)\quad\text{and}\quad\Phi\left(t\right)=\int_{-\infty}^{t}v\left(\tau\right)\,\mathrm{d}\tau+\Phi\left(-\infty\right)\;.$$ (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 $\vec{J}$ and the electric-field vector $\vec{E}$ are usually projected on the same orientation.

From charge conservation, it follows that the time-instantaneous current $i\left(t\right)$ and the charge $Q\left(t\right)$ having passed through the element, obey a relation similar to that of the voltage and flux,

$$i\left(t\right)=\frac{\mathrm{d}}{\mathrm{d}t}Q\left(t\right)\quad\text{and}\quad Q\left(t\right)=\int_{-\infty}^{t}i\left(\tau\right)\,\mathrm{d}\tau+Q\left(-\infty\right)\;.$$ (A.2)

All four variables have their support on a non-compact set; i.e., $i,\,v,\,Q,\,\Phi\in\left[-\infty,\infty\right]$.

Reference state and initial conditions.

In Eqs. (A.1) and (A.2), we now assume zero-valued initial conditions, $\Phi\left(-\infty\right)=Q\left(-\infty\right)=0$. In the case of a circuit frustrated by sources such that the equilibrium system state is non-zero, we can define the variables $i,\,v,\,Q$, and $\Phi$ 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 $p\left(t\right)$ delivered to an element and the total energy absorbed by the element $\mathcal{E}\left(t\right)$ are

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{E}\left(t\right)=p\left(t\right)=v\left(t\right)i\left(t\right)\quad\text{and}\quad\mathcal{E}\left(t\right)=\int_{-\infty}^{t}p\left(t\right)\,\mathrm{d}\tau\;.$$ (A.3)

Given our convention for the orientation of $i$ and $v$, power delivered to the element is positive if $p\left(t\right)$ is positive. Passive elements (i.e., non-source elements) obey $\mathcal{E}\left(t\right)-\mathcal{E}\left(-\infty\right)\geq 0$ for all $t$, and lossless elements convert all of their energy into stored electric or magnetic energy.

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 $Q\left(t\right)=Cv\left(t\right)$, where $C$ is its capacitance—a positive, real constant. The dual relationship $\Phi\left(t\right)=Li\left(t\right)$ defines the simplest inductor, where $L$ 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

$$\mathcal{E}_{\mathrm{cap}}\left(\dot{\Phi}\right)=\frac{1}{2}C\dot{\Phi}^{2}\quad\text{and}\quad\mathcal{E}_{\mathrm{ind}}\left(\Phi\right)=\frac{1}{2L}\Phi^{2}\;,$$ (A.4)

respectively (and assuming $\Phi\left(-\infty\right)=Q\left(-\infty\right)=0$). Since $C,L\geq 0$ and $\Phi\left(t\right)\in\mathbb{R}$ for all $t$, 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 (Josephson1962, ; Devoret1995, ; Girvin2014, ), characterized by the flux-controlled inductive relationship $i\left(t\right)=I_{0}\sin\left(\Phi\left(t\right)/\phi_{0}\right)$, where $\phi_{0}\coloneqq\hbar/2e$ is the reduced magnetic flux quantum and $I_{0}$ is the junction critical current. The energy function of the junction in terms of flux is

$$\mathcal{E}_{J}\left(\Phi\right)=E_{J}\left(1-\cos\left(\Phi/\phi_{0}\right)\right)\,,$$ (A.5)

where $E_{J}\coloneqq I_{0}\phi_{0}$ denotes the Josephson energy. For small deviations about $\Phi=0$, ignoring the constant energy offset,

$$\mathcal{E}_{J}\left(\Phi\right)\approx\frac{1}{2}E_{J}\left(\Phi/\phi_{0}\right)^{2}-\frac{1}{4!}E_{J}\left(\Phi/\phi_{0}\right)^{4}+\mathcal{O}\left(\Phi^{6}\right)\;.$$ (A.6)

To lowest order, the junction responds as a linear inductor with inductance $L_{J}\coloneqq E_{J}/\phi_{0}^{2}$ [compare this to Eq. (A.4)] It is useful to introduce the reduced magnetic flux $\varphi\coloneqq\Phi/\phi_{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 multi-valued function with an infinite number of branches; i.e., $\varphi=\sin^{-1}\left(i\left(t\right)/I_{0}\right)+2\pi k$ or $\varphi=\pi-\sin^{-1}\left(i\left(t\right)/I_{0}\right)+2\pi k$, where $k$ 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 $\varphi$. For deviations $\varphi-\varphi_{\mathrm{eq}}$ away from the equilibrium value $\varphi_{\mathrm{eq}}$, see also Eq. (A.64),

$$\mathcal{E}_{J}\left(\varphi\right)\approx E_{J}\bigg{[}\sin\left(\varphi_{\mathrm{eq}}\right)\left(\varphi-\varphi_{\mathrm{eq}}\right)+\frac{1}{2}\cos\left(\varphi_{\mathrm{eq}}\right)\left(\varphi-\varphi_{\mathrm{eq}}\right)^{2}\\ -\frac{1}{6}\sin\left(\varphi_{\mathrm{eq}}\right)\left(\varphi-\varphi_{\mathrm{eq}}\right)^{3}\bigg{]}+\mathcal{O}\left(\varphi^{4}\right)\;.$$ (A.7)

We can identify the differential inductance of the junction at $\varphi=\varphi_{\mathrm{eq}}$ to be $L_{J}\left(\varphi_{\mathrm{eq}}\right)=L_{J}/\cos\left(\varphi_{\mathrm{eq}}\right)$. In Sec. A8, we discuss how $\varphi$ can be taken as a deviation away from the equilibrium value $\varphi_{\mathrm{eq}}$, in effect mapping $\varphi-\varphi_{\mathrm{eq}}\mapsto\varphi$, see Eq. (A.63). Also in Sec. A8, we discuss sources terms such as the one presented by $E_{J}\sin\left(\varphi_{\mathrm{eq}}\right)\left(\varphi-\varphi_{\mathrm{eq}}\right)$. Such energy terms linear in $\varphi$ 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 $I\left(t\right)=\frac{\partial\mathcal{E}_{J}}{\partial\Phi}=I_{0}\sin\left(\varphi_{\mathrm{eq}}\right)+L_{J}^{-1}\left(\varphi_{\mathrm{eq}}\right)\Phi\left(t\right)+\mathcal{O}\left(\Phi^{2}\right)$ contains the constant term $I_{0}\sin\left(\varphi_{\mathrm{eq}}\right)$, which acts like a current source.

Relationship between the gauge-invariant superconducting phase difference of a tunnel junction and the reduced magnetic flux $\varphi$ (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 $\theta_{j}$ of the superconducting phases of the two islands (tinkham2004-book, ; Clarke2004, ). This macroscopic variable is a phase angle—a compact variable in the half-open interval $\theta\in[0,2\pi[$; in contrast with the non-compact variable $\varphi\in\left[-\infty,+\infty\right]$, 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 $\varphi$. The relationship between these two collective, macroscopic variables is $\theta=\Phi/\phi_{0}\mod 2\pi$. Although the variable $\varphi$ is non-compact, the associated wavefunction $\psi\left(\varphi\right)$ 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 $\psi_{n}\left(\varphi\right)$, with $n\in\mathbb{N}$. Quantum-mechanically, the representation of $\psi\left(\varphi\right)$ and $\psi\left(\theta\right)$ is therefore not very different since $\psi\left(\theta\right)$ is represented by the discrete rotor basis, $\psi_{k}\left(\theta\right)$, with $k\in\mathbb{Z}$.

A broken symmetry.

The compact support of $\theta$ corresponds to a symmetry that is usually broken in most circuits—that of the impossibility to distinguish between different values of $\varphi$ differing by $2\pi$. Losses associated with the junction, or coupling to other elements, such as in the RF-SQUID or fluxonium qubit (Manucharyan2009, ), render $2\pi$ turns of $\varphi$ 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 non-compact version of $\varphi$ and recover the compact version by a limit procedure (Koch2007, ). Thus, in this article, we use only gauge-invariant phases with non-compact support; i.e., $\varphi\in\left[-\infty,+\infty\right]$.

Flux-controlled inductor.

To generalize from the Josephson junction and introduce more general non-linear elements, consider the flux-controlled inductor. It is defined by an algebraic relationship $i\left(t\right)=h\left(\Phi\left(t\right),t\right)$, where $h$ is a single-valued function. If $h$ is time-invariant, such as for the case of the Josephson tunnel junction, the energy function [see Eq. (A.3)] is

$$\mathcal{E}_{\mathrm{ind}}\left(\Phi\right)=\int_{\Phi\left(-\infty\right)}^{\Phi\left(t\right)}h\left(\Phi^{\prime}\right)\,\mathrm{d}\Phi^{\prime}\;.$$ (A.8)

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 $b$-th element is associated with the instantaneous voltage $v_{b}\left(t\right)$, current $i_{b}\left(t\right)$, charge $Q_{b}\left(t\right)$, flux $\Phi_{b}\left(t\right)$, and reduced flux $\varphi_{b}\left(t\right)$. The universal relationships Eqs. (A.1)–(A.3) link the variables at each branch. Variables across different branches are linked by Kirchhoff’s laws.

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’s voltage law (KVL) is the lumped-element manifestation of the Maxwell-Faraday equation, $\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}$. By applying Stokes’ theorem to the Maxwell-Faraday equation along an oriented loop of lumped elements (and a surface associated with the closed loop), one finds a relationship valid for any closed circuit loop: the oriented sum of the fluxes along the $l$-th loop is equal to the external applied flux $\Phi_{l}^{\mathrm{ext}}\left(t\right)$ threading the loop,

$$\sum_{b\in\mathrm{loop}_{l}}\pm\Phi_{b}(t)=\Phi_{l}^{\mathrm{ext}}\left(t\right)\;,$$ (A.9)

where the sum runs over all branches $b$ that form the $l$-th loop. For a given branch $b$, 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 $\Phi_{b}$, the KVL conditions express a set of holonomic constraints that need to be eliminated in order to obtain a Lagrangian of the second kind (Landau1982, ).

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 $b$ connected to node $n$,

$$\sum_{b\in\mathrm{node}_{n}}\pm i_{b}(t)=0\;.$$ (A.10)

The negative sign is chosen for branches whose current reference direction points toward the $n$-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 $\Phi_{b}$. We can systematically choose a minimal set of independent branch fluxes using the minimum spanning-tree graph method (Devoret1995, ; Burkard2004, ; Girvin2014, ; Vool2017, ). 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 $\mathrm{t}_{1},\,\mathrm{t}_{2},\,\ldots$ The flux of the $k$-th spanning-tree branch is denoted $\Phi_{\mathrm{t}_{k}}\left(t\right)$. In subscripts, roman (resp., italic) symbols denote labels (resp., variables). The spanning-tree can be organized in a column vector

$$\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\left(t\right)\coloneqq\begin{pmatrix}\Phi_{t_{1}}\left(t\right)\\ \Phi_{t_{2}}\left(t\right)\\ \vdots\end{pmatrix}\;,$$ (A.11)

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 $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$. We define the energy functions and Lagrangian of the system in terms of $\mathrm{\mathbf{\Phi}}_{\mathrm{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 finite-element (FE) method. The FE method subdivides the physical layout of the system. Using a set of basis functions, the electromagnetic circuit is discretized (Louisell1973, ; Jin2014, ). The discretized circuit can be represented by a lumped-element 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, $Q_{m}\gg 1$, where $m$ 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 $j$-th Josepson dipole in the system is fully specified by its energy function $\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j,\mathrm{ext}}\right)$, see Eq. (A.8). The energy depends only on the magnetic flux $\Phi_{j}\left(t\right)$ across the terminals of the dipole and on any external parameters $\Phi_{j,\mathrm{ext}}$ that control the energy landscape; these can include a voltage bias, a current bias, and an external magnetic field bias. The index $j$ runs from unit to the total number $J$ 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 sub-circuits. 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 $\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j,\mathrm{ext}}\right)$ into a linear and non-linear part,

$$\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j,\mathrm{ext}}\right)=\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j};\Phi_{j,\mathrm{ext}}\right)+\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j};\Phi_{j,\mathrm{ext}}\right)\;.$$ (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 $\Phi_{j}$ as deviations away from the equilibrium [discussed in more detail in Sec. (A8)], the partitions take the concrete form

$$\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right) \coloneqq \frac{1}{2}E_{j}\left(\Phi_{j}/\phi_{0}\right)^{2}\;,$$ (A.13a)

$$\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right) \coloneqq \mathcal{E}_{j}\left(\Phi_{j}\right)-\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)$$ (A.13b)

$$= E_{j}\sum_{p=3}^{\infty}c_{jp}\left(\Phi_{j}/\phi_{0}\right)^{p}\;,$$ (A.13c)

where $E_{j}$ 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 $\mathcal{E}_{j}^{\mathrm{nl}}$ around the equilibrium state of the circuit and have introduced the dimensionless expansion coefficients $c_{jp}$. 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 $\mathcal{E}_{j}^{\mathrm{lin}},\mathcal{E}_{j}^{\mathrm{nl}},E_{j}$, and $c_{jp}$ on the external bias parameter $\Phi_{j,\mathrm{ext}}$ is made implicit, and we will continue to do so henceforth; in other words, keep in mind that $E_{j}\coloneqq E_{j}\left(\Phi_{j,\mathrm{ext}}\right)$ and $c_{jp}\coloneqq c_{jp}\left(\Phi_{j,\mathrm{ext}}\right)$.

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

$$\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right) \coloneqq \frac{1}{2}E_{j}\left(\Phi_{j}/\phi_{0}\right)^{2}\;,$$ (A.14a)

$$\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right) \coloneqq -E_{j}\left[\cos\left(\Phi_{j}/\phi_{0}\right)+\frac{1}{2}\left(\Phi_{j}/\phi_{0}\right)^{2}\right]\;,$$ (A.14b)

where $E_{j}$ is the Josephson energy. The energy function $\mathcal{E}_{j}^{\mathrm{lin}}$ is associated with the linear response of the junction. It presents the inductance $L_{j}=\phi_{0}^{2}/E_{j}$. The energy $\mathcal{E}_{j}^{\mathrm{nl}}$ is associated with the response of non-linear inductor. The expansion coefficient of $\mathcal{E}_{j}^{\mathrm{nl}}$ as defined in Eq. (A.13c) are

$$c_{jp}=\begin{cases}\frac{\left(-1\right)^{p/2+1}}{p!}&\text{for even }p\;,\\ 0&\text{for odd }p\;.\end{cases}$$ (A.15)

In partitioning $\mathcal{E}_{j}$, we included all of the linear response of the junction in $\mathcal{E}_{j}^{\mathrm{lin}}$, leaving none for $\mathcal{E}_{j}^{\mathrm{nl}}$; i.e., $\mathcal{E}_{j}^{\mathrm{nl}}$ lacks quadratic terms in $\Phi_{j}$. However, this is not required. There are certain cases for which retaining some part of the linear response in $\mathcal{E}_{j}^{\mathrm{nl}}$ is advantageous.

For notational ease, we now introduce the reduced flux $\varphi_{j}\left(t\right)\coloneqq\Phi\left(t\right)/\phi_{0}$.

Example of a Josephson dipole in a frustrated circuit.

Imagine a Josephson tunnel junction incorporated in a closed loop of several circuit elements. Suppose the loop supports the flow of a direct current. A current source in the path of the loop (or perhaps an external magnetic flux threading the loop) establishes a persistent current. The equilibrium flux $\varphi_{j}$ of the junction shifts to a non-zero equilibrium value $\varphi_{\mathrm{eq},j}$, 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 $\varphi_{j}$,

$$\mathcal{E}_{j}\left(\varphi_{j};\varphi_{\mathrm{eq},j}\right)=\frac{1}{2}E_{j}\left(\varphi_{\mathrm{eq},j}\right)\varphi_{j}^{2}\\ +E_{j}\left(\varphi_{\mathrm{eq},j}\right)c_{j3}\left(\varphi_{\mathrm{eq},j}\right)\varphi_{j}^{3}+\cdots\;,$$ (A.16)

where $E_{j}\left(\varphi_{\mathrm{eq},j}\right)=E_{J}/\cos\left(\varphi_{\mathrm{eq},j}\right)$, $E_{J}$ is the $j$-th junction Josephson energy, and $c_{j3}\left(\varphi_{\mathrm{eq},j}\right)=-\frac{1}{6}\sin\left(\varphi_{\mathrm{eq},j}\right)$. We emphasize that $\varphi_{j}$ denotes deviations away from the equilibrium; compare Eq. (A.16) to Eq. (A.7).

A3 Energy of the Josephson circuit and its Lagrangian

Capacitive energy.

The total capacitive energy $\mathcal{E}_{\mathrm{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, $\mathcal{E}_{\mathrm{cap}}\coloneqq\sum_{b\in\mathrm{cap.}}\frac{1}{2}C_{b}\dot{\Phi}_{b}^{2}$, where $C_{b}$ is the capacitance of branch $b$. Each of these fluxes can be expressed in terms of the linearly-independent spanning-tree fluxes $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$. The energy function is quadratic in $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$,

$$\mathcal{E}_{\mathrm{cap}}\left(\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right)=\frac{1}{2}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{C}}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\;,$$ (A.17)

where $\mathrm{\mathbf{C}}$ is the capacitance matrix of the circuit (Yurke1984, ; Devoret1995, ; Girvin2014, ). It follows from KVL and the constitutive relationships of the capacitors that $\mathrm{\mathbf{C}}$ is a positive-definite, real, symmetric (PDRS) matrix. In the continuous limit of space, the total capacitive energy $\mathcal{E}_{\mathrm{cap}}$ can be found using Eq. (C.3).

Inductive energy.

The total inductive energy in the circuit $\mathcal{E}_{\mathrm{ind}}$ is similarly the algebraic sum of the total energy of all circuit inductive branches; i.e., $\mathcal{E}_{\mathrm{ind}}\coloneqq\sum_{b\in\mathrm{ind.}}\mathcal{E}_{b}\left(\Phi_{b}\right)$. In the Josephson system, inductive branches come in two distinct flavors: linear and non-linear. Physically, linear inductive branches are associated with the geometry and magnetic fields. We denote their total energy $\mathcal{E}_{\mathrm{mag}}$. The non-linear inductive branches (Josephson dipoles) are generally associated with the kinetic inductance of electrons. We denote their total energy $\mathcal{E}_{\mathrm{kin}}$. Hence,

$$\mathcal{E}_{\mathrm{ind}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)=\mathcal{E}_{\mathrm{mag}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)+\mathcal{E}_{\mathrm{kin}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)\;.$$ (A.18)

The energy of the linear branches $\mathcal{E}_{\mathrm{mag}}$ is the dual of Eq. (A.17). It is also a quadratic form,

$$\mathcal{E}_{\mathrm{mag}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)=\frac{1}{2}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{L}}{}_{\text{mag}}^{-1}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\;,$$ (A.19)

where the inductance matrix $\mathrm{\mathbf{L}}{}_{\text{mag}}^{-1}$ completely describes all linear, magnetic-in-origin inductances in the circuit (Yurke1984, ; Devoret1995, ; Girvin2014, ). Due to its nature, $\mathrm{\mathbf{L}}{}_{\text{mag}}^{-1}$ is PDRS. In the continuous limit of space, the total magnetic inductive energy $\mathcal{E}_{\mathrm{mag}}$ can be found using Eq. (C.4).

The total inductive kinetic energy of the circuit, associated with the non-linear dipoles, is

$$\mathcal{E}_{\mathrm{kin}}=\sum_{j=1}^{J}\mathcal{E}_{j}\left(\Phi_{j}\right)=\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{lin}}\left(\Phi_{j}\right)+\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)\;.$$ (A.20)

This energy is not stored in the magnetic fields. However, we can group the magnetic energy and the linear part of $\mathcal{E}_{\mathrm{kin}}$ together to express the inductive energy as a partition of linear and non-linear contributions

$$\mathcal{E}_{\mathrm{ind}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right) =\frac{1}{2}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{L}}{}^{-1}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}+\mathcal{E}_{\mathrm{kin}}(\mathrm{\mathbf{\Phi}}_{\mathrm{t}})\;,$$ (A.21a)

$$\mathrm{\mathbf{L}}{}^{-1} \coloneqq\mathrm{\mathbf{L}}{}_{\text{mag}}^{-1}+\frac{1}{2}\sum_{j=1}^{J}E_{j}\left(\Phi_{j}/\phi_{0}\right)^{2}\;,$$ (A.21b)

$$\mathcal{E}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right) \coloneqq\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)\;,$$ (A.21c)

where we have introduced the total inductance matrix of the circuit $\mathrm{\mathbf{L}}{}^{-1}$ and the total non-linear energy function of the circuit $\mathcal{E}_{\mathrm{nl}}$. For later use, we can obtain the series expansion of $\mathcal{E}_{\mathrm{nl}}$ in terms of that of $\mathcal{E}_{j}^{\mathrm{nl}}$, defined in Eq. (A.13b),

$$\mathcal{E}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)=\sum_{j=1}^{J}\sum_{p=3}^{\infty}E_{j}c_{jp}\left(\Phi_{j}/\phi_{0}\right)^{p}\;.$$ (A.22)

Generalized coordinates for the Lagrangian.

We have expressed $\mathcal{E}_{\mathrm{cap}}$ and $\mathcal{E}_{\mathrm{ind}}$ in terms of the independent set of spanning-tree fluxes $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$. We could equivalently have expressed $\mathcal{E}_{\mathrm{cap}}$ and $\mathcal{E}_{\mathrm{ind}}$ in terms of the charge variables $Q_{b}$. However, as discussed in Sec. A1, the flux-controlled Josephson dipoles present an asymmetry which favors treatment in the flux basis. We hence employ $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$ as the generalized position coordinates in the Lagrangian description of the circuit, and $\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}$ as the generalized velocity.

Lagrangian of the Josephson circuit.

The Lagrangian function of the Josephson circuit follows from KCL. Energy functions with $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$ as their argument (resp., $\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}$) play the role of potential (resp., kinetic) energies. The system Lagrangian is the difference of the total kinetic and potential energy functions,

$$\mathcal{L}_{\mathrm{full}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}},\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right) \coloneqq\mathcal{E}_{\mathrm{cap}}\left(\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right)-\mathcal{E}_{\mathrm{ind}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)$$ (A.23a)

$$=\mathcal{\mathcal{L}_{\mathrm{lin}}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}},\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right)+\mathcal{L}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)\;,$$ (A.23b)

where we have partitioned the Lagrangian into a linear $\mathcal{L}_{\mathrm{lin}}$ and nonlinear $\mathcal{L}_{\mathrm{nl}}$ part. Substituting in Eqs. (A.17) and (A.21),

$$\mathcal{L}_{\mathrm{lin}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}},\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right) = \frac{1}{2}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{C}}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}-\frac{1}{2}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{L}}^{-1}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\;,$$ (A.24a)

$$\mathcal{L}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right) = -\mathcal{E}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)=-\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)\;.$$ (A.24b)

Equation (A.24) explicitly constructs the Josephson system Lagrangian. It partitions it into a linear and nonlinear part. In Sec. A4, we diagonalize $\mathcal{L}_{\mathrm{lin}}$ to find the eigenmodes of the linearized system. In Sec. B, we treat the effect of $\mathcal{L}_{\mathrm{nl}}$.

A4 Eigenmodes of the linearized Josephson circuit

In this sub-sectino, we diagonalize $\mathcal{L}_{\mathrm{lin}}$ and find its eigenmodes, eigenfrequencies $\omega_{m}$ 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 $\mathcal{L}_{\mathrm{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 $\mathcal{L}_{\mathrm{lin}}$ is the sum of two quadratic forms, see Eq. (A.24). We use the standard method for their simultaneous diagonalization (Landau1982, ), 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 $\mathrm{\mathbf{L}}^{-1}$ is a PDRS matrix, we can diagonalize it with a real orthogonal matrix $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}$, obeying $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\intercal}=\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\intercal}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}=\mathrm{\mathbf{I}}$, where $\mathrm{\mathbf{I}}$ is the identity matrix,

$$\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\mathrm{T}}\mathrm{\mathbf{L}}^{-1}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}=\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{-1}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{L}}}^{-1}\;,$$ (A.25)

where $\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{-1}$ is a diagonal matrix comprising the (dimensionless) eigenvalue magnitudes, and $\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathrm{\mathbf{L}}}}$ is the identity matrix with physical dimensions of inductance. The eigenvectors of $\mathrm{\mathbf{L}}^{-1}$ form the columns of $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}$.

Employing Eq. (A.25) with Eq. (A.24), the system Lagrangian takes the suggestive form

$$\mathcal{L}_{\mathrm{full}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}},\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}\right)=\frac{1}{2}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{C}}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{t}}+\mathcal{L}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)\\ -\frac{1}{2}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}^{\mathrm{T}}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{-1/2}\right)\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathrm{\mathbf{L}}}}^{-1}\left(\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{-1/2}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\mathrm{T}}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)\;,$$ (A.26)

which motivates the principle-axis transformation of the magnetic flux defined by

$$\breve{\mathrm{\mathbf{\Phi}}}\coloneqq\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{-1/2}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\mathrm{T}}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\;,$$ (A.27)

where $\breve{\mathrm{\mathbf{\Phi}}}$ is a rotated and then scaled version of $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$. Under this transformation, the transformed Lagrangian function becomes

$$\begin{split}\breve{\mathcal{L}}_{\mathrm{full}}\left(\breve{\mathrm{\mathbf{\Phi}}},\dot{\breve{\mathrm{\mathbf{\Phi}}}}\right)\coloneqq\frac{1}{2}\dot{\breve{\mathrm{\mathbf{\Phi}}}}^{\mathrm{T}}\mathrm{\mathbf{\breve{C}}}\dot{\breve{\mathrm{\mathbf{\Phi}}}}-\frac{1}{2}\breve{\mathrm{\mathbf{\Phi}}}^{\mathrm{T}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathrm{\mathbf{L}}}}^{-1}\breve{\mathrm{\mathbf{\Phi}}}+\breve{\mathcal{L}}_{\mathrm{nl}}\left(\breve{\mathrm{\mathbf{\Phi}}}\right)\;,\end{split}$$ (A.28)

where $\breve{\mathcal{L}}_{\mathrm{nl}}\left(\breve{\mathrm{\mathbf{\Phi}}}\right)\coloneqq\mathcal{L}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}\left(\breve{\mathrm{\mathbf{\Phi}}}\right)\right)$ and

$$\mathrm{\mathbf{\breve{C}}}\coloneqq\left(\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{1/2}\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}^{\mathrm{T}}\right)\mathrm{\mathbf{C}}\left(\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{1/2}\right)\;.$$ (A.29)

Since the capacitance matrix $\mathrm{\mathbf{C}}$ is PDRS and it is transformed by a rotation and then a dilation, it follows that $\mathrm{\mathbf{\breve{C}}}$ 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.

Since $\mathrm{\mathbf{\breve{C}}}$ is PDRS, we diagonalize it with a real, orthogonal transformation $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{\breve{C}}}}$, such that

$$\mathrm{\mathbf{O}}_{\mathrm{\mathbf{\breve{C}}}}^{\mathrm{T}}\mathrm{\mathbf{\mathrm{\mathbf{\breve{C}}}}}\mathrm{\mathbf{O}}_{\breve{\mathrm{\mathbf{C}}}}=\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{C}}}\;,$$ (A.30)

where $\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}$ is the dimensionless, diagonal matrix constructed from the eigenvalues of $\mathrm{\mathbf{\mathrm{\mathbf{\breve{C}}}}}$ and $\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{C}}}$ is the identity matrix with physical dimensions of capacitance.

Employing the orthogonality transformation $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{\breve{C}}}}^{\mathrm{T}}$ in a manner similar to the one used with $\mathrm{\mathbf{O}}_{\mathrm{\mathbf{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

$$\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\coloneqq\mathrm{\mathbf{O}}_{\breve{\mathrm{\mathbf{C}}}}^{\mathrm{T}}\breve{\mathrm{\mathbf{\Phi}}}\;,$$ (A.31)

in terms of which the Lagrangian $\mathcal{L}_{\mathrm{lin}}$ is diagonal,

$$\tilde{\mathcal{L}}_{\mathrm{full}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}\right)\coloneqq\frac{1}{2}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}^{\mathrm{T}}\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{C}}}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}-\frac{1}{2}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}^{\mathrm{T}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{L}}}^{-1}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\\ +\tilde{\mathcal{L}}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\right)\;,$$ (A.32)

where the nonlinear part under the transformation is

$$\tilde{\mathcal{L}}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\right)\coloneqq\mathcal{L}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\right)\right)\;.$$ (A.33)

Equations of motion.

The Lagrangian equations of motion, $\frac{\partial\tilde{\mathcal{L}}}{\partial\mathrm{\mathbf{\Phi}}_{\mathrm{m}}}-\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial\tilde{\mathcal{L}}}{\partial\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}}=\mathrm{\mathbf{0}}$, yield the harmonic eigenvalue equation $\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}+\mathrm{\mathbf{\Omega}}^{2}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}=\mathrm{\mathbf{0}}$, where $\mathrm{\mathbf{0}}$ is the column vector of all zero elements and $\mathrm{\mathbf{\Omega}}^{2}\coloneqq\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}^{-1}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{\omega}}}^{2}$. The identity matrix $\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{\omega}}}$ has physical dimensions of circular frequency. The generalized momentum canonical to $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$ is the vector of charge variables $\mathrm{\mathbf{Q}}_{\mathrm{m}}\coloneqq\frac{\partial\tilde{\mathcal{L}}}{\partial\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}}=\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{C}}}\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}$.

Diagonalized Hamiltonian of the Josephson system.

The system Hamiltonian follows from the Legendre transform on $\tilde{\mathcal{L}}_{\mathrm{full}}$, $\mathcal{H}_{\mathrm{full}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right)=\left(\dot{\mathrm{\mathbf{\Phi}}}_{\mathrm{m}}(\mathrm{\mathbf{Q}}_{\mathrm{m}})\right)^{\mathrm{T}}\mathrm{\mathbf{Q}}_{\mathrm{m}}-\tilde{\mathcal{L}}_{\mathrm{full}}$, which we can partition into

$$\mathcal{H}_{\mathrm{full}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right)=\mathcal{H}_{\mathrm{lin}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right)+\mathcal{H}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right)\;,$$ (A.34)

where the linear and nonlinear parts of the Hamiltonian expressed in the eigenmode coordinates are

$$\mathcal{H}_{\mathrm{lin}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right) \coloneqq\frac{1}{2}\mathrm{\mathbf{Q}}_{\mathrm{m}}^{\mathrm{T}}\mathrm{\mathbf{\Omega}}^{2}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{L}}}\mathrm{\mathbf{Q}}_{\mathrm{m}}+\frac{1}{2}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}^{\mathrm{T}}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{L}}}^{-1}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$$ (A.35a)

$$=\sum_{m=1}^{M}\frac{1_{L}}{2}\omega_{m}^{2}Q_{m}^{2}+\frac{1}{2}1_{L}\Phi_{m}^{2}\;,$$ (A.35b)

$$\mathcal{H}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}},\mathrm{\mathbf{Q}}_{\mathrm{m}}\right) \coloneqq-\tilde{\mathcal{L}}_{\mathrm{nl}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\right)\;,$$ (A.35c)

where the diagonal eigenfrequency matrix $\mathrm{\mathbf{\Omega}}$ of the Hamiltonian $\mathcal{H}_{\mathrm{lin}}$ is

$$\mathrm{\mathbf{\Omega}} \coloneqq\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{\breve{C}}}}^{-1/2}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{\mathbf{\omega}}}=\begin{pmatrix}\omega_{1}\\ &\ddots\\ &&\omega_{M}\end{pmatrix}\;,$$ (A.36)

and $1_{L}$ is unity carrying physical dimensions of inductance. The entries of $\mathrm{\mathbf{\Omega}}$ are the eigenmode frequencies $\omega_{m}$ of the linearized circuit. These correspond to the eigenfrequencies solved for by the FE eigenanalysis.

Eigenvectors of the Josephson system.

The eigenvector matrix $\mathrm{\mathbf{E}}$ relates the spanning-tree fluxes $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$ to the eigenmode ones $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}$,

$$\mathrm{\mathbf{\Phi}}_{\mathrm{t}}=\mathrm{\mathbf{E}}\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\;.$$ (A.37)

It is found by concatenating the principle-axis transformations defined in Eqs. (A.27) and (A.31),

$$\mathrm{\mathbf{E}}\coloneqq\mathrm{\mathbf{O}}_{\mathrm{\mathbf{L}}}\mathrm{\mathbf{\Lambda}}_{\mathrm{\mathbf{L}}}^{1/2}\mathrm{\mathbf{O}}_{{\breve{\mathrm{\mathbf{C}}}}}\;.$$ (A.38)

The eigenvector matrix $\mathrm{\mathbf{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, $\mathrm{\mathbf{E}}\mathrm{\mathbf{E}}^{\intercal}=\mathrm{\mathbf{L}}\mathrm{\mathbf{I}}_{\mathrm{H}}^{-1}$ and $\left(\mathrm{\mathbf{E}}^{-1}\right)^{\intercal}\mathrm{\mathbf{E}}^{-1}=\mathrm{\mathbf{L}}^{-1}\mathrm{\mathbf{I}}_{\mathrm{H}}$. The eigenvector matrix $\mathrm{\mathbf{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 $\mathcal{H}_{\mathrm{full}}$ using Dirac’s canonical approach (Dirac1982-Book, ). Before taking passage from classical to quantum, we introduce the complex mode amplitude operator $\alpha_{m}$—the classical analog of the bosonic amplitude operator $\hat{a}_{m}$. This provides a direct path to second quantization in the eigenmode basis of $\mathcal{H}_{\mathrm{lin}}$.

Complex action-angle variables.

We define the vector of action-angle variables $\bm{\alpha}=\left(\alpha_{1},\ldots,\alpha_{M}\right)^{\mathrm{T}}$ by the non-canonical, complex transformation

$$\bm{\alpha}\left(t\right)\coloneqq\frac{1}{\sqrt{2\hbar\mathrm{\mathbf{\Omega}}}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{m}}\left(t\right)1_{\mathrm{H}}^{-1/2}+i\mathrm{\mathbf{\Omega}}\mathbf{Q}_{\mathrm{m}}\left(t\right)1_{\mathrm{H}}^{1/2}\right)\;,$$ (A.39)

where $1_{\mathrm{H}}$ is unity with dimensions of inductance. The normalization $1/\sqrt{2\hbar\mathrm{\mathbf{\Omega}}}$ is chosen so that the Poisson bracket of the action-angles is $\left\{\alpha_{m},\alpha_{m^{\prime}}^{*}\right\}_{\mathrm{P}}=1/(i\hbar)\delta_{mm^{\prime}}$. In terms of the action-angles, the Hamiltonian remain diagonal,

$$\mathcal{H}_{\mathrm{lin}}=\frac{\hbar}{2}\left(\bm{\alpha}^{T}\mathrm{\mathbf{\Omega}}\bm{\alpha}^{*}+\bm{\alpha}^{*T}\mathrm{\mathbf{\Omega}}\bm{\alpha}\right)\;.$$ (A.40)

We have symmetrized $\mathcal{H}_{\mathrm{lin}}$ to avoid operator order ambiguity. The flux is $\mathrm{\mathbf{\Phi}}_{\mathrm{m}}=\sqrt{\frac{\hbar\mathrm{\mathbf{\Omega}}1_{\mathrm{H}}}{2}}\left(\bm{\alpha}^{*}+\bm{\alpha}\right)$.

Quantizing the action-angle variables.

Following Dirac’s prescription (Dirac1982-Book, ), we supplant Poisson brackets by commutators and promote observables to operators. The action-angles $\alpha_{m}$ and $\alpha_{m}^{*}$ promote into the ladder annihilation $\hat{a}_{m}$ and creation $\hat{a}_{m}^{\dagger}$ operators. Their commutator follows from the Poisson bracket (Louisell1961, ; Yurke1984, ; Devoret1995, ), $i\hbar\times\{\alpha_{m},\alpha_{m}^{*}\}\mapsto\left[\hat{a}_{m},\hat{a}_{m}^{\dagger}\right]$,

$$\left[\hat{a}_{m},\hat{a}_{k}^{\dagger}\right]=\delta_{mk}\hat{I}\;,$$ (A.41)

where $\left[\hat{A},\hat{B}\right]\coloneqq\hat{A}\hat{B}-\hat{B}\hat{A}$ is the commutator, $\delta_{mk}$ is the Kronecker delta function, and $\hat{I}$ is the identity operator. Particles of the electromagnetic field (photons) are distinguishable bosons (Devoret1995, ); symmetrization of the many-body wavefunction is not required. Operators are in the Schrödinger picture, unless otherwise indicated by an explicit time argument.

Hamiltonian.

The quantized form of $\mathcal{H}_{\mathrm{lin}}$, see Eq. (A.40), is

$$\boxed{\hat{H}_{\mathrm{lin}}=\sum_{m=1}^{M}\hbar\omega_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}\;.}$$ (A.42)

This is Eq. (13) of the main text. The quantized form of $\mathcal{H}_{\mathrm{nl}}$ follows from combining Eqs. (A.21), (A.22), (A.24b) and (A.35c),

$$\boxed{\hat{H}_{\mathrm{nl}}=\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}\left(\hat{\varphi}_{j}\right)=\sum_{j=1}^{J}\sum_{p=3}^{\infty}E_{j}c_{jp}\hat{\varphi}_{j}^{p}\;.}$$ (A.43)

where $\hat{\varphi}_{j}\coloneqq\hat{\Phi}_{j}/\phi_{0}$ is the reduced magnetic-flux operator for Josephson dipole $j$. This is Eq. (14) of the main text. In the following, we decompose $\hat{\varphi}_{j}$ in terms of $\hat{a}_{m}$; i.e., in second quantization with respect to the eigenmodes of $\mathcal{H}_{\mathrm{lin}}$. In writing Eq. (A.42), we have omitted the $\frac{1}{2}\hbar\omega_{m}$ 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 $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$, denotes a vector or matrix, whose elements are constants or variables. Since variable, such as $\Phi_{\mathrm{t}_{1}}$ are promoted to quantum operators, $\hat{\Phi}_{\mathrm{t}_{1}}$, we can accommodate vectors of operators in our notation with a hat symbol; e.g.,

$$\mathrm{\mathbf{\hat{\Phi}}}_{\mathrm{t}}=\begin{pmatrix}\hat{\Phi}_{\mathrm{t}_{1}}\\ \hat{\Phi}_{\mathrm{t}_{2}}\\ \vdots\end{pmatrix}\quad\text{and}\quad\mathrm{\mathbf{\hat{\Phi}}}_{\mathrm{m}}=\begin{pmatrix}\hat{\Phi}_{\mathrm{m}_{1}}\\ \vdots\\ \hat{\Phi}_{\mathrm{m}_{M}}\end{pmatrix}\;.$$ (A.44)

The spanning-tree flux operator $\hat{\Phi}_{\mathrm{t}_{k}}$ corresponding to the $k$-th spanning-tree-branch flux variable $\Phi_{\mathrm{t}_{k}}$. Similarly, the eigenmode flux operator $\hat{\Phi}_{\mathrm{m}_{k}}$ corresponding to the $k$-th eigenflux variable $\Phi_{\mathrm{m}_{k}}$.

Zero-point fluctuations of the eigenoperators.

Inverting Eq. (A.39), one finds the vectors of eigenflux and eigencharge operators,

$$\bm{\hat{\mathrm{\Phi}}}_{\mathrm{m}} \coloneqq \mathbf{\Phi_{\mathrm{m}}^{ZPF}}\left(\bm{\hat{\mathrm{a}}}^{\dagger}+\bm{\hat{\mathrm{a}}}\right)\;,$$ (A.45a)

$$\bm{\hat{\mathrm{Q}}}_{\mathrm{m}} \coloneqq i\mathbf{Q_{\mathrm{m}}^{ZPF}}\left(\bm{\hat{\mathrm{a}}}^{\dagger}-\bm{\hat{\mathrm{a}}}\right)\;,$$ (A.45b)

respectively, where the diagonal matrices of the quantum ZPF of the operators are

$$\mathbf{\Phi_{\mathrm{m}}^{ZPF}} \coloneqq \sqrt{\frac{\hbar}{2}}\mathrm{\mathbf{\Omega}}^{1/2}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{H}^{1/2}}$$

$$= \begin{pmatrix}\sqrt{\frac{\hbar\omega_{1}}{2}1_{\mathrm{H}}}\\ &\ddots\\ &&\sqrt{\frac{\hbar\omega_{1}}{2}1_{\mathrm{H}}}\end{pmatrix}\;,$$

$$\mathbf{Q_{\mathrm{m}}^{ZPF}} \coloneqq \sqrt{\frac{\hbar}{2}}\mathrm{\mathbf{\Omega}}^{-1/2}\mathbb{\mathrm{\mathbf{I}}}_{\mathrm{H}^{-1/2}}\;.$$ (A.46b)

We recall that the operator $\hat{\Phi}_{m}=\Phi_{m}^{\mathrm{ZPF}}\left(\hat{a}_{m}^{\dagger}+\hat{a}_{m}\right)$ has a zero-mean, gaussian-distributed distribution in the ground state; i.e., its mean is $\left\langle 0\middle|\hat{\Phi}_{m}\middle|0\right\rangle=0$ and its variance is $\left\langle 0\middle|\left(\hat{\Phi}_{m}\right)^{2}\middle|0\right\rangle=\left(\Phi_{m}^{\mathrm{ZPF}}\right)^{2}$. The non-zero variance is representative of the ground state energy and the quantum zero-point fluctuations of the flux. The ZPF saturate the Heisenberg uncertainty bound, $\mathbf{\Phi_{\mathrm{m}}^{ZPF}}\mathbf{Q_{\mathrm{m}}^{ZPF}}=\frac{\hbar}{2}\mathrm{\mathbf{I}}$. All mode ZPFs are positive, $\Phi_{m}^{\mathrm{ZPF}},Q_{m}^{\mathrm{ZPF}}\in\mathbb{R}_{>0}$.

Interpretation of the eigenmode ZPFs and effective mode inductances, and impedances.

The mode impedance $Z_{m}=\Phi_{m}^{\mathrm{ZPF}}/Q_{m}^{\mathrm{ZPF}}=\omega_{m}1_{\mathrm{H}}$ and ZPF $\Phi_{m}^{\mathrm{ZPF}}=\sqrt{\hbar\omega_{m}/2\times 1_{\mathrm{H}}}$ depend only on $\omega_{m}$. 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.

From abstract to physical ZPFs.

The zero-point fluctuations of the spanning-tree fluxes follow from Eqs. (A.38) and (A.46),

$$\mathbf{\Phi_{\mathrm{t}}^{ZPF}}=\mathrm{\mathbf{E}}\mathbf{\Phi_{\mathrm{m}}^{ZPF}}\;.$$ (A.47)

The $\left(k,m\right)$ element of $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}^{\mathrm{ZPF}}$ is the ZPF of the $k$-th spanning-tree branch due to mode $m$ and is either positive or negative. The overall sign is arbitrary. Hover, the relative sign between two branches $k$ and $k^{\prime}$ in the same mode determines if they are excited by the mode in-phase or out-of phase with each other. The eigenvalue matrix $\mathrm{\mathbf{E}}$ is constructed canonically so as to correctly relate the fluctuations. The eigenvector matrix of the standard product $\mathrm{\mathbf{C}}^{-1}\mathrm{\mathbf{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 $\mathrm{\mathbf{E}}$ is real, the reduced-magnetic-flux ZPF amplitudes $\varphi_{mj}$ can be chosen to be real-valued numbers; Eq. 14(c) of the main text follows,

$$\boxed{\hat{\varphi}_{j}\coloneqq\hat{\Phi}_{j}/\phi_{0}=\sum_{m=1}^{M}\varphi_{mj}\left(\hat{a}_{m}+\hat{a}_{m}^{\dagger}\right)\;.}$$ (A.48)

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 $\varphi_{mj}$.

Definition of EPR in plain english.

The EPR $p_{mj}$ of Josephson dipole $j$ in eigenmode $m$ is the fraction of inductive energy allocated to the Josephson dipole when mode $m$ 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, $\varphi_{q}^{\mathrm{ZPF}}\sim 1$, and the therefore much smaller fluctuations of the junction in the cavity mode $\varphi_{q}^{\mathrm{ZPF}}\sim 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 $\hat{a}$ and $\hat{a}^{\dagger}$, 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 $\pi/\omega_{m}$, 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,

$$\overline{\left\langle\hat{\mathcal{E}}_{\mathrm{ind}}\right\rangle}=\frac{1}{2}\overline{\left\langle\hat{H}_{\mathrm{lin}}\right\rangle}=\frac{1}{2}\sum_{m}\hbar\omega_{m}\overline{\left\langle\hat{a}_{m}^{\dagger}\hat{a}_{m}\right\rangle}\;.$$ (A.49)

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,

$$p_{mj}\coloneqq\overline{\left<\hat{\mathcal{E}}_{j\mathrm{,lin}}\right>}/\overline{\left<\hat{\mathcal{E}}_{\mathrm{ind}}\right>}\;,$$ (A.50)

where the over-line denotes a time average and the expectation value is taken over a state with an excitation only in mode $m$. 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. (A.49) and (A.42), we simplify the denominator, $\overline{\left<\hat{\mathcal{E}}_{\mathrm{ind}}\right>}=\sum_{m}\frac{1}{2}\hat{a}_{m}^{\dagger}\hat{a}_{m}$. The operator for the junction energy follows from Eq. (A.13a), $\hat{\mathcal{E}}_{j\mathrm{,lin}}=\frac{1}{2}E_{j}\left(\hat{\varphi}_{j}^{2}-\left\langle\hat{\varphi}_{j}^{2}\right\rangle_{0}\right)$, where $\left\langle\,\right\rangle_{0}$ indicates an expectation value over the ground state. Thus,

$$p_{mj}=\frac{\overline{\left\langle\frac{1}{2}E_{j}\hat{\varphi}_{j}^{2}\right\rangle}-\overline{\left\langle\frac{1}{2}E_{j}\hat{\varphi}_{j}^{2}\right\rangle_{0}}}{\frac{1}{2}\sum_{m}\hbar\omega_{m}\left\langle\hat{a}_{m}^{\dagger}\hat{a}_{m}\right\rangle}\;.$$ (A.51)

Using Eq. (A.48),

$$\hat{\varphi}_{j}^{2}-\left\langle\hat{\varphi}_{j}^{2}\right\rangle_{0}=\sum_{m=1}^{M}\varphi_{mj}^{2}\left(2\hat{a}_{m}^{\dagger}\hat{a}_{m}+\hat{a}_{m}^{2}+\hat{a}_{m}^{\dagger 2}\right)+\\ \sum_{m^{\prime}\neq m}^{M}\varphi_{mj}\varphi_{m^{\prime}j}\left(\hat{a}_{m}\hat{a}_{m^{\prime}}+\hat{a}_{m}\hat{a}_{m^{\prime}}^{\dagger}+\hat{a}_{m}^{\dagger}\hat{a}_{m^{\prime}}+\hat{a}_{m}^{\dagger}\hat{a}_{m^{\prime}}^{\dagger}\right)\;.$$ (A.52)

EPR for a Fock state excitation.

We can take expectation value in Eq. (A.51) for a Fock excitation of $n$ photons in mode $m$, denoted $\left|n_{m}\right\rangle$. We recall that the many-body vacuum state is $\left|0\right\rangle\coloneqq\left|\mathrm{vac}\right\rangle\coloneqq\prod_{k=1}^{M}\otimes\left|0\right\rangle_{k}$, where $\left|0\right\rangle_{k}$ denotes the single-particle, zero Fock state of mode $k$; thus, $\left|n_{m}\right\rangle=\frac{1}{\sqrt{n_{m}!}}\left(\hat{a}_{m}^{\dagger}\right)^{n_{m}}\left|\mathrm{vac}\right\rangle$. To take the expectation value of the Josephson dipole flux $\left\langle\hat{\varphi}_{j}^{2}\right\rangle\leavevmode\nobreak\$, we recall that Fock states have trivial time dynamic, and use Eq. (A.52). We observe that the expectation value of any product of two amplitude operators from different modes is zero. The only non-zero term in Eq. (A.52) is the $\hat{a}_{m}^{\dagger}\hat{a}_{m}$ term. Thus, the EPR of the Josephson dipole for the Fock state $\left|n_{m}\right\rangle$ is

$$p_{mj}=\frac{E_{j}\varphi_{mj}^{2}}{\frac{1}{2}\hbar\omega_{m}}\;.$$ (A.53)

The EPR is independent of the excitation amplitude $n_{m}$.

EPR for a coherent state excitation.

A coherent state excitation of mode $m$ is denoted $\left|\beta_{m}\right\rangle\coloneqq\hat{D}\left(\beta_{m},\hat{a}_{m}\right)\left|\mathrm{vac}\right\rangle$, where $\hat{D}$ is the displacement operator of the $m$th mode, $\hat{D}\left(\beta_{m},\hat{a}_{m}\right)\coloneqq\exp\left(\beta_{m}\hat{a}_{m}^{\dagger}-\beta_{m}^{\ast}\hat{a}_{m}\right)$, and $\beta$ is a non-zero complex number. The coherent state time-evolution is a rotation, $\left|\beta_{m}\left(t\right)\right\rangle=\left|\beta_{m}\left(0\right)e^{-i\omega_{m}t}\right\rangle$. 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,

$$\boxed{\varphi_{mj}^{2}=p_{mj}\frac{\hbar\omega_{m}}{2E_{j}}\;.}$$ (A.54)

In the case of a single Josephson dipole in the circuit, Eq. (A.54) reduces to Eq. (6) of the main text. The left-hand side of Eq. (A.54) gives the root-mean-square deviation of the quantum fluctuations of the reduced magnetic flux of Josephson dipole $j$, referenced to its equilibrium value. The right-hand side of Eq. (A.54) comprises three classically-known parameters: the eigenmode frequency $\omega_{m}$, the Josephson dipole energy scale $E_{j}$, and the EPR $p_{mj}$.

Sign of the EPR.

Solving Eq. (A.54) explicitly,

$$\varphi_{mj}=s_{mj}\sqrt{p_{mj}\frac{\hbar\omega_{m}}{2E_{j}}}\;,$$ (A.55)

where $s_{mj}$ is the EPR sign, $s_{mj}\in\left\{-1,+1\right\}$. In practice, the sign is calculated using Eq. (C.8). Algebraically, the sign can be found from the eigenvector matrix $\mathrm{\mathbf{E}}$. If the $j$th row of the eigenvector matrix corresponds to the $j$th spanning-tree branch flux and the $m$th column corresponds to the $m$th eigenmode, then  $s_{mj}=\operatorname{sign}\left(\left[\mathrm{\mathbf{E}}\right]_{mj}\right)$, where $\left[\mathrm{\mathbf{E}}\right]_{mj}$ is the $\left(m,j\right)$th entry of the matrix.

EPR sign: sign freedom.

The value of an individual sign $s_{mj}$ 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 $\pi$-phase shift. The sign $s_{mj}$ derives its physical meaning in relationship to other signs $s_{mj^{\prime}}$ 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 $p_{mj}$ 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 $\hat{O}$, one takes the (Wick) normal-ordered form of the operator (Gerry2005, ), denoted by a colon on either side of the operator, ${:\mathrel{\mspace{1.0mu}O\mspace{1.0mu}}:}$. Thus, ${:\mathrel{\mspace{1.0mu}\hat{a}\hat{a}^{\dagger}\mspace{1.0mu}}:}=\hat{a}^{\dagger}\hat{a}$ and ${:\mathrel{\mspace{1.0mu}\left(\hat{a}^{\dagger}\hat{a}\right)^{2}\mspace{1.0mu}}:}=\hat{a}^{\dagger 2}\hat{a}^{2}$. The normal-ordered form of the operators should not be confused with the normal-ordering transformation $\mathscr{N}$ applied to the operator, which takes the commutation relations into account; e.g., $\text{$\text{$\mathscr{N}$}$}\left[\hat{a}\hat{a}^{\dagger}\right]=\hat{a}^{\dagger}\hat{a}+\hat{I}$ and $\text{$\mathscr{N}$}\left[\left(\hat{a}^{\dagger}\hat{a}\right)^{2}\right]=\hat{a}^{\dagger 2}a^{2}+\hat{a}^{\dagger}\hat{a}$. 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:

$$p_{mj}\coloneqq\overline{\left<{:\mathrel{\mspace{1.0mu}\hat{\mathcal{E}}_{j\mathrm{,lin}}\mspace{1.0mu}}:}\right>}/\overline{\left<{:\mathrel{\mspace{1.0mu}\hat{\mathcal{E}}_{\mathrm{ind}}\mspace{1.0mu}}:}\right>}\;.$$ (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 $H_{\mathrm{lin}}$, obtained in Sec. A4 and are linked to the ZPF, as discussed in Sec. A6. In the following, $M$ 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,

$$0\leq p_{mj}\leq 1\;.$$ (A.57)

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,

$$\sum_{m=1}^{M}p_{mj}=1\quad\text{for}\ j\in\left\{1,\ldots,J\right\}\;,$$ (A.58)

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 $m$ can at most have a total EPR of unity, and no less than zero,

$$0\leq\sum_{j=1}^{J}p_{mj}\leq 1\quad\text{for}\ m\in\left\{1,\ldots,M\right\}\;.$$ (A.59)

The upper bound is saturated when there are no linear inductors excited in the mode.

The vector EPRs of two dipoles are orthogonal.

For each dipole, we can define a vector EPR by the components $\left\{s_{mj}\sqrt{p_{mj}}\,|\,m=1,\ldots,M\right\}$. The vector EPRs of two dipoles are orthogonal in the following sense

$$\sum_{m=1}^{M}s_{mj}s_{mj^{\prime}}\sqrt{p_{mj}p_{mj^{\prime}}}\, =0\quad\text{for}\ j\neq j^{\prime}\;,$$ (A.60)

where $s_{mj}$ 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

Supplementary Figure S1: Example of a Josephson dipole embedded in frustrated loop. (a) Illustration of a Josephson tunnel junction (with Josephson energy $E_{j}$) embedded in a distributed, superconducting ring subjected to an external magnetic-flux bias $\Phi_{\mathrm{ext}}$. (b) Lumped-element model a junction in an inductive ring (with inductance $L$) frustrated by an external flux bias $\Phi_{\mathrm{ext}}$, a voltage source $v_{s}$, and a current source $i_{s}$.

The magnetic flux across a Josephson dipole $\Phi_{j}$, 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 $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}$ from equilibrium guarantees that the total inductive energy $\mathcal{E}_{\mathrm{ind}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t}}\right)$ is at a minimum for $\mathrm{\mathbf{\Phi}}_{\mathrm{t}}=\mathrm{\mathbf{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, $\frac{\mathrm{d}^{2}\mathrm{\mathbf{\Phi}}_{\mathrm{t}}}{\mathrm{d}t^{2}}=\mathrm{\mathbf{0}}$. The corresponding magnetic flux of the spanning-tree branches in static equilibrium $\mathrm{\mathbf{\Phi}}_{\mathrm{t,eq}}$ can be found by extremizing the Lagrangian with respect to the generalized position variables (Landau1982, ),

$$-\frac{\partial\mathcal{L}}{\partial\mathrm{\mathbf{\Phi}}_{\mathrm{t}}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t,eq.}}\right)=\frac{\partial\mathcal{E}_{\mathrm{ind}}}{\partial\mathrm{\mathbf{\Phi}}_{\mathrm{t}}}\left(\mathrm{\mathbf{\Phi}}_{\mathrm{t,eq}}\right)=\mathrm{\mathbf{0}}\;.$$ (A.61)

Solving these conditions for the equilibrium flux $\mathrm{\mathbf{\Phi}}_{\mathrm{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 $\frac{\partial\mathcal{E}_{j}}{\partial\Phi_{j}}=0$, where $\mathcal{E}_{j}$ is the energy function of the Josephson dipole and $\Phi_{j}$ 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.

Supplementary Figure S2: Mechanical analogy of a Josephson dipole in isolation vs. in a frustrated system. (a) Depiction of a mechanical spring in isolation. Without applied forces on the spring, the native equilibrium length of the spring is $L_{\mathrm{eq}}^{0}$. Stretching or compressing the spring from equilibrium is measured by the deviation $x$ away from the native equilibrium. The spring constant is $k$. (b) Depiction of two connected springs, with spring constant $k^{\prime}$ and $k^{\prime\prime}$, constrained between two walls separated by a fixed distance $L_{\mathrm{ext}}$. The equilibrium lengths of the two springs in the system are $L_{\mathrm{eq}}^{\prime}$ and $L_{\mathrm{eq}}^{\prime\prime}$. These differ from their equilibrium lengths in isolation; i.e., $L_{\mathrm{eq}}\neq L_{\mathrm{eq}}^{\prime}$. Deviations from the system equilibrium are denoted by $x$.

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 $\mathcal{E}_{i}\left(\Phi_{i}\right)=-E_{J}\cos\left(\Phi_{i}/\phi_{0}\right)$, where $\Phi_{i}$ is the total magnetic flux across the junction in isolation. (Here, the subscript $i$ denotes the junction in isolation.) The energy is minimized at the native equilibrium-flux value $\Phi_{\mathrm{eq}}^{\mathrm{nat}}=0$. Now, consider embedding this junction in a superconducting ring frustrated by an external magnetic flux $\Phi_{\mathrm{ext}}$, as depicted in Supplementary Figure S1(a). Suppose the geometric ring inductance is $L$, see Supplementary Figure S1(b). The equilibrium condition of the circuit, given by Eq. (A.61), reduces to Kepler’s transcendental equation

$$L^{-1}\Phi_{i}+\phi_{0}^{-1}E_{J}\sin\left(\Phi_{i}/\phi_{0}\right)=i_{s}^{\mathrm{eff}}\;,$$ (A.62)

where the effective dc current sourced to the loop is $i_{s}^{\mathrm{eff}}=L^{-1}\Phi_{\mathrm{ext}}$. In the additional presence of a voltage source $v_{s}$ and current source $i_{s}$ as depicted in see Supplementary Figure S1(b), the net effective current is $i_{s}^{\mathrm{eff}}=L^{-1}\Phi_{\mathrm{ext}}+i_{s}$.

A solution $\Phi_{\mathrm{eq}}\left(\Phi_{\mathrm{ext}},\Phi_{s},i_{s}\right)$ 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 $\Phi_{\mathrm{eq}}$ is then

$$\Phi_{j}\coloneqq\Phi_{i}-\Phi_{\mathrm{eq}}\;.$$ (A.63)

The energy function of the junction now in terms of the deviation $\Phi_{j}$ is

$$\mathcal{E}_{j}\left(\Phi_{j}\right)\coloneqq-E_{J}\cos\left(\left(\Phi_{\mathrm{eq}}+\Phi_{j}\right)/\phi_{0}\right)\;.$$ (A.64)

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 $E_{J}\left(\Phi_{\mathrm{ext}},\Phi_{s},i_{s}\right)\coloneqq E_{J}\cos\left(\Phi_{\mathrm{eq}}/\phi_{0}\right)$ and the leading-order non-linear coefficient is $c_{j3}\left(\Phi_{\mathrm{ext}},\Phi_{s},i_{s}\right)=-\frac{1}{6}\tan\left(\Phi_{\mathrm{eq}}/\phi_{0}\right)$.

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 $L_{\mathrm{eq}}^{0}$. Stretching the spring to length $L$ results in the force $F=-k\left(x_{i}\right)$ on its ends, where  $x_{i}\coloneqq L-L_{\mathrm{eq}}^{0}$ denotes a deviation away from equilibrium. In general, $k\left(x_{i}\right)$ is a non-linear function; however, in equilibrium, $k\left(0\right)=0$. For a linear spring, the energy, $\mathcal{E}_{i}\left(x_{i}\right)=\frac{1}{2}kx_{i}^{2}$, is minimized for $x=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 $L^{\prime}+L^{\prime\prime}=L_{\mathrm{ext}}$, where $L^{\prime}$ is the length of the second spring, and $L_{\mathrm{ext}}$ is the total distance between the two walls. For linear springs, the length of the spring in the equilibrium of the system is $L_{\mathrm{eq}}=\left[kL_{\mathrm{eq}}^{\prime}+k^{\prime}\left(L_{\mathrm{ext}}-L_{\mathrm{eq}}^{\prime\prime}\right)\right]/\left(k+k^{\prime}\right)$, where $k^{\prime}$ and $L_{\mathrm{eq}}^{\prime\prime}$ are the spring constant and native equilibrium length of the second spring, respectively. The energy of the spring is $\mathcal{E}\left(x\right)=\frac{1}{2}kx^{2}$, where we omit constant energy terms and use $x\coloneqq L-L_{\mathrm{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 $v_{s}$ does not appear in the effective current bias $i_{s}$ 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, $\Phi_{s}=v_{s}t$, 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

Supplementary Figure S3: Schematic of five types of elementary flux-biased Josephson circuits. (a) Example of a fully linear circuit (no Josephson dipoles). A superconducting ring is frustrated by an external magnetic flux $\Phi_{\mathrm{ext}}$ and, through a mutual inductive coupling, a current source $i_{s}$. (b) Example of system incorporating a Josephson dipole (here, a Josephson tunnel junction) in a dc-conducting distributed loop frustrated by $\Phi_{\mathrm{ext}}$ and $i_{s}$. (c) Same as panel (b), but with a Superconducting Nonlinear Asymmetric Inductive eLement (SNAIL) instead of a junction in the loop. The SNAIL, a composite Josephson dipole, contains internal nodes and loops. Its internal loop is frustrated by the external magnetic flux $\Phi_{\mathrm{ext}}^{\prime}$. (d) Example of Josephson tunnel junction embedded in a distributed non-dc-conducting loop. In the region of the gap, we define the segment of the line-contour for the magnetic flux $\Phi_{\mathrm{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.

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.

A frustrated loop not incorporating Josephson dipoles.

If an external bias sets up a persistent current involving only linear parts of the circuit [no Josephson dipoles; see Supplementary Figure S3(a)], the equilibrium of the circuit does not need to be calculated, since the dc current will not affect the eigenmodes of the circuit. Put another way, the alternating-current (ac) response of a linear sub-circuit is independent of its dc configuration.

An open loop incorporating a Josephson dipole without internal loops.

If a Josephson dipole is only a member of loops in the circuit that cannot not support a dc current [see Supplementary Figure S3(d)], any external frustration or bias of the loop is impotent. Thus, the equilibrium flux of the Josephson dipole in isolation $\Phi_{\mathrm{eq}}^{\mathrm{nat}}$ can be taken as the equilibrium flux of the Josephson dipole in the circuit. In other words, we can treat the dipole locally in the steady-state analysis—i.e., we can forget about the rest of the distributed circuit in which it is embedded. This is one of the most common cases of practical interest, and is the one of the transmon and SQUID-based qubit.

An open loop incorporating a Josephson dipole with internal loops.

A Josephson dipole may have frustrated loops internal to its structure [see Supplementary Figure S3(e)]. For example, a lumped-element SQUID, RF-SQUID, or SNAIL is subjected to an external magnetic flux $\Phi_{\mathrm{ext}}^{\prime}$. However, since a Josephson dipole is purely inductive and if it is not part of a global conducting loop, the equilibrium conditions of its internal nodes are solved locally.

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 $\mathcal{E}_{\mathrm{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 (Matveev2002, ; Manucharyan2012, ; Pop2012, ; Masluk2012, ; Astafiev2012, ; Rastelli2013, ).

Supplementary Section B Nonlinear interactions, effective Hamiltonians, and the EPR

In this section, we find the amplitude of a general nonlinear interaction due to $\hat{H}_{\mathrm{nl}}$ explicitly. In Sec. B1, we find the effective interaction Hamiltonian of weakly non-linear systems in the dispersive regime to leading-order. In Sec. B2, to more systematically handle large systems, we introduce the EPR matrix $\mathrm{\mathbf{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 $\hat{H}_{\mathrm{nl}}$ 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 $\varphi_{mj}\ll 1$, for all modes $m$ and junctions j, the Hamiltonian $\hat{H}_{\mathrm{nl}}$ exerts only a perturbative effect on the eigenspectrum of $\hat{H}_{\mathrm{lin}}$. Hence, we treat its effect order-by-order using the expansion $\hat{H}_{\mathrm{nl}}=\sum_{j=1}^{J}\sum_{p=3}^{\infty}E_{j}c_{jp}\hat{\varphi}_{j}^{p}$, 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, $\hbar\left(\omega_{k}-\omega_{m}\right)\gg E_{j}c_{jp}\left|\hat{\varphi}_{j}^{p}\right|$ for all $k,m,j$, and $p$.

Example: transmon coupled to a readout cavity mode.

Recall the simple example circuit quantized at the beginning of the main text. A qubit ($q$) is coupled to a cavity ($c$). To leading order in $p$, the circuit Hamiltonian $\hat{H}_{\mathrm{nl}}$ is, using Eq. (A.43),

$$\hat{H}_{\mathrm{nl}} = -\frac{E_{J}}{24}\left(\varphi_{c}\hat{a}_{c}+\varphi_{q}\hat{a}_{q}+\mathrm{H.c.}\right)^{4}+\mathcal{O}\left(\hat{\varphi}_{J}^{6}\right)\;,$$ (B.1)

where $E_{J}$ and $\hat{\varphi}_{J}$ are the Josephson energy and flux operator of the tunnel junction, and $\varphi_{q}$ and $\varphi_{c}$ are the qubit and cavity mode ZPF, respectively.

Expanding the multinomial.

Expanding the $p=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 $\hat{a}_{q}\left(\hat{a}_{q}^{\dagger}\right)^{3}$ and $\hat{a}_{c}\hat{a}_{c}^{\dagger}\hat{a}_{q}\hat{a}_{q}^{\dagger}$. Using the commutation relations [Eq. (A.41)], we normal order the polynomial. For example, the qubit-only operators, this yields

$$\hat{a}_{q}^{4}+4\hat{a}_{q}^{\dagger}\hat{a}_{q}^{3}+6\hat{a}_{q}^{\dagger 2}\hat{a}_{q}^{2}+4\hat{a}_{q}^{\dagger 3}\hat{a}_{q}+\hat{a}_{q}^{\dagger 4}\\ +\text{$6\hat{a}_{q}^{2}+12\hat{a}_{q}^{\dagger}\hat{a}_{q}+6\hat{a}_{q}^{\dagger 2}$}+3\hat{I}\;.$$

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 $\hat{a}_{q}^{\dagger}\hat{a}_{q}$ and $3\hat{I}$. Such quadratic terms (e.g., $\hat{a}_{q}^{\dagger}\hat{a}_{q}$, $\hat{a}_{c}^{\dagger}\hat{a}_{c}$, and $\hat{a}_{c}^{\dagger}\hat{a}_{q}$) dress the modes of $\hat{H}_{\mathrm{lin}}$ in a linear manner—both renormalizing their frequencies and hybridizing them. These new linear couplings cannot be a-priori straightforwardly included in $\hat{H}_{\mathrm{lin}}$, 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 $H_{\mathrm{lin}}$ is possible under certain assumptions (Solgun2017, ).

Rotating-wave approximation in the dispersive regime.

To leading-order, under the rotating-wave approximation (RWA), only excitation-number conserving interactions in $\hat{H}_{\mathrm{nl}}$ (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 Hamiltonian

$$\begin{split}\hat{\overline{H}}_{p=\mathrm{4}}=\sum_{m}-\hbar\Delta_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}-\frac{1}{2}\hbar\alpha_{m}\hat{a}_{m}^{\dagger 2}\hat{a}_{m}^{2}\\ -\sum_{n\neq m}\frac{1}{2}\hbar\chi_{mn}\hat{a}_{m}^{\dagger}\hat{a}_{m}\hat{a}_{n}^{\dagger}\hat{a}_{n}\;,\end{split}$$ (B.2)

where $m\in\left\{q,c\right\}$ is the mode label, $\Delta_{m}$ is the effective Lamb shift, $\alpha_{m}$ is the anharmonicity, and $\chi_{mn}$ is the total cross-Kerr frequency shift induced between modes $m$ and $n$. The bar over $H$ denotes the RWA; the subscript $p=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),

$$\chi_{mn} =\sum_{j=1}^{J}\hbar^{-1}E_{j}\varphi_{mj}^{2}\varphi_{nj}^{2}=\sum_{j=1}^{J}\frac{\hbar\omega_{m}\omega_{n}}{4E_{j}}p_{mj}p_{nj}\;,$$ (B.3a)

$$\Delta_{m} =\frac{1}{2}\sum_{n=1}^{M}\chi_{mn}\;,$$ (B.3b)

$$\alpha_{m} =\frac{1}{2}\chi_{mm}\;.$$ (B.3c)

Constraints in the Hamiltonian.

First, the structure of the non-linearity imposes certain constraints. For example, $\chi_{mn}\leq 2\sqrt{\alpha_{m}\alpha_{n}}$, where the equality holds only for the one-junction case, $J=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, $\omega_{k}\approx z\times\omega_{m}$, where $z\in\mathbb{Z}$. In the case of such a degeneracy, additional terms survive the RWA; for example, if $\omega_{c}=3\omega_{q}$, then the non-excitation conserving term $\hat{a}_{q}^{3}\hat{a}_{c}^{\dagger}$ 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), $\left|\chi_{\mathrm{BC}}\right|\gg 0$, and the antisymmetric modes does not, $\left|\chi_{\mathrm{DC}}\right|=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 ($p_{D1}=p_{D2}=p_{B1}=p_{B2}$), the couplings to the cavity are essentially equal, $\chi_{\mathrm{DC}}\approx\chi_{\mathrm{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 $M$ modes under the simple extension $m\in\left\{1,\ldots,M\right\}$.

Supplementary Figure S4: Excitation spectrum (not-to-scale) of a transmon-cavity circuit, corresponding to the effective dispersive Hamiltonian $\hat{H}_{\mathrm{lin}}+\hat{\overline{H}}_{p=4}$, see Eq. (B.2). (a) The cavity and qubit frequencies in the ground state are $\omega_{c}^{\prime}$ and $\omega_{q}^{\prime}$, respectively. The qubit anharmonicity is $\alpha_{q}$, and the qubit-cavity dispersive cross-Kerr shift is $\chi_{qc}$. (b) Zoom-in on the qubit spectrum in the vicinity of $\omega_{q}^{\prime}$. Each photon in the resonator loads the qubit frequency down by $\chi_{qc}$. (c) Zoom-in on the cavity spectrum in the vicinity of $\omega_{c}^{\prime}$.

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 $\mathrm{\mathbf{P}}$, comprising the $M\times J$ EPRs $p_{mj}$ as its entries,

$$\mathrm{\mathbf{P}}\coloneqq\left(\begin{array}[]{ccc}p_{11}&\cdots&p_{1J}\\ \vdots&\ddots&\vdots\\ p_{M1}&\cdots&p_{MJ}\end{array}\right)\;.$$ (B.4)

The cross-Kerr interaction amplitudes $\chi_{mn}^{(p)}$ due to the $p$-th order terms of $\hat{H}_{\mathrm{nl}}$ [for $p=4$, see Eq. (B.3a)] can similarly be organized in an $M\times M$ matrix,

$$\bm{\chi}_{p}\coloneqq\left(\begin{array}[]{ccc}\chi_{11}^{(p)}&\cdots&\chi_{1M}^{(p)}\\ \vdots&\ddots&\vdots\\ \chi_{M1}^{(p)}&\cdots&\chi_{MM}^{(p)}\end{array}\right)\;.$$ (B.5)

It follows from Eq. (B.3a) that the leading-order Kerr matrix is

$$\boxed{\bm{\chi}_{4}=\frac{\hbar}{4}\left(\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}\right)\mathrm{\mathbf{E}}_{j}^{-1}\left(\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}\right)^{\mathrm{T}}\;,}$$ (B.6)

where $\mathrm{\mathbf{\Omega}}$ is the diagonal eigenfrequency matrix and

$$\mathrm{\mathbf{E}}_{j}^{-1}\coloneqq\begin{pmatrix}E_{1}^{-1}\\ &\ddots\\ &&E_{J}^{-1}\end{pmatrix}$$ (B.7)

is the diagonal matrix of inverse Josephson energies.

To leading order, $p=4$, the vector of anharmonicities is the diagonal of $\bm{\chi}_{4}$,

$$\bm{\alpha}_{4}\coloneqq\left(\alpha_{1}^{(4)},\ldots,\alpha_{M}^{(4)}\right)^{\mathrm{T}}=\frac{1}{2}\mathrm{diag}\,\bm{\chi}_{4}\;,$$ (B.8)

and the Lamb shift of mode $m$ is the $m$-th row sum of the Kerr matrix.

$$\mathrm{\mathbf{\Delta}}_{4}\coloneqq\left(\Delta_{1}^{(4)},\ldots,\Delta_{M}^{(4)}\right)^{\mathrm{T}}=\frac{1}{2}\bm{\chi}_{4}\mathrm{\mathbf{1}}_{M}\;,$$ (B.9)

where $\mathrm{\mathbf{1}}_{M}$ denotes a column vector of length $M$ 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 $\mathrm{\mathbf{E_{J}}}^{-1}$ are diluted by $\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}$ through a congruence transform to the Kerr coefficients. While the eigenfrequencies $\mathrm{\mathbf{\Omega}}$ weighs the contribution to each mode, the participation matrix $\mathrm{\mathbf{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 $p$-th other correction to the Kerr matrix is

$$\bm{\chi}_{p}=\hbar c_{p}\left(\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}\right)\mathrm{\mathbf{E}}_{j}^{-1}\bm{\varphi}_{\mathrm{tot}}^{p-4}\left(\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}\right)^{\mathrm{T}}\;,$$ (B.10)

where $\bm{\varphi}_{\mathrm{tot}}=\operatorname{Diag}\left(\varphi_{1,\mathrm{tot}},\ldots,\varphi_{J,\mathrm{tot}}\right)$ is the diagonal matrix of the total ZPF fluctuation of the Josephson dipole reduced fluxes, defined in Eq. (B.18). The Kerr matrix incorporating corrections due to all orders of the nonlinearity, see Eq. (B.24), is $\bm{\chi}\coloneqq\sum_{p=3}^{\infty}\bm{\chi}_{p}$. The congruence-transformation form of Eq. (B.10) is identical to that of Eq. (B.6); it governs the dilution of the nonlinearity in the same manner, subject to $\mathrm{\mathbf{\Omega}}\mathrm{\mathbf{P}}$.

Similarly to the results of Eqs. (B.8) and (B.9), the $p$-th order corrections to the anharmonicity and Lamb-shift vectors are

$$\bm{\alpha}_{p}=\frac{1}{2}\mathrm{diag}\,\bm{\chi}_{p}\quad\text{and}\quad\mathrm{\mathbf{\Delta}}_{p}=c_{p}\left(\bm{\Omega}\mathrm{\mathbf{P}}\right)\bm{\varphi}_{\mathrm{tot}}^{p-2}\bm{1}_{M}\;.$$ (B.11)

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 $\hat{H}_{\mathrm{lin}}$ due to $\hat{H}_{\mathrm{nl}}$. We expressed the eigenmode interactions using normal-ordered many-body terms, such as $\chi\hat{a}_{q}^{\dagger}\hat{a}_{c}^{\dagger}\hat{a}_{q}\hat{a}_{c}$. Here, we extend the analysis and compute the amplitude of any general normal-ordered term in $\hat{H}_{\mathrm{nl}}$.

The form a general many-body interaction.

A general normal-ordered many-body interaction has the form

$$C_{\bm{\alpha},\bm{\beta}}^{p}\hat{a}_{1}^{\dagger\beta_{1}}\cdots\hat{a}_{M}^{\dagger\beta_{M}}\hat{a}_{1}^{\alpha_{1}}\cdots\hat{a}_{M}^{\alpha_{M}}\eqqcolon C_{\bm{\alpha},\bm{\beta}}^{p}\bm{\hat{a}^{\dagger\beta}}\bm{\hat{a}^{\alpha}}\;,$$ (B.12)

where $C_{\bm{\alpha},\bm{\beta}}^{p}$ is an energy-dimensioned amplitude, determined by the $p$-th-order of $\hat{H}_{\mathrm{nl}}$ [see Eqs. (A.43) and (A.48)], and the multi-index $M$-tuples

$$\bm{\alpha}\coloneqq\left(\alpha_{1},\ldots,\alpha_{M}\right)\quad\text{and}\quad\bm{\beta}\coloneqq\left(\beta_{1},\ldots,\beta_{M}\right)$$ (B.13)

describe the distribution of annihilation and creation operators involved in the interaction among the $M$ modes, respectively. The entires of the multi-index tuples $\bm{\alpha}$ and $\bm{\beta}$ are non-negative integers, $\alpha_{m},\beta_{m}\in\mathbb{Z}_{\geq 0}$. The right-hand side of Eq. (B.12) introduces the multi-index shorthand notation for the interaction terms.

Multi-index shorthand and operator powers.

The total number of lowering and raising operators in the expression of Eq. (B.12) is equal to the 1-norm of $\bm{\alpha}$ and $\bm{\beta}$,

$$\left|\bm{\alpha}\right|\coloneqq\sum_{m=1}^{M}\alpha_{m}\quad\text{and}\quad\left|\bm{\beta}\right|\coloneqq\sum_{m=1}^{M}\beta_{m}\;,$$ (B.14)

respectively. For a given power $p$ of $\hat{H}_{\mathrm{nl}}$, the total number of resulting operators is bounded, $\left|\bm{\alpha}\right|+\left|\bm{\beta}\right|\leq p$.

Expanding the $\hat{H}_{\mathrm{nl}}$ multinomials.

To arrive at Eq. (B.12) from Eq. (A.43), we first group the sum of the annihilation operators for a Josephson dipole, see Eq. (A.48), and define

$$\hat{A}_{j} \coloneqq\sum_{m=1}^{M}\varphi_{mj}\hat{a}_{m}\;;$$ (B.15)

$$\hat{H}_{\mathrm{nl}} =\sum_{p=3}^{\infty}\sum_{j=1}^{J}E_{j}c_{jp}\left(\hat{A}_{j}+\hat{A}_{j}^{\dagger}\right)^{p}\;.$$ (B.16)

Importantly, the commutator $\left[\hat{A}_{j},\hat{A}_{j}^{\dagger}\right]=\sum_{m=1}^{M}\left|\varphi_{mj}\right|^{2}$ is scalar-valued, which allows us to use the normal-ordering non-commutative binomial theorem to expand the $p$-th power term of $\hat{H}_{\mathrm{nl}}$ (Cohen1966, ; Blasiak2005, ; Reagor2015, ). Using the non-commuting expansion,

$$\begin{split}\left(\hat{A}_{j}+\hat{A}_{j}^{\dagger}\right)^{p}=\sum_{k=0}^{\mathrm{floor}\frac{p}{2}}\sum_{i=0}^{p-2k}\frac{p!}{2^{k}k!i!(p-2k-i)!}\times\\ \left(\varphi_{j,\mathrm{tot}}^{2}\right)^{k}\left(\hat{A}_{j}^{\dagger}\right)^{i}\left(\hat{A}_{j}\right)^{p-2k-i}\;,\end{split}$$ (B.17)

where $\mathrm{floor}\frac{p}{2}$ gives the greatest integer that is less than or equal to $\frac{p}{2}$ and the total variance of the ZPF of the $j$-th dipole is

$$\varphi_{j,\mathrm{tot}}^{2}\coloneqq \left[\hat{A},\hat{A}^{\dagger}\right]=\frac{\hbar}{2}E_{j}^{-1}\sum_{m=1}^{M}p_{mj}\omega_{m}\;.$$ (B.18)

Since $\hat{A}_{j}$ is the sum of operators that commute, see Eq. (B.15), we can now expand the powers of $\hat{A}_{j}$ using the classical multinomial theorem,

$$\left(\hat{A}_{j}\right)^{i} = \sum_{\left|\bm{\alpha}\right|=i}\left(\begin{array}[]{c}\left|\bm{\alpha}\right|\\ \bm{\alpha}\end{array}\right)\varphi_{mj}^{\bm{\alpha}}\bm{\hat{a}}^{\bm{\alpha}}\;,$$ (B.21)

where the multi-index shorthand $\left(\varphi_{mj}\right)^{\bm{\alpha}}\coloneqq\prod_{m=1}^{M}\varphi_{mj}^{\alpha_{m}}$, the multinomial coefficient is

$$\left(\begin{array}[]{c}\left|\bm{\alpha}\right|\\ \bm{\alpha}\end{array}\right)\coloneqq\frac{\left|\bm{\alpha}\right|!}{\alpha_{1}!\cdots\alpha_{M}!}\;,$$ (B.22)

and the sum condition $\left|\bm{\alpha}\right|=i$ means that the sum includes terms with all possible tuples $\bm{\alpha}$ such that their 1-norm has the value $i$.

The general normal-ordered form.

Combining Eqs. (B.16)–(B.21), we find the general normal-ordered many-body form of $\hat{H}_{\mathrm{nl}}$ to all orders in $p$ and without approximations,

$$\boxed{\begin{split}\hat{H}_{\mathrm{nl}}=\sum_{p=3}^{\infty}\sum_{k=0}^{\mathrm{floor}\frac{p}{2}}\sum_{i=0}^{p-2k}\sum_{\begin{subarray}{c}\left|\bm{\beta}\right|=i,\\ \left|\bm{\alpha}\right|=p-2k-i\end{subarray}}C_{\bm{\alpha},\bm{\beta}}^{p}\bm{\hat{a}^{\dagger\bm{\beta}}}\bm{\hat{a}}^{\bm{\alpha}}\;.\end{split}}$$ (B.23)

The energy-dimensioned amplitude of an interaction due to the $p$-th power of $\hat{H}_{\mathrm{nl}}$ for a general many body mode-interaction term is the sum of the individual junction contributions,

$$\boxed{C_{\bm{\alpha},\bm{\beta}}^{p}\coloneqq\frac{p!}{\bm{\alpha}!\bm{\beta}!k!2^{k}}\,\sum_{j=1}^{J}E_{j}c_{jp}\varphi_{mj}^{\bm{\alpha}}\varphi_{mj}^{\bm{\beta}}\varphi_{j,\mathrm{tot}}^{2k}\;,}$$ (B.24)

where $k\coloneqq\frac{1}{2}\left(p-\left|\bm{\alpha}\right|-\left|\bm{\beta}\right|\right)$.

Equation (B.24) provides the amplitude for any mode interaction. Its exact value is calculated using the EPR to obtain the ZPF $\varphi_{mj}$, using Eq. (A.55). Thus, we have analytically fully constructed $\hat{H}_{\mathrm{nl}}$, and, individually, every term contained within, from the FE simulations through the EPR.

Use cases.

Equation (B.24) can be used to explicitly calculate higher-order corrections to an effective Hamiltonian; for example, see Eqs. (B.10) and (B.11). Using Eq. (B.24), we can calculate mode parameters even when the numerical diagonalization of $\hat{H}_{\mathrm{nl}}$ becomes intractable—an issue that occurs at even moderate number of modes. Moreover, Eq. (B.24) can be used to engineer pumped multi-photon drive processes and to activate non-RWA interactions (Mundhada2018, ), as discussed in Sec. B4.

B4 The driven Josephson system: parametrically-activated interactions

Supplementary Figure S5: Depiction of four-wave mixing in a Josephson junction used to parametrically activate a bilinear interaction between two modes. Off-resonant pumping of the Josephson circuit at $\omega_{p1}$ and $\omega_{p2}$ with effective straights $\xi_{1}$and $\xi_{2}$, respectively, can activate a specific four-wave-mixing process present in $\hat{H}_{\mathrm{nl}}$. The pump frequencies $\omega_{p1}$ and $\omega_{p2}$ could be degenerate. The activation condition is determined by the storage- and readout-mode frequencies $\omega_{c}$ and $\omega_{r}$, respectively; their corresponding mode operators are $\hat{a}_{c}$ and $\hat{a}_{q}$. (a), (b) The activated interaction is a beam-splitter-like conversion process, containing exactly one mode annihilation and one creation operator. The resonance condition is determined by $\omega_{c}-\omega_{r}$. (c), (d) Two-mode squeezing interaction, contains exactly two creation operators. The resonance condition is set by $\omega_{c}+\omega_{r}$. In the first column, panels (a) and (c), the resonance condition is set by $\omega_{p1}-\omega_{p2}$, whereas in the second column, panels (b) and (d), it is set by $\omega_{p1}+\omega_{p2}$. Note, for each of the four diagram, the conjugate process (all arrow directions flipped) is also activated by the pumps.

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.24) 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 (Leghtas2015, ; Pfaff2017, ). The system has a high-quality, storage-cavity mode ($c$) and another low-quality, readout-cavity mode ($r$). The two modes are far detuned and obey the conditions outlined at the start of Sec. B1. From the system $\hat{H}_{\mathrm{nl}}$, we aim to obtain the effective BS Hamiltonian

$$\hat{H}_{\mathrm{eff}}=\hbar g\hat{a}_{c}^{\dagger}\hat{a}_{r}+\hbar g^{*}\hat{a}_{r}^{\dagger}\hat{a}_{c}\;,$$ (B.25)

where $g$ is the rate excitation exchange.

Interaction in the rotating frame.

In the rotating frame with respect to $\hat{H}_{\mathrm{lin}}$, defined by the transform $\hat{U}\left(t\right)\coloneqq\exp\left[it\left(\omega_{c}\hat{a}_{c}^{\dagger}\hat{a}_{c}+\omega_{r}\hat{a}_{r}^{\dagger}\hat{a}_{r}\right)\right]$, the mode operators $\hat{a}_{m}$ acquire a harmonic time dependence, $\hat{a}_{m}\mapsto\hat{a}_{m}\left(t\right)\coloneqq\hat{U}^{\dagger}\left(t\right)\hat{a}_{m}\hat{U}\left(t\right)=\hat{a}_{m}e^{-i\omega_{m}t}$, where $m\in\{c,r\}$. While a term in $\hat{H}_{\mathrm{nl}}$ of the form $\hat{a}_{c}^{\dagger}\hat{a}_{r}$ exists, this terms is non-stationary, $C_{\left(1,0\right),\left(0,1\right)}^{p=4}e^{-it\left(\omega_{r}-\omega_{c}\right)}\hat{a}_{c}^{\dagger}\hat{a}_{r}+\text{H.c.}\;$and is eliminated in the RWA in deriving the effective Hamiltonian.

Parametric activation and interaction rate.

A third mode of the system, indexed by $p$, is off-resonantly driven at frequencies $\omega_{P}$ and $\omega_{P}^{\prime}$ with amplitudes $\epsilon_{1}$ and $\epsilon_{2}$, respectively; within the RWA, the drive Hamiltonian is

$$\hat{H}_{p}\coloneqq\epsilon_{p}e^{-i\omega_{P}t}\hat{a}_{p}+\epsilon_{p}^{*}e^{+i\omega_{P}t}\hat{a}_{p}^{\dagger}\;.$$ (B.26)

Consider the following term contained in $\hat{H}_{\mathrm{nl}}$, see Eq. (B.23),

$$C_{\left(0,1,2\right),\left(1,0,0\right)}^{p=4}\hat{a}_{r}^{\dagger}\left(t\right)\hat{a}_{c}\left(t\right)\hat{a}_{p}^{2}\;,$$ (B.27)

where the mode label order is $\left(r,c,p\right)$. Moving into a displaced and rotating frame defined by $\hat{D}\left(\xi_{p}\right)\coloneqq\exp\left(\xi_{p}\hat{a}_{p}^{\dagger}-\xi_{p}^{*}\hat{a}_{p}\right)\exp\left(-i\hbar\omega_{P}t\hat{a}_{p}^{\dagger}\hat{a}_{p}\right)$, where $\xi_{P}\coloneqq\epsilon/\left(\omega_{P}-\omega_{p}\right)$, the $p$ mode operator is shifted and rotated, $\hat{a}_{p}\mapsto\hat{a}_{p}\left(t\right)\coloneqq\hat{D}\left(\xi_{p}\right)^{-1}\hat{a}_{p}\hat{D}\left(\xi_{p}\right)=\left(\hat{a}_{p}+\xi_{p}\right)e^{-i\omega_{P}t}$ (Leghtas2015, ; Touzard2019, ). In this frame, expanding Eq. (B.27) yields the suggestive interaction term

$$C_{\left(0,1,2\right),\left(1,0,0\right)}^{p=4}\hat{a}_{r}^{\dagger}\hat{a}_{c}\xi_{p}^{2}e^{-it\left(\omega_{r}-\omega_{c}-2\omega_{P}\right)}\;,$$ (B.28)

which is resonantly activated when $\omega_{r}-\omega_{c}-2\omega_{P}=0$, see Supplementary Figure S5. Under this condition, this term survives the RWA when deriving the effective Hamiltonian. We find the effective BS rate by casting Eq. (B.28) in the canonical BS form, Eq. (B.25),

$$\hbar g=C_{\left(0,1,2\right),\left(1,0,0\right)}^{p=4}\xi_{p}^{2}\;.$$ (B.29)

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.24). The procedure is useful for more complex and even cascaded processes (Mundhada2017, ; Mundhada2018, ) used in dissipation engineering (Kapit2017, ). It is reported that the limit of procedure is typically reached when $\xi_{P}$ approaches unity, and other activated nonlinear processes become non-negligible (Verney2019, ; Lescanne2019, ; Tripathi2019, ; Malekakhlagh2020, ).

Supplementary Section 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 $L_{j}$. In Secs. C2 and C3, we extract the EPR $p_{mj}$ and EPR signs $s_{mj}$ from the result of an eigenanalysis simulation of the linearized Josephson circuit, corresponding to $\hat{H}_{\mathrm{lin}}$. 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 , 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 (Nigg2012, ), 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 $L_{j}$, 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 $j$-th Josephson dipole is modeled as a two-dimensional sheet $S_{j}$, see Supplementary Figure S6. The sheet is assigned a surface-impedance boundary condition, which imposes $\vec{E}_{\parallel}=Z_{s}(\hat{n}\times\vec{H}_{\parallel})$ across the sheet, where $\vec{E}_{\parallel}$ and $\vec{H}_{\parallel}$ are the tangential electric and magnetic fields of the sheet, respectively, $\hat{n}$ is the sheet normal vector, and $Z_{s}$ is the complex-valued surface impedance corresponding to a total sheet inductance of $L_{j}$. The hat symbol over $n$ 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 $L_{j}$ 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 non-negligible 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.

Supplementary Figure S6: Illustration (not-to-scale) of the representation of a Josephson dipole in a finite-element (FE) simulation model. The two large grey rectangles on the edges of the illustration depict metal pads of a transmon qubit. The light-grey sheet $S_{j}$ in the center depicts the sheet used to model the Josephson dipole (and potentially its leads) in the FE model. The sheet is assigned a lumped-element inductive boundary condition, with inductance $L_{j}$ (sheet inductance symbolized by floating inductor-element symbol), corresponding to the Josephson dipole inductance with respect to the operating equilibrium point, see Sec. A9. The structural details of the Josephson dipole (location marked by brown cross) and its lead wires (brown wires connected to cross) can be abstracted away.

Performance tip: mesh operations.

The sheet $S_{j}$ 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 $S_{j}$ 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 $p_{mj}$ as a function of the simulation adaptive pass number and by verifying the conditions detailed in Sec. A7.

C2 Calculating the EPR $p_{m}$ in the case of a single Josephson dipole

If a Josephson circuit incorporates exactly one Josephson dipole, then we use the global electric and magnetic eigenmode energies to directly calculate the EPR $p_{m}$ of the dipole in the mode.

Energy balance.

The time-averaged electromagnetic energy in a resonantly excited mode is equally split into an inductive $\mathcal{E}_{\mathrm{ind}}$ and capacitive $\mathcal{E}_{\mathrm{cap}}$ contribution (Pozar, ). This detailed balance, $\mathcal{E}_{\mathrm{ind}}=\mathcal{E}_{\mathrm{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 $\mathcal{E}_{\mathrm{mag}}$ and a kinetic $\mathcal{E}_{\mathrm{kin}}$ contribution; $\mathcal{E}_{\mathrm{ind}}=\mathcal{E}_{\mathrm{mag}}+\mathcal{E}_{\mathrm{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 $\vec{H}_{m}$ and the kinetic energy is stored in the lumped-element boundary condition on $S_{j}$. If lumped-element capacitive boundary conditions are absent from the model, then the capacitive eigenmode energy is stored entirely in electric eigenfields $\vec{E}_{m}$, $\mathcal{E}_{\mathrm{cap}}=\mathcal{E}_{\mathrm{elec}}$; hence,

$$\mathcal{E}_{\mathrm{elec}}=\mathcal{E}_{\mathrm{cap}}=\mathcal{E}_{\mathrm{ind}}=\mathcal{E}_{\mathrm{mag}}+\mathcal{E}_{\mathrm{kin}}\;.$$ (C.1)

Calculating EPR from the energy balance.

Using Eq. (C.1) and Eq. (5) of the main text, $p_{m}=\mathcal{E}_{\mathrm{kin}}/\mathcal{E}_{\mathrm{ind}}$, the EPR is calculated from the ratio of global energy quantities,

$$\boxed{p_{m}=\frac{\mathcal{E}_{\mathrm{elec}}-\mathcal{E}_{\mathrm{mag}}}{\mathcal{E}_{\mathrm{elec}}}\;,}$$ (C.2)

where the total magnetic- and electric-field energies are computed from the eigenfields phasors,

$$\mathcal{E}_{\mathrm{elec}} = \frac{1}{4}\mathrm{Re}\int_{V}\vec{E}_{\text{max}}^{*}\overleftrightarrow{\epsilon}\vec{E}_{\text{max}}\,\mathrm{d}v\;,$$ (C.3)

$$\mathcal{E}_{\mathrm{mag}} = \frac{1}{4}\mathrm{Re}\int_{V}\vec{H}_{\text{max}}^{*}\overleftrightarrow{\mu}\vec{H}_{\text{max}}\,\mathrm{d}v\;,$$ (C.4)

where $\vec{E}_{\mathrm{max}}(x,y,z)$ [resp., $\vec{H}_{\mathrm{max}}(x,y,z)$] is the eigenmode electric (resp., magnetic) phasor, and $\overleftrightarrow{\epsilon}$ (resp., $\overleftrightarrow{\mu}$) denotes the electric-permittivity (resp., magnetic-permeability) tensor. The spatial integrals are performed over total volume $V$ of the device. The electric eigenfield is related to the phasor by

$$\vec{E}\left(x,y,z,t\right)=\mathrm{Re}\left[\vec{E}_{\mathrm{max}}\left(x,y,z\right)e^{i\omega_{m}t}\right]\;.$$ (C.5)

Supplementary Figure S7: Model of a transmon qubit overlaid with the surface-current eigendensity $\vec{J}_{s}$$\left(\vec{r}\right)$ of the qubit eigenmode, obtained using a FE simulation. Transmon pads depicted by the two large gray rectangles, separated by distance $l_{j}$. Small center rectangle represents the sheet model $S_{j}$ of the Josephson junction.

C3 Calculating the EPR $p_{mj}$ in the case of multiple Josephson dipoles

EPR $p_{mj}$ of the $j$-th Josephson dipole.

In the case of multiple Josephson dipoles, the total kinetic energy $\mathcal{E}_{\mathrm{kin}}$ is itself split among the $J$ dipoles, see Eq. (A.20). We calculate the EPR $p_{mj}$ of junction $j$ in mode $m$ using the EPR definition, Eq. (A.50),

$$p_{mj}=\frac{1}{2}L_{j}I_{mj}{}^{2}/\mathcal{E}_{\mathrm{ind}}\;,$$ (C.6)

where $I_{mj}$ is the peak value of the Josephson dipole current in mode $m$. The current $I_{mj}$ is calculated from the integral of the mode surface-current density $\vec{J}_{s,m}\left(x,y,z\right)$ over the dipole sheet $S_{j}$,

$$\left|I_{mj}\right|=l_{j}^{-1}\int_{S_{j}}\left|\vec{J}_{s,m}\right|\,\mathrm{d}s\;,$$ (C.7)

where $l_{j}$ is the length of the sheet, see Supplementary Figure S7.

EPR sign.

The EPR sign is calculated from the orientation of the current $I_{mj}$. The absolute orientation is relative, as described in the main text, and is determined by defining a directed line $\mathrm{DL}_{j}$ across $S_{j}$. If the phasor $\vec{J}_{s,m}$ is aligned with $\mathrm{DL}_{j}$ then $s_{mj}=+1$; otherwise, $s_{mj}=-1$. Hence, we extract the sign using

$$s_{mj}=\operatorname{sign}\int_{\mathrm{DL}_{j}}\vec{J}_{s,m}\cdot\mathrm{d}\vec{l}\;.$$ (C.8)

The direction of the line merely establishes a convention for a positive EPR sign.

Enforcing energy balance.

The convergence of $p_{mj}$, 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 $\sum_{j=1}^{J}p_{mj}=\mathcal{E}_{\mathrm{kin}}/\mathcal{E}_{\mathrm{ind}}$, which is a globally calculated quantity, and is therefore expected to converge quicker.

The above calculations are automated by pyEPR (Note1, ).

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 $M$ 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 impedance-response 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.

Supplementary Section 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 loss-contribution bound of each lossy and radiative element, is extracted from the same eigensolution used to calculate the EPR $p_{mj}$, in essentially the same way—by calculating the fraction of the $m$-th mode energy stored in the $l$-th lossy element—the lossy EPR $p_{ml}$. While lossy EPR signs $s_{ml}$ 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 (Pozar, ; Gao2008, ; Wang2009, ; Wenner2011, ; Zmuidzinas2012, ; Geerlings2013, ; Wang2015, ; Reagor2016, ; Gambetta2017a, ; McRae2020, ). We classily losses as capacitive, inductive, or radiative, and summarize their calculation here. Each loss mechanism is described by a corresponding lossy EPR $p_{ml}$ and an intrinsic quality factor $Q$. 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 $m$-th mode is (Zmuidzinas2012, ; Geerlings2013, ; Reagor2016, )

$$\frac{1}{Q_{\text{total}}}=\frac{1}{Q_{\text{cap}}}+\frac{1}{Q_{\text{ind}}}+\frac{1}{Q_{\text{rad}}}\;,$$ (D.1)

where $Q_{\mathrm{cap}}$ and $Q_{\mathrm{ind}}$ are the total mode quality factors due to capacitive and inductive losses (i.e., those proportional to the intensity of the electric $\left|\vec{E}\right|^{2}$ and magnetic field $\left|\vec{H}\right|^{2}$, respectively, see Secs. D2 and D3), and $Q_{\mathrm{rad}}$ is the total radiative mode quality factor, see Sec. D4. In this section, the mode index $m$ 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 (Vepsalainen2020, ; Cardani2020, ).

D2 Capacitive loss

The total capacitive mode quality $Q_{\mathrm{cap}}$ is the weighted sum of intrinsic quality factors $Q_{l}^{\text{cap}}$ (i.e., the inverse of the dielectric loss tangent) of all lossy dielectrics $l$ in the Josephson circuit (Zmuidzinas2012, ; Geerlings2013, ),

$$\frac{1}{Q_{\text{cap}}}=\sum_{l}\frac{p_{ml}^{\text{cap}}}{Q_{l}^{\text{cap}}}\;,$$ (D.2)

where $p_{ml}^{\text{cap}}$ is the lossy energy-participation ratio of the $l$-th dielectric in the mode—i.e., $p_{ml}^{\text{cap}}$ is the fraction of capacitive energy stored in the dielectric element $l$. 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 (Martinis2014, ; Dial2016, ; Kamal2016-anneal, ). On the other hand, the surface of the substrate, which may be a surface-dielectric layer is classified as a surface-loss element (Martinis2014, ; Wang2015, ). The lossy EPR for bulk capacitive loss is calculated from the eigenfield solutions,

$$p_{ml}^{\text{cap}}=\frac{1}{\mathcal{E}_{\mathrm{elec}}}\frac{1}{4}\operatorname{Re}\int_{V_{l}}\vec{E}_{\text{max}}^{*}\overleftrightarrow{\epsilon}\vec{E}_{\text{max}}\,\mathrm{d}v\;,$$ (D.3)

where the integral is carried over the volume $V_{l}$ of the $l$-th bulk dielectric object; the total electric energy $\mathcal{E}_{\mathrm{elec}}$ is defined in Eq. (C.3). The lossy EPR for a surface dielectric is approximated by

$$p_{ml}^{\text{cap,surf}}=\frac{1}{\mathcal{E}_{\mathrm{elec}}}\frac{t_{l}\epsilon_{l}}{4}\operatorname{Re}\int_{\text{surf}_{l}}\left|\vec{E}_{\text{max}}\right|^{2}\,\mathrm{d}s\;,$$ (D.4)

where the surface-dielectric layer thickness is $t_{l}$ and its permittivity is $\epsilon_{l}$.

D3 Inductive loss

Physically, inductive losses originate from the dissipative flow of electrical current in metals or through metal-metal 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 $Q_{l}^{\mathrm{ind}}$. The bound on the total inductive-loss quality factor $Q_{\text{ind}}$ of mode $m$ is a weighted sum of $Q_{l}^{\mathrm{ind}}$,

$$\frac{1}{Q_{\text{ind}}}=\sum_{l}{\frac{p_{ml}^{\mathrm{ind}}}{Q_{l}^{\mathrm{ind}}}}\;,$$ (D.5)

where $p_{ml}^{\mathrm{ind}}$ is the lossy EPR for inductive element $l$. 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 $\lambda_{0}$ of metal surface $l$, denoted $\mathrm{surf}_{l}$, if the lossy EPR of the surface, and is obtained from the eigenfield solutions,

$$p_{ml}^{\text{ind,surf}}=\frac{1}{\mathcal{E}_{\mathrm{mag}}}\frac{\lambda_{0}\mu_{l}}{4}\operatorname{Re}\int_{\text{surf}_{l}}\left|\vec{H}_{\text{max},\parallel}\right|^{2}\,\mathrm{d}s\;,$$ (D.6)

where $\mu_{l}$ is the magnetic permeability of the surface; typically, $\mu_{l}=\mu_{0}$, and $\vec{H}_{\text{max},\parallel}$ is the magnetic field phasor parallel to the surface. For superconductors, $p_{ml}^{\text{ind,surf}}$ is the kinetic inductance fraction (Gao2008, ; Zmuidzinas2012, ), commonly denoted $\alpha$. The total magnetic energy $\mathcal{E}_{\mathrm{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 (Pozar, ), $Q_{l}^{\text{ind,surf}}\approx 1$. However, in the superconducting state, the lower bound on the quality factor is typically found to be in the range of several thousand (Reagor2013, ), $Q_{l}^{\text{ind,surf}}\gtrsim 10^{3}$. The lower bound for thin-film superconducting aluminum has been measured to exceed $Q_{l}^{\text{ind,surf}}>10^{5}$ (Minev2013, ).

Bulk magnetic loss.

The mode magnetic field can couple to bulk magnetic impurities present in the volume $V_{l}$ of lossy object $l$. The lossy EPR for the bulk-inductive-loss mechanism of volume $V_{l}$ is

$$p_{ml}^{\mathrm{ind,bulk}}=\frac{1}{\mathcal{E}_{\mathrm{mag}}}\frac{1}{4}\operatorname{Re}\int_{V_{l}}\vec{H}_{\text{max}}^{*}\overleftrightarrow{\mu}\vec{H}_{\text{max}}\,\mathrm{d}v\;.$$ (D.7)

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 (Brecht2015, ). 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 $\text{seam}_{l}$ tracing out the location of the seam at the surface of the two mating walls. The inductive lossy EPR participation for the seam is

$$p_{ml}^{\text{ind,seam}}=\frac{1}{\mathcal{E}_{\mathrm{mag}}}\frac{\lambda_{0}t_{l}\mu_{l}}{4}\operatorname{Re}\int_{\text{seam}_{l}}\left|\vec{H}_{\text{max},\perp}\right|^{2}\,\mathrm{d}l\;,$$ (D.8)

where the seam thickness is denoted $t_{l}$, its magnetic permeability $\mu_{l}$, and its the penetration depth $\lambda_{0}$. It is convenient to rewrite the total seam loss due to seam $l$ in terms of an effective seam admittance $g_{\text{seam}}$, defined in Ref. Brecht2015, ,

$$\frac{p_{ml}^{\mathrm{ind,seam}}}{Q_{\text{seam}}}=\frac{1}{g_{\text{seam}}}\frac{\int_{\text{seam}_{l}}\left|\vec{J}_{s}\times\vec{l}\right|^{2}\,\mathrm{d}l}{\omega\mu_{0}\int_{\text{all}}\left|H_{\text{max}}\right|^{2}\,\mathrm{d}V}\;.$$ (D.9)

D4 Radiative loss and input-output coupling

Supplementary Figure S8: Schematic of a transmon-qubit circuit (mode operator $\hat{a}$), comprising a Josephson tunnel junction (flux operator $\hat{\Phi}_{j}$) and a capacitor, coupled to two input-output ports. The input and output fields of the left (resp., right) transmission line are $\hat{a}_{\mathrm{in},1}$ and $\hat{a}_{\mathrm{out},1}$ (resp., $\hat{a}_{\mathrm{in},2}$ and $\hat{a}_{\mathrm{out},2}$). (a) Ports 1 and 2 both coupled to the top qubit node; hence, both port EPR signs $s_{mp}$ are equal, $s_{q1}=s_{q2}$. (b) Port 1 couples to the top node, but port 2 couples to the bottom qubit node. The port EPR signs $s_{mp}$ have opposite signs, $s_{q1}=-s_{q2}$.

The Josephson circuit incorporates radiative boundaries, typically purposefully introduced to serve as input-output ports of the circuit (Pozar, ), 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 co-planar waveguide. The total radiative quality factor $Q_{\text{rad}}$ of mode $m$ is the sum of the individual port contributions $Q_{\text{rad}}^{-1}=\sum_{p=1}^{P}Q_{mp}^{-1}$, where $P$ is the total number of ports and $Q_{mp}$ is the quality factor due to port $p$. Below, we describe how port $p$ is modeled in the FE simulation to extract $Q_{mp}$ from the eigensolutions.

Radiative energy loss.

Energy stored in mode $m$ can leak at a rate $\kappa_{mp}$ through port $p$ 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 (Houck2008-Purcell, ); while, for the readout mode of the same structure, the energy loss to the readout port sets the rate of information gain (Clerk2010, ). The rate $\kappa_{mp}$ is calculated from port EPR $p_{mp}$, as detailed in the following for both wanted and spurious terms. The port EPR sign $s_{mp}$ 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 $S_{p}$ 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 $R_{p}$. In the case of a 50 $\Omega$ port line, $R_{p}=50\,\Omega$. 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 $\kappa_{mp}$ from the port EPR $p_{mp}$.

From the lossless eigensolutions of mode $m$, the energy lost to the effective resistor $R_{p}$ of port $p$ during one mode oscillation period $T_{m}=2\pi/\omega_{m}$ is

$$P_{mp}=\frac{1}{2}R_{p}I_{mp}^{2}T_{m}\;,$$ (D.10)

where $I_{mp}$ is the peak current across the port due to the excitation of mode $m$; see Eq. (C.7). The total mode energy $\mathcal{E}_{m}\left(t\right)$ at time $t$ decays at rate

$$\frac{d}{dt}\mathcal{E}_{m}=-\sum_{p}{\kappa_{mp}}\mathcal{E}_{m}\;.$$ (D.11)

Hence, assuming a high quality mode, the energy loss during one oscillation period, between times $t=0$ and $T_{m}$, is

$$\mathcal{E}_{m}(T_{m})=\mathcal{E}_{m}(0)-\sum_{p}{T_{m}\kappa_{mp}\mathcal{E}_{m}(0)}\;.$$ (D.12)

This shows that the loss to port $p$ during one period is $P_{mp}=T_{m}\kappa_{mp}\mathcal{E}_{m}(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,

$$\kappa_{mp}=\frac{\frac{1}{2}RI_{mp}^{2}}{\mathcal{E}_{m}(0)}\;.$$ (D.13)

Thus, the mode coupling quality factor is

$$Q_{mp}\coloneqq\omega_{m}/\kappa_{mp}=\frac{\omega_{m}\mathcal{E}_{m}(0)}{\frac{1}{2}RI_{mp}^{2}}\;.$$ (D.14)

Sign of the I-O participation.

If there are multiple ports, the port EPR sign $s_{mp}$, calculated using Eq. (C.8), is important in a manner similar to that explained for the case of the EPR sign $s_{mj}$ 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 $S_{p}$ 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 $R_{p}=50\,\Omega$ (Malekakhlagh2016-A2, ; Malekakhlagh2017-Cutoff-Free, ). 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 $Q_{m}$ approaches unity. For lossless treatment of $S_{p}$, suitable when $Q_{m}\gg 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 $R_{p}$.

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 (Burkhart2020, ). The reactive part of the port structure thus houses energy in its internal modes. There can be treated by modifying the boundary condition on $S_{p}$ 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 (Wang2019-cav-atten, ).