Effective Hamiltonians for interacting superconducting qubits -- local basis reduction and the Schrieffer-Wolff transformation

An open question in designing superconducting quantum circuits is how best to reduce the full circuit Hamiltonian which describes their dynamics to an effective two-level qubit Hamiltonian which is appropriate for manipulation of quantum information. Despite advances in numerical methods to simulate the spectral properties of multi-element superconducting circuits, the literature lacks a consistent and effective method of determining the effective qubit Hamiltonian. Here we address this problem by introducing a novel local basis reduction method. This method does not require any ad hoc assumption on the structure of the Hamiltonian such as its linear response to applied fields. We numerically benchmark the local basis reduction method against other Hamiltonian reduction methods in the literature and show that it is applicable over a wider parameter range, particularly for superconducting qubits with reduced anharmonicity, including the capacitively-shunted flux qubit. By combining the local basis reduction method with the Schrieffer-Wolff transformation we further extend its applicability to systems of interacting qubits and use it to extract both non-stoquastic two-qubit Hamiltonians and three-local interaction terms in three-qubit Hamiltonians.


Introduction
Since their first appearance, superconducting (SC) circuits including Josephson junctions have proved to be one of the most promising platforms for quantum information processing applications [4,5,6,7]. The lithographic fabrication process allows fine tuning of the physical properties of each superconducting circuit, thus resulting in qubits with different spectral properties. Individual qubits can be manufactured in large arrays, with electrostatic and magnetic interactions coupling pairs of them. The strength of the local fields on each qubit and of the two-qubit interactions can further be adjusted dynamically by applying external electrostatic and magnetic fields, making for a flexible and scalable architecture for both gate-based quantum computation (GBQC) and quantum annealing (QA) [8,5,4,9].
A two decades quest to improve the coherence metrics of superconducting qubits, by materials and circuit engineering, has led to a number of SC qubit designs, such as capacitively-shunted flux qubits and transmons, having T 1 and T 2 times in the 100 μs range [10,11]. These circuits, as much as the earlier designs, including rf-SQUID qubits [12], persistent-current qubits [8] and single-Cooper-pair boxes [13] are, by construction, char-acterised by the fact that, under specific operation conditions, they can be regarded as two-level systems (in the sense that any additional stationary state of the system has a substantially higher energy and a small probability of being populated) [14].
The fundamental theory describing SC circuits, i.e. quantum network theory, is well established and can be used, at least in some approximate form, to numerically determine the energy spectrum of an arbitrary SC qubit circuit [1,2,3]. The literature seems, however, to be missing an agreed and consistent way of connecting the electromagnetic HamiltonianĤ e.m. of an arbitrary system of n SC qubits to the corresponding effective qubit HamiltonianĤ 1 , or, equivalently, of numerically determining the parameters of an n-spin Hamiltonian which reproduces the low-energy spectrum ofĤ e.m. , as well as the expectation values of some set of the system observables. As we will see below, where such mapping methods do exist (see, for instance, supplementary materials of Ref. [15,16]), they are not guaranteed to reproduce the correct low-energy spectrum of the circuit.
A general scheme for reducing the circuit Hamiltonian of an arbitrary SC qubit system to the correct effective qubit Hamiltonian would serve several purposes. Firstly it would improve our intuitive understand-ing of the system in the language of qubits, allowing us to determine an economical description of the system, which still retains all the information about the computational subspace and which can be easily used for comparison with experiment, or as a tool for further numerical simulations. Secondly, specifically in the context of adiabatic quantum computing (AQC), identification of non-stoquastic and multi-local terms in the qubit Hamiltonian could help the engineering of such terms, which are fundamental to implement nonstoquastic AQC (which is thought to be more powerful than its stoquastic counterpart [17]) and error suppression protocols based on stabiliser codes [18], respectively. Lastly, in the case of single qubits, it would provide an improved way of calculating tunnelling amplitudes between semi-classical potential minima that is an alternative to instanton-based approaches and therefore potentially more accurate, especially in the limit of large tunnelling amplitudes [19,20].
In this paper we propose a method of implementing Hamiltonian reduction based on a natural local definition of the computational basis. Our method does not require any ad hoc assumption on the structure of the Hamiltonian, such as its linear response to the applied electrostatic and magnetic fields, which is at the core of standard perturbative reduction methods [15]. Additionally the scheme can be applied to individual SC qubits of any kind, as well as to systems of qubits and coupler circuits, interacting magnetically or electrostatically. In the interacting case the scheme makes use of the Schrieffer-Wolff transformation to separate the low energy subspace of the Hilbert space from its complement [21].
The article is structured as follows: in the next section we revise how to write a general electromagnetic Hamiltonian for isolated and coupled superconducting circuits. In section 3 we introduce some of the state-of-the-art reduction methods in the literature and then present our novel approach to the problem, in the context of both single and interacting qubits. In section 4 we present the numerical results of Hamiltonian reduction applied to systems of superconducting qubits, with a specific reference to recent publications. Finally we summarise our conclusions.

Circuit Hamiltonians from Quantum Circuit Analysis
Since we want to establish a way to numerically derive an effective qubit Hamiltonian from the full Hamiltonian describing the superconducting circuit, we begin this paper by revising how to write down the circuit Hamiltonian for a generic non-dissipative circuit. We start with isolated circuits and later consider the presence of interactions. The framework which we use is that of quantum network theory, which is the quantum version of Lagrangian mechanics applied to electrical circuits [1,22]. Following the standard procedure we will first write the classical Hamiltonian and then quantise it by replacing the variables with the corresponding operators. The reader who is familiar with these concepts may wish to skip to the next section.

Isolated circuits
The first key assumption we make in order to apply quantum network theory is that the size of our qubit is sufficiently small relative to microwave wavelengths, such that it is appropriate to use a lumped-element description of the circuit [23]. Because the system is superconductive, this circuit will consist of nodes connected by branches containing only non-dissipative elements, namely inductors, capacitors and Josephson junctions. An example representing the equivalent circuit of an rf-SQUID flux qubit is shown in figure 1. Then, without loss of generality, we can arbitrarily assign one of the circuit nodes to ground. (For a floating qubit there will be a capacitor between the ground node and the rest of the circuit.) At this point, in order to later take into account the effect of external magnetic fields, we need to choose a spanning tree, i.e. a path of connected branches going from the ground node to every other node, without generating loops. The specific choice of the spanning tree will not affect our final results [22]. A possible choice of the spanning tree for the rf-SQUID flux qubit in figure 1 is highlighted in red. We will indicate the set of branches in the spanning tree by T and the complementary set of closure branches by C. Every closure branch is associated with an irreducible loop in the circuit, which is the smallest loop formed by that closure branch and by other Figure 1: Equivalent lumped-element circuit of an rf-SQUID flux qubit [12]. One of the two possible choices of the spanning tree is highlighted in red. branches in the spanning tree. For instance, in the flux qubit in Fig. 1 the closure branch b 01 is associated with the single loop in the circuit [22].
Every state of our circuit is defined by specifying the instantaneous voltages at each of the nodes. Alternatively, we can define, for every node j (excluding ground), a node flux variable Φ j , representing the integral over time of its voltage, i.e.
The ground node acts as the voltage reference, so its associated voltage and flux are set to be identically equal to 0 [22]. The node fluxes can be used, together with the voltages, to write down the circuit Lagrangian L e.m. ({Φ i }, {Φ i }), which in turn allows to define the variables canonically conjugate to the node fluxes, i.e. the node charges [22]: For brevity we omit here the derivation of the system Lagrangian (which can be found, for instance, in [22]) and we simply report the final form we obtain for the circuit Hamiltonian, If we take care to define the spanning tree so as not to leave any inductive branch in the closure set C, this takes a particularly simple and general form: where is its linear part, with N the number of circuit nodes, C and L are the (N × N ) capacitance and inductance matrices of the circuit, respectively, (see appendix A.1 for their definition) and where is the Josephson energy component. Here E J,bij is the Josephson energy of the Josephson junction in the branch b ij connecting nodes i and j (the index 0 refers to the ground node here) and Φ 0 = h/(2e) 2.0678 · 10 −15 Wb is the magnetic flux quantum. The branch fluxes {Φ bij } j>i=0,...,N appearing inside the expression are defined as where Φ ext ij is the external magnetic flux threading the irreducible loop associated with b ij .
Our definition of the branch fluxes includes the effect of external magnetic fields on the energy of the system. Current and voltage biases, however, may also be applied to the circuit and each of them will contribute with its own term to the Hamiltonian. In the case of current bias, this is applied through a dangling inductive branch. Let a be the origin node of this branch, L a its inductance and I ext the bias current; the corresponding Hamiltonian term is [22]: In order to apply a voltage bias, a voltage source V g is connected to the desired circuit node a through a gate capacitor C g . The resulting effect on the Hamiltonian is to change the capacitance matrix C →C (to take into account that the total capacitance attached to node a has increased by C g ) and to introduce the additional term [22]: Now that we have put together all the necessary Hamiltonian terms, we can finally obtain the quantum Hamiltonian of the circuitĤ e.m. by simply replacing the variables {Φ j , Q j } i=1,...,N with the corresponding Hermitian operators. These will obey the canonical commutation relations [1]:

