CircuitQ: An open-source toolbox for superconducting circuits

We introduce CircuitQ, an open-source toolbox for the analysis of superconducting circuits implemented in Python. It features the automated construction of a symbolic Hamiltonian of the input circuit and a dynamic numerical representation of the Hamiltonian with a variable basis choice. Additional features include the estimation of the T1 lifetimes of the circuit states under various noise mechanisms. We review previously established circuit quantization methods and formulate them in a way that facilitates the software implementation. The toolbox is then showcased by applying it to practically relevant qubit circuits and comparing it to specialized circuit solvers. Our circuit quantization is applicable to circuit inputs from a large design space, and the software is open-sourced. We thereby add an important resource for the design of new quantum circuits for quantum information processing applications.


I. INTRODUCTION
Superconducting circuits are one of the most versatile and promising platforms in the development of chipbased quantum processors [1,2]. The development of cutting-edge qubit designs is currently driven by the effort to realize quantum computers with long coherence times as well as high-fidelity control and readout, all in a scalable design [3]. Combining these requirements is a grand challenge for quantum hardware design. Therefore, considerable effort is invested into the study of new and improved qubit circuits for quantum information processing applications [4][5][6][7].
A major part of the analysis of superconducting circuits is the construction of a quantum model to describe the system theoretically. Such a model is obtained by using general methods to construct the corresponding Hamiltonian [8][9][10]. A numerical implementation of the algebraic description is then needed to analyze the quantum properties of the circuit. A number of open source software packages has been developed for this purpose. The library scQubits [11], for example, simulates qubits from a specific set of circuits. The package QuCAT [12] offers a more general circuit input, as it permits a combination of Josephson junctions, inductances, capacitances and resonators. The quantization is performed in the basis of normal modes, which is suitable for weakly anharmonic systems. Another useful toolbox is provided by the SuperQuantPackage [13], which includes an algorithm to provide the user with a numerical representation * philipp.aumann@uibk.ac.at † wolfgang.lechner@uibk.ac.at of the Hamiltonian for a given input circuit. It performs a coordinate transformation prior to quantization. Qiskit Metal [14] and KQCircuits [15] enable the analysis of superconducting circuits based on their physical layout on the chip. While such software packages have been proven to be useful for specific circuit design tasks, we expand on prior work by presenting a toolbox that works for a generic variety of circuits, determines a symbolic and numerical Hamiltonian, provides an automated choice of implementation basis and includes a measure for several T 1 contributions.
In this work, we provide a structured review of the superconducting circuit quantization procedure implemented in the software toolbox CircuitQ. This provides an insight into the software implementation, but also serves as a more general review of the quantization process of superconducting circuits by constructing the Hamiltonian. CircuitQ is written in Python and can be used to analyze superconducting circuits that are a userdefined combination of Josephson junctions, linear inductances, and capacitances. It takes the circuit and optionally the circuit component parameters as an input and returns the quantum physical properties of the circuit, particularly the corresponding Hamiltonian in symbolic and numerical form. Figure 1 provides an overview of this conceptual procedure for three example circuits. Depending on the shape of the inductive potential, CircuitQ can dynamically perform the numerical implementation in the charge basis, the flux basis, or in a mixture of both. Therefore, it provides a toolbox for circuits comprising different parameter regimes. As detailed in section IV, this implementation works well for few-node circuits, where the direct implementation of the node variables with the flux and charge basis is a natural choice, FIG. 1. Overview of the functionality of CircuitQ for three example circuits, as discussed in section IV. The general procedure to analyze a superconducting circuit (top row) is to interpret it as a graph whose edges correspond to the circuit elements (middle row). The graph serves as the input of the toolbox, and a symbolic and a numerical Hamiltonian are constructed from it. The potential and eigenspectrum of the Hamiltonian are shown in the bottom row. The symbolic Hamiltonians and parameter values can be found in section IV, and the code to generate these instances is listed in Appendix B. (a) Transmon circuit. The last row depicts the scaled absolute value squared of the wavefunctions of the lowest five eigenstates as a function of the flux variable Φ1, which is the node flux variable connected to node 1 of the respective graph in the middle row, while the scaling factor Φ0 is the flux quantum. Each wave function is offset by its corresponding eigenvalues. The potential of the Hamiltonian is shown as a black line. (b) Fluxonium circuit. The scaled wavefunctions of the eigenstates and potential energy are shown in the same fashion as in (a). (c) Persistent-current flux qubit circuit. Due to the additional node, the problem becomes two-dimensional. Similar to before, Φ1 and Φ2 are the node variables of node 1 and 2 in the circuit graph. We show an overlap of two contour plots, where one depicts the scaled absolute value squared of the wavefunction of the lowest eigenstate and the other visualizes the potential.
whereas the limits of this implementations are reached for more complex circuits. The interface offers the possibility for a general circuit input. A variety of features and degrees of freedom can be adjusted by the user, such as the circuit composition, the component parameters, ground nodes, offset charges, and loop fluxes. In order to analyze the circuit in view of noisy environments, we implemented an estimation of the T 1 lifetime of an eigenstate by considering three common relaxation mechanisms. Finally, we showcase the software by comparing its accuracy to specialized circuit solvers for several prominent qubit circuits. We provide community access to CircuitQ by making the code and documentation publicly available on GitHub [16].
The article is organized as follows. In section II, we present the procedure to generate the symbolic Hamiltonian from a given input circuit. Subsequently, the numerical implementation of the Hamiltonian and features of the toolbox are presented in section III. Lastly, applications of the toolbox to practically relevant circuit examples are provided in section IV.

