Variational tight-binding method for simulating large superconducting circuits

We generalize solid-state tight-binding techniques for the spectral analysis of large superconducting circuits. We find that tight-binding states can be better suited for approximating the low-energy excitations than charge-basis states, as illustrated for the interesting example of the current-mirror circuit. The use of tight binding can dramatically lower the Hilbert space dimension required for convergence to the true spectrum, and allows for the accurate simulation of larger circuits that are out of reach of charge basis diagonalization.


I. INTRODUCTION
Increasing coherence and noise resilience in superconducting qubits is a key requirement on the roadmap for developing the next generation of error-corrected quantum processors surpassing the NISQ era. Intrinsic noise protection in superconducting circuits has therefore become an important focus of research [1][2][3][4][5][6][7][8][9][10][11][12][13]. However, achieving simultaneous protection from depolarization and dephasing is impossible for small circuits like the transmon, and instead necessitates circuits with two or more degrees of freedom. Such larger circuits, especially of the size considered for the current-mirror circuit [3] or rhombi lattice [4,5], pose significant challenges for the quantitative analysis of energy spectra and prediction of coherence times. Consequently, the development of more efficient numerical tools capable of solving for eigenstates and eigenenergies of large superconducting circuits has emerged as a vital imperative. Strategies recently introduced for that purpose include hierarchical diagonalization [14], adaptive mode decoupling [15], and DMRG methods [16][17][18]. Here, we propose variational tight binding as another strategy complementing the former ones and illustrate its application.
Since the Hilbert space dimension d of even a single transmon circuit is infinite, it is not fully accurate to blame the "growth" of d for the challenges encountered with circuits of larger size. Nonetheless, when representing the Hamiltonian in a basis not specifically tailored for the problem at hand, the dimension of the truncated Hilbert space typically grows exponentially when choosing the truncation level such that a particular level of convergence is reached. This turns the numerical diagonalization of the circuit Hamiltonian into a hard problem. An approach to address this challenge, which is also implicitly represented by the strategies mentioned above, consists of constructing basis states which more closely approximate the desired low-energy eigenstates from the very beginning. As long as construction of the tailored basis and decomposition of the Hamiltonian in that basis can be accomplished efficiently, this approach will allow for reduced truncation levels and hence enable coverage of circuit sizes otherwise inaccessible numerically.
Our construction of such tailored basis states is based on the observation that low-lying eigenstates of superconducting circuits are often localized in the vicinity of minima of the potential energy, when expressed in terms of appropriate generalized-flux variables. If the potential energy is periodic or at least periodic along certain axes, then the situation resembles the setting of a particle in a periodic potential, as commonly encountered in solid-state physics when considering electrons inside a crystal lattice. In the regime where tunneling between atomic orbitals of different atoms is weak, tight-binding methods are appropriate for band structure calculations [19,20]. An analogous treatment has previously been applied to small circuits; see, for example, the discussions of tunneling between minima in the flux qubit [21][22][23], the derivation of an asymptotic expression for the charge dispersion in the transmon qubit [24], or the analysis of charge noise in the fluxonium circuit [25]. Chirolli and Burkard carry out a full tight-binding description of the low-energy physics of the flux qubit, considering Bloch sums of harmonic oscillator ground-state wavefunctions localized in each minimum at the half-flux sweet spot [22]. Motivated by the new interest in circuits of increased size and complexity, we build upon this research in two specific ways. First, we consider multiple basis states in each minimum, to both improve ground-state energy estimates and extract excited-state energies. Second, we consider minima that are not necessarily identical, and introduce an efficient means of calculating matrix elements between states localized in such inequivalent minima. These techniques allow us to demonstrate that tight-binding methods can be adapted for efficient computation of energy spectra of large circuits.
Our paper is organized as follows. In Sec. II we review the tight-binding treatment in solid-state physics and develop its adaptation for superconducting circuits. We then detail our numerical implementation, involving the calculation of tight-binding matrix elements using ladder operator algebra. In Sec. III we first apply our numerical method to the simple example of a flux qubit and illustrate that tight binding can accurately reproduce the well-known results for the spectrum. In Sec. IV we then apply tight binding to the current-mirror circuit with up to nine degrees of freedom, and demonstrate that tight binding outperforms the diagonalization of the Hamiltonian in charge-basis representation in terms of accuracy, convergence behavior, and memory efficiency. We close with our conclusions in Sec. V.