Interacting circuits
Let us now consider a system of N superconducting circuits of the kind just considered which are interacting with each other. The total electromagnetic Hamiltonian of the system will have the general formĤ e.m.
is the unperturbed part, withĤ i the Hamiltonian of the i -th circuit, in the form of Eq. (4), and describes the interactions between pairs of different circuits. Here, {Ô i k } k=1,2,... is a set of operators (either node or branch operators) acting on the i -th circuit and the α i k ,j l 's are the interaction constants.
In practice the interactions can be electrostatic, mediated by the charge operators, and magnetostatic, involving the flux operators. (In principle, there could also be additional interactions mediated by Josephson junctions shared between two circuits, but, for simplicity, we will not consider these here.) The electrostatic interaction is achieved by connecting the k-th node of circuit i with the l-th node of circuit j = i with a coupling capacitor C i k ,j l . This has two effects on the system Hamiltonian: it rescales the inverse capacitance matrices of the two circuits (known as capacitive loading), as shown explicitly in appendix A.2, and introduces the interaction term where C −1 m is a suitable inverse mutual capacitance matrix (see appendix A.2) [24].
The magnetostatic interactions are the result of the mutual inductive coupling between pairs of branches belonging to two different circuits, say b i k and b j l . The effect of this mutual inductance is again twofold: it rescales the inverse inductance matrices of the circuits (inductive loading), and introduces in the Hamiltonian the interaction term whereΦ bi is the branch-flux operator associated with the branch b i (see appendix A.2 for the definitions of L −1 i , L −1 j and M −1 ) [24]. Notice that the uncoupled Hamiltonians {Ĥ i } in equation (11) are intended to be corrected for capacitive and inductive loading.

Hamiltonian reduction methods
In this section we review some of the state-of-the-art numerical Hamiltonian reduction approaches and successively introduce two novel protocols, one for single qubits (subsection 3.1) and one for multiple interacting qubits (subsection 3.2). We also point out the key differences between the standard methods and our new method and demonstrate how the latter improves the range of applicability of the reduction. The standard reduction protocols described here will be used in numerical simulations (section 4) for a comparison against the new protocols.

Single qubits
Let us begin by introducing a formal definition of the reduction process. In the case of one isolated qubit, this amounts to finding an effective single-spin Hamiltonian, that is: Definition 3.1 (Effective Single-Qubit Hamiltonian:) A Hermitian operatorĤ q acting on a Hilbert space with dimension 2, whose spectrum matches the two lowest energy eigenstates (E 0 and E 1 ) of the SC qubit circuit HamiltonianĤ e.m. .
Assuming that the SC qubit is at thermal equilibrium with an environment at temperature T, then, if k B T is small compared to the transition energy to the second excited state, E 2 − E 0 , the probability that this state, or any further excited state, is occupied at any given time is exponentially small. In fact, in the absence of any resonant drive term in the Hamiltonian, the higher excited states of the qubit circuit can only be occupied as a result of environment-induced relaxation. The stationary probability that the system occupies a state with energy E i at the end of this process is [25]. Under this hypothesis, the dynamics of the qubit are effectively restricted to the eigenspace associated with the two lowest energy eigenstates of the (potentially timedependent) circuit HamiltonianĤ e.m. (t) (i.e. the qubit subspace H q = Span{|E 0 (t) , |E 1 (t) }) and can be described in terms of an (instantaneous) effective single qubit Hamiltonian [14].
Let us now consider the spectral decomposition of the circuit Hamiltonian, (17) where we have sorted the energy eigenvalues in increasing order. By considering the definition of the qubit Hamiltonian, we see immediately that a good candidate forĤ q is the restriction ofĤ e.m. to the qubit subspace, that is: This expression, however, is not particularly useful to describe the evolution of the qubit in a quantum computation process. In fact, the computational basis used to encode the information on the quantum computer does not correspond, in general, to the system energy eigenbasis. (Note that, in this basis, the Hamiltonian is diagonal, and therefore classical [26].) It is therefore necessary to define the two computational states and their relationship to the energy eigenstates [14].
The computational basis for a superconducting qubit is defined in terms of two eigenstates of an observable which is used in practice to measure the qubit state. This operational definition distinguishes, therefore, between the two main categories of SC qubit design. For circuits of the flux-qubit type (including rf-SQUID qubits[27], three and four-Josephson-junction persistent current qubits [8,28] and C-shunt flux qubits [11]), the computational states are identified with two states with opposite and well-defined values of persistent current in the qubit loop. For charge-qubit-type designs (including single Cooper-pair box qubits [13] and transmons [10]), |0 and |1 are instead identified with states with a different number of Cooper pairs on the superconducting island [4].

Perturbative reduction (PR) method
The usual approach to identifying the computational basis states for theory and simulations, which is extensively used in the literature (cf. for example [15,11,13,29]), is based on a series expansion of the circuit Hamiltonian around a fixed value of one of its bias parameters (voltage or magnetic flux bias). For clarity, let us consider the specific case of the rf-SQUID qubit, whose circuit is shown in figure 1. Following the method introduced in section 2, we can write its circuit Hamiltonian (up to an additive constant) as [12]: where f z = Φ z /Φ 0 := Φ ext 01 /Φ 0 is the magnetic flux applied externally to the rf-SQUID loop, in units of Φ 0 . When f z 0.5, we can rewrite the previous equation as: where δf z = f z − 0.5 and withÎ the loop current operator, which will define our computational basis. Notice that we used Kirchhoff's current law to go from the first to the second line in the last equation [1]. At this point, we can invoke stationary perturbation theory to write the n-th eigenstate ofĤ e.m. (f z ), up to first order in δf z as [30]: rf-SQUID qubit), it is useful to consider the following equivalent derivation, which has a straightforward extension to the interacting qubit case. Once we have found the computational states according to (24), we can use the homomorphism between C 2 and qubit subspace (H q = Span{|0 , |1 }) to introduce the following four operators, which represent the action of the Pauli matrices on H q : Then, using the following property of the Pauli matrices, we find that: Notice that here, bothĤ e.m. andσ I,x,y,z are conveniently expressed in whatever basis we initially choose forĤ e.m. . The perturbative reduction approach has a clear disadvantage: the effective Hamiltonian (25) reproduces the two lowest energy levels of the full circuit Hamiltonian only in the limit in which the first order perturbative expansion (22) holds. This entails two requirements. Firstly that the spectrum of the unperturbed Hamiltonian (in other words, the circuit Hamiltonian at the point of the expansion) is highly anharmonic, which is only true for some SC qubit designs and not for others (such as the capacitively-shunted flux qubit and the transmon) [11,10]. Secondly, the perturbation to the bias parameter must be small, for instance |δf z | 1 for the rf-SQUID qubit [15].

Instanton approach
A second common approach to the numerical calculation of the Pauli coefficients is the use of semi-classical theory. In this case the quantum state of the system is approximated by one that minimises its semi-classical potential, which is the part of the classical Hamiltonian depending on the coordinate variable (i.e. the flux in a flux qubit and the charge in a charge qubit). At the operational point the semi-classical potential of qubit circuits assumes a general double-well shape (or, more generally, that of a system of wells in more than one dimension), with two local minima very close in energy, such that quantum tunnelling can occur between them.
In this picture, the longitudinal Pauli coefficient h z is identified with the difference in energy between the two potential minima, whereas the effective transverse field h x corresponds to the tunnelling energy. This is calculated using the semi-classical instanton method (or equivalently the WKB approximation) [31]. These calculations are only accurate in the limit in which the tunnelling action across the potential barrier is very large, which implies that the tunnelling energy has to be exponentially small [32]. The instanton calculation of the transverse field for the rf-SQUID qubit is described in detail in appendix A.4.