II. FROM THE CIRCUIT TO THE SYMBOLIC HAMILTONIAN
The first step in the analysis of the quantum properties of a superconducting circuit is the construction of the circuit Hamiltonian. Our software implementation automates the process to derive the symbolic Hamiltonian. A multitude of techniques for circuit quantization has been developed, including several based on the method of nodes [8][9][10]17], black box quantization [18,19] and others [20,21]. Here we follow the method of nodesbased approach, which is generally applicable to any superconducting circuit that includes capacitances, inductances, and Josephson junctions, given the realistic condition that spurious capacitances exist between all circuit nodes.
The starting point of circuit quantization is a circuit diagram such as the one shown in figure 2, which is a lumped-element representation of a prospective on-chip microfabricated device. A similar version of this circuit example can be found in reference [22]. The circuit diagram can be seen as a graph with the circuit elements on its edges. The nodes are then home to the conjugate charge and flux variables, which represent the charge stored on capacitances connected to the node, and the flux along a specific path from the node to ground, respectively. They are denoted as the node charge q i and node flux Φ i , while q and Φ are the vectors of all charge and flux variables ordered by node index. One or more ground nodes can be specified by the user when initializing an instance of the CircuitQ class. Additionally, all active nodes with only one neighbouring node are added to the ground nodes if the neighbouring node is ungrounded. An active node is a node that is connected to a capacitance as well as an inductive element. If no ground node could be specified, an active node is chosen to be the ground node. Node 0 is the ground node of the circuit depicted in figure 2. Ground nodes do not appear in the constructed Hamiltonian and the associated variables are removed from vector q and Φ.
In the software, the circuit graph is implemented as a MultiGraph instance of the NetworkX package [23]. A code example to initialize the circuit in figure 2 is given in Appendix C. We note that a graph is automatically simplified using the common rules for parallel and series capacitors.
The definition of the node fluxes requires the choice of a unique path from each node in the circuit graph to ground. The set of such paths for all nodes is termed the spanning tree of the graph, which cannot contain loops. There may be multiple ways to choose it for a given circuit. One such choice is highlighted as a green sub-graph in figure 2. The choice of spanning tree is equivalent to setting a gauge and therefore does not change the physics of a circuit [9]. However, it does affect the circuit quantization procedure and the form of the Hamiltonian. In this work, we follow common practice and route the spanning tree through capacitive circuit elements only. Since each node pair is connected by spurious capacitances, a spanning tree can always be defined in such a way. For the implementation, CircuitQ makes use of the spanning tree functionality of the NetworkX package.
In order to determine the contribution of an inductive circuit element to the Hamiltonian, one needs to evaluate the flux difference across the respective edge. An inductive loop consisting of inductances and junctions can encircle an external flux, and the boundary condition has to be fulfilled that all the fluxes around a loop sum to a multiple of the flux quantum Φ 0 = 2e . We highlight one such loop in orange in figure 2. The external flux enters the circuit Hamiltonian by being added to the flux difference across an inductive element between two of the nodes that are part of the loop.
In CircuitQ, the automated evaluation of flux differences is started by defining the set of inductive edges L and splitting it into the subset S that is in parallel to spanning tree edges and into the remaining edges B = L − S. As we allow for multiple parallel inductive elements between two nodes, a restriction needs to be introduced that only the first parallel inductive edge is added to S. The routine then iteratively steps through the edges in B. If the current edge does not close an inductive loop, we add it to a subset B o ⊆ B and do not assign an external flux to it. In case it does close an inductive loop, we add it to B c ⊆ B and do assign an external flux to it. Using the fluxoid quantization condition, the directionality of the edges in S can be chosen such that we obtain the following relation for the edge flux Φ eijn of the n-th inductive edge e ijn between nodes i and j: For the circuit in figure 2, this procedure identifies the loop fluxΦ 120 and applies it to the flux difference between nodes 1 and 2. We note that the loop fluxes are degrees of freedom that can be used to tune the circuit properties. At the same time, they open a path for undesired fluctuations from the environment to couple to the circuit. CircuitQ determines a set of loop fluxes automatically when initializing an instance for a given circuit. These fluxes are treated as conventional circuit parameters whose numerical values can be specified by the user.
Given the relationship between branch and node fluxes, we can explicitly state the inductive potential of the Hamiltonian, which is the sum of linear inductive and Josephson potentials: where N ij L is the number of linear inductors between node i and j and N ij J is the number of Josephson junctions, while I L and I J represent the set of node pairs connected by inductors and Josephson junctions. Our formulation of the potential expands upon prior work on general circuit quantization formulations in that it allows for multiple parallel inductive elements per node pair. While this is relevant for Josephson junctions, for example in modelling a frequency tunable transmon, where two Josephson junctions form a SQUID loop, multiple parallel linear inductors are of limited practical relevance for quantum information processing applications because they lead to sensitivity of the circuit energy to static flux offsets.
The kinetic energy T cap of the Hamiltonian is given by the capacitive energy of the circuit. In general, T cap takes the form Here, C is the node capacitance matrix, which contains the sum of all capacitances connected to a node as the respective diagonal entry and the negative capacitance between two nodes on the off-diagonals. That is, the i-th diagonal element of the matrix is given by by summing over all capacitances C in connected to node i, while the off-diagonal elements are given by for i = j and capacitance C ij linking nodes i and j. The capacitances are invariant under permutation of the indices: C ij = C ji . We note that the rows and columns of C that are associated to ground nodes are removed from the matrix. The form of the kinetic energy arises in the Legendre transformation of the circuit Lagrangian.
Charge offsets on a node are taken into account by directly adding the offset to the respective charge operator. The circuit Hamiltonian in terms of the conjugate coordinates Φ and q is given by the sum of the kinetic and potential terms: In CircuitQ, it is returned as a symbolic SymPy object [24]. The Hamiltonian that is constructed for the example circuit in figure 2 can be found in Appendix E.
For the quantum physical treatment of the system described by the Hamiltonian in equation 6, we perform the usual quantization procedure by promoting the conjugate variables to operators: Those operators fulfill the canonical commutation relations Generating the Hamiltonian in a symbolic form is the first step towards the automated description of a circuit. The second step is to express the quantized Hamiltonian numerically.

