Computer-aided quantization and numerical analysis of superconducting circuits

The development of new superconducting circuits and the improvement of existing ones rely on the accurate modeling of spectral properties which are key to achieving the needed advances in qubit performance. Systematic circuit analysis at the lumped-element level, starting from a circuit network and culminating in a Hamiltonian appropriately describing the quantum properties of the circuit, is a well-established procedure, yet cumbersome to carry out manually for larger circuits. We present work utilizing symbolic computer algebra and numerical diagonalization routines versatile enough to tackle a variety of circuits. Results from this work are accessible through a newly released module of the scqubits package.


Introduction
The past two decades have seen dramatic progress in the design and control of quantum devices based on superconducting circuits [1][2][3]. A central part of this effort is the design of qubits which are robust against decoherence and can be controlled and measured with high fidelity [4][5][6]. The systematic discovery of new qubits requires exploring the wealth of possible circuits formed by combining inductors, capacitors, and Josephson junctions as the basic circuit elements. This has led to the design of a growing number of qubit types, each with its own advantages [4,[7][8][9][10][11][12]. The task of exploring new circuit designs, however, can be tedious and become numerically challenging since the Hilbert space dimension needed to describe the low-energy physics of the system tends to grow exponentially with the number of elements in the circuit. In order to make progress in the design of devices of increasing complexity, it is essential to promote tools that facilitate the efficient analysis of spectra and related key properties of superconducting circuits.
Several prescriptions have been proposed to approximate the low-energy Hamiltonian of a superconducting circuit, each depending on the level of detail needed for an accurate description [13][14][15][16][17][18]. For example, several approaches and software tools model the full geometry and distributed nature of the circuit [13,15,18]. This is computationally costly, and often works most readily for circuits with only weak anharmonicity. The most basic, yet quite widely applicable quantization framework is based on the lumped-element description of a circuit [17,[19][20][21][22][23]. Despite the fact that this approach is the simplest one, it nonetheless presents challenges once the number of circuit elements is increased beyond the bare minimum of two (found, for example, in the transmon qubit). This has led to the development of software packages that aim to automate the construction of the Hamiltonian [24][25][26][27].
In this work, we present a streamlined computer-aided approach to quantize lumped-element systems. By examining the circuit topology, we show how to classify and systematically detect four different types of degrees of freedom in a user-specified circuit: free, frozen, periodic and extended degrees of freedom. The identification of these variables can help reduce the number of degrees of freedom needed to describe the system and makes manifest the basis that should be used for an efficient representation of the Hamiltonian. This functionality has recently been incorporated into the scqubits Python package [28], leveraging the full array of capabilities already available. We complement our approach with a hierarchical-diagonalization algorithm. This further mitigates the circuit-size limitations governed generally by the size of the Hamiltonian matrix which describes the low-energy spectrum of interest.
For illustration, consider the task of finding the circuit Hamiltonian and numerical eigenenergies of a KITE circuit [29]. Using a simple input file KITE.yaml that encodes the circuit graph and associated circuit-element parameters (see Appendix A.1), generating the circuit Hamiltonian and performing numerical diagonalization is readily achieved by a few lines of code: The outline of the paper is as follows. In Section 2 we set the stage for the conceptual discussion of our approach by reviewing the lumped-element description of superconducting circuits. This includes a simple discussion of constructing the circuit Lagrangian in terms of node variables. In Section 2.2, we show how to systematically identify the four types of degrees of freedom in circuits consisting of capacitances, inductances, * and Josephson junctions. In Section 4, we describe how the identification of these variables is implemented in the scqubits package to efficiently represent the Hamiltonian of a user-provided circuit. Finally, Section 5 concludes with a summary and outlook of future developments.

Brief review of circuit Lagrangians
The common starting point for the quantization of a superconducting circuit is the construction of a circuit Lagrangian [20,21,[30][31][32][33]. In this paper, we exclusively focus on lumped-element circuits which are composed of capacitances, inductances, and Josephson tunnel junctions. (For more general discussions including additional elements such as gyrators, or extending the discussion beyond the lumped-element regime, consult Refs. [34][35][36][37].) Although seemingly restrictive, this small set of circuit elements spans a large family of different circuits, and includes the majority of superconducting qubits studied in the past, and used today. Each circuit element forms a branch in the connected circuit network.
Formulation of the Lagrangian depends on the definition of an appropriate set of generalized coordinates (which we also refer to as variables). In circuits including Josephson junctions, generalized-flux variables are the most convenient choice. For each branch b, one defines where t 0 is a reference time † and V b (t) is the voltage across branch b. The energy stored in a branch b is given by the for a capacitance C b , inductance L b , and junction with Josephson energy E Jb , respectively.
The collection of all branch variables Φ b (t) is in general not a set of independent variables: for each elementary loop L in the circuit, Kirchhoff's loop rule in conjunction with fluxoid quantization [21,38] imposes a constraint on the corresponding branch fluxes: Here, Φ ext is the external magnetic flux through the -th loop, and the branch sign (±) b is determined according to the directionality of the branch relative to the chosen loop orientation. The construction of a spanning tree provides a systematic way to determine a set of independent branch variables consistent with the constraints Eq. (2). A spanning tree T is a subset of connected branches such that every circuit node is connected to ground by one, and only one, path within the spanning tree [21] [see the example in Fig. 1(a)]. * Each branch not part of the spanning tree is a closure branch. Iterative inclusion of closure branches, i.e., starting from the spanning tree, adding the first closure branch, then the second, etc., leads to the closing of one circuit loop at a time. † Once the constraints have been resolved in this way, the circuit Lagrangian can be expressed in terms of the branch variables Φ b∈T belonging to the spanning tree, with Φ b / ∈T eliminated by means of Eq. (2), . (3) Here, notation follows a simple scheme assigning the capacitance C b to branch b = (n , n) if a capacitance indeed connects nodes n and n . Otherwise, the term is omitted from the summation. Analogous statements hold for inductance and junction contributions to the potential energy. We note in passing that time-dependent external flux renders the holonomic constraint (2) among the coordinates time-dependent, and hence also yields a constraint among the generalized velocitiesΦ b [37,40]; for simplicity we limit ourselves to the time-independent case here. Node variables present a useful alternative to branch variables. The node variable φ n at node n corresponds to the generalized flux obtained by adding up branch fluxes along the spanning tree; specifically, the spanning tree defines one path P n connecting the ground node to node n. Then, the node variable φ n is obtained by summing the branch * Spanning trees are not unique. This reflects the fact that the constraints (2) leave open the choice of which branch variable in a loop to eliminate. † Careful analysis of the circuit in Fig. 1 illustrates that there are cases in which the inclusion of the next closure branch appears to close multiple loops at a time. However, given the motivating purpose of resolving the loop constraints, it is always the newly closed loop containing exactly one closure branch which determines the relevant external loop flux. A more extensive discussion of this point has been presented in [39] variables along that path * : as illustrated in Fig. 1(a). Branch variables within the spanning tree as well as the generalized flux associated with each closure branch can then be expressed in terms of node variables via Here n, n are the nodes connected by the branch b = (n, n ), and Φ ext is the external flux threading the loop formed by the closure branch b and the spanning tree. The circuit Lagrangian (2) can then be re-expressed in terms of node variables in the general form where φ = (φ 1 , . . . , φ N ), and the N × N capacitance matrix C has non-zero off-diagonal entries (C) nn = −C b whenever a capacitance C b connects nodes n = n ; diagonal entries are given by (C) nn = − m =n (C) nm . The potential energy in node variables is obtained via Eq. (3) by substitution, where cl b is the indicator for closure branches, i.e. cl b = 1 for b / ∈ T and cl b = 0 otherwise. In principle, one might now proceed with either the branch-or node-variable Lagrangian, perform a Legendre transformation to obtain the Hamiltonian, and subsequently impose canonical quantization rules. However, it is generally worthwhile (especially where the ultimate aim involves numerical treatment) to consider a further transformation of variables which will achieve two important goals: (1) we may reduce the number of relevant variables, for example by identifying conserved quantities, and (2) we cleanly distinguish variables associated with different boundary conditions upon quantization. A suitable variable transformation can thus help make the quantum description of the circuit both more transparent and numerically efficient.

Classification of variables
Any variable θ, no matter whether branch, node, or collective variable of neither type, can be classified based on four categories: • Extended: a variable θ e is extended if it describes a direction in configuration space in which the potential energy is confining, expressed as lim θe→±∞ V (θ) = ∞. • Periodic: a variable θ p is periodic if it aligns with a direction in configuration space in which the potential energy is periodic, i.e., V (θ + Φ 0 e p ) = V (θ). Here, e p is the unit vector pointing along the axis associated with the θ p coordinate. We exclude the case of a variable in which the potential energy is constant; this situation is covered by the following category.
This means that θ c does not appear in the potential energy, and the corresponding conjugate charge Q c = ∂θ c L(θ, θ) is conserved (the subscript c refers to constant charge).
In this caseθ f is missing from the kinetic energy. This somewhat pathological case only arises when too many parasitic capacitances are neglected in the lumped element model. In the simplest setting, a frozen variable leads to an additional constraint given by ∂ θ f L(θ, θ) = 0. * Note: while the voltage across branch b is well-defined, one cannot generally fall back to a description in terms of electric-potential differences.
In the presence of non-zero external magnetic flux, even if constant in the present, the magnetic field had to be switched on after t 0 and be ramped up to its present value. As a result, the integral in equation (1) necessarily extends over a time interval during which the electric field was not curl-free, and could not be represented as a gradient field of the electric potential. † Another common, though potentially misleading, name for this variable type is "cyclic."

free periodic extended frozen
The diagram on the right shows the relation among these four categories. The free, periodic, and extended categories are mutually exclusive properties of the potential energy. The frozen category regards the kinetic energy and as such can occur for each of the former three categories.
This classification holds several advantages when analyzing the Hamiltonian of a circuit. Recognizing free and frozen variables helps reduce the number of degrees of freedom needed to describe the system. Furthermore, knowledge of periodic and extended variables guides a consistent and efficient choice of basis respecting the appropriate boundary conditions. Neither branch nor node variables are generally the suitable choice which recognizes the full number of periodic, free, and frozen degrees of freedom. Failure to find all the variables of each type makes diagonalization of the quantum Hamiltonian usually more difficult, and can potentially lead to wrong conclusions about offset-charge dependence and charge-noise sensitivity of circuit eigenstates.

Adaptive variable transformation
In the following, we show how to systematically construct a linear variable transformation which will pinpoint all periodic, free, and frozen degrees of freedom. Against naive instinct, it is helpful to define this transformation by expressing the "old" node variables in terms of the new variables, i.e., φ n = N m=1 z nm θ m , where Z = (z nm ) is an N × N invertible transformation matrix.
We determine the z nm by studying the subcircuits of the system, defined as connected sub-graphs of the circuit. We distinguish three relevant classes of subcircuits: (1) A subcircuit forms an isolated island if it is connected to the rest of the circuit exclusively by capacitances; (2) a subcircuit forms a superconducting island if it is coupled to the remainder of the circuit by capacitances and at least one junction (but not via any inductances); (3) a subcircuit is purely inductively coupled if it connects to the rest of the circuit only through inductances. Examples of each type of subcircuit are shown in Fig. 1(b).

Periodic and extended variables
For the analysis of periodic variables, consider a subcircuit S p forming a superconducting island as in the example of Fig. 2. There, S p is shaded in gray and couples to its circuit complement S p exclusively through capacitances and junctions. Any inductances present in the circuit are entirely inside S p or inside S p , but not at their interface. Given this partitioning, we see there is a simple choice of coefficients z n1 that will render the variable θ 1 periodic, namely: In other words, a change in the new flux variable θ 1 will simultaneously raise the flux of each node variable inside the superconducting island by equal amounts, but keep the flux on all nodes outside the superconducting island unchanged. Since θ 1 clearly produces no flux difference across any inductance, the potential energy must indeed be periodic with respect to θ 1 . The full set of periodic variables is obtained by identifying all superconducting islands in the circuit. This indicates that the full number of periodic variables, N p , is equal to the number of irreducible superconducting islands S p which cannot be decomposed into smaller superconducting islands. Once all periodic variables are assigned, any remaining variable θ e must now appear in the quadratic potential energy terms generated by the inductances. These terms have the usual confining property lim θe→±∞ V (θ) = ∞. As a result, all these degrees of freedom are extended.

Free variables
Free variables are treated in a manner nearly identical to the periodic case. Suppose θ c = θ 1 is one of the free variables. In this case, changes in θ c cannot produce node-flux differences across any of the inductances or junctions. Again, partitioning into subcircuits facilitates the definition of new variables. For free variables, the subcircuit S c must now be an isolated island connected to the remainder of the circuit S c via capacitances exclusively. Carrying over the definition of z n1 to the isolated island case, equips us with a free variable: since flux differences generated by θ 1 are across capacitances only, θ 1 will be absent from all potential energy terms associated with inductances and junctions. As a result, ∂ θc V (θ) = 0, so θ c is indeed free. The full number of free variables, N c , corresponds to the total number of irreducible isolated islands in the given circuit.

Frozen variables
As the last category, we discuss the case of frozen variables. For a frozen variable, the relevant subcircuit is of type (3): a subcircuit coupled to the remaining circuit through inductances only. Once such a subcircuit is identified, a frozen variable θ f is readily defined following the now familiar way: z-coefficients for nodes inside and outside the subcircuit are set to 1 and 0, respectively. As a result, the only flux differences that occur under a change in θ f are across inductance terms in the Lagrangian. Kinetic energy terms due to capacitances remain unchanged such that ∂θ f L = 0, and θ f is, hence, indeed a frozen variable. The full number of frozen variables corresponds to the total number of irreducible, purely inductively coupled subcircuits.
A few remarks about frozen variables are in order. Variables only appear frozen in situations where too many parasitic capacitances have been dropped from the description and the capacitance matrix has become singular. A foolproof, physically accurate but not necessarily efficient way of handling this situation is to re-admit the omitted parasitic capacitances which physically are always present. This is not always desirable as it generally increases the number of degrees of freedom retained in the modeling of the circuit.
Dropping a parasitic capacitance from the charging-energy contributions is akin to setting this capacitance to zero. In the mechanically analogous system, this corresponds to setting the mass of a particle to zero. The C → 0 limit is generally a singular limit [41] that exhibits subtleties particularly significant when junction capacitances are omitted [34]. We exclude this case here and focus on the more benign scenario in which one of the Euler-Lagrange equations turns into a linear holonomic constraint for C = 0, see Subsection 3.5.

Full transformation
As an example, consider the circuit in Fig. 1(b), which has three subcircuits. The corresponding transformation is This transformation identifies a free, a frozen and a periodic variable. The remaining variables are extended. The matrix elements corresponding to the extended variables {θ e,i } 4 i=1 are chosen such that the transformation is invertible. The Lagrangian in terms of the new variables is where C = Z C Z.

Elimination of non-dynamical variables
To obtain the circuit Hamiltonian, we first eliminate the frozen and free variables. Eliminating frozen variables is, in fact, a necessary step, since the capacitance matrix [Eq. (11)] is singular in the presence of frozen variables. To carry out the elimination of frozen variables from the Lagrangian, we use that ∂ θ f,i L = d dt ∂θ f,i L = 0, where i = 1, . . . , N f labels the inductively-coupled subcircuits. In general, this set of N f equations are nonlinear due to the cosine terms in the Lagrangian. As a crucial simplification, we assume that none of the junction capacitances are omitted in the description. In that case, only subcircuits with purely linear-inductive coupling to the remainder of the circuit lead to frozen variables. Thus, the variables {θ f,i } can be solved for using the constraints emerging from the Euler-Lagrange equations. Upon eliminating the {θ f,i } from the Lagrangian, the kinetic energy is expressed in terms of the invertible capacitance matrix C cpe obtained by removing the rows and columns in C that were associated with frozen variables. At this point, we perform a Legendre transform to obtain the classical Hamiltonian where Q cpe = ∂θ cpe L and θ cpe corresponds to free, periodic and extended variables. For every free variable θ c , the conjugate charge Q c is a conserved quantity. As shown in Appendix B, these conserved charges have no impact on the dynamics of the remaining degrees of freedom, and they may hence be omitted. The kinetic energy of the Hamiltonian can thus be written in the form T = 1 2 Q pe C −1 pe Q pe , where C −1 pe is the matrix block of C −1 cpe that corresponds to the extended and periodic variables only. Similarly, the vector Q pe only includes the non-conserved charges conjugate to the variables θ pe .

Quantization
Finally, canonical quantization yields the quantum mechanical Hamiltonian. Quantization is achieved by promoting the coordinate and conjugate momenta to operatorsθ n andQ n = −i ∂ θn satisfying canonical commutation relations. The Hamiltonian finally takes the form where we have introduced a vector Q g representing off-set charges. Since periodic and extended degrees of freedom are identified explicitly, the boundary conditions required for the wavefunction Ψ(θ pe ) can be stated transparently: The basis choice is further simplified by the distinction between periodic and extended variables. Periodic degrees of freedom can be represented in the discrete charge basis, extended variables in the phase basis or an appropriately chosen oscillator basis. The process of constructing the HamiltonianĤ of a given circuit and obtaining its spectrum has been implemented in the scqubits package [28], as we explain in the following.

Variable transformation
The variable transform discussed in Section 2.2 is implemented programmatically in order to automate the quantization of custom superconducting circuits. Several steps described in terms of subcircuits are equivalently formulated in more compact linear-algebra language applied to the vector space of node variables. The algorithm proceeds by dividing the vector space of node variables into independent subspaces associated with the different variable categories. For instance, the subspace corresponding to periodic variables is spanned by those vectors z p = (z 1p , z 2p , . . . , z N p ) which are compatible with the condition of vanishing flux difference across all inductances. Similar subspaces are constructed for free and frozen degrees of freedom present in the circuit. The full transformation matrix * is obtained by combining the vectors spanning the above three subspaces, and adding additional linearly independent vectors, if needed, to complete the basis. Fig. 3 summarizes the steps involved in this algorithm which has been implemented as a new module extending the scqubits package [28]. The Circuit class added for this purpose facilitates the generation and manipulation of symbolic circuit Hamiltonians by leveraging the Python package SymPy.

Numerical diagonalization
To obtain the spectrum of the circuit numerically, its Hamiltonian must be represented as a matrix with respect to a suitable basis. We choose bases {|µ i } (i = 1, 2 . . .) for all individual degrees of freedom. For the periodic degrees of freedom we employ the discrete charge eigenstates as a convenient (though not necessarily ideal) basis. For extended degrees of freedom, we either apply simple discretization with multi-point stencils turning derivatives into finite differences, or employ the harmonic oscillator basis associated with the potential energy obtained when omitting junction terms.
In the ordinary diagonalization scheme, the basis of the full Hilbert space is chosen as the Kronecker product |µ 1 ⊗ |µ 2 ⊗ . . . For each degree of freedom, individual truncation levels must be monitored to ensure convergence of the procedure. By contrast, hierarchical diagonalization aims to mitigate the fact that the dimension of the Hilbert space grows exponentially with the number of nodes † . The Circuit class offers an implementation of the hierarchical diagonalization procedure following ideas inspired by Kerman [42]. In the first stage, subsystems are diagonalized individually to obtain low-lying energy eigenstates {|E i }. The basis of the full Hilbert space used in this case is the Kronecker product |E 1 ⊗ |E 2 ⊗ . . . Under favorable conditions of moderate hybridization among subsystems, this procedure lowers the overall Hilbert space dimension. Once either of the two diagonalization schemes has been chosen, the full functionality of scqubits for visualizing spectra and wavefunctions, performing parameter sweeps, etc. is available to explore the behavior of the custom circuit.
As a concrete illustration of the hierarchical-diagonalization scheme, consider the fluxonium coupler circuit (see Appendix A.2), composed of five nodes and four degrees of freedom. Simulation of the full circuit Hamiltonian generally calls for diagonalization of matrices of size 10 6 × 10 6 . In suitable parameter regimes of practical interest, the system can be meaningfully divided into weakly coupled sub-systems, and hierarchical diagonalization proves to be a key step that decreases matrix sizes to 10 3 × 10 3 and dramatically reduces the required runtime by two orders of magnitude.

Conclusions
We described a computer-aided framework to quantize general lumped-element superconducting circuit systems. We showed how to systematically identify periodic, extended, and non-dynamical (free and frozen) degrees of freedom in a system. Eliminating the non-dynamical degrees of freedom leads to an efficient representation of the Hamiltonian. This framework is implemented programmatically in the new module circuit of scqubits Python library, which can simulate a user-defined circuit described as a graph. Hierarchical diagonalization was also implemented, which made numerical analysis for larger circuits accessible.
Future work will further extend the scope of the module by incorporating coherence-time calculations for custom circuits, handle the quantization and numerical analysis of circuits in the presence of time-dependent magnetic flux, and utilize the computation of perturbative corrections for hierarchical diagonalization discussed in Ref. [42]. Additional opportunities for enhancing the performance of our code lie in optimizing the choice of variable transformation which is generally not completely fixed by the protocol reported here. The remaining freedom may offer routes to increasing numerical efficiency, especially if variable transformations can reduce explicit coupling terms [17]. * Note that the transformation matrix generated can have any real valued coefficients. The only condition is that z ip -z jp ∈ Z such that i ∈ Sp and j / ∈ Sp, which ensures that the associated conjugate charge is periodic in a single Cooper pair. In some instances, having non-integer elements in the transformation matrix can help to increase the sparsity of the resultant matrix representation of the Hamiltonian, which can significantly reduce the diagonalization time [17]. † Details of convergence depend not only on the node number, but also critically on the specific parameter regime determined by the circuit parameters.

Note added
During preparation of this manuscript, we were made aware of ongoing work in Amir Safavi-Naeini's group with a similar scope. Our two groups did not discuss details of our projects, but agreed to post preprints on the arXiv simultaneously.
The parameters are (in units of GHz): E C = 2.5, E L1 = 0.23, E L2 = 0.36, E J = 5.9, E CJ = 6.6. The following snippet sets the truncation for the Hilbert space, chooses a more efficient transformation following Ref. [29] and also chooses the closure branches.

A.2. Hierarchical diagonalization for two fluxonia with a tunable coupler
The following string for the input file fluxonium-coupler.yaml encodes the circuit graph for the fluxonium coupler [43,44]: The circuit can be diagonalized using hierarchical diagonalization by the following lines of code: In: import scqubits as scq fluxonium_coupler scq.Circuit.from_yaml("fluxonium-coupler.yaml", ext_basis "harmonic", initiate_sym_calc False, basis_completion "canonical") We have verified that these eigenvalues agree with results obtained independently by the authors of [43].

Appendix B. Role of conserved charges in the dynamics of superconducting circuits
In this appendix, we show that conserved charges do not impact the physics of the dynamical variables. To simplify the analysis, we will assume that there is a single variable θ 1 that is free. We perform a partial Legendre transform with respect to the free variable to obtain the Routhian function R(θ pe , θ pe ; Q 1 ) = Q 1θ1 − L(θ, θ), (B.1) where Q 1 = ∂θ pe L(θ, θ) and θ pe represents the variables excluding θ 1 . The Routhian can be written in the form R(θ pe , θ pe ; Q 1 ) = −L 0 (θ pe , θ pe ) + 1 where L 0 (θ pe , θ pe ) is the Lagrangian of the system without the free degree of freedom, and 3) The dynamics of the variables θ pe is determined by the Euler-Lagrange equations d dt ∂θ pe R = ∂ θpe R. Since the conserved charge Q 1 enters through G(Q 1 , t) as a total time derivative in the Routhian, the value of Q 1 has no impact on the dynamics. As a result, after obtaining the classical Hamiltonian, the charge Q 1 can be set to zero. The generalization to more than one free variable can be carried out straightforwardly.