Local basis reduction (LR) method
In order to overcome the difficulties of the standard reduction approaches outlined above, we propose an alternative reduction method which relies on a local definition of the computational basis, i.e. one that explicitly depends on all of the circuit bias parameters. In other words, in this case the computational basis states are built as a linear combination of the two local circuit low-energy states: whereĤ e.m. |E i = E i |E i andĤ e.m. is the local circuit Hamiltonian. In order for these two states to be appropriately orthonormal, the u ij 's have to be the elements of a unitary matrix, which we will have to find. The unitarity condition ensures that when we transform from the energy eigenbasis {|E 0 , |E 1 } to the local computational basis {|0 , |1 } the spectrum of the effective qubit Hamiltonian (18) is unchanged and the two lowest-energy levels of the circuit Hamiltonian are preserved.
Owing to the orthonormality of U columns, we can always rewrite U, up to an irrelevant global phase multiplication factor, as: where θ = acos|u 00 | = acos|u 11 |, so that θ ∈ [0, π] and ϕ 1 , ϕ 2 ∈ [0, π/2]. Now we consider again the operational definition of the computational states. This specifies that these should be eigenstates of a certain observableÔ. For a flux qubit O =Î, the current operator associated with the qubit SC loop, whereas for a charge qubitÔ =Q represents the charge on the qubit SC island. One can easily see that imposing this condition on the states (30) is equivalent to finding the two eigenstates of the operator associated with a non-zero eigenvalue 1 , that is Definition 3.3 (Computational basis states (local)) |0 and |1 such that Notice that this definition coincides with the one used in the perturbative method at the specific bias point at which the Hamiltonian expansion is performed (for instance at f z = 0.5 in the the rf-SQUID qubit case).
By identifying |E 0 with the vector (1, 0) and |E 1 with (0, 1), we can rewriteÔ p in the matrix form Finding the eigenvalues and the eigenvectors of this 2×2 matrix is straightforward. In particular, for the eigenvalues, we have: where t = Tr(O p ) and d = det(O p ). In accordance with the operational definitions given above, we need to enforce one condition on these eigenvalues. For a flux-qubit type circuit, we need to have u 1 < 0 < u 0 , which implies det(I p ) < 0, or, more explicitly For a qubit of the charge type, instead, we will require u 1 = u 0 ± 2e (up to some suitably small numerical error). If this condition is not satisfied, then the circuit cannot be operated as a qubit with the desired computational states and the reduction protocol fails. Note, however, that since we are not making use of a perturbative expansion or a semi-classical approximation here, the range of applicability of this local reduction method should be wider than that of the standard methods presented before. 1 One can easily show thatÔp achieves its maximum rank of two as long asÔ|E 0 andÔ|E 1 are linearly independent.
If we now write the eigenvectors of O p as u 0 = (u 00 , u 01 ) and u 1 = (u 10 , u 11 ), then equation (30) returns our desired computational basis states, which makeÔ p diagonal. Armed with u 0 and u 1 , we can easily calculate the general expression of the effective qubit Hamiltonian in the computational basis. If we keep working with 2×2 matrices, the qubit Hamiltonian in the energy eigenbasis (18) takes the obvious diagonal form Going from this basis to the computational basis amounts to applying the unitary transformation U defined above; this gives the effective qubit Hamiltonian in the computational basis as: We observe that, by rescaling the computational states u 0 and u 1 by two phase factors, say e iφ0 and e iφ1 , i.e. by applying some local gauge transformation in the qubit subspace, we can always remove the imaginary component h y σ σ σ y of H q . In fact such a gauge transformation G(φ 0 , φ 1 ) corresponds to a spin rotation around the z axis, multiplied by a global phase: (43) Hence G(−ϕ 1 , −ϕ 2 ), which represents a rotation around z by the angle ϕ 2 − ϕ 1 = −ϕ (followed by a rescaling by e −i(ϕ1+ϕ2)/2 ), transforms cos ϕσ σ σ x + sin ϕσ σ σ y into σ σ σ x , and makes the effective qubit Hamiltonian real, that is: Notice that this gauge transformation can equivalently be written as: Expression (44) for the effective qubit Hamiltonian, is the one adopted by most of the literature on SC qubits [23,15,33]. (Note that the coefficient ∆ is usually further assumed to be positive, a condition which can also always be achieved with a π rotation about z.) An equivalent and more convenient way of calculating the four Pauli coefficients h i , i = I, x, y, z than using equations (33), (42) and (26) together is again to use the computational states to build the Pauli operators and then to apply equation (29).
In section 4 we will present numerical simulations which benchmark the performance of the local reduction method against the standard methods and demonstrate the increased accuracy of the former relative to the latter ones.

Multiple qubits
Let us now consider the Hamiltonian reduction process in the case of multiple interacting superconducting qubits. Given a system of N qubits and M additional coupling circuits, coupled inductively and/or capacitively, its effective qubit Hamiltonian is one that reproduces the lowest 2 N energy levels of the total system Hamiltonian, as well as the expectation values of the qubit operators. Notice that any such Hamiltonian can be written in the general form where η = (η 1 , . . . , η N ), η i ∈ {I, x, y, z} and σ σ σ η = σ σ σ η1 ⊗ · · · ⊗ σ σ σ η N is a 2 N × 2 N matrix in the Pauli group G N . Recalling the equality (28) and using the following property of the trace we can see that the real Pauli coefficients h η obey the equation According to section 2.2, the circuit Hamiltonian of the system can be written aŝ withĤ i (Ĥ c,i ) the unperturbed Hamiltonian of the i -th qubit (coupler) circuit and whereĤ int includes all the interaction terms. Notice that the unperturbed Hamiltonians are assumed to be corrected for capacitive and inductive loading (cf. section 2.2). Now we can define the qubit subspace, in analogy with the single-qubit case, to be the one spanned by the lowest two eigenstates the unperturbed Hamiltonian of each qubit. Since the couplers are designed to be classical elements which always remain in their ground state, while adiabatically following the qubits, the qubit subspace will at the same time be the one spanned by the ground state of each coupler circuit Hamiltonian. We therefore have, in symbolic form: One could then think of defining the qubit Hamiltonian for this N -qubit system simply as in Eq. (18): where againP 0 is the projector on H q . This operatorĤ q , however, does not have the correct spectrum, matching the lowest 2 N energy levels ofĤ e.m. , and therefore does not satisfy our initial definition of qubit Hamiltonian. The reason for this is that the interaction described byĤ int mixes the states in H q with those outside it, i.e. the higher excited states of the individual circuits. Such mixed states become the new low-energy eigenstates ofĤ e.m. [21,16].
Contrary to the single-qubit case, the literature concerning Hamiltonian reduction for multiple interacting SC qubits is relatively scarce. In the following subsections we present two protocols adopted in recent publications and later present a new alternative reduction method, which overcomes some of their limitations and explicitly addresses the problem of the mixing of the qubit subspace with the rest of the Hilbert space by using the Schrieffer-Wolff transformation theory [21].

Approximate rotation method
In this subsection we briefly review the reduction method outlined in a recent work by Ozfidan et al. [16]. This method starts by writing the low-energy part of the total circuit HamiltonianĤ e.m. , i.e. the component associated with its lowest 2 N eigenvalues, in its diagonal form: Since orthogonal operations do not change the spectrum of an operator, this protocol guarantees by construction that the spectrum of H q matches the low-energy spectrum of the circuit Hamiltonian.
The first rotation applied in this protocol, R 1 , maps from the low-energy eigenbasis of the total Hamiltonian H e.m. , {|E 0 , . . . , |E 2 N −1 } to that of the unperturbed 2 N −1 }, and is initially calculated as However, as we pointed out before, |E i also has components outside of the subspace Span{E 2 N −1 }, which implies that this matrix is not orthogonal. R 1 must therefore be explicitly orthonormalised, for instance using the Gram-Schmidt procedure. This step is only justified if the columns of R 1 are already approximately orthonormal [16]. Since in our case orthogonality follows from normalisation, it suffices to check that before we apply the Gram-Schmidt procedure.
To obtain the qubit Hamiltonian we now need the second rotation R 2 to map from the basis of the energy eigenstates |E (0) 0 , . . . to the computational basis. We then take where |i = |i 2 N −1 ⊗ · · · ⊗ |i 0 is an outer product of single qubit computational states, with These computational states are found from the reduction of the unperturbed single-qubit Hamiltonians. If the local reduction method is used for this, the rotation matrix R 2 is guaranteed to be orthogonal. Note that, although the effective Hamiltonian calculated with this method has the correct spectrum, the procedure is based on the approximate equality (52), which is not often satisfied, particularly in the case of relatively large interactions. (This can be seen by considering, once again, the perturbative expansion (22).) Additionally, the previous derivation implicitly assumes that the circuit Hamiltonian is real, so that all the eigenstates and computational states can be chosen to have only real components. This ensures that R 1 , R 2 ∈ SO(2 N ). Some circuits, however, may have an efficient matrix representation of the Hamiltonian which is complex. In this case the definition of the two rotations would lead to the presence of arbitrary complex phases in their elements, which would need to be somehow taken care of. (Notice that even in the real case the scalar products defining the elements of R 1 and R 2 are only defined up to an arbitrary sign.)