III. FROM THE SYMBOLIC HAMILTONIAN TO ITS NUMERICAL IMPLEMENTATION
The numerical implementation of the symbolic Hamiltonian is crucial for the analysis of the quantum properties of the input circuit. In this section, we first present important steps of the automated numerical implementation that is part of the toolbox. Subsequently, we describe numerical analysis tools that are implemented in CircuitQ.

A. Implementation
Numerical values for the circuit parameters can be specified individually. This includes capacitances, inductances, Josephson energies, external fluxes and charge offsets. If the values are not specified by the user, they are set to a default value. In addition to the numerical parameter values, the charge and flux operators that appear in the Hamiltonian have to be implemented as numerical matrices. These matrices can be either formulated in the flux basis or in the charge basis, which is analogous to choosing a position or momentum space representation. If the potential is periodic along the direction of a node flux variable, the charge basis is the preferred choice for implementing the variables of this node. The connection of a linear inductance to a node leads to a non-periodic harmonic contribution of the potential, which makes a flux basis implementation for the corresponding node variables more desirable. To distinguish these two cases, we label a node as periodic if there is no linear inductance connected to it and if neighbouring nodes that are connected via a Josephson junction are periodic as well. Periodic node variables are automatically implemented in the charge basis and non-periodic variables in the flux basis.