II. TIGHT BINDING FOR SUPERCONDUCTING CIRCUITS
Research on superconducting qubits has repeatedly encountered physics familiar from models and phenomena in solidstate physics. Examples include the close connection between the Cooper pair box and a particle in a one-dimensional crystal, or the interpretation of the fluxonium Hamiltonian in terms of Bloch states subject to interband coupling [26]. Another analogy, which points to the computational technique applied to circuits in this article, is the consideration of crystal electrons in the tight-binding limit. In this regime, tunneling between electronic orbitals of different atoms is weak, and linear combinations of atomic orbitals constructed by "periodically repeating" localized wavefunctions serve as a meaningful basis. The tight-binding method then employs this basis in an approximate solution to the Schrödinger equation. An analogous scenario can be encountered for superconducting circuits as shown in Fig. 3. Minima of the potential energy may give rise to localized states that are only weakly connected by tunneling to partner states in other potential minima. The "atomic orbitals" which we will refer to as "local wavefunctions" in this case can be identified with the harmonic oscillator states associated with a local Taylor expansion around each minimum.

A. Local-wavefunction construction
The starting point for this treatment is the full circuit Hamiltonian H = T + V . To stress the analogy with the setting of an infinite crystal, we first focus on a purely periodic potential V ( φ), as realized by a circuit that does not include any inductors. (Including inductors is possible, which we comment on further in Sec. II B). In terms of the node variables φ = (φ 1 , . . . , φ N ) t , the potential energy obeys the periodicity condition V ( φ + 2π j) = V ( φ) with j ∈ Z N and thus forms a (hyper-)cubic Bravais lattice. Within the central unit cell defined by φ ∈ [−π, π) ×N , the potential energy will exhibit a set of M minima located at positions θ m where m = 0, 1, . . . , M . In the language of solid-state physics, this set of minima corresponds to the multi-atomic basis associated with the Bravais lattice.
The analogy with solid-state physics is further strengthened by considering a gauge where the offset charge dependence is shifted from the Hamiltonian to the wavefunctions [22,27]. In this representation, solutions |ψ to the full Hamiltonian H obey quasiperiodic boundary conditions for every θ in the Bravais lattice, where T is the translation operator and n g = (n g1 , · · · , n gN ) t is the vector of offset charges. We recognize Eq. (1) as an expression of Bloch's theorem with wavevector − n g (typically denoted as k in a solid-state context). The construction of the local wavefunctions now proceeds by considering the individual harmonic oscillator Hamiltoni-ans H m = T + V m where each local potential is obtained by Taylor expansion around the respective m th minimum, Here, ϕ 0 = /2e is the reduced flux quantum, Γ (m) ij = ϕ −2 0 ∂ φi ∂ φj V | θm the inverse of the inductance matrix and φ (m) = φ − θ m the "position" relative to the minimum location. The local Hamiltonian then takes the form where n i is the charge number operator for node i and C is the capacitance matrix. Hereafter, explicit references to m will be omitted for notational simplicity. To obtain the eigenstates of the coupled oscillator Hamiltonian Eq. (3) we first determine its normal modes. This is accomplished most efficiently based on the corresponding classical Lagrangian Using the usual oscillatory solution ansatz φ = ξ µ e −iωµt reduces the equations of motion to the generalized eigenvalue problem Γ ξ µ = ω 2 µ C ξ µ [28]. Here, Latin indices refer to node variables, and Greek indices to normal mode variables. The eigenmode vectors ξ µ are only determined up to normal- where c µ is undetermined. This normalization will be fixed when we return to the quantum-mechanical description, in such a way that the Hamiltonian for each mode takes the standard form Here, ζ = (ζ µ ) collects the normal-mode variables related to the original generalized fluxes via φ = Ξ ζ, where Ξ is the matrix of column vectors ξ µ . 1 In these new variables, both bilinear forms in L are diagonal. Legendre transform and quantization thus readily yield To cast H into the form suggested by Eq. (5) we now choose c µ = (2e) 2 / ω µ as our normalization constants. We denote the eigenstates of H by | s, m . Here, s collects the excitation numbers s µ = 0, 1, . . . for each mode µ and m specifies the minimum of interest.
FIG. 1. Comparison between tight binding as applied to solids and superconducting circuits. On the left is an example two-dimensional lattice with a two-atom basis, signified by the grey and blue atoms. On the far left are example atomic wavefunctions. On the right is the potential of the flux qubit, which has two inequivalent minima in each unit cell at the chosen value of flux. We only color the potential below a cutoff value to draw the eye to the potential minima locations. Near the minima the potential is approximately harmonic, therefore the local wavefunctions take the form of harmonic oscillator states. Example wavefunctions are shown on the right. Here, N is the number of unit cells and | s, m; j the wavefunction localized in minimum m in the unit cell located at 2π j. (Note that these kets | s, m; j are implictly offset-charge dependent). It is straightforward to show that |ψ ng, s,m satisfies the quasiperiodicity condition (1). We now represent the Schrödinger equation in terms of these basis states. Due to their lack of orthogonality, this transforms the Schrödinger equation into the generalized eigenvalue problem Eq. (8) can be simplified by performing one of the sums over lattice vectors [19]. This can be done by expressing the kets explicitly in terms of the translation operators | s, m; j = e −i ng·(2π j+ θm) T 2π j | s, m; 0 and noting that the operator T 2π j commutes with the Hamiltonian. The summation yields a factor of N and we obtain Formally, Eq. (10) now has the standard form of a generalized eigenvalue problem with two semidefinite positive Hermitean matrices and can be handled numerically by an appropriate solver. To accomplish this, the crucial remaining task consists of the efficient evaluation of the matrix elements and state overlaps in Eq. (10). Note that an alternative route to this equation is application of the variational principle to ψ ng |H|ψ ng = E ψ ng |ψ ng [29]; the benefit of this viewpoint is that the eigenenergies thus obtained represent upper bounds to the true eigenenergies of the system [29,30]. Our analysis thus far has assumed a purely periodic potential, allowing for a direct analogy with the theory of tight binding as applied to solids. Including inductive terms in the potential immediately implies that associated degrees of freedom are no longer subject to (quasi-)periodic boundary conditions. Alternatively, we can say that the unit cell no longer has finite volume, but must extend along the relevant axes. To include such inductive potential terms, we therefore do not perform periodic summation in Eq. (7) along these non-periodic directions. We have successfully implemented the tight-binding method for the symmetric 0 − π qubit [6,8], a circuit with one periodic and one extended degree of freedom. The low-energy spectra thus obtained are in excellent agreement with exact results over a wide range of circuit parameters. In this paper we will continue to focus on circuits with purely periodic potentials, the natural setting for the tight-binding method. We defer a detailed discussion of our results for the 0 − π qubit to a future publication.