Diagonal Hamiltonian method
A second method of determining the effective Hamiltonian of a multi-qubit system is presented in a recent work by Melanson et al. [34]. This method works under the more restrictive assumption that the effective Hamiltonian is diagonal in the computational basis. In this case the lowest 2 N eigenstates of the circuit are also eigenstates of the single-qubit operatorsÔ i specifying the computational basis and the corresponding eigenvalues can be calculated numerically as the expectation values E n |Ô i |E n . Additionally the 2 N non-zero Pauli coefficients of the system can be expressed as a linear combination of its low-energy eigenvalues [34]. For instance, in the two-qubit case one has: 1} is the eigenvalue of the circuit Hamiltonian corresponding to the computational state |i |j . We can therefore determine the Pauli coefficients of the two-qubit system by finding the lowest four energy eigenvalues of its circuit Hamiltonian, calculating the expectation value of the operatorsÔ 1,2 on the each eigenstate to identify its corresponding computational state and by inverting the previous equation to get The same procedure can be applied to systems with three or more qubits (plus eventual additional couplers).
In practice, an effective Hamiltonian diagonal in the computational basis is verified when the qubit tunnelling barriers are high (negligible transverse field h x ) and the qubits are coupled only through their z degree of freedom (that is when the coupling is inductive between flux qubits or capacitive between charge qubits). A Hamiltonian of this form is however classical and cannot be sufficient for universal quantum computation [26]. This method can nevertheless still be useful when it is reasonable to assume that the different non-commuting terms of the qubit Hamiltonian can be turned on and off independently.

Schrieffer-Wolff transformation method
In this final subsection we introduce a new reduction protocol for multi-qubit systems which overcomes some of the limitations of the methods described above. In particular, this method does not require the mixing between the qubit subspace (Eq. (50)) and its complement, resulting from the interactions, to be negligible, which is a crucial assumption of the approximate rotation reduction. Secondly, unlike the approximate rotation reduction, it can be applied directly to circuit Hamiltonians with complex elements, since the arbitrary phase choices made when numerically evaluating the Hamiltonian eigenvectors cancel out in all the necessary expressions. Thirdly, the reduction method introduced here can be applied to find arbitrary non-diagonal effective Hamiltonians. This is all made possible by the Schrieffer-Wolff transformation (SWT), which by construction maps the total circuit HamiltonianĤ e.m. to a new Hermitian operator acting on the qubit subspace H q and whose spectrum matches the low-energy spectrum ofĤ e.m. , which is precisely what we expect from the effective qubit Hamiltonian [21].
The SWT relies on a single assumption regarding the form of the full system Hamiltonian, namely that the spectrum of the unperturbed part of the Hamiltonian (excluding the interactions) has a sufficiently large gap, as we will see below. For the purpose of this reduction method, we will replace this assumption with an equivalent pair of two distinct conditions. In order to state the first one, let us rewrite the unperturbed part of the N -qubit M -coupler system Hamiltonian (49) aŝ is the projector on the low-energy eigenspace H low , spanned by the eigenstates corresponding to the lowest 2 N eigenvalues ofĤ 0 , andQ 0 =Î −P 0 projects on the complementary subspace H \Ĥ The first assumption of our reduction is that H q ≡Ĥ (0) low , i.e. that no additional excited state of the independent circuits is mixed in the low energy subspace ofĤ 0 , and that the two sets S 2 N −1 | ≥ ∆. This composite condition can be written, more explicitly, in the following form: where we have introduced the notation ∆E i,j = E i,j − E i,0 and ∆E ci,j = E ci,j − E ci,0 (with E i,j (E ci,j ) again the j -th eigenstate of the i -th qubit (coupler) unperturbed Hamiltonian). Since the summations above grow linearly with the number of qubits in the system, this condition limits the size of the systems to which we can apply our reduction method. Intuitively this limit reflects the impossibility of finding any coherent description of the low-energy spectrum of a composite system in terms of interacting two-level subsystems, whenever the second excited state of one of these subsystems appears in the spectrum. Therefore if we are interested in characterising a very large circuit, we should first subdivide it into smaller connected subsystems for which the inequalities (58) hold.
The second requirement is simply that the strength of the interaction Hamiltonian should be small compared to the spectral gap ofĤ 0 , ∆. Namely: where · op is the operator norm: with · the 2-norm · | · . Since the addition of the interaction termĤ int can shift the eigenvalues ofĤ 0 by at most Ĥ int op , this second inequality implies that the spectrum ofĤ e.m. remains gapped. This in turn allows us to rewrite the total Hamiltonian in the block-diagonal formĤ e.m. = PĤ e.m.P +QĤ e.m.Q , whereP is the projector on the 2 N -dimensional low-energy eigenspace ofĤ e.m. , H low , andQ =Î −P [21].
Additionally, according to [21], since H low and H q have the same dimension, they are connected by a direct rotationÛ such that U is called the Schrieffer-Wolff transformation and can be written, in terms of the projectors, as [21]: The principal square root √ · above is well-defined as long as P −P 0 op < 1, which in our case can be shown to be equivalent to (59) [21]. Now the action of the SWT onĤ e.m. is given bŷ where we used the identitiesÛP =P 0Û andÛQ =Q 0Û . According to equation (64),ÛĤ e.m.Û † is block-diagonal with respect toP 0 andQ 0 . This finally leads us to the conclusion thatĤ is an Hermitian operator, acting on H q , whose 2 N nonzero eigenvalues are the same as the lowest eigenvalues of the original interacting HamiltonianĤ e.m. (because the unitaryÛ leaves the spectrum ofPĤ e.m.P unchanged) [21].Ĥ q therefore represents our effective qubit Hamiltonian, from which we can directly extract the Pauli coefficients by rewriting equation (48) as In this case, the Pauli operatorσ η =σ η1 ⊗· · ·⊗σ η N ⊗P c is built from the single-qubit Pauli operators {σ ηi }, which, in turn, are obtained as in the single-qubit case, starting from the unperturbed HamiltonianĤ i of each qubit and the appropriate operatorÔ p,i . The operator represents the required identities acting on each of the ground-state energy subspaces of the coupler circuits. Finally note that since both the approximate rotation reduction and the SWT reduction method determine an effective qubit Hamiltonian with the correct spectrum, the two results must be equivalent up to a unitary transformation. However, as mentioned before, the SWT reduction extends the range of applicability of the method to Hamiltonians with complex elements and does not involve the restrictive assumption (52).

Numerical results
In this section we present some numerical examples of Hamiltonian reduction for different SC qubit designs and interacting systems. For concreteness, we will focus on qubits of the flux-type ( [8,11,27,28]) and we will consider circuits and physical parameters from works in the recent literature.
For these simulations the circuit Hamiltonians and all other circuit operators were represented in matrix form by projection on a truncated orthonormal basis. The approximate Krylov-Schur method, implemented by the MATLAB c function eigs [35], was used to determine the relevant subsets of the operator eigenvalue-eigenvector pairs. This approach can be much faster than the complete diagonalisation of the operator, especially when it is very large and sparse, as is usually true for SC qubit Hamiltonians [36].