Flux Basis
To implement a node variable in the flux basis, we confine the numerical flux values to a finite grid with grid spacing δ. The grid length can be decided by the user or is set to a default value otherwise. Consequently, we can assign a diagonal matrix to the flux variable: As the conjugate momentum of the flux, the charge variable can be associated with the derivative with respect To implement the derivative as a Hermitian operator, we use the finite difference method: To generate Hermitian matrices, we have chosen a different discretization for the first and second derivative. The cosine terms in equation 2, referring to the energy contribution of Josephson junctions, can be represented by diagonal matrices in the flux basis, where the diagonal elements are the cosine of the corresponding numerical edge flux value. CircuitQ is capable of working with charge and flux offsets, where the flux offsets Φ are associated with loop fluxes and charge offsetsq with node charges. They are implemented by multiplying them with the identity matrix: A Hilbert space is assigned to every node which is not set to ground. To obtain a numerical description of the full Hamiltonian, these subspaces are combined into a composite space using the tensor product by substitutinĝ Here, the variable corresponding to node i is implemented by placing the respective matrix at the position of the composite space which corresponds to node i. We note that node numbering may change due to the elimination of the ground nodes. The sequence of the nodes is deduced by the algorithm and kept consistent throughout the evaluation of an instance.

Charge Basis
Similar to the flux basis, we restrict the charge variables to a finite grid when using the charge basis. The charge is truncated at a cutoff number of Cooper pairs n cutoff , which leads to the charge grid 2e[−n cutoff , . . . , n cutoff ] = [−q max , . . . , q max ]. In this setting, we can express the charge variable as a diagonal matrix: Flux variables that correspond to periodic nodes appear exclusively in the arguments of the cosine terms in the Hamiltonian. The cosine acts as a hopping operator in the Cooper pair number basis [25]: with |n being the charge state corresponding to n Cooper pairs. We can then use the decomposition of the cosine into complex exponentials, to represent these terms numerically. This procedure is used in scQubits [11]. The exponential of a flux operator describes the tunneling process of a Cooper pair through a Josephson junction, and it can be described as a jump operator in the charge basis [25]: The composite space has to be considered if there are multiple flux variables in the argument of the cosine. In this case, we again make use of the tensor product: The full cosine function can thus be implemented as To account for a flux offsetΦ, the exponential function in equation 19 can be multiplied by the complex scalar e −iΦ . The choice of an appropriate numerical value for the discretization parameters Φ max , δ and n cutoff , which determine the numerical representation of the flux and charge variables, depend on the particular circuit. An appropriate regime can be found by increasing (for the case of Φ max and n cutoff ) or decreasing (for the case of δ) the numerical value of the parameter until convergence is reached, such that the resulting spectrum of the circuit becomes almost invariant under a slight modification of those values.
The numerical grid for the numerical Hamiltonian is generated using the lambdify function of SymPy [24] with the parameters and matrices that have been described in this section as inputs. The final implementation is returned as a sparse matrix in SciPy format [26].
For some analyses, it may be helpful to visualize the eigenstates as a function of the flux variable even when an implementation in the charge basis has been used. For this purpose, CircuitQ provides a method which transforms the eigenvectors from the charge to the flux basis. To transform a state vector, given in the charge basis B q = {|q i } i , to a representation in the flux basis B Φ = {|Φ i } i , the transformation matrix T can be defined, which maps the state vector from the charge to the flux basis. The coefficients of this matrix read: with d being the number of basis vectors. We follow the same procedure in our numerical implementation, however we use a modified transformation matrix which respects the construction of the composite Hilbert space, which, in general, consists of subspaces that are either implemented in the charge or the flux basis. We note that depending on the size of the numerical matrices, this transformation can be numerically demanding and consequently may lead to a bottleneck in computation time.