C. Efficient computation of matrix elements and overlaps
The relevant matrix elements involve harmonic-oscillator states at different locations and, possibly, with different normal-mode orientations and oscillator lengths. The calculation of these quantities proceeds either via use of ladder operators or by explicit integration within the position representation. Even though integration can in principle be accomplished analytically, the expressions become increasingly tedious in higher dimensions. (The integrals are generally twocenter integrals that lead to two-variable Hermite polynomials [31]). By contrast, the ladder-operator formalism is more readily adapted for the numerical calculations of the matrix elements in question. Therefore, we focus on this approach.
The matrix elements and overlaps to be evaluated have the form where O is either the Hamiltonian H or the identity. To facilitate the use of the ladder-operator formalism, we next reexpress operators and states in terms of the creation and annihilation operators associated with the m = 0 minimum in the central unit cell. Since inequivalent minima differ in locations and curvatures, local wavefunctions are shifted and possibly squeezed relative to each other, where | s ≡ | s, 0 and we have taken the location of the m = 0 minimum to be the origin, θ m=0 = 0. The intuitive interpretation of Eq. (12) is based on a two-step process: first the harmonic oscillator states for m = 0 are deformed to match the local curvature of the m th minimum and are then shifted over to the appropriate location of that minimum. According to Eq. (12), the matrix elements take the form The expression for the states is readily obtained, Likewise the treatment of the translation operator T θ = e −i θ· n is straightforward: making use of the relation between the number operators and ladder operators the translation operator can be expressed as Here we use the compact notation a = (a 1 , · · · , a N ) t and a † = (a 1 † , · · · , a N † ) t and denote The expression for the squeezing operator can be found by considering a simplified situation of two harmonic Hamiltonians H a , H c of the form of Eq. (3), but defined at the same center point. The Hamiltonian H a is diagonalized by the ladder operators a, a † and H c is diagonalized by c, c † . The respective eigenfunctions | s a and | s c are related by a unitary squeezing transformation, which is equivalent to S aS † = c. To obtain a concrete expression for S we note that the bosonic ladder operators are related by a Bogoliubov transformation [32,33] u where u, v are N × N matrices. These can be found by considering the two decompositions of the phase and number operators in terms of the differing sets of ladder operators where the matrix Ξ is defined for H a and Ξ for H c as in Sec. II A. Solving Eq. (19) for the ladder operators c, c † yields the real-valued Bogoliubov matrices As shown in Ref. [34], the multimode squeezing operator can now be expressed in terms of u, v as follows: where Returning now to our original notation, we identify S m with S in Eq. (21), where the m dependence carries forward to the Bogoliubov u, v matrices. With Eqs. (14), (16), (21) we have, in principle, collected all ingredients necessary for the evaluation of the matrix elements and overlaps [Eq. (13)]. However, numerical implementation necessarily involves truncation, and we will show in the following that normal-ordering operator expressions is essential for maximizing accuracy. A standard approach for truncating the infinite-dimensional operators a µ , a † µ consists of excitation cutoffs s max applied to each mode. This is a fine strategy for moderate sized systems, but quickly becomes intractable for larger systems. To mitigate this bottleneck, one can instead use a global excitation number cutoff Σ max , which institutes a maximum Manhattan length of the excitation number vector, s 1 ≤ Σ max [35].
Given a specific truncation level, it makes a difference whether operator expressions are normal ordered or not. Denoting the truncated operators asã µ ,ã † µ , the nominally identical expressionsã µã † µ andã † µãµ + δ µµ in fact give different results as seen, for instance, in Here, the "wrong" result of the first expression can be circumvented by using the normal-ordered version in the second expression. This example is indicative of a general result, that it is beneficial to normal order ladder-operator expressions before further numerical evaluation. The translation operator T θ can be normal ordered via the the Baker-Campbell-Hausdorff (BCH) formula [36], which takes the form e X e Y = e X+Y + 1 2 [X,Y ] when X and Y are operators that commute with their commutator [X, Y ]. This yields where Expressions for commuting V operators past operators such as n j and e iφj (which enter O in Eq. (11)) can be easily obtained: where in Eq. (26) we have used the identity [37] e X Y = (Y + [X, Y ])e X , again valid when X and Y commute with [X, Y ].
The normal-ordering procedure for the squeezing operators is more involved and we defer it to Appendix A. In the following sections where we apply the tight-binding method to several example systems, we find that the sets of basis states constructed with (proper tight binding) or without (improper tight binding) squeezing may yield similar numerical performance. This is naturally the case if minima contributing to the low-energy spectrum have similar curvatures. Whenever possible, omitting squeezing from the construction of basis states significantly simplifies the numerical treatment.
The final step in setting up the generalized eigenvalue problem (10), is truncating the sum over vectors j. A typical truncation scheme is the nearest-neighbor approximation which selects only those unit cells that have the minimal Euclidean distance from the central unit cell. This strategy however does not account for any anisotropy in the harmonic lengths, which results in local wavefunctions whose Gaussian tails extend further in some directions than in others. We therefore use a different criterion based on the overlap of local wavefunctions. Whether the unit cell centered at 2π j is a nearest neighbor to the central unit cell now generally depends on the minima under consideration. Specifically, given a minimum m in the central unit cell, and a minimum m in the unit cell at vector 2π j, we determine the nearest-neighbor character by computing the overlap of the two harmonic oscillator groundstate wavefunctions. For a given overlap threshold value , we call the two unit cells nearest neighbors with respect to m and m if Here, we have defined δθ = 2π j + θ m − θ m and ∆ m = Ξ −t m Ξ −1 m , where Ξ m is defined relative to minimum m. With this definition in place, we truncate the sum over j by selecting neighbors up to a certain degree. (Note that the overlap threshold must ultimately be adjusted adaptively in order to ensure convergence.) A possible challenge for the numerical treatment, which we have observed in several cases, is that the overlap matrix ψ ng, s ,m |ψ ng, s,m may approach singularity (and possibly become indefinite due to rounding errors). This is a familiar problem in quantum chemistry calculations [38,39] and arises when the set of "basis" states {|ψ 1 , |ψ 2 , . . . , |ψ h } is approximately linearly dependent. A common technique for resolving this issue which we have implemented here is the canonical orthogonalization procedure of Löwdin [40]. One diagonalizes the inner product matrix to obtain the eigenvalues {∆ 1 , ∆ 2 , . . . , ∆ h } and matrix of column eigenvectors U . The orthonormalized states are [40] Choosing a cutoff ∆ min allows for the rejection of states |ψ k where ∆ k < ∆ min . The Hamiltonian H is then projected onto the deflated basis and we are left with a standard eigenvalue problem.