rf-SQUID flux qubit
We start by considering the simplest example of a flux qubit, i.e. the rf-SQUID circuit. As shown in figure  1, this consists of a Josephson junction, with tunnelling energy E J , shunted by a superconducting inductive loop with self-inductance L and in parallel with its intrinsic capacitance C J [4]. Figure 2 shows the lowest five energy eigenvalues of the circuit, calculated as a function of the dimensionless external magnetic flux f z = Φ z /Φ 0 ≡ Φ ext 01 /Φ 0 . (Note that the constant offset E 0 (f z = 0.49) has been subtracted from all the energies.) The parameters used for the simulations are E J = 125 GHz, C J = 5 fF and L = 2.5 nH, which are typical for this type of device [12]. In this case, the Hamiltonian was represented in a basis of harmonic oscillator occupation number states, truncated at a maximum occupation number of 40, which ensured the convergence of the low energy spectrum (cf. appendix A.3) [37].
As we can see from the graph in figure 2, the lowest two energy levels of the system (i.e. the qubit states), vary approximately linearly with the flux f z , except around the symmetry point f z = 0.5, where they show a characteristic avoided crossing. In fact, as we saw previously, for small values of |δf z | = |f z −0.5| the rf-SQUID Hamiltonian is well approximated by its first order expansion in δf z . This maps to an effective qubit Hamiltonian of the form (see Eq. (25)) where we have neglected the term proportional to the identity, ∆ = (E 1 (f z ) − E 0 (f z )) fz=0. 5 and ε(f z ) = 2Φ 0 I p δf z . The lowest two energy levels of the circuit are therefore approximately E 0,1 = const. ∓ ∆ 2 + 4Φ 2 0 I 2 p |δf z | 2 , which become linear in f z for larger values of |δf z |. Figure 3a shows the values of the system Pauli coefficients as a function of f z , calculated using equation (29). The solid lines correspond to values obtained by defining the Pauli operators according to the local reduction (LR) method introduced here (subsection 3.1.3). These are compared with the result of the perturbative (PR, empty circles) and instanton (crosses) methods. As we can see, the three reduction methods produce largely compatible results for this circuit. In particular, away from the symmetry point the LR method finds a 10% increase in the transverse field h x at the boundary of the flux interval considered, compared to its centre. The  result of PR is instead independent of f z , in agreement with Eq. (25). The values of h z and h I calculated with the LR and PR methods are compatible to 1% over the whole flux bias range. This implies that the definition of the computational basis in the LR method coincides, as it should, with that of the standard PR method in the limit in which the series expansion (20) and perturbation theory apply. As for the semi-classical calculations, these appear to over-estimate both the longitudinal field (by 40%) and the transverse field (by up to 6%), compared to the other two reduction methods.
Since the semi-classical approximation applies in the limit where is much smaller than the actions at play in the system, i.e. S , and since the tunnelling energy h x decreases exponentially with the tunnelling action, h x ∝ e −S/ (see appendix A.4), we expect the result of the instanton calculations to be more accurate in the limit where h x is small [31,32]. To verify this, we determined the qubit transverse field in the case of biasing at the symmetry point f z = 0.5 for increasing values of the loop inductance L. As we can see in figure 3b, increasing L causes the barrier between the two semi-classical potential wells (blue line) to rise, therefore suppressing the tunnelling h x (data in red). Since f z = 0.5, the perturbative and local reduction methods coincide, and they both determine the correct value of the tunnelling energy: h x = −∆E/2, where ∆E is the energy separation between the ground and first excited state of the circuit (cf. dots and solid line in Fig. 2). As expected, the instanton method result (crosses in Fig. 2) closely approaches that of the Hamiltonian reduction only as L increases and |h x | becomes smaller. At this point, as a consistency check, we can calculate the spectrum of the reduced qubit Hamiltonian simply as E 0,1 = h I ∓ h 2 x + h 2 y + h 2 z and compare it with that obtained from the full circuit model. PR and LR do a good job in reproducing the low-energy spectrum of the rf-SQUID qubit, as we can see from the plot in figure  3c. This also shows the spectrum derived from the semiclassical model (crosses), which does not agree with the correct circuit spectrum as well.
Notice that LR is guaranteed to exactly reproduce the circuit levels as long as f z 0.5. As mentioned in the previous section, the LR protocol only fails when, as |f z − 0.5| increases, the two eigenvalues ofÎ p (f z ) = P 0 (f z )ÎP 0 (f z ) begin to have the same sign, meaning that no measurement distinguishing two qubit states with opposite persistent current is possible at the given bias. For the particular rf-SQUID circuit considered here, the local reduction method breaks down for |f z − 0.5| 0.035, as shown in figure 3d (region shaded in red). As we can see in this plot, as we approach this region the behaviour of the Pauli coefficients starts changing. In particular the transverse field increases considerably in magnitude, while the longitudinal field saturates. The green dotted lines in figure 3d show the circuit energy levels. We see that at the boundary of the unshaded region the second excited state starts mixing with the first, leading to an avoided crossing. This mixing means that, at this point, the two-level approximation does not hold any more, which leads to the failure of the LR.
Finally we might want to consider how well the reduced Hamiltonians are able to reproduce the correct expectation values of some circuit operatorÔ, i.e. whether the following relationship holds where {|E i } i=0,1 are energy eigenstates ofĤ e.m. and {|E i } i=0,1 are the eigenstates of the corresponding effective qubit Hamiltonian.Ô p is defined locally aŝ P 0 (f z )ÔP 0 (f z ) (whereP 0 (f z ) is the projector on the twodimensional low-energy subspace ofĤ e.m. (f z )) in the LR method case, and is defined globally asP 0 (0.5)ÔP 0 (0.5) in the PR case. Figure 4 shows the matrix elements of the loop current operatorÎ between qubit states, calculated with both the full and the reduced operators. We observe that LR is ensured to give the exact result, while PR produces a reasonable result.
We have shown here that the approximations inherent in the perturbative reduction method are valid and sufficient for for determining the reduced Hamiltonian in the case of the simple rf-SQUID qubit of Fig. 1. We will see in the next subsection, however, that this is not true in the general case and that the local reduction method has a wider range of validity.

C-shunt flux qubit
The accuracy of the perturbative reduction method deteriorates when we consider other flux qubit designs, particularly those with reduced anharmonicity like the capacitively-shunted flux qubit shown in figure 5. This consists of a superconducting loop interrupted by three Josephson junctions. The area of one junction is a factor α < 1 smaller than that of the other two and is shunted by a relatively large capacitor C sh C JT . The capacitive shunt reduces the qubit sensitivity to charge noise, while improving the device reproducibility (by compen- Figure 5: Equivalent lumped-element circuit of a capacitively-shunted flux qubit (as described in Ref. [11]). A possible choice of the spanning tree is highlighted in red.
sating for the fabrication variability of the junction size, which affects C JT ). At the same time, the effect of flux noise is mitigated by choosing small values of α (typically 0.125 < α < 0.5), which reduce the magnitude of the persistent current and therefore the magnetic dipole moment of the circuit [11]. The result is superconducting qubits with typical measured relaxation times T 1 in excess of 40μs (three orders of magnitude longer than the standard rf-SQUID T 1 ) and decoherence times approaching the relaxation limit T 2 = 2T 1 [11].
This substantial coherence enhancement comes at the cost of a decrease in the spectrum anharmonicity. We can see this by looking at figure 6a, which shows the calculated low energy spectrum of a C-shunt qubit circuit as a function of f z = Φ z /Φ 0 = Φ ext 23 /Φ 0 , and comparing it with Fig. 2. The physical parameters used for the simulation are shown in table 1 (cf. Fig. 5 for the mean-ing of the symbols). For the two lower junctions we used E JL = E JR = E JT /α and C JL = C JR = C JT /α. These parameters are compatible with those reported in the experiments in Ref. [11]. In this case, the Hamiltonian Parameter Value E JT 45 GHz C JT 1.8 fF α 0.43 C sh 50 fF L 100 pH Table 1 was represented numerically by projecting on a finite basis consisting of harmonic oscillator states for the mode associated with the circuit node 1 and charge number  states for the modes associated with nodes 2 and 3 (cf. appendix A.3) [37,33,38]. As we can see from figure 6a, the two dispersion relations E 0,1 (f z ) have first derivatives with the same sign everywhere. Since this means that the average persistent currents in the two energy eigenstates have equal sign (cf. Fig. 6d). This is in contrast with the rf-SQUID flux qubit [11], but does not preclude the possibility to find two current eigenstates with opposite sign in the qubit subspace. Figure 6b shows the Pauli coefficients obtained by the perturbative (circles) and local (lines) reduction methods. As anticipated, there is a clear discrepancy between the two results. In fact, owing to the much smaller anharmonicity of this circuit compared to the rf-SQUID, the two low-energy eigenstates of the circuit Hamiltonian at f z = 0.5 are not a good approximation for those away from f z = 0.5. This implies that projectingĤ e.m. (f z ) on the states (24) does not preserve its low-energy spectrum and does not lead to the correct reduction. From the numerical results we see that the slope of h z (f z ) in the local reduction case is smaller than in the perturbative reduction and further decreases away from f z = 0.5. Additionally, the transverse field h x (f z ) shows a clear negative curvature in the LR results, whereas it is roughly constant in f z in the PR case (as in the rf-SQUID). The strong dependence of the transverse field on f z is a known distinguishing feature of the C-shunt flux qubit design when compared to more standard flux qubit circuits like the rf-SQUID [8,11].
Calculating the spectra of the two reduced Hamiltonians leads to the result shown in figure 6c. The local reduction result (filled dots) again reproduces the circuit ground and first excited states (lines) exactly, while the perturbative reduction fails to accurately predict the first excited state. Finally figure 6d shows the matrix elements of the current operator between the qubit energy eigenstates, calculated using the full circuit model (lines) and the two reduced two-level models (circles). The PR (empty circles) gives incorrect expectation values, which are opposite in sign for the two states.