B. Features for circuit analysis
Since CircuitQ constructs a numerical implementation of the circuit Hamiltonian, it can be used as a tool for the analysis of the quantum properties of superconducting circuits. This includes the energy spectrum of the Hamiltonian and relaxation times of the energy eigenstates.

Spectrum
To calculate the energy spectrum of the numerical Hamiltonian, the toolbox provides a method that returns the lowest eigenstates and eigenenergies of the numerical Hamiltonian matrix. We use the SciPy library for the (partial) diagonalization, which in turn makes use of efficient ARPACK routines [33]. With this functionality, the toolbox can be used to investigate how a change of parameter values -for example an external flux -or a change in the circuit composition affects the energy spectrum and eigenstates.
For the description of the superconducting circuit as a qubit, we associate the lowest two energy levels that have a nonvanishing energy difference with the qubit states |0 and |1 by default. However, it is possible to declare a different state as the excited qubit state manually. The corresponding energy levels should not be degenerate. To operate a circuit as a qubit, a high degree of anharmonicity of its spectrum is desired. CircuitQ provides a method which gives an estimate of the harmonicity of the spectrum in a quantified form.

Relaxation time T1
In order to determine the performance of the circuit as a qubit, it is crucial to study its sensitivity to various noise sources. The qubit can decay to its ground state as a result of its interactions with the environment. The sensitivity to this relaxation process is quantified by the T 1 time. We note that an undesired excitation of the qubit state may also result from such interactions. Table I provides an overview of the noise contributions that can be estimated with CircuitQ. We included relaxation due to quasiparticle tunneling, dielectric loss and flux noise.
a. Quasiparticle tunneling In experiments, significant non-vanishing densities of unpaired electrons could be observed, which are referred to as quasiparticles in this context [34]. Tunneling of such quasiparticles through the junction barrier can lead to a relaxation of the qubit. This effect is separated into two contributions as given in equation 22. The first contribution concerns the junctions in the circuit (sum over j), while the second contribution is associated with every linear inductance in the circuit (sum over l) [27]. The implementation of the  Tc = 1.2 K: Critical temperature of aluminum [29] Dielectric loss [6] T1 SQ i (ωq) = Qcap(ωq )C i 1 + coth ωq 2k B T : Noise spectral density [6,30]   T co witĥ Finally, to evaluate the transition element M eg in equation 26 in the charge basis, we compute b. Dielectric Loss Another noise channel present in superconducting qubits is relaxation due to the fact that the electrical field, which stores capacitive energy, couples to charged fluctuators [36]. The resulting effect on the T 1 time is calculated with equation 23. Here we sum over all capacitors in the circuit.
c. Flux noise The fluctuation of spins on the superconducting material are suspected to be the origin of flux noise [36]. Those fluctuations perturb the magnetic field, which stores inductive energy, effectively leading to fluctuations of the electrical current of the inductive elements. We estimate the corresponding contribution to the relaxation time with equation 24. As described in section II, we assign an external flux to a subset of edges B c of the circuit graph. We therefore sum over all such edges to estimate the flux noise. If the associated element is a Josephson junction, the current operator is given bŷ I = I C sinΦ Φ0 , while the expressionÎ =Φ L is used for linear inductors. We included the possibility to obtain a lower bound on the T 1 estimate by summing not only over the edges in B c but including all inductive edges.
The focus of this toolbox is on pure qubit design without considering qubit control such as state preparation. Therefore, we do not consider noise due to the Purcell effect for now. We also did not include pure dephasing mechanisms explicitly, which can be attributed to the fluctuation of the qubit frequency due to various noise channels. Estimating those dephasing processes would entail the calculation of the derivative of the qubit frequency with respect to the particular noise source. This derivative could be either calculated numerically or even symbolically, depending on the efficiency of those approaches. Adding dephasing processes to the analysis is an important part of our outlook, as it is an essential part of an extensive and general study of superconducting circuits.
In comparison to the related open-source software toolbox scQubits [11], we follow a similar strategy by estimating the coherence times using Fermi's golden rule combined with the specific noise spectral densities from literature. However, our expressions for the noise spectral densities differ in the case of noise due to quasiparticle tunneling and flux noise, where we follow Nguyen et al.'s study of the fluxonium qubit [6]. We also add the contribution due to linear inductances to calculate noise due to quasiparticle tunneling.