D. Optimization and anharmonicity correction of the ansatz wavefunctions
One of the main goals of this work is the construction of basis states that closely approximate the low-energy eigenstates of superconducting circuits. We can optimize the tightbinding wavefunctions (7) for this purpose by recognizing that sufficiently far from each minimum location, the potential ceases to be strictly harmonic. The low-energy eigenfunctions typically have spatial spreads that are broader if the leadingorder anharmonic term is negative and narrower if it is positive. We take this effect into account and improve the tightbinding wavefunctions by treating the harmonic length of each mode as a variational parameter. Specifically, we modify the matrix Ξ by optimizing the magnitude of the eigenvectors ξ µ , leaving the directions unchanged, ξ µ → λ µ ξ µ , where λ µ is optimized. We perform this optimization procedure for the The resulting harmonic lengths are then used for all other states defined in the same minimum m. 2 We term this optimization scheme "anharmonicity correction," which combined with improper and proper tight binding leads to additional choices for constructing tight-binding states: (IPAC) improper with anharmonic correction and (PAC) proper with anharmonic correction of the m = 0 minimum, see Fig. 2. We could further envision the construction of states according to: proper with anharmonicity correction of all minima. However, we find that this scheme yields no numerical benefit over the other methods in all cases considered here, thus we do not discuss it further.