ZZ plus XX coupling
We begin this subsection on coupled SC qubit systems by considering a simple two-qubit system, without any nonlinear coupling element. As one such example we consider the system which Ozfidan et al. characterised experimentally in [16]. This is composed of two compound-Josephson-junction rf-SQUID qubits (where the single Josephson junction is replaced by two junctions in paral-lel, forming a dc-SQUID) coupled both inductively and capacitively, as shown in figure 7. Assuming that the dc-SQUID loop is very small (such that its inductance is much smaller than the both main loop inductance and the Josephson inductance (Φ 0 /2π) 2 /E J ), we can effectively describe it as a single junction whose Josephson energy depends on the flux Φ x threading the dc-SQUID [39]: here is the sum of the energies of the two junctions in parallel, which we are assuming to be equal. Within this approximation, the Hamiltonian describing our circuit is: where C 1(2) = C 1(2) +C 12 C 2(1) /(C 2(1) +C 12 ) and L 1(2) = L 1(2) − M 2 12 /L 2(1) [16]. Using the physical parameters given in Ref. [16], i.e. C 12 = 132fF and those in table 2, and calculating the lowest four eigenvalues of our Hamiltonian for different values of mutual inductance in the range −2pH < M 12 < 2pH, we obtained the graph shown in figure 8a. This graph matches well with the corresponding one present in Fig. 3c of Ref. [16]. The avoided levelcrossing at M 12 0.7pH is proportional to the capacitive coupling C 12 and only occurs at finite longitudinal fields, Figure 7: Circuit diagram of the system of two interacting qubits studied in [16]. Highlighted in different colours are the coupling elements and the magnetic bias fluxes.
i.e. Φ z,i = 0 [16]. (Notice that when −1 < Φ x,i /Φ 0 < 0, the effective Josephson energy E J,i (Φ x,i ) is negative and the symmetry point where h z = 0 is displaced from Φ z,i = Φ 0 /2 to Φ z,i = 0 [39].) It is worth noting that, in order to efficiently represent a composite circuit Hamiltonian like (72), we cannot retain the representation of the circuit operators that we used for single circuits. In that case, the size of the total Hamiltonian matrix would equal the product of the sizes of all the individual circuit Hamiltonians, and would rapidly become unmanageable. Since we are, once again, only interested in the low-energy properties of the system, a good alternative basis choice is that of the outer products of some small number N i of low-energy eigenstates of each unperturbed (i.e. non-interacting) circuit HamiltonianĤ i /Ĥ c,i . In this case, for example, we can write: with the meanings of the symbols introduced before. To ensure the convergence of our results, we first used 40 1.603 · 10 3 119.5 231.9 −0.6538 1 · 10 −4 Q2 1.568 · 10 3 116.4 239 −0.6526 1 · 10 −4 Table 2 harmonic oscillator number states to represent the single qubit Hamiltonians and then projected onto their N 1 = N 2 = 10 lowest-energy eigenstates. Now that we have determined the low energy spectrum of the system, we can apply some reduction method to calculate the effective qubit Hamiltonian. We begin with the Schrieffer-Wolff transformation method, introduced in section 3.2.3. After verifying that the hypotheses of its construction are satisfied, in particular observing that P −P 0 op 0.5 in the whole range of M 12 , we extracted the Pauli coefficients. These were calculated by defining the computational states and the Pauli operators locally for each qubit and then projecting the effective qubit Hamiltonian on them, as shown in section 3.2.3. The six one-local coefficients are shown in figure 8b by solid lines. (We do not consider the coefficient h II = Tr(Ĥ q ) here since we are focusing on relative energies.) The dashed lines represent the same coefficients obtained by applying the SWT reduction to the non-interacting part of the circuit Hamiltonian, i.e. to the sum of the Hamiltonians of the isolated qubits (corrected for the static inductive and capacitive loading). Since for h zI and h Iz the solid and the dashed lines overlap, the values of the longitudinal fields of the coupled system are completely determined by the static loading of the unperturbed Hamiltonians. This effect appears approximately linear in M 12 . The values of the transverse fields for the coupled system, instead, are ∼ 25% lower in magnitude than those resulting from the loaded single-qubit Hamiltonians. The interaction with the other qubit, then, has an additional effect, which we call dynamic loading. The change in (a) Low energy circuit spectrum, relative to the ground state, of the circuit in figure 7, as a function of the mutual inductance M 12 .
(b) One-local Pauli coefficients calculated, as a function of M 12 , by applying the Schrieffer-Wolff transformation reduction method to the full (solid lines) and the unperturbed Hamiltonian (dashed lines) of the circuit in Fig. 7 (notice that the solid and dashed lines for h zI and h Iz all overlap at this scale). Circles: same coefficients, calculated using the approximate rotation reduction of [16].  (b) One and two-local Pauli coefficients determined with the approximate rotation method, after the application of the local rotation removing XZ and ZX terms (circles), compared against the ones calculated with the SWT method (solid lines). Note that the local Pauli coefficients for the two qubits overlap almost completely at this scale. Figure 9 transverse field appears approximately quadratic in M 12 and is not centred around M 12 = 0 due to the presence of the capacitive coupling (as we verified by comparing against the case C 12 = 0). As usual, the components of the local field along the y direction have been removed by making the appropriate local gauge transformation. (Actually, the circuit Hamiltonian in this case is completely real, so that no imaginary terms can appear in the reduced Hamiltonian; the gauge transformation only ensures that the signs of different coefficients are consistent across the range of M 12 .) The empty circles in figure 8b are the one-local Pauli coefficients determined with the approximate rotation method, introduced in [16] and reviewed in section 3.2.1. Comparing with the previous results, we can see that we obtain qualitatively similar, but quantitatively different results. In particular the values for the transverse fields are close to those obtained with the SWT reduction, while the new longitudinal fields are everywhere smaller in magnitude, and, in this case, do not agree with their unperturbed values (dashed lines). Figure 9a shows the coefficients of the nine effective qubit Hamiltonian two-local terms. According to the reduction based on the SWT (lines), the only nonnegligible terms in the Hamiltonian are those proportional to σ σ σ z,1 σ σ σ z,2 , σ σ σ x,1 σ σ σ x,2 and σ σ σ y,1 σ σ σ y,2 . The first term represents the inductive interaction,Û M ∝ M 12Φ1Φ2 , the flux being our z degree of freedom, and it indeed scales linearly with M 12 . Since we have chosen to identify a flux degree of freedom with the real operator σ σ σ z , the canonically conjugate charge operator must be com-plex (since [Φ,Q] = i ), and therefore must be identified with σ σ σ y . The YY term, then, describes the capacitive interaction and, in fact, appears to be largely independent of M 12 . Finally the XX term is a result of the presence of the higher excited states of the system [16]. It is related to both the inductive and the capacitive Hamiltonian terms and appears to scale linearly with M 12 .
According to reference [40], a two-local two-qubit Hamiltonian of the form H = h xI σ σ σ x,1 + h Ix σ σ σ x,2 + h zI σ σ σ z,1 + h Iz σ σ σ z,2 + +h xx σ σ σ x,1 σ σ σ x,2 + h yy σ σ σ y,1 σ σ σ y,2 + h zz σ σ σ z,1 σ σ σ z,2 is non-stoquastic, and remains such after arbitrary local rotations, as long as h xI , h Ix , h zI , h Iz = 0 and |h yy | > |h xx |, |h zz |. The region where this condition is satisfied is highlighted in green in figure 9a. Non-stoquastic twolocal catalyst Hamiltonians are know to provide an exponential speed-up to the convergence of quantum adiabatic optimisation, at least with specific problem classes, including the ferromagnetic p-spin model [41]. For this reason, they might be key to establish a quantum advantage over classical optimisation routines such as Quantum Monte Carlo [16,42]. Again, our implementation of the approximate rotation reduction produces qualitatively similar results to the SWT reduction for the two-local Pauli coefficients (see hollow circles in figure 9a), except for h xz h zx (purple circles), which are now of the same order of magnitude as the other coefficients. As we mentioned in section 3.2.3, the approximate rotation and the SWT re-duction methods actually find equivalent effective qubit Hamiltonians, modulo a unitary. This was in fact verified by showing that both sets of coefficients lead to qubit Hamiltonians with the same spectrum.
Notice that Ref. [16] actually reports the two h xz h zx coefficients to be negligible, which we ascribe to the fact that the authors used a different form for the circuit Hamiltonian, and potentially a different definition of the computational basis, and hence of R 2 , as defined in section 3.2.1 [16]. (In our case the computational basis was defined locally as shown in section 3.1.3.) In fact, any mixed two-local term, involving different Pauli operators acting on the two qubits, can be eliminated from a two-qubit Hamiltonian by performing a local change of basis [40]. Applying this transformation produces a new set of coefficients which are within 5% of those found by the SWT reduction method (see Fig. 9b). In this case, then, the unitary mapping between the two is a local transformation.