IV. DEMONSTRATION AND BENCHMARK
To demonstrate the capabilities of CircuitQ, we use three well known circuits from the literature, i.e., the fixed-frequency transmon [37], the fluxonium [38] and the persistent-current flux qubit [39]. The examples are initialized with the corresponding input graph, from which CircuitQ computes the symbolic and numerical Hamiltonian. The latter can be diagonalized to analyze the spectrum and eigenstates of the system. Figure 1 gives an overview of this procedure while figure 3 depicts the T 1 contributions that are estimated by the toolbox for these circuits. For the noise estimates, we considered all three depolarization channels introduced in section III B. For the fixed frequency transmon circuit, we consider two noise channels: quasiparticle tunneling and dielectric loss. For the chosen circuit parameters, the second contribution is observed to be the limiting factor. The lifetime of planar (2D) transmon fabrications are reported to be limited by dielectric loss [40]. As we simulate the fluxonium at the sweet spotΦ ext = π · Φ 0 here, the quasiparticle noise of the small junction (first term in equation 22) is suppressed and we can ascribe the T 1 decay mostly We calculated the spectrum for those instances with reference implementations (see main text) and compared the transition energy between the ground state and excited state to the corresponding output of CircuitQ. The deviation is scaled by the qubit energy of the respective reference implementation. The test has been repeated for various dimensions of the subsystem matrices, which are the numerical representation of the node fluxes and charges. We note that the total dimension scales exponentially with the number of ungrounded nodes, which leads to longer evaluation times for larger circuits.
to the linear inductance (second term in equation 22). Similar to the transmon, the flux qubit is limited by the dielectric loss when comparing the three relaxation processes. The parameter values used for estimating the lifetime of a flux qubit refer to Qubit B in reference [32], where the lifetime at the sweet spotΦ ext = π · Φ 0 seems to be limited by flux noise. Our findings indicate that our specifications chosen for the flux and charge noise estimates do not resemble this particular experimental setup accurately. However, our estimate for the effective T 1 time lies in the same order of magnitude compared to the findings in this reference. The effective T 1 time for the transmon and fluxonium range in between 10 2 − 10 3 µs, while this value is reduced by one order of magnitude for the flux qubit. These numbers are in accordance with values from literature [3]. We note that the exact estimate of the T 1 time depends on parameters like the density of quasiparticles or the dielectric quality factor. Those parameters depend on the specific realization of the circuits and will vary from experiment to experiment. Although we have chosen representative values, the computed T 1 times should not be understood to be exact for a specific circuit layout but should serve as estimates to classify the sensitivity of a circuit to certain noise channels. In Appendix B, code samples are provided to initialize the instances for the example circuits. The symbolic Hamiltonian that is generated for the transmon by the toolbox is Here, Φ 1 and q 1 are the flux and charge variables of node 1, andq 1 is a charge offset that can be introduced upon initialization (see Appendix B). We associate E J010 with the Josephson energy of the 0-th junction between node 0 and 1, which is shunted by the capacitance C 01 . This notation, which assigns circuit elements like a Josephson junction to Hamiltonian parameters like a Joesphson energy by providing the corresponding edge nodes in the index of the symbols, is kept consistent. For better readability, we do not define all the symbols in the following Hamiltonians individually. The flux quantum Φ 0 will be displayed as Φ o in the toolbox, to distinguish it from the node flux of node 0. The numerical values for the corresponding spectrum plot in figure 1, which shows the lowest eigenstates of the weakly anharmonic cosinepotential, are the default values, i.e. C 01 = 100 fF and E J010 ≈ 9.69 GHz·h.
For the fluxonium qubit, the toolbox determines the symbolic Hamiltonian whereΦ 010 labels the offset flux of the inductive loop. An excerpt of the spectrum of this Hamiltonian is shown in figure 1 with C 01 = 10 fF, L 010 = 0.5 µH, E J010 ≈ 48.43 GHz·h andΦ 010 = πΦ 0 . It shows the typical low-lying 0 and 1 states that are localized in the wells, with higher plasma states several GHz above.
Finally, the persistent-current flux qubit Hamiltonian constructed by the toolbox can be written as We depict the ground state of this Hamiltonian in figure 1 for α = 0.7, C 01 = C 02 = C12 α = 50 fF, E J010 = E J020 = E J020 α ≈ 9.69 GHz·h andΦ 020 = πΦ 0 . The ground state is localized in the double well potential, which is repeated periodically throughout the chosen flux grid. As for the transmon circuit, due to the periodicity of the potential, the Hamiltonian is implemented in the charge basis. For figure 1, we use the transformation method of the toolbox to visualize the eigenstates in the flux basis.
In order to test the software and to check the accuracy of our numerical implementation, we perform benchmark tests for a variety of circuits. In addition to the three example circuits that have been discussed in this section, we complete the benchmark by adding the 0-π-Qubit [41], the 4-body-coupler which is referred to as Circuit C in reference [42] and the circuit of a transmon that is capacitively coupled to a resonator to the list of test circuits. The latter circuit is called the QuCAT circuit here, as a similar version is discussed in the corresponding reference [12]. The code to construct the CircuitQ instances can be found in Appendix D. We use existing software implementations to calculate the spectrum of the test circuits as a benchmark and compare the results of CircuitQ to it. As a reference, we used the toolbox sc-Qubits [11] for the transmon, fluxonium, 0-π-qubit and persistent-current flux qubit circuit. The QuCAT circuit has been compared to its implementation in the QuCAT toolbox [12]. The 4-body-coupler has been tested against a direct and individualized software implementation. As the test outcome depends on the size of the numerical matrices which represent the charge and flux variables, we vary the dimension of those matrices. The result is shown in figure 4. CircuitQ automatically implements the transmon and flux qubit in the charge basis, and the fluxonium, 0-π-qubit and 4-body-coupler in the flux basis. The QuCAT circuit is implemented in a mixture of both bases. For the transmon, fluxonium and flux qubit as well as for the QuCAT circuit, we find a good agreement between CircuitQ and the benchmark implementation. As detailed in appendix F, it is still possible to observe numerical limitations on less complex circuits like the fluxonium qubit. We observe larger deviations for the 4-body-coupler and the 0-π-qubit, which are more complex circuits, even for large numerical matrices. For some circuits, a quantization of the node variables is not the most natural choice, as characteristic modes of the system might be a combination of several node variables.
To find a more natural quantization, a coordinate transformation of the node variables can be performed prior to quantization. Therefore, the deviation for the 4-bodycoupler and the 0-π-qubit can be attributed to the lack of an appropriate coordinate transformation.