E. Applicability of tight binding
A natural question to ask is whether the tight-binding method is appropriate for obtaining the eigenspectrum of a given superconducting circuit. A general and systematic answer to this question is difficult to obtain and we do not aim to give a comprehensive answer here. Instead we seek to motivate a "rule-of-thumb" criteria that serves as an indicator of whether the tight-binding method can produce meaningful results.
If the spatial spread of the localized harmonic-oscillator states [eigenfunctions of Eq. (6)] is small compared to distances between minima then the tight-binding approach is physically well motivated and we expect tight-binding wavefunctions to serve as good approximations to low-energy eigenstates. If, on the other hand, the wavefunctions have large spatial spread and significant overlap, then the weakperiodic-potential approximation is more appropriate for describing the low-energy excitations.
To quantify this discussion, we define length scales to compare the spatial spread of wavefunctions with the distance between minima. Examining the exponential dependence exp(− 1 2 φ t Ξ −t Ξ −1 φ) of the local harmonic wavefunctions, we can extract the effective harmonic length mm along the unit vectorû mm separating two minima m and m It is natural to compare mm to d mm /2, half the distance between the minima. Our rule-of-thumb for application of the tight-binding method is based on the smallness of the localization ratios compared to unity. This provides a rough threshold for judging whether the tight-binding method might be appropriate.

III. TIGHT BINDING APPLIED TO THE FLUX QUBIT
In order to evaluate the accuracy of the tight-binding method, we first apply it to the familiar case of the threejunction flux qubit. The spectrum of the flux qubit is well understood [21,41], but applying the method in this context is of interest and nontrivial because the flux qubit has multiple degrees of freedom and multiple inequivalent minima in the central unit cell. Additionally, the flux qubit is typically operated in a parameter regime where tight-binding techniques are applicable. Indeed, many authors have used tight-binding techniques to get analytical estimates of tunneling rates and lowenergy eigenvalues [21][22][23]. We extend this previous research by using multiple tight-binding basis states in each inequivalent minimum to obtain improved low-energy eigenvalue estimates.
We consider the case where two of the junctions are identical with junction energy E J and capacitance C J , while the third has junction energy and capacitance reduced by a factor of α. The Hamiltonian is [21] where ϕ ext = 2πΦ ext /Φ 0 , Φ 0 = h/2e is the flux quantum, E C = e 2 2 C −1 is the charging energy matrix and the constant term is included to ensure that the spectrum of H flux is positive. The capacitance matrix C is where C g is the capacitance to ground of each island. (See Refs. [21,41] for details on the derivation of Eq. (33)).
In order to demonstrate quantitative accuracy of the tightbinding method, we calculate the flux and offset-charge dependence of the spectrum, see Fig. 3. For the parameters considered, the localization ratios are large compared with unity, indicating that the parameter regimes are amenable to tight binding. Figs. 3(a-b) show the spectrum as a function of flux and offset charge, respectively, with improper-tightbinding results overlaying the exact spectrum obtained via charge-basis diagonalization. While the spectra from the two different methods are indistinguishable in the upper panels of Figs. 3(a-b), we explicitly visualize the residuals for the four lowest eigenergies in the lower panels of Figs. 3(a-b). For Σ max = 5, the residuals are all below 1 MHz for flux and offset-charge variation. Further suppression of the absolute error below 1 kHz is possible by increasing the global excitation number cutoff to Σ max = 10.
Even for relatively greedy cutoffs of the global excitation number Σ max , the improper-tight-binding method can provide accurate estimates of the eigenspectrum. To compare results obtained using tight-binding methods with results from exact diagonalization, we compute the relative deviation from the exact low-energy spectrum, averaged over the four lowestenergy eigenvalues Here, i is the exact eigenenergy of the state indexed by i and E i is the approximate eigenenergy. We also define the minimum and maximum relative deviations To monitor convergence and assess the memory requirements for reaching a desired accuracy, we plot in Fig. 4 η avg as a function of nonzero Hamiltonian matrix elements (n H ). We use n H rather than Hilbert-space dimension as a proxy for memory usage to account for the different cases of sparse vs. dense matrix numerics encountered for diagonalization in the charge basis vs. tight binding. For a cutoff as greedy as Σ max = 1, we find η avg < 7 · 10 −3 using improper tight binding. Note that for the flux qubit in the parameter regimes considered here, neither the proper-tight-binding technique nor anharmonicity correction provided any appreciable benefit in terms of convergence to the spectrum over improper tight binding.
To benchmark convergence of the tight-binding method, we compare against results obtained using truncated diagonalization in the charge basis. We compute the average relative deviation η avg using energy estimates obtained via a choice of charge-basis cutoff n cut . By increasing n cut , we increase n H and can thereby perform a direct comparison with η avg values obtained via tight binding, see Fig. 4. The shaded region of Fig. 4 indicates where tight binding outperforms approximate diagonalization in the charge basis for a given n H . The advantage region for tight binding is for small values of n H , indicating that when keeping few basis states, tight-binding states yield a closer approximation to the true low-energy eigenstates than charge-basis states. At larger values of n H , chargebasis diagonalization begins to outperform tight binding.