ZZZ coupling
As the final example we consider a proposed circuit implementing a three-local ZZZ interaction between three flux qubits, presented in [34]. The circuit diagram is shown in figure 10a and consists of the three flux qubits (in this case rf-SQUID qubits) and two compound-Josephson-junction rf-SQUID couplers. The main loops of the two couplers, one of which contains a twist, medi-(a) Circuit diagram of three flux qubits and the three-local ZZZ interaction circuit, described in [34], consisting of two compound-Josephson-junction rf-SQUID tunable magnetic couplers, one of which, c 2 , has a twist in the main loop.
(b) Solid lines: Pauli coefficients extracted, using the SWTbased reduction method, for the system of three qubits and coupler c 1 , as a function of f x,c1 . (Note that c 2 is absent here.) Filled circles: same coefficients, extracted using the diagonal Hamiltonian reduction method.  ate a magnetic interaction between the superconducting loops of qubits q 1 and q 2 (see Fig. 10a). If the flux applied to the coupler main loop, Φ z,ci is kept constant, the flux applied to its dc-SQUID loop, Φ x,ci , controls the effective mutual inductance between the qubits and therefore the magnitude and sign of the effective ZZ interaction [39]. By magnetically coupling the current loop of qubit q 3 to the coupler dc-SQUID loop, one can control the two local interaction between q 1 and q 2 with the current state of q 3 , therefore obtaining a three-local h zzzσz1σz2σz3 interaction [34].
The solid lines in figure 10b show the effective Hamiltonian coefficients for the system consisting of the three flux qubits and the single coupler c 1 , extracted using the SWT reduction method. The main loop of the coupler and those of the three qubits are all biased at Φ z,c1 = Φ z,i = Φ 0 /2, such that the qubit longitudinal fields are all zero. The transverse fields are also zero with the physical parameters considered (which are given below). As expected, we find a three-local interaction term ∝ h zzz , in addition to a residual two-local interaction between qubits q 1 and q 2 , ∝ h zzI and a large longitudinal field h IIz on qubit q 3 .
The parameters used in the simulations are as follows: all qubits (i = 1, 2, 3) have E J,i = 99.3GHz, L q,i = 4.5nH and a large shunting capacitance C sh,i = 45fF; the two coupler junction Josephson energies are E J1,c1 = E J2,c1 = 233.4GHz, the coupler main loop inductance is L z,c1 = 550pH, while the small loop has an inductance of L x,c1 = 170pH and is shunted by a capacitance C sh,c1 = 10fF; all mutual inductances are 50pH. As in the previous simulations, the rf-SQUID qubit Hamiltonians have been expressed in a basis of 40 occupation number states. The three degrees of freedom of the coupler are expressed using 20 occupation number states for the small plasma frequency mode and 7 for the higher plasma frequency modes. The total Hamiltonian is projected on the lowest 8 unperturbed eigenstates of each qubit and on the lowest 5 unperturbed coupler eigenstates. Since the effective Hamiltonian here is diagonal in the computational basis, its coefficients can also be calculated with the method used in [34] and reviewed in section 3.2.2. The result of this reduction is represented by the filled dots in figure 10b and matches very well with the result of the SWT reduction.
Introducing a twist in the coupler, for instance changing the mutual inductance between the coupler and qubit q 2 from 50pH to −50pH (as in coupler c 2 ), and changing the sign of the coupler x -bias, reverses the sign not only of the two-local coefficient h zzI , but also of h IIz . The three-local interaction coefficient, however, remains of the same sign. Therefore attaching both couplers c 1 and c 2 to the qubits leaves us with a purely three-local Hamiltonian. The numerical simulation of the full system agrees with this picture. The Pauli coefficients extracted, as a function of f x,c1 = −f x,c2 , are shown in figure 10c, with the solid lines and the dots being the result of the SWT and the diagonal Hamiltonian reduction method, respectively. Coupler c 2 shares the same physical parameters as c 1 and is also biased at f z,c2 = 0.5. Its lowest 5 unperturbed eigenstates are kept for representing the full system Hamiltonian. As we can see, the size of the three-local ZZZ interaction can be changed from zero to as much as 700MHz in the range of fluxes considered. Its sign can also be changed to negative by biasing at f x,c1 = −f x,c2 = 2 − f x,c1 [34].
Finally we can check that the reduced Hamiltonian has the correct spectrum. This is shown in figure 10d, where the filled dots represent the effective qubit Hamiltonian transition energies and the solid lines those of the circuit Hamiltonian. The levels are grouped in two manifolds each of four degenerate levels, separated by an energy of 2|h zzz |. In the ground state manifold the expectation value of the product of the qubit currents, Î 1Î2Î3 , and therefore σ z1σz2σz3 in the reduced model, is negative, while it is positive in the excited manifold states. At energies above 8GHz we see the additional states of the system, specifically the first excited states of the couplers. As we can see, the interaction does not close the spectral gap of the Hamiltonian, which allows us to use the Schrieffer-Wolff transformation reduction method.

Conclusions
We have developed a systematic numerical method for determining the effective spin Hamiltonian, written in the appropriate computational basis, describing a system of interacting superconducting circuits. Our starting point was a numerical representation of the circuit Hamiltonian, in which each component is described as a lumped-element circuit, with potential magnetic and electrostatic biases, and interacts with the other components through mutual inductive or electrostatic interactions.
Comparison with other reduction approaches in the literature and self-consistency checks on the system spectrum allowed us to demonstrate the validity of our reduced model. At the same time, our approach is based on more general assumptions than other reduction methods in the literature. Therefore, in the case of isolated superconducting qubits we have seen that choosing the local computational basis with explicit reference to the measurement operator improves the accuracy of the reduced Hamiltonian, in terms of both the spectrum and expectation values of circuit operators. This is especially true for qubit designs with reduced anharmonicity, such as the capacitively-shunted flux qubit. In the multiple-qubit case, the Schrieffer-Wolff transformation theory provided the basis for calculating the effective spin Hamiltonian, the only requirement for its application being that the size of the spectral gap of the un-perturbed Hamiltonian should be larger than the size of the interaction. In principle this limitation can be circumvented, as long as one is able to partition the system in smaller units, and as long as the qubits in each unit display sufficient anharmonicity. Numerical calculations of the effective multiple-qubit Hamiltonians provided results in good agreement with the existing reduction methods, when these were used within their range of applicability.
This reduction method should prove useful in different areas of applied quantum computation, where complex systems of continuous variable circuits are described in terms of interacting two-level systems. In practice one could start by fitting the parameters in the circuit model to some preliminary data, then extract the effective qubit Hamiltonian as a function of the control biases. The reduced model could then be verified with additional experiments, for instance spectroscopic or state population oscillation measurements, and successively be employed as the reference model for the operation of the system [16]. In the context of circuit design this method can be used to model the interplay between different qubit Hamiltonian terms, for instance the effect of the coupler bias on the qubit transverse fields [12] (i.e. dynamic inductive loading), or to predict the size of non-Ising terms like non-stoquastic or many-body interactions (as well as of Ising terms like the transverse fields, beyond the instanton approximation).