V. CONCLUSION
We presented the core functionalities of CircuitQ. With the ability to derive a symbolic and numerical Hamiltonian from a superconducting circuit in an automated way, CircuitQ can serve as a toolbox for the community to analyze superconducting circuits. The input circuit can be a general superconducting circuit that combines Josephson junctions, linear inductances and capacitances. An automated procedure to analyze superconducting circuits is a benefitial tool for the study of superconducting circuits within the context of quantum information. Apart from the application to computing, superconducting circuits can be also used as a platform in other areas of application like sensing [43] or studying thermodynamics [44].
While the toolbox is currently limited to the computation of few-node circuits, future work should address the optimization of speed and scalability. As an example, on a conventional personal computer, it took below 1 s to initialize an instance of the transmon circuit and calculate the lowest 10 eigenstates and eigenvalues of the numerical Hamiltonian for subsystem matrix dimensions 40 and 80, while for the 0-π-Qubit, it took around 72 s for a subsystem dimension 40 and around 105 s for a subsystem dimension 80, still with lacking accuracy as described in the section IV. A key feature of CircuitQ is its dynamic implementation in the charge and flux basis. At the moment, the variables that are quantized are always the node variables of the circuit graph. For some circuits, it is crucial to perform a variable transformation prior to quantization. Adding a suitable transformation would represent an important development step towards the goal of increasing the calculation speed. Another improvement can be made by implementing hierarchical diagonalization such as discussed in reference [10]. In addition, the toolbox is written in a modular fashion that allows for extensions towards time-dependent simulations.
Thanks to the general functionality, the possibility to include charge and flux offsets, as well as the the incorporation of noise estimates, CircuitQ can serve as a versatile tool for the design of superconducting qubits. Adding more noise channels, especially estimates for the dephasing time, would be an important future addition to the software. Moreover, adding the possibility to incorporate external impedances as circuit elements would allow to estimate noise from first principles [8].   In table II, we list the parameter values chosen for the purpose of this illustration. The numbers for the fluxonium correspond to Qubit A from reference [6], while the values for the flux qubit are associated to Qubit B from reference [32].
Appendix F: Limitations of the numerical treatment As discussed in section IV, the limitations of our numerical implementation become evident for complex circuits like the 0-π-Qubit. However, it is possible to investigate the numerical limitations on more simple circuits like the fluxonium qubit. In figure 5(a) and 5(b), we depict the lowest eigenstates together with the potential energy of the fluxonium qubit for low and high capacitive energy E C . The inductive energy E L has been kept small to avoid strong confinement. We observe the lowest eigenstates to be located within the potential wells for the case of low E C , while the eigenstates for higher E C tend to become more delocalized. In CircuitQ, the fluxonium circuit will be implemented in the flux basis, which works well for localized states that are trapped in a harmonic potential. To measure the numerical accuracy of the software implementation, figure 5(c) shows the deviation of the eigenenergies for the instances of subfigure 5(a) and 5(b) as a function of the numerical matrix dimension. In an ideal case, the spectrum is almost invariant under slight changes of the matrix dimension. Such a convergence can be observed for high values of matrix dimension. However, for lower values of matrix dimension, we observe a significant deviation of the energy, which is, besides the regime of very low matrix dimension, drastically higher for the case of high E C . This indicates that for the case of weakly localized wavefunctions, the flux basis implementation of the toolbox is reaching its numerical limitation. Moreover, another inaccuracy is introduced for higher lying states due to the cut-off of the numerical flux grid from −4π to 4π.
In figure 6, we present a similar study for the fixedfrequency transmon circuit (see figure 1(a)), which will be implemented in the charge basis. Figure 6(a) depicts the spectrum of the transmon circuit for low E C with E J /E C = 50. Here, the wavefunctions are periodically localized within the potential wells. In figure 6(b) we show the spectrum for higher E C with E J /E C ≈ 7, where higher eingestates become less confined. As before, we depict the difference in energy as a function of matrix dimension in figure 6(c). The energy difference drops of fast for both instances and fluctuates due to numerical fluctuations at small values for high matrix dimension. In comparison to the fluxonium study in figure 5(c), here, the depicted values of energy difference are small as the y-axis is scaled to smaller values. This indicates, that, in contrast to the flux basis, the charge basis is well suited to describe delocalized states. The states lying energetically above the potential can be associated with free particles, which, in the case of conventional mechanics, are efficiently described in the momentum basis. As the charge basis is analogous to the momentum basis, we find a more accurate description for the transmon circuit with high capacitive energy in comparison to the case of the fluxonium circuit. To calculate the energy difference for a given matrix dimension we sum over the differences between the lowest 5 eigenvalues and the eigenvalues calculated for the previous matrix dimension. The energy difference is in units of GHz·h and scaled logarithmically. Here, we solely use odd numbers as matrix dimension, as even numbers for the dimension will be converted to odd numbers by CircuitQ.