IV. TIGHT BINDING APPLIED TO THE CURRENT-MIRROR CIRCUIT
We expect the tight-binding method to be most useful in the study of larger circuits, where keeping a generous number of basis states is not feasible due to memory requirements. To demonstrate the tight-binding method on such a larger circuit, we apply it to the current mirror circuit [3], described by the Hamiltonian [17,18] where N B refers to the number of big capacitors. The charging energy matrix E C involves contributions from individual charging energies E C B , E C J , E Cg due to the big-shunt, junction and ground capacitances respectively [17]. An example circuit with N B = 3 without the capacitors to ground is shown in the inset of Fig. 5. The number of degrees of freedom of the circuit is given by 2N B − 1. The interest in this circuit originates from Kitaev's prediction that quantum information should be protected relaxation and dephasing in the current mirror [3]. For a representative choice of parameters, one can identify N B ≈ 12 as the ideal value of N B [17]. Circuit sizes with such large values of N B exceed our capabilities for finding eigenstates and eigenenergies via diagonalization in the charge basis; the maximum value of N B where we can achieve spectral convergence is N B = 3. 3 We show below that the The colored shaded regions encompass the range between ηmin and ηmax. The gray shaded region represents the nH values for which tight binding yields an advantage over charge-basis diagonalization, comparing ηavg for a given nH . Improper tight binding allows for an accurate estimate of the low-energy eigenspectrum already at Σmax = 1, yielding ηavg < 7 · 10 −3 and maximum absolute error of less than 25 MHz. We choose the parameters of Fig. 3(a), as well as ϕext = 0.47, ng1 = 0.2, ng2 = 0.3. We can perform the same calculation for the parameters of Fig. 3(b) and obtain similar results, with tight binding outperforming charge-basis diagonalization for small nH . The inset shows a schematic of the flux-qubit circuit.
tight-binding method is an advantageous alternative for simulating the current-mirror circuit at larger values of N B . Implementation of the tight-binding method for the currentmirror circuit proceeds in a manner analogous to the case of the flux qubit, as neither circuit contains inductors. We choose a set of protected circuit parameters given by E C B /h = 0.2 GHz, E C J /h = 35 GHz, E Cg /h = 45 GHz, E J /h = 10 GHz, ϕ ext = 0 and n gi = 0. To establish qualitatively that the current-mirror circuit with these parameters is amenable to a tight-binding treatment, we fix a value of N B and verify that the localization ratios are all of order unity or larger. We observe that the localization ratios generally increase with N B , indicating that the tight-binding method should become increasingly accurate with larger N B . For an independent quantitative assessment of the validity of tight binding, we will compare spectra obtained with tight-binding methods with exact results. For this purpose we first apply the tight-binding method to the N B = 3 current-mirror circuit.
We can obtain excellent agreement between spectra obtained via tight binding and exact results, with average relative deviations η avg below 2 · 10 −5 , see Fig. 5. The best agreement is for the energy of the ground state, for which we obtain agreement to within 16 kHz. For the first-and second-excited states these results correspond to sub-MHz agreement, while FIG. 5. Performance of the tight-binding method as applied to the NB = 3 current-mirror circuit. Similarly to the case of the flux qubit, we plot ηavg for improper tight binding (blue, solid), improper tight binding with anharmonicity correction (red, solid) and approximate diagonalization in the charge basis (green, dashed) as a function of nH . Improper tight binding with anharmonicity correction outperforms charge-basis diagonalization across approximately four orders of magnitude in nH , as indicated by the shaded region. The sharp cliff in ηmin for improper tight binding with anharmonicity correction at nH ≈ 10 5 is due to the inclusion of new basis states that contribute to the ground state, yielding ηmin ≈ 10 −4 . The inset shows a schematic of the NB = 3 current mirror circuit. We choose ϕext = 0 and ngi = 0, with circuit parameters given in the main text.
for the third-excited state agreement is on the order of a MHz. The use of the anharmonicity correction yields a substantial benefit that is critical for achieving this level of accuracy. We find that the proper-tight-binding method yields nearly identical results to those produced by improper tight binding, and therefore those results are not shown in Fig. 5. Our highest accuracy approximations are obtained with Σ max = 8, beyond which we encounter numerical instabilities. We emphasize that one can actually obtain a reasonable approximation to the spectrum based on moderate values of Σ max , as shown in Fig. 5. For example with Σ max = 1, 2, improper tight binding with anharmonicity correction yields η avg ≈ 8 · 10 −3 , corresponding to absolute errors of about 300 MHz.
We can contrast these results with those obtained using truncated diagonalization in the charge basis. Using the same metric for memory efficiency previously applied to the fluxqubit example, we find that tight binding is advantageous over a wide range of n H values, see Fig. 5. Specifically, to achieve η avg ≈ 8 · 10 −3 , truncated diagonalization in the charge basis requires about three more orders of magnitude in memory resources as compared to tight binding. The advantage region for tight binding extends over approximately four orders of magnitude 10 2 n H 10 6 , as shown in the shaded area of Fig. 5.
To extend toward the regime of ideal N B , we apply the tight-binding method to obtain the spectrum of the N B = 5 current-mirror circuit, which has 9 degrees of freedom. We compute the ground-state energy E 0 and first-excited-state energy E 1 using the four tight-binding techniques [improper, proper, (IPAC) improper with anharmonicity correction and (PAC) proper with anharmonicity correction of the m = 0 minimum], see Tab. I. By the variational principle, our computed eigenenergies are upper bounds to the true eigenenergies [29,30,40]. Therefore, lower eigenenergy values always imply higher accuracy. The proper, IPAC and PAC tight-binding schemes all perform similarly but collectively outperform the improper scheme. The lowest eigenenergies are obtained using IPAC, with bounds 0 ≤ 81.6472 GHz and 1 ≤ 82.7224 GHz. The largest cutoff we can handle is Σ max = 5, beyond which we encounter memory issues. We observe that for this circuit, ansatz states localized in minima aside from the m = 0 minimum contribute to the low-energy spectrum, and moreover the curvatures of those minima differ from those of the m = 0 minimum. Otherwise, there would be no difference between the eigenenergies computed with improper and proper tight binding. Note that schemes IPAC, PAC allow for rough estimates of the eigenspectrum with a greedy cutoff Σ max = 2, with calculated E 0 , E 1 less than 200 MHz greater than the lowest obtained respective values.
We next compare tight-binding results with those from approximate diagonalization using the truncated charge basis. For N B = 5 the maximum possible charge cutoff we can handle is n cut = 3, corresponding to a Hilbert-space dimension of (2·n cut +1) 9 = 7 9 ≈ 4·10 7 . The best estimates for E 0 and E 1 obtained using approximate diagonalization in the charge basis are in fact higher and therefore less accurate than the lowest obtained values using tight-binding methods, see Tab. I. Moreover, the tight-binding methods consistently yield lower eigenenergy approximations across all n H values. Fig. 6 illustrates this point for the ground-state energy E 0 , and similar results hold for the first-excited-state energy E 1 . We thus find that the tight-binding method is more memory efficient than charge-basis diagonalization for the N B = 5 current-mirror circuit. More broadly, this may indicate that the tight-binding method can serve as an interesting and useful method in the context of large circuits.