A.1 Capacitance and inverse inductance matrices
In this appendix we give the definition of the capacitance and inverse inductance matrices used to specify the linear part of the circuit HamiltonianĤ LC . For a circuit with N nodes (ground node excluded), these are two symmetric N × N matrices. In the capacitance matrix, each diagonal element (C) ii represents the sum of the capacitances connected to the i -th node, while, for every pair of nodes i = j, the off-diagonal element (C) ij equals minus the total capacitance between i and j. For the circuit in figure 5, for instance, the capacitance matrix is whose inverse is where C = C JT + C sh . Notice that 1/(C −1 ) ii corresponds to the effective capacitance between node i and ground.
In analogy with C, the inverse inductance matrix L −1 has, along the diagonal, the sums of the inverse inductances connected to each node and, in the off-diagonal elements, the total inverse inductance between pairs of nodes. The inverse inductance matrix for the circuit in figure 5 is, for instance,

A.2 Capacitance and inverse inductance matrices: interacting circuits case
In this appendix we show how to modify the capacitance and inverse inductance matrices of two circuits in order to take into account their interactions. The following definitions can easily be extended to the case of more than two interacting circuits. Let C 1 and C 2 be the two original capacitance matrices of the two circuits (as defined in appendix A.1), and let their sizes be N × N and M × M , respectively. Let C 12 be the N × M matrix whose elements are the capacitances between pairs of nodes belonging to different circuits. Consider then the following (N +M )×(N +M ) matrix: where the primed matrices include the additional capacitance attached to each node, i.e.: Notice that C is nothing but the capacitance matrix defined for the extended circuit including all the nodes of the two interacting circuits. By inverting it, we get: where C −1 1 and C −1 2 are the new inverse capacitance matrices of the two circuits (cf. Eq. (13)) which include the effect of the external capacitive loading, and C −1 m is the inverse mutual capacitance matrix, describing the interaction between the two circuits, which appears in equation (14).
For the inductive interactions, these involve pairs of inductive branches belonging to different circuits, coupled by their mutual inductance. Let N and M be the number of branches in the two circuits and consider the following (N + M ) × (N + M ) matrix: where L bi is the inductance matrix of circuit i in the branch representation, having along the diagonal the self-inductance of each branch (L bi ) kk = L bi k and zeros everywhere else, and M is the N × M matrix whose elements are the mutual inductances between pairs of inductive branches. Inverting L b , we obtain where M −1 is the matrix appearing in equation (16). L −1 b1 and L −1 b2 can be used to rescale the inverse inductance matrices of the two circuits (see Eq. (15)). This is accomplished by replacing each branch inductance L bi k appearing in the expression of L −1 i with 1/(L −1 bi ) kk .

A.3 Spectrum convergence
In this section we consider the convergence of the numerical spectrum of a qubit circuit as a function of the number of states included in the basis used to describe each of its modes. We refer to this number as the (mode) truncation.
The circuit examined here is that of the capacitivelyshunted flux qubit shown in Fig.5. By inspecting its circuit Hamiltonian, we find that the mode associated with node 1 (O1 ) is conveniently expressed in a basis of harmonic oscillator states below a certain occupation number N max O , while those associated with nodes 2 (C1) and 3 (C2) are better expressed in the charge number basis, keeping only integer charges lower in absolute value than N max C1 (N max C2 ) [37,33,38]. Figure 11 shows the lowest 20 eigenvalues of the approximate circuit Hamiltonian H (N ) e.m. , as a function of its linear size N = (N max O1 + 1) · (2N max C1 + 1) · (2N max C2 + 1), as well as the time required to evaluate them (shown by the pink line and indicated on the right vertical axis). The qubit is taken to be biased at the optimal point f z = Φ ext 23 /Φ 0 = 0.5 and its other physical parameters are given in section 4.1.1 of the main text. In the graph the values of the truncations N max O1 , N max C1 and N max C2 are increased sequentially going from left to right, starting from the values (N max O1 , N max C1 , N max C2 ) = (2, 3, 3). As we can see, all of the 20 lowest eigenvalues have converged for the set of truncations (9, 10, 10), corresponding to a Hamiltonian of linear size N = 4410. As it turns out, the convergence is mainly determined by the Josephson modes, and the set (3, 10, 10) (N = 1323) is already sufficient to obtain the same eigenvalues. Also notice that the lowest three eigenvalues already converge for the set of truncations (3,5,5) and N = 363.
The eigenvalue evaluation times refer to the use of MATLAB c eigs algorithm [35], run on a quad-core laptop CPU. As the pink line in the graph shows, the run time scales as a power law of the linear matrix size (notice the log-log scale), namely t run (1.1 · 10 −5 s) · N 1.4 , as results from a non-linear fit.

A.4 Tunnelling rates in the rf-SQUID qubit with the instanton method
The semi-classical description of tunnelling through a potential barrier is a very well-known subject in quantum mechanics and is routinely used in many applications of chemistry and quantum physics [31,43,44]. In order to describe the tunnelling between the two opposite persistent current states of the rf-SQUID qubit, we are going to use the formalism developed in [45], which applies to a generic, potentially asymmetric double-well potential. Let us first write the semi-classical potential of the circuit [4]: where ϕ = 2πΦ/Φ 0 is the dimensionless total flux, ϕ ext is the externally applied flux, U L = (Φ 0 /2π) 2 /L is the characteristic inductive energy and β L = E J L(2π/Φ 0 ) 2 is called the screening parameter. When β L 1 and ϕ ext /2π 0.5, this potential has three stationary points, given by the solutions of the transcendental equation Two of the solutions, say ϕ L and ϕ R , correspond to the minima of the left and right potential wells, respectively, while the third, ϕ M , is the maximum of the barrier between them (ϕ L < ϕ M < ϕ R ). For instance, when U L = 65GHz, β L = 1.9 and ϕ ext /2π = 0.49, we obtain the potential profile shown in figure 12 (solid black line, sitting below the dashed lines). According to the semi-classical theory, the low energy behaviour of the rf-SQUID system can be described in terms of the tunnelling between the lowest bound states in its two potential wells, Ψ L (ϕ) and Ψ R (ϕ) [31]. Figure 12: rf-SQUID semi-classical potential (black line) for ϕ ext /2π = 0.49, and its two symmetrised versions (dashed lines). Also shown are the energies of the lowest bound states in the two wells.
These represent the local solutions to the stationary Schrödinger equation, in the limit where the two wells are completely isolated from each other (eg. ϕ L ϕ R ). One way to approximately identify these solutions is by considering the second-order series expansion of the potential around its minima: Then Ψ L (ϕ) and Ψ R (ϕ) approximately correspond to the vacuum states of two displaced harmonic oscillators, such that with C the total capacitance across the Josephson junction, and The oscillator frequency here is Notice that these states have a phase expectation value of φ i = ϕ i and an average persistent current of Therefore, since ϕ L < ϕ ext < ϕ R , the bound states also correspond to persistent current states of opposite sign, as expected. Quantum tunnelling across the potential barrier couples the two bound states, leading to the repulsion between their energy levels. The resulting eigenstates of the system are determined by the following two-level Hamiltonian, expressed in the persistent current basis {|Ψ R , |Ψ L }: where ∆ is the tunnelling energy. This represents the effective qubit Hamiltonian of the circuit, and is again in the standard form of Eq. (26). Finally, following reference [45], we can write the tunnelling energy explicitly as: where A = 1 2 with V 0 = V (ϕ M ), and where ∆ L,R is the tunnelling energy relative to the symmetric double-wells V L (ϕ) and V R (ϕ), obtained by reflecting V (ϕ) about the local maximum ϕ M (cf. dashed lines in figure 12): V L (ϕ) = V (min(ϕ, 2ϕ M − ϕ)), The instanton result for the symmetric double-well tunnelling energies reads: with S i the tunnelling action, given by: where ϕ i,1 = 2ϕ M − ϕ i,2 are the two points at which the potential barrier intersects the energy level: V i (ϕ i,1 ) = V i (ϕ i,2 ) = E i . This semi-classical formula holds when S i and therefore in the limit of small tunnelling energies [32].