QuCAT: Quantum Circuit Analyzer Tool in Python

Quantum circuits constructed from Josephson junctions and superconducting electronics are key to many quantum computing and quantum optics applications. Designing these circuits involves calculating the Hamiltonian describing their quantum behavior. Here we present QuCAT, or"Quantum Circuit Analyzer Tool", an open-source framework to help in this task. This open-source Python library features an intuitive graphical or programmatical interface to create circuits, the ability to compute their Hamiltonian, and a set of complimentary functionalities such as calculating dissipation rates or visualizing current flow in the circuit.


I. INTRODUCTION
Quantum circuits, constructed from superconducting electronics and involving one or more Josephson junctions, have steadily gained prominence in experimental and theoretical physics over the past twenty years. Foremost, they are one of the most successful platforms in the quest to build a quantum computer (Devoret and Schoelkopf, 2013). The control that can be gained over their quantum state, and the flexibility in their design have also made these circuits an excellent test-bed to probe fundamental quantum effects (Gu et al., 2017). They can also be coupled to other systems, such as atoms, spins, acoustic vibrations or mechanical oscillators, acting as a tool to measure and manipulate these systems at a quantum level (Xiang et al., 2013).
Any application mentioned above generally translates to a desired Hamiltonian, which governs the physics of the circuit. The task of the quantum circuit designer is to determine which circuit components to use, how to inter-connect them, and calculate the corresponding Hamiltonian (Nigg et al., 2012;Vool and Devoret, 2017). Performing this task analytically can be time consuming or even challenging.
Here we present QuCAT, which stands for "Quantum Circuit Analyzer Tool", an open-source Python framework to help in analyzing and understanding quantum circuits. We provide an easy interface to create and visualize circuits, either programmatically or through a graphical user interface. A Hamiltonian can then be generated for further analysis in QuTiP (Johansson et al., 2012(Johansson et al., , 2013. The current version of QuCAT supports quantization in the basis of normal modes of the linear circuit (Nigg et al., 2012), making it suited for the analysis of weakly anharmonic circuits with small losses. The properties of these modes: their frequency, dissipation rates, anharmonicity and cross-Kerr couplings can be directly calculated. The user can also visualize the current flows in the circuit associated with each normal mode. The library covers lumped element circuits featuring an arbitrary number of Josephson junctions, inductors, capacitors and resistors. Through equivalent lumped element circuits, certain distributed elements such as waveguide resonators can also be analyzed (see Sec. A.3). The software relies on the symbolic manipulation of the circuits equations, making it reliable even for vastly different circuits and parameters. It also results in efficient parameter sweeps, as analytical manipulations need not be repeated for different circuit parameters. In a few seconds, circuits featuring 10 nodes (or degrees of freedom), corresponding to between 10 and 30 circuit elements can be simulated.
In the main section of this article, we cover the functionalities of the software. We start by showing how to create circuits, first using the graphical user interface, then programmatically. We then demonstrate how to generate the corresponding Hamiltonian. Lastly, we show how to extract the characteristics of the circuit modes: frequencies, dissipation, anharmonicity and cross-Kerr coupling and present a tool to visualize these modes. This main section will feature as an example the standard circuit of a transmon qubit coupled to a resonator (Koch et al., 2007). In the appendices, we will first use Qu-CAT to analyze some recent experiments: a tuneable coupler (Kounalakis et al., 2018), a multi-mode ultra-strong coupling circuit , a microwave optomechanics circuit (Ockeloen-Korppi et al., 2016) and a Josephson-ring based qubit (Roy et al., 2017). We then provide an overview of the circuit quantization method used and the algorithmic methods which implement it. The limitations of these methods regarding weak anharmonicity and circuit size will then be presented. Finally we will explain how to install QuCAT and we provide a summary of all its functions. More tutorials and examples are available on the QuCAT website https://qucat.org/.

II. CIRCUIT CONSTRUCTION
Any use of QuCAT will start with importing the qucat library import qucat import qucat circuit = qucat.GUI ('netlist.txt') circuit.show() 100 fF Lj 1 fF 100 fF 50 Ω 10 nH 500 aF FIG. 1 Construction of a circuit: code and output. The circuit used as an example in this section comprises of a transmon qubit on the left, coupled through a 1 fF capacitor to an LC-oscillator. Dissipation arises from the capacitive coupling of the LC-oscillator to a 50 Ω resistor on the right. After importing the qucat package, the circuit object is created manually through a graphical user interface (GUI) opened after calling qucat . GUI ( " netlist . txt " ). All information necessary to construct the circuit is stored in the text file netlist . txt. After closing the GUI, this information is also stored in the variable circuit. The show method finally displays the circuit.
One should then create a circuit. These are named Qcircuit, short for "quantum circuit" in QuCAT. There are two ways of creating a Qcircuit: using the graphical user interface (GUI), or programmatically.

A. Creating a circuit with the GUI
We first cover how to create a circuit with the GUI. This is done through this command circuit = qucat . GUI ( ' netlist . txt ') which opens the GUI. The GUI will appear as a separate window, which will block the execution of the rest of the Python script until the window is closed. The user can drag-in and drop capacitors, inductors, resistors or Josephson junctions, or grounds. These components can then be inter-connected with wires. Each change made to the circuit will be automatically be saved in the ' netlist . txt ' file. After closing the GUI, the Qcircuit object will be stored in the variable named circuit which we will use for further analysis.

B. Creating a circuit programmatically
Alternatively, one can create a circuit with only Python code. This is done by creating a list of circuit components with the functions J, L, C and R for junctions, inductors, capacitors and resistors respectively. For the circuit of Fig. 1: circuit_components = [ qucat . C (0 ,1 ,100 e -15) , # transmon qucat . J (0 ,1 , ' Lj ') , qucat . C (0 ,2 ,100 e -15) , # resonator qucat . L (0 ,2 ,10 e -9) , qucat . C (1 ,2 ,1 e -15) , # coupling capacitor qucat . C (2 ,3 ,0.5 e -15) , # ext . coupl . cap . qucat . R (3 ,0 ,50) # 50 Ohm load ] All circuit components take as first two argument integers referring to the negative and positive node of the circuit components. Here 0 corresponds to the ground node for example. The third argument is either a float giving the component a value, or a string which labels the component parameter to be specified later. Doing the latter avoids performing the computationally expensive initialization of the Qcircuit object multiple times when sweeping a parameter. By default, junctions are parametrized by their Josephson inductance L j = φ 2 0 /E j where φ 0 = /2e is the reduced flux quantum, and E j (in Joules) is the Josephson energy.
Once the list of components is built, we can create a Qcircuit object via the Network function circuit = Network ( circuit_components ) as with a construction via the GUI, the Qcircuit object will be stored in the variable named circuit which we will use for further analysis.

III. GENERATING A HAMILTONIAN
The Hamiltonian of a Josephson circuit is given bŷ It is written in the basis of its normal modes. These have an angular frequency ω m and we write the operator which creates (annihilates) photons in the modê a † m (â m ). The cosine potential of each Josephson junction j with Josephson energy E j has been Taylor expanded to order n for small values of its phase fluctuationsφ j across it. The phase fluctuations are a function of the annihilation and creation operators of the modeŝ ϕ j = m ϕ zpf,m,j (â † m +â m ). For a detailed derivation of this Hamiltonian, and the method used to obtain its parameters, see Sec. B.
There are three different parameters that the user should fix 1. the set of modes to include 2. for each of these modes, the number of excitations to consider 3. the order of the Taylor expansion.
The more modes and excitations are included, and the higher Taylor expansion order, the more faithful the Hamiltonian will be to physical reality. The resulting increase in Hilbert space size will however make it more  Fig. 1. The Hamiltonian is expressed in the basis of circuit normal modes m with frequencies fm = ωm/2π, annihilation operatorsâm, and zero-point phase fluctuations ϕ zpf,m,j across junction j with Josephson energy Ej. The junction non-linearities are expressed through a Taylor expansion of the cosine potentials, where the user chooses the degree of Taylor expansion. The other arguments are the list of modes to include, the number of excitations to consider for each of these modes, and any unspecified component value, here Lj. The returned Hamiltonian is a QuTiP object, giving the user access to an extensive set of tools for further analysis (Johansson et al., 2012(Johansson et al., , 2013. As an example, we compute the eigenenergies of the Hamiltonian, and plot the two first transition frequencies, as a function of Lj.
computationally expensive to perform further calculations. Typically, larger degrees of anharmonicity require a larger Hilbert space, with a fundamental limitation on the maximum anharmonicity due to the choice of basis. We expand on these topics in Sec. D.2.
Such a Hamiltonian is generated through the method hamiltonian. More specifically, this function returns a QuTiP object (Johansson et al., 2012(Johansson et al., , 2013, enabling an easy treatment of the Hamiltonian. All QuCAT functions use units of Hertz, so the function is actually returninĝ H/h. As an example, we generate a Hamiltonian for the circuit of Fig. 1  With modes = [0 ,1], we are specifying that we wish to consider the first and second modes of the circuit. Modes are numbered with increasing frequency, so here we are selecting the two lowest frequency modes of the circuit. With excitations = [10 ,12], we specify that for mode 0 (1) we wish to consider 10 (12) excitations. With taylor = 4, we are specifying that we wish to expand the cosine potential to fourth order, this is the lowest order which will give an anharmonic behavior. The unspecified Josephson inductance must now be fixed through a keyword argument Lj = 8e -9. Doing so avoids initializing the Qcircuit objects multiple times during parameter sweeps, as initialization is the most computationally expensive task. We calculate these energies with different values of the Josephson inductance, and the first two transition frequencies are plotted in Fig. 2, showing the typical avoided crossing seen in a coupled qubit-resonator system. f,k,A,chi = circuit.f_k_A_chi( Lj = numpy.linspace(11e-9,9e-9)) FIG. 3 Extracting eigenfrequencies, loss-rates, anharmonicities, and cross-Kerr couplings. We apply the f_k_A_chi method to circuit defined in Fig. 1 to obtain a list of eigenfrequencies (f), loss-rates (k), anharmonicities (A), and cross-Kerr couplings (chi), for all the normal modes of the circuit. There is one unspecified variable in the circuit, the Josephson inductance Lj, which is here specified with a list of values. In QuCAT can also return the parameters of the (already diagonal) Hamiltonian in first-order perturbation theorŷ valid for weak anharmonicity χ mn , A m ω m . The physics of this Hamiltonian can be understood by considering that an excitation of one of the circuit modes may lead to current traversing a Josephson junction. This will change the effective inductance of the junction, hence changing its own mode frequency, as well as the mode frequencies of all other modes. This is quantified through the anharmonicity or self-Kerr A m and cross-Kerr χ mn respectively. When no mode is excited, vacuum-fluctuations in current through the junction give rise to shifted mode energies ω m − A m − n χ mn /2.
In a circuit featuring resistors, these anharmonic modes will be dissipative. A mode m will lose energy at a rate κ m . If these rates are specified in angular frequencies, the relaxation time T 1,m of mode m is given by T 1,m = 1/κ m . A standard method to include the loss rates in a mathematical description of the circuit is through the Lindblad equation (Johansson et al., 2012), where the losses would be included as collapse operators √ κ mâm The frequencies, dissipation rates, and Kerr parameters can all be obtained via methods of the Qcircuit object. These methods will return numerical values, and we should always specify the values of symbolically defined circuit parameters as keyword arguments. Lists, or Numpy arrays, can be provided here making it easy to perform parameter sweeps. Additionally, initializing the circuit is the most computationally expensive operation, so this will be by far the fastest method to perform parameter sweeps.
We will assume that we want to determine the parameters of the Hamiltonian (2) for the circuit of Fig. 1 at different values of L j . The values for L j are stored as a Numpy array 101) We can assign the frequency, dissipation rates, self-Kerr, and cross-Kerr parameters to the variables f, k, A and chi respectively, by calling f = circuit . eigenfrequencies ( Lj = Lj_list ) k = circuit . loss_rates ( Lj = Lj_list ) A = circuit . anharmonicities ( Lj = Lj_list ) chi = circuit . kerr ( Lj = Lj_list ) or alternatively through a single function call: f ,k ,A , chi = circuit . f_k_A_chi ( Lj = Lj_list ) All values returned by these methods are given in Hertz, not in angular frequency. With respect to the conventional way of writing the Hamiltonian, which we have also adopted in (2), we thus return the frequencies as ω m /2π, the loss rates as κ m /2π and the Kerr parameters as A m /h and χ mn /h. Note that f, k, A, are arrays, where the index m corresponds to mode m, and modes are ordered with increasing frequencies.
For example, f [0] will be an array of length 101, which stores the frequencies of the lowest frequency mode as Lj is swept from 11 to 9 nH. The variable chi has an extra dimension, such that chi [m , n ] corresponds to the cross-Kerr between modes m and n, and chi [m , m ] is the self-Kerr of mode m, which has the same value as A [ m ]. These generated values are plotted in Fig. 3.
We can also print these parameters in a visually pleasing way to get an overview of the circuit characteristics for a given set of circuit parameters. For a Josephson inductance of 9 nH, this is done through the command circuit . f_k_A_chi ( Lj = 10 e -9 , pretty_print = True ) which will print We see that mode 1 is significantly more anharmonic than mode 0, whereas mode 0 has however a higher dissipation. We would expect that mode 1 is thus the resonance which has current fluctuations mostly located in the junction, whilst mode 0 is located on the other side to the coupling capacitor, where it can couple more strongly to the resistor.
Such interpretations can be verified by plotting a visual representation of the normal modes on top of the circuit as explained below. This can be done by plotting either the current, voltage, charge or flux distribution, overlaid on top of the circuit schematic. As shown in Fig. 4, this is done by adding arrows, representing one of these quantities at each circuit component and annotating it with the value of that component. The annotation corresponds to the complex amplitude, or phasor, of a quantity across the component, if the mode was populated with a single photon amplitude coherent state. The absolute value of this annotation corresponds to the contribution of a mode to the zero-point fluctuations of the given quantity across the component. The direction of the arrows indicates what direction we take for 0 phase for that component.
We note that an independantly developped Julia platform also allows the calculation of normal mode frequencies and dissipation rates for circuits (?).  show_normal_mode takes as argument any unspecified circuit parameter, here we specify Lj =10 e -9 where the two modes undergo an avoided crossing. We plot each mode by specifying mode = 0 or mode = 1 and see that for mode 0, the anti-symmetric mode, the voltage has opposite signs on each side of the coupling capacitor, leading to a larger voltage across the coupler (and hence a larger effective capacitance and lower frequency) than the symmetric mode.

V. OUTLOOK
We have presented QuCAT, a Python library to automatize and speed up the design process and analysis of superconducting circuits. By facilitating quick tests of different circuit designs, and helping develop an intuition for the physics of quantum circuits, we also hope that Qu-CAT will enable users to develop even more innovative circuits.
Possible extensions of the QuCAT features could include black-box impedance components to model distributed components (Nigg et al., 2012), more precisely modeling lossy circuits (Solgun et al., 2014;Solgun and DiVincenzo, 2015), handling static offsets in flux or charge through DC sources, additional elements such as coupled inductors or superconducting quantum interfer-ence devices (SQUIDS) and different quantization methods, enabling for example quantization in the charge or flux basis. The latter would extend QuCAT beyond the scope of weakly-anharmonic circuits.
In terms of performance, QuCAT would benefit from delegating analytical calculations to a more efficient, compiled language, with the exciting prospect of simulating large scale circuits (Kelly et al., 2019). Note however that there is a strong limitation on the maximum Hilbert space size that one can simulate after extracting the Hamiltonian. A 50 Ω load, representing a cable, is connected to an LC resonator through two LC band pass filters. (b) The dissipation rate of the resonator is plotted as a function of inductance and capacitance of the filter using the loss_rates method.
In this application we show how QuCAT can be used to design classical microwave components. We study here a band pass filter made from two LC oscillators with the inductor inline and a capacitive shunt to ground. Such a filter can be used to stop a DC bias line from inducing losses, whilst being galvanically connected to a resonator, see for example Ref. (Viennot et al., 2018). In this case we are interested in the loss rate κ of a LC resonator connected through this filter to a 50 Ω load, which could emulate a typical microwave transmission line. We want to study how κ varies as a function of the inductance L and capacitance C of its components.
The QuCAT GUI function can be used to open the GUI, the user will manually create the circuit, and upon closing the GUI a Qcircuit object is stored in the variable filtered_cavity. By calling the method show, we display the circuit as shown in Fig. 5(a). These steps are accomplished with the code # Open the GUI and manually build the circuit filtered_cavity = qucat . GUI ( ' netlist . txt ') # Display the circuit filtered_cavity . show () We can then access the loss rates of the different circuit modes through the method loss_rates. Since the values of C and L were not specified in the construction of the circuit, their values have to be passed as keyword arguments upon calling loss_rates. For example, the loss rate for a 1 pF capacitor and 100 nH inductor is obtained through # Loss rates of all modes k_all = filtered_cavity . loss_rates ( C = 1e -12 , L = 100 e -9) # Resonator loss rate k = k_all [ -1] Since the filter capacitance and inductance is large relative to the capacitance and inductance of the resonator, the modes associated with the filter will have a much lower frequency. We can thus access the loss rate of the resonator by always selecting the last element of the array of loss rates with the command k_all [ -1] The dissipation rates for different values of the capacitance and inductance are plotted in Fig. 5(b).

Computing optomechanical coupling
In this application, we show how QuCAT can be used for analyzing another classical system, that of microwave optomechanics. One common implementation of microwave optomechanics involves a mechanically compliant capacitor, or drum, embedded in one or many microwave resonators (Teufel et al., 2011). One quantity of interest is the single-photon optomechanical coupling. This quantity is the change in mode frequency ω m that occurs for a displacement x zpf of the drum (the zero-point fluctuations in displacement) The change in mode frequency as the drum head moves ∂ω m /∂x is not straightforward to compute for complicated circuits. One such example is that of Ref. (Ockeloen-Korppi et al., 2016), where two microwave resonators are coupled to a drum via a network of capacitances as shown in Fig. 6(a). Here, we will use QuCAT QuCAT reconstruction of the circuit. By specifying a label for the mechanically compliant capacitances, we have the possibility to compute the eigenfrequencies ωm with the method eigenfrequencies for small variation in C d (x). This enables an easy computation of the optomechanical coupling ∝ dωm/dx. to calculate the optomechanical coupling of the drums to both resonator modes of the circuit.
We start by reproducing the circuit of Fig. 6(a), excluding the capacitive connections on the far left and right. This is done via the graphical user interface opened with the qucat . GUI function. Upon closing the graphical user interface, the resulting Qcircuit is stored in the variable OM, and the show method is used to display the schematic of Fig. 6(a). These steps are accomplished with the code below # Open the GUI and manually build the circuit OM = qucat . GUI ( ' netlist . txt ') # Display the circuit OM . show () We use realistic values for the circuit components without trying to be faithful to Ref. (Ockeloen-Korppi et al., 2016), the aim of this section is to illustrate a method to obtain g 0 . Crucially, the mechanically compliant capacitors have been parametrized by the symbolic variable C d . We can now calculate the resonance frequencies of the circuit with the method eigenfrequencies as a function of a keyword argument C d .
The next step is to define an expression for C d as a function of the mechanical displacement x of the drum head with respect to the immobile capacitive plate below it.
def Cd ( x ) : # Radius of the drumhead radius = 10 e -6 # Formula for half a circular parallel plate capacitor return eps * pi * radius **2/ x /2 where pi and eps have been set to the values of π and the vacuum permittivity respectively. We have divided the usual formula for parallel plate capacitance by 2 since, as shown in Fig. 6(a), the capacitive plate below the drum head is split in two electrodes. We are now ready to compute g 0 . Following Ref. (Teufel et al., 2011), we assume the rest position of the drum to be D = 50 nm above the capacitive plate below. And we assume the zero-point fluctuations in displacement to be x zpf = 4 fm. We start by differentiating the mode frequencies with respect to drum displacement using a finite differences formula and we find that the current is majoritarily located in the inductor of resonator 1.

Convergence in multi-mode cQED
In this section we use QuCAT to study the convergence of parameters in the first order Hamiltonian (Eq. 2) of an ultra-strongly coupled multi-mode circuit QED system.
Using a length of coplanar waveguide terminated with engineered boundary conditions is a common way of building a high quality factor microwave resonator. One CPW resonator Transmon implementation is a λ/4 resonator terminated on one end by a large shunt capacitor, acting as a near-perfect short circuit for microwaves such that only a small amount of radiation may enter or leave the resonator. On the other end one places a small capacitance to ground: an open circuit. The shunt capacitor creates a voltage node, and at the open end the voltage is free to oscillate. This resonator hosts a number of normal modes, justifying its lumped element equivalent circuit: a series of LC oscillators with increasing resonance frequency . Here, we study such a resonator with a transmon circuit capacitively coupled to the open end. In particular we consider this coupling to be strong enough for the circuit to be in the multi-mode ultra-strong coupling regime as studied experimentally in Ref.  and theoretically in Ref. . The particularity of this regime is that the transmon has a considerable coupling to multiple modes of the resonator. It then becomes unclear how many of these modes to consider for a realistic modeling of the system. This regime is reached by maximizing the coupling capacitance of the transmon to the resonator and minimizing the ca-pacitance of the transmon to ground. The experimental device accomplishing this is shown in Fig. 7(a), with its schematic equivalent in Fig. 7(b), and the lumpedelement model in Fig. 7(c). We will use QuCAT to track the evolution of different characteristics of the system as the number of considered modes N increases. For this application, programmatically building the circuit is more appropriate than using the GUI. We start by defining some constants Note that 12 is the index of the ground node.
We can now access some parameters of the system. Only the first mode of the resonator has a lower frequency than the transmon. The transmon-like mode is thus indexed as mode 1. Its frequency is given by mmusc . eigenfrequencies () [1] and the anharmonicity of the transmon, computed from first order perturbation theory (see Eq. 2) with mmusc . anharmonicities () [1] Finally the Lamb shift, or shift in the transmon frequency resulting from the zero-point fluctuations of the resonator modes, is given following Eq. (2) by the sum of half the cross-Kerr couplings between the transmon mode and the others These parameters for different total number of modes are plotted in Figs 7(d-f).
From this analysis, we find that as we reach 10, the plotted parameters are converging. Surprisingly, adding even the highest modes significantly modifies the total Lamb shift of the Transmon despite large frequency detunings.

Modeling a tuneable coupler
In this section, we study the circuit of Ref. (Kounalakis et al., 2018) where two transmon qubits are coupled through a tuneable coupler. This tuneable coupler is built from a capacitor and a Superconducting Quantum Interference Device, or SQUID. By flux biasing the SQUID, we change the effective Josephson energy of the coupler, which modifies the coupling between the two transmons. We will present how the normal mode visualization tool helps in understanding the physics of the device. Secondly, we will show how a Hamiltonian generated with QuCAT accurately reproduces experimental measurements of the device.
We start by building the device shown in Fig. 8(a). More specifically, we are interested in the part of the device in the dashed box, consisting of the two transmons and the tuneable coupler. The other circuitry, the flux line, drive line and readout resonator could be included to determine external losses, or the dispersive coupling of the transmons to their readout resonator. We will omit these features for simplicity here. After opening the GUI with the qucat . GUI function, manually constructing the circuit, then closing the GUI, the resulting Qcircuit is stored in a variable TC.
The inductance L j of the junction which models the SQUID is given symbolically, and will have to be specified when calling Qcircuit functions. Since L j is controlled through flux φ in the experiment, we define a function which translates φ (in units of the flux quantum) to L j def Lj ( phi ) : # maximum Josephson energy Ejmax = 6.5 e9 # junction asymmetry d = 0.0769 # flux to Josephson energy Ej = Ejmax * numpy . cos ( pi * phi ) * numpy . sqrt (1+ d **2 * numpy . tan ( pi * phi ) **2) # Josephson energy to inductance return ( hbar /2/ e ) **2/( Ej * h ) where pi, h, hbar, e were assigned the value of π, Plancks constant, Plancks reduced constant and the electron charge respectively.
By visualizing the normal modes of the circuit, we can understand the mechanism behind the tuneable coupler. We plot the highest frequency mode at φ = 0, as shown in Fig. 8  that the value of L j influences the oscillation frequency of the mode. Conversely, if we plot the anti-symmetric mode instead, where currents are flowing away from the coupler in each transmon, we find a current through the coupler junction and capacitor on the order of 10 −21 A. This mode frequency should not vary as a function of L j . When the bare frequency of the coupler matches the coupled transmon frequencies, the coupler acts as a band-stop filter, and lets no current traverse. At this point, both symmetric and anti-symmetric modes should have identical frequencies.
In Fig. 8(c) this effect is shown experimentally through a measure of the first transitions of the two non-linear modes. One is tuned with flux (symmetric mode), the other barely changes (anti-symmetric mode). We can reproduce this experiment by generating a Hamiltonian with QuCAT and diagonalizing it with QuTiP for different values of the flux. For example, at 0 flux, the two first two transition frequencies f1 and f2 can be generated from f1 and f2 is plotted in Fig. 8(d) for different vales of flux and closely matches the experimental data. Note that we have constructed a Hamiltonian with modes 1 and 2, excluding mode 0, which corresponds to oscillations of current majoritarily located in the tuneable coupler. One can verify this fact by plotting the distribution of currents for mode 0 using the show_normal_mode method.
This experiment can be viewed as two "bare" transmon qubits coupled by the interaction where left and right transmons are labeled L and R and  FIG. 9 Trimon device and Purcell-decay-protected mode visualization. (a) Schematic of the cross-cut of a 3D microwave cavity. Dark gray shows metal whilst light gray show the hollowed out section forming the cavity. Arrows represent the electric field of the T E101, or "vertical" cavity mode. In the cavity is placed a chip hosting the trimon circuit shown in the optical micrograph (b). The circuit has 4 capacitive pads labeled from 1 to 4. These pads are connected by the Josephson junction ring shown in the scanning electron microscope image (c). Scale bars correspond 200 and 2 µm for panels (b) and (c) respectively. (d) Lumped-element equivalent circuit of the device constructed using the QuCAT GUI and displayed with show_normal_mode. The four pads of the trimon are color-coded to match (b). The capacitor Ca formed by pads 1 and 3 forms an electrical dipole which couples to a vertical cavity mode, and the capacitor C b formed by pads 2 and 4 forms an electrical dipole which couples to modes with horizontal electric fields. The show_normal_mode overlays the voltage across different components if a single-photon amplitude coherent state was populating mode 2. This mode has a particularity that the voltage is concentrated across the junctions and their parallel capacitors without leading to a buildup of voltage across the capacitors Ca or C b . This decouples mode 2 from the cavity mode decay (no Purcell effect) whilst the presence of voltage fluctuations across the junctions will ensure cross-Kerr coupling to the other modes of the system. Concerning panels (a-c): reprinted figures with permission from T. Roy et al., Phys. Rev. Appl. 7 (5), 054025 (2017). Copyright 2017 by the American Physical Society. σ x is the x Pauli operator. The coupling strength g reflects the rate at which the two transmons can exchange quanta of energy. If the transmons are resonant a spectroscopy experiment reveals a hybridization of the two qubits, which manifests as two spectroscopic absorption peaks separated in frequency by 2g. From this point of view, this experiment thus implements a coupling which is tuneable from an appreciable value to near 0 coupling.

Studying a Josephson-ring-based qubit
In this section, we demonstrate the ability for QuCAT to analyze more complex circuits. The experiment of Ref. (Roy et al., 2017) features a Josephson ring geom-etry, which is a Wheatstone-bridge-like circuit, typically difficult to analyze as it cannot be decomposed in series and parallel connections. We consider the coupling of this ring to two lossy modes of a cavity, bringing the total number of modes in the circuit to 5. We aim to understand the key feature of this circuit: that one qubit-like mode acts as a quadrupole with little coupling to the resonator modes.
The studied device consists of a 3D cavity ( Fig. 9(a)) hosting a number of microwave modes, in which is positioned a chip patterned with the trimon circuit. The trimon circuit has four capacitive pads in a cross shape ( Fig. 9(b)) which have an appreciable coupling between each other making up the capacitance of the trimon qubit   10 Other modes of the Trimon.
Using the f_k_A_chi method together with the pretty_print option gives the user an overview of the different modes frequencies, dissipations rates and levels of anharmonicity. Here we have overlaid the output of the method with schematics of the corresponding trimon and cavity modes adapted from Ref. (Roy et al., 2017). One can identify a mode to the schematic by observing where the currents or voltages are mostly located in the circuit using the show_normal_mode method as in Fig. 9(d). Since the only resistors of the circuit are located in the cavity modes, all dissipation in transmon modes 0 through 2 are due to the Purcell effect. Mode 2 is better protected from this effect by 3 orders of magnitude with respect to the two other transmon modes. modes. The two vertically (horizontally) positioned pads will couple to modes of the 3D cavity featuring vertical (horizontal) electric fields. We will consider both a vertical and a horizontal cavity mode in our model. We number these pads from 1 to 4 as displayed in Fig. 9(b). Each pad is connected to its two nearest neighbors by a Josephson junction (Fig. 9(c)), forming a Josephson ring.
Using the QuCAT GUI, we build a lumped element model of this device, generating a Qcircuit object we store in the variable trimon.

trimon = qucat . GUI ( ' netlist . txt ')
The cavity modes are modeled as RLC oscillators with each plate of their capacitors capacitively coupled to a pad of the trimon circuit. The junction inductances are assigned different values, first to reflect experimental reality, but also to avoid infinities arising in the QuCAT analysis. Indeed, the voltage transfer function of this Josephson ring between nodes 1,3 and nodes 2,4 will be exactly 0, which will cause errors when initializing the Qcircuit object. Component parameters are chosen to only approximatively match the experimental results of Ref. (Roy et al., 2017), the objective here is to demonstrate QuCAT features rather than accurately model the experiment.
The particularity of this circuit is that it hosts a quadrupole mode. It corresponds here to the second highest frequency mode and can be visualized by calling trimon . show_normal_mode ( mode = 2 , quantity = ' voltage ') the result of which is displayed in Fig. 9(d). The voltage oscillations are majoritarily located in the junctions, indicating this is not a cavity mode, but a mode of the trimon circuit. Crucially, the polarity of voltages across the junctions is such that the total voltage between pads 1 and 3 and the total voltage across pads 2 and 4 is 0, warranting the name of "quadrupole mode". Due to the orientation of the chip in the cavity, the vertically and horizontally orientated cavity modes will only be sensitive to voltage oscillations across pads 1 and 3 or 2 and 4. This ensures that the mode displayed here is decoupled from the cavity modes, and from any loss channels they may incur. We can verify this fact by computing the losses of the different modes, and comparing the losses of mode 2 to the other qubit-like modes of the circuit. We perform this calculation by calling trimon . f_k_A_chi ( pretty_print = True ) which will calculate and return the loss rates of the modes, along with their eigenfrequencies, anharmonicities and Kerr parameters. Setting the keyword argument pretty_print to True prints a table containing all this information, which is shown in Fig. 10. To be succinct, we have not shown the table providing the cross-Kerr couplings. By using the show_normal_mode method to plot all the other modes of the circuit, and noting where currents or voltages are majoritarily located, we can identify each mode with the schematics provided in Fig. 10. The three lowest frequency modes are located in the trimon chip, and we notice that as expected the quadrupole mode 2 has a loss rate (due to resistive losses in the cavity modes) which is three orders of magnitude below the other two. Despite this apparent decoupling, the quadrupole mode will still be coupled to both cavity mode through the cross-Kerr coupling, given by twice the square-root of the product of the quadrupole and cavity mode anharmonicities.

Appendix B: Circuit quantization overview
In this section we summarize the quantization method used in QuCAT, which is an expansion on the work of Ref. (Nigg et al., 2012). This approach is only valid in the weak anharmonic limit, where charge dispersion is negligible. See (Nigg et al., 2012) or Sec. D.2 for a detailed discussion of this condition.
The idea behind the quantization method is as follows. We first consider the "linearized" circuit. This is a circuit where the junctions are replaced by their Josephson . 11 Example of equivalent circuit construction to prepare for quantization. We use the same example as used in Fig 1-4. (a) The circuit is linearized by replacing the junction with an inductance Lj. The circuit is characterized at the nodes of the junction by its admittance Yj. (b) In the limit of small dissipation, this circuit is equivalent to a series combination of RLC resonators.
inductances L j = φ 2 0 /E j , where E j is the Josephson energy and the reduced flux quantum is given by φ 0 = /2e. We determine the oscillation frequencies and dissipation rates of the different normal modes of this linearized circuit. Then, we calculate the amplitude of phase oscillations across each junction when a given mode is excited. This will determine how non-linear each mode is. All this information will finally allow us to build a Hamiltonian for the circuit.

Circuit simplification to series of RLC resonators (Foster circuit)
The eigenfrequencies and non-linearity of each mode is obtained by transforming the linearized circuit to a geometry we can easily analyze. We will first describe this process assuming there is only a single junction in the circuit, the case of multiple junctions will follow. We consider the example circuit of Fig. 1. After replacing the junction with its Josephson inductance, we determine the admittance Y j (ω) = I j (ω)/V j (ω) evaluated at the nodes of the junction. This admittance is the inverse of the impedance measured at the nodes of the junction. It relates the amplitude |V j | and phase θ(V j ) of the voltage oscillating at frequency ω that would build up across the junction if one would feed a current oscillating at ω with amplitude |I j | and phase θ(I j ) to one of its nodes through a infinite impedance current source. In Fig. 11(a) we show a schematic describing this quantity. In the case where all normal modes of the circuit have small dissipation rates, this circuit has an approximate equivalent shown in Fig. 11(b), consisting of a series of RLC resonators (Solgun et al., 2014). By equivalent, we mean that the admittance Y j of the circuit is approximatively equal to that of a series combination of RLC resonators which each have an admittance Each RLC resonator represents a normal mode of the circuit, with resonance frequency ω m = 1/ √ L m C m , and dissipation rate κ m = 1/R m C m . Since this equivalent circuit comes from an extension of Foster's reactance theorem (Foster, 1924) to lossy circuits, we call this the Foster circuit.

Hamiltonian of the Foster circuit
The advantage of this circuit form, is that it is easy to write its corresponding Hamiltonian following standard quantization methods (see Ref. (Vool and Devoret, 2017)). In the absence of junction non-linearity, it is given by the sum of the Hamiltonians of the independent harmonic RLC oscillators: The annihilation operatorâ m for photons in mode m is related to the expression of the phase difference between the two nodes of the oscillator where ϕ zpf,m are the zero-point fluctuations in phase of mode m. The total phase difference across the Josephson junctionφ j is then the sum of these phase differencesφ j = mφ m,j , and we can add the Junction nonlinearity to the Hamiltonian Since the linear part of the Hamiltonian corresponds to the circuit with junctions replaced by inductors, the linear part already contains the quadratic contribution of the junction potential ∝φ 2 j , and it is subtracted from the cosine junction potential.

Calculating Foster circuit parameters
Both ω m and κ m can be determined from Y (ω) since we have Y (ω m +iκ m /2) = 0 for low loss circuits. This can be proven by noticing that the admittance Y m of mode m has two zeros at and ζ * m . The approximate equality holds in the limit of large quality factor Q m = R m / L m /C m 1. From Eq. (B1) we see that the zeros of Y are exactly the zeros of the admittances Y k . The solutions of Y (ω) = 0, which come in conjugate pairs ζ m and ζ * m , thus provide us with both resonance frequencies ω m = Re[ζ m ] and dissipation rates κ m = 2Im[ζ m ].
Additionally, we need to determine the effective capacitances C m in order to obtain the zero-point fluctuations in phase of each mode. We focus on one mode k, and start by rewriting the admittance in Eq. (B1) as Its derivative with respect to ω is The capacitance is thus approximatively given by

Multiple junctions
When more than a single junction is present, we start by choosing a single reference junction, labeled r. All junctions will be again replaced by their inductances, and by using the admittance Y r across the reference junction, we can determine the Hamiltonian including the nonlinearity of the reference junction through the procedure described above.
In this section, we will describe how to obtain the Hamiltonian including the non-linearity of all other junctions toô whereφ j is the phase across the j-th junction. This phase is determined by first calculating the zero-point fluctuations in phase ϕ zpf,m,r through the reference junction r for each mode m given by Eq. (B4). For each junction j, we then calculate the (complex) transfer function T jr (ω) which converts phase in the reference junction to phase in junction j. We can then calculate the total phase across a junction j with respect to the reference phase of junction r, summing the contributions of all modes and both quadratures of the phasê The definition of phase (Vool and Devoret, 2017) where v j is the voltage across junction j translates in the frequency domain to ϕ j (ω) = iωφ −1 0 V j (ω). Finding the transfer function T jr for phase is thus equivalent to finding a transfer function for voltage T jr (ω) = V j (ω)/V r (ω). This is a standard task in microwave network analysis (see Sec. C.3 for more details).

Further treatment of the Hamiltonian
The cosine potential in Eq. (B11) can be expressed in the Fock basis by Taylor expanding it around small values of the phase. This yieldŝ which is the form returned by the QuCAT hamiltonian method. By keeping only the fourth power in the Taylor expansion and performing first order perturbation theory, we obtain Where the anharmonicity or self-Kerr of mode m is as returned by the anharmonicites method, where is the contribution of junction j to the total anharmonicity of a mode m. The cross-Kerr coupling between mode m and n is Both self and cross-Kerr parameters are computed by the kerr method. Note in Eq. B14 that the harmonic frequency of the Hamiltonian is shifted by A m and n =m χ nm /2. The former comes from the change in Josephson inductance induced by phase fluctuations of mode m. The latter is called the Lamb shift (Gely et al., 2018) and is induced by phase fluctuations of the other modes of the circuit.

Appendix C: Algorithmic methods
There are three calculations to accomplish in order to obtain all the parameters necessary to write the circuit Hamiltonian. We need: • the eigen-frequencies ω m and loss rates κ m fulfilling Y r (ζ m = ω m + iκ m /2) = 0 where Y r is the admittance across a reference junction • the derivative of this admittance evaluated at ζ m • the transfer functions T jr between junctions j and the reference junction r In this section, we cover the algorithmic methods used to calculate these three quantities 1. Resonance frequency and dissipation rate

a. Theoretical background
In order to obtain an expression for the admittance across the reference junction, we start by writing the set of equations governing the physics of the circuit. We first determine a list of nodes, which are points at which circuit components connect. Each node, labeled n, is assigned a voltage v n . We name r ± the positive and negative nodes of the reference junction.
We are interested in the steady-state oscillatory behavior of the system. We can thus move to the frequency domain, with complex node voltages |V n (ω)|e i(ωt+θ(Vn(ωn))) , fully described by their phasors, the complex numbers V n = |V n (ω)|e iθ(Vn(ωn)) . In this mathematical construct, the real-part of the complex voltages describes the voltage one would measure at the node in reality. Current conservation dictates that the sum of all currents arriving at any node n, from the other nodes k of the circuit should be equal to the oscillatory current injected at node n by a hypothetical, infinite impedance current source. This current is also characterized by a phasor I n . This can be compactly written as where k label the other nodes of the circuit and Y nk is the admittance directly connecting nodes k and n. Note that in this notation, if a node k 1 can only reach node n through another node k 2 , then Y nk1 = 0. Inductors (with inductance L), capacitors (with capacitance C) and resistors (with resistance R) then have admittances 1/iLω, iCω and 1/R respectively. Expanding Eq. C1 yields which can be written in matrix form as Since voltage is the electric potential of a node relative to another, we still have the freedom of choosing a ground node. Equivalently, conservation of currents imposes that current exciting that node is equal to the sum of currents entering the others, there is thus a redundant degree of freedom in Eq.(C3). For simplicity, we will choose node 0 as ground. Since we are only interested in the admittance across the reference junction, we set all currents to zero, except the currents entering the positive and negative reference junction nodes: I r+ and I r− = −I r+ respectively. The admittance is defined by For Y r = 0, Eq C4 has a solution for only specific values of ω = ζ m . These are the values which make the admittance matrix singular, i.e. which make its determinant zero The determinant is a polynomial in ω, so the problem of finding ζ m = ω m + iκ m /2 reduces to finding the roots of this polynomial. Note that plugging ζ m into the frequency domain expression for the node voltages yields V k (ζ m )e iωmt e −κmt/2 , such that the energy ∝ v k (t) * v k (t) ∝ e −κmt decays at a rate κ m , which explains the division by two in the expression of ζ m . Also note, that we would have obtained equation Eq. (C6) regardless of the choice of reference element.

b. Algorithm
We now describe the algorithm used to determine the solutions ζ m = ω m + iκ m /2 of Eq. C6. As an example, we consider the circuit of Fig. 12(a) that a user would have built with the GUI. The algorithm is as follows 12 Example circuit to illustrate the mode frequency finder algorithm. (a) Example of a circuit built through the GUI by a user. (b) Application of the first step of the algorithm which removes the wires and grounds to obtain a minimal number of nodes without removing any components.
1. Eliminate wires and grounds. In this case, nodes N 0 , N 5 would be grouped under a single node labeled 0 and nodes N 1 , N 2 , N 3 would be grouped under node 1, we label node N 4 node 2, as shown in Fig. 12(b).
2. Compute the un-grounded admittance matrix. For each component present between the different couples of nodes, we append the admittance matrix with the components admittance. The matrix is then multiplied by ω such that all components are polynomials in ω, ensuring that the determinant is also a polynomial. In this example, the matrix is 3. Choose a ground node. The node which has a corresponding column with the most components is chosen as the ground node (to reduce computation time). These rows and columns are erased from the matrix, yielding the final form of the admittance matrix 4. Compute the determinant. Even if the capacitance, inductance and resistance were specified numerically, the admittance matrix would still be a function of the symbolic variable ω. We thus rely on a symbolic Berkowitz determinant calculation algorithm (Berkowitz, 1984;Kerber, 2009) implemented in the Sympy library through the berkowitz_det function. In this example, one would obtain 5. Find the roots of the polynomial. Whilst the above steps have to be performed only once for a given circuit, this one should be performed each time the user edits the value of a component. The rootfinding is divided in the following steps as prescribed by Ref. (Press et al., 2007).
Diagonalize the polynomials companion matrix (Horn and Johnson, 1985) to obtain an exhaustive list of all roots of the polynomial. This is implemented in the NumPy library through the roots function.
Refine the precision of the roots using multiple iterations of Halley's gradient based root finder (Press et al., 2007) until iterations do not improve the root value beyond a predefined tolerance given by the Qcircuit argument r oo t _r el a ti ve _ t o l e r a n c e . The maximum number of iterations that may be carried out is determined by the Qcircuit argument root_max_iterations . If the imaginary part relative to the real part of the root is lower than the relative tolerance, the imaginary part will be set to zero. The relative tolerance thus sets the highest quality factor that QuCAT can detect.
Remove identical roots (equal up to the relative tolerance), roots with negative imaginary or real parts, 0-frequency roots, roots for which Y l (ω m ) < 0 for all l, where Y l is the admittance evaluated at the nodes of an inductive element l, and roots for which Q m <Qcircuit . Q_min. The user is warned of a root being discarded when one of these cases is unexpected.
The roots ζ m obtained through this algorithm are accessed through the method eigenfrequencies which returns the oscillatory frequency in Hertz of all the modes Re[ζ m ]/2π or loss_rates which returns 2Im[ζ m ]/2π.

Derivative of the admittance
The zero-point fluctuations in phase ϕ zpf,m,r for each mode m across a reference junction r is the starting point to computing a Hamiltonian for the non-linear potential of the Junctions. As expressed in Eq. (B4), this quantity depends on the derivative Y r of the admittance Y r calculated at the nodes of the reference element. In this section we first cover the algorithm used to obtain the admittance at the nodes of an arbitrary component. From this admittance we then describe the method to obtain the derivative of the admittance on which ϕ zpf,m,r depends Finally we describe how to choose a (mode-dependent) reference element.

a. Computing the admittance
Here we describe a method to compute the admittance of a network between two arbitrary nodes. We will continue using the example circuit of Fig. 12, assuming we want to compute the admittance at the nodes of the inductor.
1. Eliminate wires and grounds as in the resonance finding algorithm, nodes N 0 , N 5 would be grouped under a single node labeled 0 and nodes N 1 , N 2 , N 3 would be grouped under node 1, we label node N 4 node 2. We thus obtain Fig. 14(a) 2. Group parallel connections. Group all components connected in parallel as a single "admittance component" equal to the sum of admittances of its parts. In this way two nodes are either disconnected, connected by a single inductor, capacitor, junction or resistor, or connected by a single "admittance component".
3. Reduce the network through star-mesh transformations. Excluding the nodes across which we want to evaluate the admittance, we utilize the star-mesh transformation described in Fig. 13 to reduce the number of nodes in the network to two. If following a star-mesh transformation, two components are found in parallel, they are grouped under a single "admittance component" as described previously. For a node connected to more than 3 other nodes the star-mesh transformation will increase the total number of components in the circuit. So we start with the least-connected nodes to maintain the total number of components in the network to a minimum. In this example, we want to keep nodes 0 and 2, but remove node 1, a start-mesh transform leads to the circuit of Fig. 14(b) then grouping parallel componnets leads to (c).
4. The admittance is that of the remaining "admittance component" once the network has been com-pletely reduced to two nodes.
The symbolic variables at this stage (Sympy Symbols) are ω, and the variables corresponding to any component with un-specified values. 13 Star-mesh transform. A node N connected to nodes A, B, C, .. through admittances YA, YB, ... can be eliminated if we interconnect nodes A, B, C, .. with impedances YAB, YAC , YBC , ... given by YXY = YX YY / M YM . We show the 5 node case, the initial network on the left is called the "star", which is then transformed to the "mesh" on the right, reducing the total number of nodes by 1.

(a) (b) (c)
14 Example to illustrate the admittance calculation algorithm. (a) Example of a circuit built through the GUI by a user, after removal of all wires and grounds. (b) Application of the star-mesh transformation to remove node 1. (c) After each application of the star-mesh transformation, parallel connections are grouped together. Only the two nodes across which we want to compute the admittance remain, the admittance is that of the remaining "admittance component".

b. Differentiating the admittance
The expression for the admittance obtained from the above algorithm will necessarily be in the form of multiple multiplication, divisions or additions of the admittance of capacitors, inductors or resistors. It is thus possible to transform Y to a rational function of ω Y (ω) = P (ω) Q(ω) = p 0 + p 1 ω + p 2 ω 2 + ... q 0 + q 1 ω + q 2 ω 2 + ...

c. Choice of reference element
For each mode m, we use as reference element r the inductor or junction which maximizes ϕ zpf,m,r as specified by Eq. B4. This corresponds to the element where the phase fluctuations are majoritarily located. We find that doing so considerably increases the success of evaluating Y (ω m ). As an example, we plot in Fig. 15 the zero-point fluctuations in phase ϕ zpf,m,r of the transmonlike mode, calculated for the circuit of Fig. 1, with the junction or inductor as reference element. What we find is that if the coupling capacitor becomes too small, resulting in modes which are nearly totally localized in either inductor or junction, choosing the wrong reference element combined with numerical inaccuracies leads to unreliable values of ϕ zpf,m,r .

Transfer functions
In this section, we describe the method used to determine the transfer function T jr between a junction j and the reference junction r. This quantity can be computed from the ABCD matrix (Pozar, 2009). The ABCD matrix relates the voltages and currents in a two port network where the convention for current direction is described in Fig. 16. By constructing the network as in Fig. 16, with the reference junction on the left and junction j on the right, the transfer function is given by To determine A, we first reduce the circuit using starmesh transformations (see Fig.13), and group parallel connections as described in the previous section, until only the nodes of junctions r and j are left. If the junctions initially shared a node, the resulting circuit will be equivalent to the network shown in Fig.17 (a). In this case, If the junctions do not share nodes, the resulting circuit will be equivalent to the network shown in Fig.17 where some admittances may be equal to 0 to represent open circuits. To compute the ABCD matrix of this resulting circuit, we make use of the property illustrated in Fig. 16: the ABCD matrix of a cascade connection of two-port networks is equal to the product of the ABCD matrices of the individual networks. We first determine the ABCD matrix of three parts of the network (separated by dashed line in Fig.17) such that the ABCD matrix of the total network reads where the A and B coefficients of the middle part of the network arẽ The ABCD matrix for the middle part of the circuit is derived in Sec. 10.11 of Ref. (Arshad, 2010), and the ABCD matrices for the circuits on either sides are provided in Ref. (Pozar, 2009).
FIG. 17 Networks after star-mesh reduction. The two non-trivial situations reached after applying star-mesh transformations to a network to obtain Trj.
This method is also applied to calculate the transfer function to capacitors, inductors and resistors, notably to visualize the normal mode with the show_normal_mode function.

Alternative algorithmic methods
Since symbolic calculations are the most computationally expensive steps in a typical use of QuCAT, we cover in this section some alternatives to the methods previously described, and the reasons why they were not chosen. a. Eigen-frequencies from the zeros of admittance One could solve Y r (ω) = 0 where Y r is the admittance computed as explained in Sec. C.2. Providing good initial guesses for all values of the zeros ζ m can be provided, a number of root-finding algorithms can then be used to obtain final values of ζ m .
A set of initial guesses could be obtained by noticing that Y r is a rational function of ω. Roots of its numerator are potentially zeros of Y, and a complete set of them is easy to obtain through a diagonalization of the companion matrix as discussed before. Note that if these roots are roots of the denominator with equal or higher multiplicity, then they are not zeros of Y. They can, however, make good initial guesses of a root-finding algorithm run on Y r . This requires a simplification of Y r , as computed through star-mesh transforms, to its rational function form. We find this last step to be as computationally expensive as obtaining a determinant.
A different approach, which does not require using a root-finding algorithm on Y r , is to simplify the rationalfunction form of Y such that the numerator and denominator share no roots. This can be done by using the extended Euclidian algorithm to find the greatest common polynomial divisor (GCD) of the numerator and denominator. However, the numerical inaccuracies in the numerator and denominator coefficients may make this method unreliable.
The success of both of these approaches is dependent on determining a good reference component r, which may be mode-dependent (see Fig. 15). This reference component is difficult to pick at this stage, when the mode frequencies are unknown. b. Finite difference estimation of the admittance derivative Rather than symbolically differentiating the admittance, one could use a numerical finite difference approximation, for example Y r can be obtained through star-mesh reductions, or from a resolution of Eq. C4. But finding a good value of δω is no easy task. As an example, we consider the circuit of Fig. 1, where we have taken as reference element the junction. As shown in Fig. 18, when the resonator and transmon decouple through a reduction of the coupling capacitor, a smaller and smaller δω is required to obtain Y r evaluated at the ζ 1 of the resonator-like mode. We have tried making use of Ridders method of polynomial extrapolation to try and reliably approach the limit δx → 0 (Press et al., 2007). However at small coupling capacitance, it always converges to the slower varying background slope of Y r , without any way of detecting the error.  Fig. 1 for different values of the coupling capacitor Cc and for Lj = 12 nH. Resonances correspond to the frequencies ω at which the admittance crosses 0, and the calculation of zero-point fluctuations depends on the derivative Im[Y ] at that point. As the resonator and transmon parts of the circuit decouple, Im[Y ] becomes larger, requiring a lower δω if the admittance is to be determined through Eq. (C18). In extreme cases (see lower panel), when the derivative is very large, the smaller variation in the background slope may be mistaken for the slope at a resonance.

c. Transfer functions from the admittance matrix
Calculating the transfer function T ij could alternatively be carried out through the resolution of the system of equations (C5). The difference in voltage of a reference elements nodes would first have to be fixed to the zero-point fluctuations computed with the method of Sec. (C.2). These equations would have to be resolved at each change of system parameters and for each mode, with ω replaced in the admittance matrix by its corresponding value for a given mode. This is to be balanced against a single symbolic derivation of T ij through starmesh transformations, and fast evaluations of the symbolic expression for different parameters. In this section we ask the question: how big a circuit can QuCAT analyze? To address this, we first consider the circuit of Fig. 7(c), and secondly the same circuit with resistors added in parallel to each capacitor. As the number of (R)LC oscillators representing the modes of a CPW resonator is increased, we measure the time necessary for the initialization of the Qcircuit object. This is typically the most computationally expensive part of a QuCAT usage, limited by the speed of symbolic manipulations in Sympy.
These symbolic manipulations include: • Calculating the determinant of the admittance matrix • Converting that determinant to a polynomial • Reducing networks through star-mesh transformations both for admittance and transfer function calculations • Rational function manipulations to prepare the admittance for differentiation Once these operations have been performed, the most computationally expensive step in a Qcircuit method is finding the root of a polynomial (the determinant of the admittance matrix) which typically takes a few milliseconds.
The results of this test are reported in Fig. 19. We find that relatively long computation times above 10 seconds are required as one goes beyond 10 circuit nodes. Due to an increased complexity of symbolic expressions, the computation time increases when resistors are included. For example, the admittance matrix of a non-resistive circuit will have no coefficients proportional to ω, only ω 2 and only real parts, translating to a polynomial in Ω = ω 2 which will have half the number of terms as a resistive circuit. However, we find that this initialization time is also greatly dependent on the circuit connectivity, and this test should be taken as only a rough guideline. On the vertical axis, we show the time necessary to initialize the Qcircuit object, which is the computationally expensive part of a typical QuCAT user case. This is plotted as a function of the number of nodes in the circuit. The test circuit used here is the multi-mode circuit of Fig. 7(c), optionally with a resistor in parallel of each capacitor. The number of nodes are increased by adding modes to the circuit. Most of the computational time is spent in the symbolic manipulations performed with the sympy library.
Making QuCAT compatible with the analysis of larger circuits will inevitably require the development of more efficient open-source symbolic manipulation tools. The development of the open-source C++ library SymEngine https://github.com/symengine/symengine, together with its Python wrappers, the symengine.py project https://github.com/symengine/symengine.py, could lead to rapid progress in this direction. An enticing prospect would then be able to analyze the large scale cQED systems underlying modern transmon-qubit-based quantum processors (Kelly et al., 2019). One should keep in mind that an increase in circuit size translates to an increase in the number of degrees of freedom of the circuit and hence of the Hilbert space size needed for further analysis once a Hamiltonian has been extracted from Qu-CAT. On the x-axis, we vary the approximate anharmonicity EC = e 2 /2C with respect to the frequency ω0 = 1/ LjC. For each value, we plot a Hilbert space size, and order of Taylor-expansion of the junction cosine-potential. Incrementing these values produces less than a 0.1 percent change in the first two transition frequencies obtained by diagonalizing the Hamiltonian. Beyond a relative anharmonicity of 8, convergence is no longer reached, even for Hilbert space sizes and Taylor expansions up to 100. (c) Frequency of the first transition ωg−e obtained from Hamiltonian diagonalization, relative to the value expected from first order perturbation theory: ω0 − EC / . (d) Anharmonicity ω e−f − ωg−e obtained from Hamiltonian diagonalization relative to the value expected from first order perturbation theory EC / . Black lines correspond to a diagonalization in the harmonic Fock basis (Eq. (B13)), blue and orange dashed lines correspond to a diagonalization of the Cooper-pair box Hamiltonian with gate charges of 0 and 1/2 respectively. The harmonic Fock basis provides reliable results up to approximatively 6 percent anharmonicity.
In this section we study the limits of the current quan-tization method used in QuCAT. More specifically, we study the applicability of the basis used to express the Hamiltonian, that of Fock-states of harmonic normal modes of the linearized circuit. To do so, we use the simplest circuit possible (Fig. 20(a)), the parallel connection of a Josephson junction and a capacitor. As the anharmonicity of this circuit becomes a greater fraction of its linearized circuit resonance, the physics of the circuit goes from that of a Transmon to that of a Cooper-pair box (Koch et al., 2007), and the Fock-state basis becomes inadequate. This test should be viewed as a guideline for the maximum acceptable amount for anharmonicity. We find that when the anharmonicity exceeds 6 percent of the eigenfrequency, a QuCAT generated Hamiltonian will not reliably describe the system. In this test, we vary the ratio of Josephson inductance L j to capacitance C, increasing the anharmonicity expected from first-order perturbation theory (see Eq. B14), called charging energy E C = e 2 /2C. The resonance frequency of the linearized circuit ω 0 = 1/ L j C is maintained constant. For each different charging energy, we use the hamiltonian method to generate a Hamiltonian of the system. We are interested in the order of the Taylor expansion of the cosine potential, and the size of the Hilbert space, necessary to obtain realistic first and second transition frequencies of the circuit, named ω g−e and ω e−f respectively. To do so, we increase the order of Taylor expansion, and for each order we sweep through increasing Hilbert space sizes. In Fig. 20(b), we show the values of these parameters at which incrementing them would not change ω g−e and ω e−f by more than 0.1 percent. Beyond a relative anharmonicity E C / ω 0 of 8 percent, such convergence is no longer reached, even for cosine expansion orders and Hilbert space sizes up to 100.
Up to the point of no convergence, we compare the results obtained from the diagonalization in the harmonic Fock basis (Hamiltonian generated by QuCAT), with a diagonalization of the Cooper-pair box Hamiltonian. In regimes of higher anharmonicity, the system becomes sensitive to the preferred charge offset between the two plates of the capacitor N g (expressed in units of Cooperpair charge 2e) imposed by the electric environment of the system. The Cooper-pair box Hamiltonian takes this into accountĤ where |N is the quantum state of the system where N Cooper-pairs have tunneled across the junction to the node indicated in Fig. 20(a). For more details on Cooperpair box physics and the derivation of this Hamiltonian, refer to Ref. (Schuster, 2007). We diagonalize this Hamiltonian in a basis of 41 |N states.
We find that beyond 6 percent anharmonicity, the Cooper-pair box Hamiltonian becomes appreciably sensitive to N g and diverges from the results obtained in the Fock basis. This corresponds to E j /E C 35 at which the charge dispersion (the difference in frequency between 0 and 0.5 charge offset) is 4 × 10 −5 and 1 × 10 −3 for the first two transitions respectively.
Beyond 8 percent anharmonicity, one cannot reach convergence with the Fock basis and just before results diverge considerably from that of the Cooper-pair box Hamiltonian. This corresponds to E j /E C 20 at which the charge dispersion is 1.5 × 10 −3 and 3 × 10 −2 for the first two transitions respectively. A possible extension of the QuCAT Hamiltonian could thus include handling static offsets in charge and different quantization methods, for example quantization in the charge basis to extend QuCAT beyond the scope of weakly-anharmonic circuits.