V. CONCLUSION
We have generalized the well-known method of tight binding for the purpose of efficiently and accurately obtaining the low-energy spectra of superconducting circuits. Construction of the Hamiltonian proceeds by using ansatz Bloch states that localize in minima of the potential. The method can handle many degrees of freedom, multiple inequivalent minima, and periodic or extended potentials. In terms of these states, the Schrödinger equation turns into a generalized eigenvalue problem. Solving it yields a spectrum that provides upper bounds to the true eigenenergies. To establish the accuracy of the tight-binding method we apply it to the flux qubit and achieve agreement with exact results at the kHz level.
Because the method is expected to be of use for larger circuits, we apply it to the N B = 3, 5 current-mirror circuits, which have 5 and 9 degrees of freedom, respectively. We TABLE I. Eigenenergies for the ground state and first-excited state of the NB = 5 current-mirror circuit. Energies were computed using tightbinding schemes (IP) improper, (P) proper, (IPAC) improper with anharmonicity correction, (PAC) proper with anharmonicity correction, as well as (AD) approximate diagonalization in the charge basis. The energies are color coded from least accurate (darkest) to most accurate (lightest). The three tight-binding flavors proper, IPAC and PAC all perform similarly and outperform improper. The most accurate results for E0 and E1 were obtained with tight binding rather than with approximate diagonalization in the charge basis (circuit parameters used are the same as in Fig. 5).

MHz
FIG. 6. Comparison of computed ground-state energies for the NB = 5 current-mirror circuit. Individual curves correspond to results obtained with tight-binding schemes improper (blue, solid), improper with anharmonicity correction (red, solid), proper (orange, solid) and proper with anharmonicity correction of the global minimum (pink, solid), as well as approximate diagonalization in the charge basis (green, dashed). Tight-binding techniques consistently yield lower and hence more accurate eigenenergies as compared to charge basis diagonalization, with a difference of 163 MHz between the best results obtained with tight binding and approximate diagonalization in the charge basis, see the inset. We used the same current-mirror circuit parameters here as in Fig. 5. find excellent agreement with exact results in the case of the N B = 3 circuit. Moreover, across multiple orders of magnitude in memory usage (as quantified by n H ), eigenenergies computed using tight binding are found to be more accurate than those calculated using the charge basis. For the N B = 5 circuit, the tight-binding method allows for the extraction of eigenenergies that are lower than any obtainable using the charge basis, given our computational resources. This work supplements recent research also focused on the efficient simulation of large superconducting circuits [14,15].
To extend and improve the tight-binding method beyond what is described here, we envision developing an improved state-optimization procedure beyond optimizing the harmonic lengths of the ansatz ground state, as well as devising a hybrid method including both tight binding and charge-basis diagonalization to accommodate circuits with both localized and delocalized degrees of freedom.
exp(a µ Z µν a ν ) exp(a † µ X µν a † ν ) = exp(a µ Z µν a ν ) exp(λ µ a † µ ) = exp(λ µ Z µν λ ν ) exp(λ µ a † µ ) exp(a µ Z µν a ν ) exp(λ µ (Z µν + Z t µν )a ν ), Here, X, Y, Y , Z and λ are arbitrary, except for the requirement of 1 1 − 4XZ and 1 1 − 4ZX to be nonsingular in Eq. (A4). We note that it is relatively straightforward to obtain Eqs. (A6)-(A8) from standard applications of the BCH formula [37]. Obtaining Eqs. (A4)-(A5) is slightly more difficult, and requires using either Lie algebra techniques [44,46] or the so-called IWOP procedure [45]. An instance of Eq. (A1) relevant for computing wavefunction overlaps is S † m exp( λ t a † ) exp(− λ t a)S m , identifying and neglecting the overall multiplicative factor [c.f. Eq. (24)]. To simplify notation, we have suppressed the dependence of λ on m, m and j. We will continue to likewise suppress the m dependence of the various matrices and distinguish between X m and X m by using the notation X and X , etc. Applying each of the relations Eq. (A4)-(A8) in a few steps of algebra leads to the normal-ordered result where P = (1 1 − XX ) −1 , P t X = 1 2 (P t X + X P ) and the matrices X, X , etc. can be taken to be symmetric. Similar expressions can be obtained when the operator O is an explicit function of the ladder operators a, a † .