Native three-body interaction in superconducting circuits

We show how a superconducting circuit consisting of three identical, non-linear oscillators in series considered in terms of its electrical modes can implement a strong, native three-body interaction among qubits. Because of strong interactions, part of the qubit-subspace is coupled to higher levels. The remaining qubit states can be used to implement a restricted Fredkin gate, which in turn implements a CNOT-gate or a spin transistor. Including non-symmetric contributions from couplings to ground and external control we alter the circuit slightly to compensate, and find average fidelities for our implementation of the above gates above 99.5% with operation times on the order of a nanosecond. Additionally we show how to analytically include all orders of the cosine contributions from Josephson junctions to the Hamiltonian of a superconducting circuit.

Here we show how the electrical modes of an idealized, symmetric superconducting circuit can implement a native three-body interaction among qubits. We call these three qubits input, output and control. This is similar to how the electrical modes are used to implement desired interactions in Refs. [25][26][27], and in Ref. [28] a Schrieffer-Wolff transformation is used similarly. Examples of quantum computation using three-body interactions can be seen in [29,30]. Being native, the couplings are strong and do not require any ancillary qubits. This contrasts perturbative interactions, which by their very nature are generally weak, and decomposition schemes, where ancilla qubits are necessary. Because the couplings are strong in our set up, interactions between the qubit subspace and higher excited states of the anharmonic oscillators are not completely suppressed. Two of the eight qubit states immediately evolve out of the qubit subspace. We therefore effectively work in a six-state Hilbert space. We refer to the degrees of freedom as qubits to give a clearer intuitive picture and underline the intention of the system as a component in quantum technology or a quantum simulator. The remaining six states can be used to implement a CNOT-gate or a spin transistor, where the control is a quantum degree of freedom. Including non-symmetric contributions from couplings to external control, readout and ground, thus making the circuit more realistic, undesired detuning and couplings are induced, and the circuit must be altered slightly to compensate. After this alteration there still remains an undesired XX-coupling, resulting in uncontrolled hopping between input and output. Even so the average fidelity of our implementation of the above gates is nonetheless above 99.5% with operation times on the order of a nanosecond. These extremely fast operation times are also a consequence of the strong couplings. This time scale is comparable to recent experimental results where a two-qubit gate with an operation time of less than a nanosecond was realised with electron spin qubits bound to phosphorous donors in silicon [31]. Our numerical results for a (restricted) three-qubit gate in the context of superconducting circuits are thus comparable to this two-qubit gate in electron spins, which to the best of our knowledge is one of the fastest realised two-qubit gates.
Furthermore, in this paper we show how to include all orders of the cosine contributions from Josephson junctions to the qubit. We do this with two levels for each qubit in the main text but in the Appendix we show how the same method can be used to include more levels. This is useful for analytical study of suppression of couplings to higher levels and for faster simulation of multiple levels.
This article is organized as follows: In Section II A we arXiv:1910.08741v2 [cond-mat.supr-con] 22 Nov 2019 go through the Hamiltonian for the idealized circuit. We discuss what the resulting gate is in Section II B. Then in Section III A we move on to the more realistic circuit by introducing couplings to external control and ground. The alteration to the circuit and the subsequent Hamiltonian is discussed in Section III B. In Section IV we go through the simulation of the system. Spin model parameters are presented in Section IV A, the state population evolution in Section IV A, and calculations of average fidelity is presented in Section IV C. Finally in Section V we discuss results and conclude the paper in Section VI.

II. IDEALIZED SYSTEM
We consider a very symmetric circuit consisting essentially of three transmon [32] blocks in series with identical parameters. This circuit is unrealistic in the sense that it is not connected to external control and readout, or to ground. These external couplings induce an asymmetry in the system, which we will compensate for in the next section. Changing basis to the eigenmodes of the capacitance network, the circuit implements three qubits interacting via an XXC-coupling which flips the input and output states conditioned on the state of the control qubit.

A. Circuit & Hamiltonian
The circuit we consider can be seen in black in Fig. 1. It is similar to that of [25]. We relate a flux coordinate, ψ j for j = 1, 2, 3, 4, to each node in the circuit [33]. Gathering these in a vector ψ, we change coordinates to the eigenmodes of the capacitance network φ = (φ CM , φ 1 is a matrix whose columns are the eigenvectors of the capacitance matrix In these coordinates the capacitance matrix becomes diagonal, thus avoiding interactions induced by the capacitances. These eigenmodes or electrical modes are associated with circuit charge oscillations between the superconducting islands of the circuit or groups of these [25]. For example φ c corresponds to charge oscillating between the islands of ψ 1 and ψ 4 and those of ψ 2 and ψ 3 . That is, in this mode charge oscillates between the two ends of circuit and the middle of it. All our modes are quadrupolar as they involve all four of the original coordinates and thus all four of the superconducting islands. In other circuits there might be dipolar modes, involving only two of the original flux coordinates, or higher than fourth order poles, involving more than four coordinates. All but φ CM will be anharmonic oscillators and give rise to transmonlike qubits [32]. The Hamiltonian in these coordinates is (see Appendix 1) where the p i , i = 1, 2, c, are the momenta conjugate to φ i , and E C = 1/8C is the charging energy of a capacitor with capacitance C (in units = 2e = 1). The mode φ CM does not enter at all and corresponds to a center of mass mode, and so we ignore it from this point on. Recasting the modes in terms of harmonic oscillators and truncating each to its two lowest levels, we can express this Hamiltonian exactly in terms of Pauli matrices. Usually one first expands the cosines to fourth order in the fluxes when working with transmon qubits [32], but by writing the cosines in terms of complex exponentials we can derive the following identity (see Appendix 2), where (·) 2 indicates truncation to two levels, a † /a are bosonic creation/annihilation operators, k is a real number, and σ x,y,z are the Pauli matrices. For a four level version of this identity see Appendix 7. With this identity the contributions from the cosines to the two-level Hamiltonian can be calculated exactly, i.e. without any Taylor approximation. The kinetic terms are easily truncated to two levels with standard methods. We find the following qubit Hamiltonian where the Ω's and J's are positive spin model parameters.
Their explicit form can be found in Appendix 2. The circuit has thus resulted in two resonant qubits and a detuned third, coupled via pairwise ZZ-couplings, a threebody ZZC-coupling and a three-body XXC-coupling. C refers to control and to the operator (1 − σ z )/2, whose eigenvectors are the ground and excited state of the qubit with eigenvalues 0 and 1 respectively. This shows that this type of coupling is turned off and on depending on whether the c-qubit is in its ground or excited state. The In black: the idealized, symmetric circuit without external coupling. In red: external couplings to control, readout and ground represented by one joint capacitance for each node. In blue: the additional inductors with parasitic capacitance, introduced to counter the asymmetry of the external couplings.
ZZC-coupling is small, being sixth order in the fluxes and would thus not have survived a usual expansion to fourth order. Notice how there are no other XX-couplings except the controlled one. The couplings from the capacitances are gone because of the choice of coordinates, and while the Josephson junctions do generally contribute such couplings, they have cancelled because of the symmetry of the circuit parameters. This shows the power of our method. The desired three-body interaction is native in these coordinates and there are no undesired couplings, which we need to deal with. The Z-type couplings merely shift the energy levels of the system, which can be compensated. As we shall see, the XXC-coupling has a strength comparable to the anharmonicities as a consequence of these two quantities arising from the same circuit elements and the same type of terms in the Hamiltonian. More explicitly, it is the same circuit parameters that determine the strength of both these two quantities and they both arise from terms that are fourth order in the node fluxes, meaning they will naturally tend to be similar in size. Hence, the coupling is additionally quite strong, making the time scale very fast, and certain couplings with states outside the qubit subspace will not be suppressed by the anharmonicities. As mentioned we nonetheless consider the degrees of freedom as qubits, though it is important to keep in mind that some of the qubit states are "forbidden", as we shall see, if we want to the system to remain in the qubit subspace. Naming qubits 1 and 2 input and output, and the c-qubit control, we see that in a rotating wave approximation the effect of the XXC-coupling is that the input and output will flip each other conditioned on the state of the control qubit. The ZZ-and ZZC-couplings effectively only induce shifts of the qubit energies depending on their states, and since the Hamiltonian is symmetric under the exchange of qubits 1 and 2, the states |0 1 1 2 and |1 1 0 2 for any state of the c-qubit will still be resonant. These are the states which will be flipped into one another by the XXC-coupling.

B. The implemented gate
The qubit transition frequencies are about two orders of magnitude larger than the coupling strengths, as we shall detail later. Therefore interactions that do not preserve the number of excitations will be highly suppressed. In the same way we found exact expressions for spin model parameters in H 2 , we can find exact expressions for the anharmonicities, α i , which can be seen in Appendix 3. To fourth order they have the following simple relations to each other and the coupling strengths Hence, the coupling strengths are of same order of magnitude as the anharmonicities, and so interactions with higher energy levels that preserve the number of excitations will not be suppressed, for example |1 1 1 2 1 c ↔ |0 1 2 2 1 c . To justify reducing our local degrees of freedom to two-level qubits, we must make sure to not populate the states which couple to outside of the qubit subspace. These "forbidden" states are |1 1 1 2 0 c and |1 1 1 2 1 c , i.e. those where both input and output qubits are excited.
States where at most one of qubits 1 and 2 is excited are protected from interaction with higher states by the conservation of excitations and the fact that the c-qubit only changes its excitation in even numbers. This is a consequence of φ c only entering in even powers in the Hamiltonian Eq. (3). Therefore, if we initialize in a state where there is at most one excitation among qubits 1 and 2, we will never populate the |1 1 1 2 states. Hence, we may consider our system to implement a gate operating on a six-state Hilbert space consisting of the states |0 1 0 2 , |0 1 1 2 and |1 1 0 2 for any state of the c-qubit. The matrix of this gate for the c-qubit in its ground state and in the excited state can be written as These two correspond to the interaction being off or on, and likewise we may consider the gate to be closed or open. We have included phases β, γ and δ (and ignored an irrelevant overall phase), which are a result of different dynamical phases on the different states. These phases depend on the exact energies of the system and therefore on the exact circuit parameters. This means they can be predicted by simulation and even tuned to some extent. Furthermore, working in an appropriately rotating frame, these phases would be zero. From these matrices we see we can implement a CNOT-gate using the states |0 1 1 2 and |1 1 0 2 as our logical |0 and |1 . We can also implement a transistor which transfers quantum information with a gate that is conditioned on our control qubit. This has been discussed recently in Refs. [16,34,35]. While [34] finds that a minimum of four qubits is required for transistor functionality, the threebody interactions in our system reduces this requirement ot just input, output and gate qubits. The transfer process can also be used to generate entanglement between input and output in the case where the control is in a superposition state. In a frame where the phases vanish, this gate is similar to a Fredkin gate [19], also called CSWAP, except we must restrict the initial state of the output qubit to always be |0 2 . We call the resulting gate a restricted Fredkin gate. In Ref. [36] they show a decomposition of the Fredkin gate into five two-qubit gate, which has been proven to be the minimal number of gates necessary [37]. The decomposition of the gate can be seen in Fig. 2. We have introduced the two-qubit gate W, which is the conjunction of a controlled-NOT and controlled-V, where V satisfies V = X 2 , as indicated in the figure. This shows the advantage of using a direct implementation of a three-qubit gate. An experimental implementation using no fewer than five gates would likely induce more noise, would have a longer operation time, and would simply require more hardware. Furthermore, the W gate is not a standard two-qubit gate and would require additional work to implement directly. By rethinking quantum circuits and algorithms in terms of many-qubit gates, we can greatly decrease complexity and the implementational issues that follow. However, it should be noted that it has been shown that some multi-qubit gates can have shorter operation times, when decomposed, than the sum of the operation times of their components [38]. Thus, though it appears that the decomposed Fredkin gate might have an operation time around five times that of a standard two-qubit gate, like the CNOT-gate, it could be significantly less. Nonetheless, it seems intuitive that a direct multi-qubit gate would be superior, and this possibility demands study.

A. External coupling
We now add capacitive couplings to external control and ground. The external control could consist of resonators for readout, flux lines for tuning, and drive lines for preparation [39]. Each external coupling adds to the diagonal of the capacitance matrix, and we assume that the total contribution for each node is identical and equal to K. The external couplings are drawn in red in Fig. 1 as a single capacitor from each node coupling to outside of the circuit. Hence, we assume that these external couplings are identical for each node in the circuit. If they are not identical it will induce interactions between the modes after the coordinate transformation. The interaction strengths will be determined by the discrepancies between the capacitances of the external couplings. However, whether or not they are identical, we can drive each of our modes individually by essentially performing the same transformation of the drives as the coordinates. The general term in the Hamiltonian for some external capacitive coupling will have the form −K j D jψj , where K j is the capacitance of the coupling and D j is the externally controlled variable. This variable could be a voltage, flux or current. Gathering these variables in a vector D ψ , we can write all couplings together as where K D is a diagonal matrix with K j as entries. If these are all identical and equal to K, and we write D ψ = V D φ for some other drives D φ , then after the coordinate transformation we have Hence, the drives D φ couple individually to the modes φ. Even if the K j are not identical we can choose D ψ to achieve the same individual coupling, (see Appendix 5), but we will as mentioned have additional XX-couplings between the modes.
Our gate U is, ignoring the phases, a Fredkin gate, except we must restrict the output qubit to always be initialized in the ground state |02 , while the other two qubits can be in any state |φ1ψc . The Fredkin gate is decomposed into five two-qubit gates [36], where the W-gate is the conjunction of a controlled-NOT and a controlled-V, where V is the square root of a Pauli X-gate.

B. Circuit & Hamiltonian
These couplings result in an asymmetry between φ 1 and φ 2 , because the entries in the capacitance matrix (after the coordinate transformation) pertaining to these two modes will no longer be identical (see Appendix 6). This creates a detuning and the XX-coupling between them no longer cancels completely. While this coupling is small, the detuning is not negligible even for K 1 fF. Neither can be tuned to zero without removing the external couplings again. To counter this we add an additional branch to our circuit connecting the two ends via an inductor with inductance L 0 and a parasitic capacitance C 0 . To cancel XX-couplings from this new inductor and to ensure that the capacitance network retains the same modes, we must furthermore add an identical inductor to the middle transmon block. This is drawn in blue in Fig. 1. The circuit is now similar to that of [26], with some non-linear inductors replaced with linear ones. Notice, that these new inductors do not add any new dynamical degrees of freedom to the physical system but simply alter the parameters of the existing qubits. So while they do add to the complexity of the circuit they should not be considered an overhead or as ancillary, but as an equal and integrated part of the altered circuit. The new capacitances actually add to the asymmetric terms and should therefore be as small as possible. The contributions from the inductors are also asymmetric, and can be used to tune the states |0 1 1 2 1 c and |1 1 0 2 1 c back into resonance, balancing out the asymmetry from the capacitances. The undesired XX-coupling cannot be removed completely, though it can be tuned to be relatively small. We will therefore not be able to completely turn off the XX-coupling between the input and output qubits, in other words we can not completely close the gate. The leakage through the closed gate, when the control qubit in its ground state and the XXC-coupling is thus off, is however quite small as we shall see.

The resulting Hamiltonian is
where C −1 i are matrix elements of the inverted capacitance matrix (which are no longer as simple as in Eq. (3)). Their exact expression can be found in Appendix 6. E L0 = 1 2L0 is the inductive energy of the newly introduced inductors. We can see that there is no interaction pertaining to E L0 , as this has cancelled between the two inductors as intended. However, when we truncate the system to three qubits, there is still an uncontrolled XX-coupling between input and output, because of the introduced asymmetry. The qubit Hamiltonian is which is the same as with the idealized circuit except for the additional XX-coupling and the fact that the explicit symmetry between qubits 1 and 2 is broken. The spin model parameters are calculated analytically in Appendix 6. Hence, we expect the same dynamics as before, except the closed case will not be completely closed. Furthermore, we must tune the circuit parameters to make |0 1 1 2 1 c and |1 1 0 2 1 c resonant, compensating partially for non-identical qubit frequencies and for the Z-type couplings which shift the energy levels.

IV. SIMULATION
Let us now move on to numerical calculations and simulations of this circuit. We first show an example of the most important spin model parameters, showing how the couplings are comparable and even larger than the anharmonicities, resulting in extremely fast dynamics. Afterwards we show how the population of the qubit states evolve depending on the initial state, which show that states with zero or one excitation are essentially stationary, |0 1 1 2 1 c and |1 1 0 2 1 c exchange their populations as we want, and |1 1 1 2 0 c and |1 1 1 2 1 c quickly evolve out of the qubit subspace. Finally, we use the average fidelity [40] as a measure of the quality of our implementation of the gate described in Eq. (7). With the altered circuit, we can however not expect the phases acquired by |0 1 1 2 0 c and |1 1 0 2 0 c to be identical. Therefore the two β's in Eq. (7) are distinguished from each other by calling them β 1 and β 2 respectively. Simulations include noise by using the Lindblad master equation where ρ is the density matrix, L k are Lindblad operators representing different noise channels with rates γ k , and the sum runs over all Lindblad operators. The Lindblad operators are |ψ ψ| for each state |ψ in the Hilbert space inducing dephasing, and annihilation operators b i for each anharmonic oscillator φ i inducing photon loss. The decoherence rates are set to the same value for both types of noise and for all states and qubits, γ k = 0.05 MHz, giving them a lifetime of 20 µs, which is well within state-of-art achievements [39,41]. We simulate four levels for each non-linear oscillator to include the (suppressed) interactions with higher levels. In addition to causing two of the qubit states to evolve out of the qubit subspace, these interactions cause effective shifts of the energies of the qubit states. This is similar to the effects of higher levels considered in Ref. [42]. These shifts can essentially be considered as additional Z-type couplings and can be compensated for as mentioned previously. All simulations make use of the Python toolbox QuTiP [43].

A. Spin model parameters comparison
In our numerical search for circuit parameters, we fix C 0 and K at 1 fF, and vary C, E, L 0 to satisfy the following conditions: the states |0 1 1 2 1 c and |1 1 0 2 1 c must be resonant; the anharmonicities must be greater than 200×2πMHz; the ratio between the kinetic and harmonic potential terms for each qubit must be small to ensure that we are in the transmon regime (see Appendix 2); the ratio J x 12 /J x 12c must be small. The detuning between |0 1 1 2 1 c and |1 1 0 2 1 c is at first calculated from the analytical expressions for their energies, and for higher precision it is then calculated from the infidelity of the process The most important spin model parameters. The ZZ-and ZZC-coupling strength in blue, the XX-coupling strength in red, the XXC-coupling strength in purple and the anharmonicities in grey. The qubit transition frequencies Ωi/2π are all close to 36 GHz.
The analytical expressions do not completely capture the detuning because of higher order energy shifts from the strong ZZ-couplings and from interactions with higher levels as mentioned. An example of the spin model parameters can be seen in Fig. 3. In this simulation we have E/2π = 200 GHz, C = 21.67 fF and L 0 = 12.67 nH. We require such a large value of E/2π in order to make the contributions from the linear inductors relatively small, as these must balance out the asymmetric contributions from the capacitors, but not become dominant themselves. The size of the asymmetric inductive contribution is given by the ratio 2E L0 /E, while (2C 0 + K)/2C determines the size of the asymmetric capacitive contribution (see Appendix 6). For the parameters given here we have E L0 /E = 0.065 and (2C 0 +K)/2C = 0.069. We can see that as expected J z 12c , which is sixth order in the fluxes, is much smaller than the others. However, both the ZZ-couplings and in particular J x 12c are comparable to the anharmonicities. This large value of J x 12c gives us very fast dynamics, but clearly the anharmonicities are not large enough to suppress interactions that have coupling strengths like these. Lastly, J x 12 is small, but unfortunately not negligible, which will reduce the average fidelity of the gate, as the closed case, where the control qubit is in the state |0 c , will not be entirely closed. In this example the three qubit transition frequencies Ω i /2π are all close to 36 GHz, so interactions that do not preserve the number of excitations are heavily suppressed as expected. As mentioned, the ZZ-and XXC-couplings are of the same order of magnitude as the anharmonicities, because they arise from the same terms of the original flux Hamiltonian Eq. (10), i.e. fourth and higher order terms from the cosines. As we must tune the anharmonicities to be large enough so that we can experimentally address the individual qubits for preparation and so on, the couplings also become large. We see that the states with zero or one excitation are nearly stationary, while |11021c and |01121c exchange their populations. Finally the states with both input and output excited quickly evolve out of the qubit subspace.

B. Population evolution
We simulate the dynamics of the system with the parameters presented in the previous section. In Fig. 4 can be seen the populations of all states as we initialise the system in one of the eight qubit states. Clearly |0 1 0 2 0 c , |1 1 0 2 0 c , |0 1 1 2 0 c and |0 1 0 2 1 c are nearly stationary, i.e. their populations do not change, though they may acquire a phase. The states |0 1 1 2 1 c and |1 1 0 2 1 c exchange their populations, as the XXC-coupling is on, corresponding to the input and output qubits flipping each other. The very large value of J x 12c has resulted in dynamics on the nanosecond scale. Lastly, |1 1 1 2 0 c and |1 1 1 2 1 c quickly evolve out of the qubit space and excite a mess of higher states. Simulating the system with the same circuit parameters but including six levels for each anharmonic oscillator makes no significant difference, indicating that we have captured the effects of the higher levels by taking just four levels into account. Additionally, it can be seen that we have properly compensated for the effective energy shifts of the energy levels caused by higher levels and the strong couplings in the system. Finally, we find in our simulations that the qubit subspace loses 0.6% or less of its population due to leakage to the higher levels, except of course when initializing in the states |1 1 1 2 0 c or |1 1 1 2 1 c , in which case most of the population goes to the higher levels after some time.

C. Average fidelity
The average fidelity can be calculated from [40] F where U is the target gate, and E is the quantum map used to implement it. The U j are a unitary basis for operators on the relevant space of states, which has dimension d. In our case the target gate is the one described in Eq. (7) (with β replaced by β 1 and β 2 ), and E(ρ(0)) = ρ(t) is the time evolution of our system. The relevant space of states is spanned by the six qubit states that do not evolve out of the qubit subspace. Thus for any set of circuit parameters and phases β 1 , β 2 , γ and δ, we can calculate the average fidelity as function of time. We choose the phases to maximise the fidelity and call the time of maximum fidelity the gate operation time T g . This time is expected to be approximately π 2(J x 12c +J x 12 ) , as we are essentially considering the oscillation between two resonant energy levels interacting with a coupling strength J x 12c + J x 12 (see Appendix 2). However, higher order interactions make T g slightly smaller than this, and so we determine it numerically. We do this for a range of circuit parameters, choosing values for E/2π in [100, 400]GHz, and tuning C and L 0 to satisfy the demands described above, in particular maximising the fidelity of the process |0 1 1 2 1 c ↔ |1 1 0 2 1 c . The results can be seen in Fig. 5. These are simulated without noise. The fidelities are above 99.5% and vary only over a short interval, increasing slightly with the value of E. There are small, noisy variations in the curve, attributed to numerical imprecision and randomness in the optimization routine. In particular the calculations are affected by there being many local minima close together in the space of phases, and in simulating dynamics we have only a finite resolution in time. The gate times are well below 1 ns and increase slightly with E. Finally the angles also vary slightly with E.

V. DISCUSSION
The fundamental reason that two of the qubit states are lost, due to interactions with higher levels, is that we have an odd interaction induced by the same circuit components that create the anharmonicities. By odd interaction we mean an interaction stemming from product terms of odd powers of different flux coordinates, like φ 1 φ 2 or φ 3 1 φ 2 , which exchange odd numbers of excitations. Such interactions with higher levels must be suppressed by the anharmonicities. Since these interactions come from the same circuit elements as the anharmonicities (and both are fourth order in the flux coordinates), they inherently are comparable in strength, and so the anharmonicities can not suppress them. This is in opposition to even interactions with higher levels, i.e. interactions stemming from terms like φ 2 1 φ 2 2 , that exchange even numbers of excitations, and which are suppressed by the approximate conservation of the number of excitations. It would require a qubit to have for example two excita-tions in order for another qubit to be twice excited, but we always initialize in states with only zero or one excitation. We could therefore employ the above method of a symmetric circuit viewed in the basis of electrical modes more successfully if the odd interactions came from different elements than the anharmonicities, such that they could possibly be tuned independently, or if we only had even interactions, as is done in Refs. [26,27].
The gate we have proposed works on the scale of a single nanosecond, which is very fast compared to other present two-or multi-qubit gates. This may pose a challenge for current implementation, but we expect our implementation to be relevant in the future as other operation times are expected to decrease. We note in passing that the one-qubit operation times can be 10 ns or shorter [44,45] and our gate times are comparable. As mentioned, recently a two-qubit gate with operation time of less than a nanosecond was experimentally achieved in the context of electron spin qubits in atoms [31]. Our proposal predicts comparable results for superconducting circuits, though for a direct three-qubit gate working in a limited part of the qubit subspace. Furthermore, joining several fast gates together in a network might result in a useful component that takes advantage of the fast operations. Alternatively to a three-qubit gate, direct threebody interactions are useful for quantum simulation of lattice gauge theories. For example an XXX-coupling is needed for the hopping of particles between lattice sites mediated by a U(1) gauge field living on the link between sites [7], and the matter-gauge interaction in Z 2 gauge theories is an XXZ-coupling [46]. Our approach could be used to derive circuits which implements such interactions directly. Indeed, we found an XXZ-coupling in the presented circuit. Likewise, three-body interactions could be used to implement three-local constraints in the context of quantum annealing schemes. Here it is most often even interactions, like a ZZZ-coupling, which are used. Such an interaction could certainly also be found with our approach, indeed we did in fact find it in the above circuit and wrote it as a ZZC-coupling. It would, however, necessarily be small because it would be sixth order in the fluxes, just as we found J z 12c to be much smaller than the other coupling strengths. If one could figure out a way for other couplings to cancel or become irrelevant, it might still be a good approach and faster than indirect implementations.
Alternatively to the above work, one could interpret the system as a qutrit, consisting of the states |0 1 0 2 , |0 1 1 2 and |1 1 0 2 , interacting with the c-qubit. The energy levels |0 1 1 2 and |1 1 0 2 of the qutrit will be degenerate and, depending on the state of the c-qubit, will either be stationary or oscillate between each other. Or one could diagonalize the two higher levels of the qutrit with the states 1 √ 2 (|0 1 1 2 ± |1 1 0 2 ), which will again either be degenerate or have an energy splitting of 2J x 12c , depending on whether the XXC-coupling is off or on. Proposals have been made on how to utilize the third level of a qutrit to reduce the number of elementary gates neces-a) b) c) Figure 6. Three different approaches for scaling of the circuit (in black), always using identical capacitative couplings to either a copy of the circuit or a mediating resonator (in red). a) Using capacitors between every node and its counterpart in the copied circuit couples each mode in the first circuit to its counterpart in the copied circuit. b) Connecting the middle two nodes of the circuit to a mediating resonator couples the modes φCM and φc to the resonator. c) Connecting the middle two nodes of the circuit to the ends of a mediating resonator couples the modes φ1 and φ2 to the resonator. sary to implement concrete quantum computations [47][48][49][50][51][52][53][54][55][56][57]. However, it is also possible to employ a set up, which actively chooses to only make use of exactly the states we have available in our system. In Refs. [16,58] two physical qubits represent a site in a lattice, which can host a mobile logical qubit. The state |00 represents a vacant site, while the states |01 and |10 represent the two internal states of a present logical qubit. In this way the logical qubits are able to move around the lattice while carrying quantum information.

A. Scalability
Finally, we consider scaling of the system by looking at the idealized circuit, and how it might be coupled to copies of itself. The challenge that needs to be addressed in regards to scaling is the fact that the flux coordinates φ are all a linear combination of all of the original node fluxes ψ. Hence, individual coupling any of these degrees of freedom, and the corresponding qubits, to external qubits will require coupling to all four nodes in our circuit. However, connecting pairs of nodes to external circuits can be used to couple pairs of the modes to the external degrees of freedom.
In Fig. 6 can be seen three approaches to scaling the idealized circuit by connecting it to a copy of itself, either directly or through a mediating resonator. The first approach is to connect each node of the circuit with identical capacitors to the corresponding nodes of the copy. This results in XX-couplings between each of the modes and their corresponding partners in the copy. The second approach connects the two middles nodes with identical capacitors to a resonator. This results in the modes φ CM and φ c coupling to the resonator. Likewise, the third approach connects the two middles nodes with identical capacitors to the two ends of a resonator. This results in the modes φ 1 and φ 2 coupling to the resonator. If might then be possible to tune the resonator and the couplings between it and the modes to mediate some interesting coupling between the main circuits.
Implementing scaling via a resonator, however, comes with the complication that the two middle nodes of the main circuits now get capacitative contributions which the end nodes do not. In order to retain the proper electrical modes of the main circuits, some capacitors would need to be added to the end nodes. Alternatively, one could also implement scaling by connecting the end nodes via resonators in the same way as the above examples which use the middle nodes. In the context of scaling the main circuit to a long chain of copies connected through mediating resonators, one might utilize an alternating scheme where the mediating resonators are alternatingly coupled to the middle nodes and to the end nodes, in order to ensure identical capacitative contributions to all nodes of the main circuits.

VI. CONCLUSION
In this article we have seen how a symmetric superconducting circuit of three transmon blocks in series can implement a direct three-body interaction between its electrical modes. The modes are non-linear oscillators and can be considered as qubits. Specifically, we found a strong controlled XX-coupling between two of these qubits conditioned on the state of the third. This interaction can be used to implement a CNOT-gate or a transistor among quantum degrees of freedom using what we call a restricted Fredkin gate. However, as this coupling is induced by the same circuit elements as the anharmonicities of the non-linear oscillators, the coupling strength is comparable to the anharmonicities. This means that some interactions with higher levels are not suppressed. Specifically, the two states where both the input and output qubits are excited will quickly evolve out of the qubit subspace. Furthermore, external couplings to control and ground induce an asymmetry in the system, which we compensate for by altering the circuit with two linear inductors. The XX-coupling is no longer completely controlled after this, but we still find average fidelities for our implementation of the restricted Fredkin gate above 99.5% with operation times below 1 ns. The extremely fast time scale is caused by the large coupling strengths.
In regard to scaling of the circuit the challenge that needs to be addressed is the fact that the flux coordinates φ are all a linear combination of all of the original node fluxes ψ. Thus, individual coupling of the qubits to external qubits require coupling to all four nodes in the circuit. This means that cross-talk may be involved and that addressing individual logical qubits in the system requires a careful tuning of couplings to different access points. We have given some examples of how scaling can be approached, either coupling the circuit directly or through a mediating resonator.

Circuit, Lagrangian and Hamiltonian
The idealized circuit we work with can be seen in black in Fig. 1 of the main text. To each node we associate a flux coordinate ψ j for j = 1, 2, 3, 4 and order these in a vector ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T . If we write up the capacitance matrix for the circuit we can write the kinetic terms of the Lagrangian as The potential part of the Lagrangian from the Josephson junctions is We now wish to change coordinates to the eigenmodes of the capacitance network, in order to remove interactions from the capacitances. This means we change to the coordinates where K becomes diagonal. The matrix used for this change of basis has the eigenvectors of K as its columns Notice that the eigenvectors are not normalized, their norms are (from left to right) 4, 4, 4(2 + √ 2) and 4(2 − √ 2). The eigenvalues of these are 0, 2C, (2 − √ 2)C and (2 + √ 2)C. If we name the new coordinates φ CM , φ 1 , φ 2 and φ c , and gather them in a vector φ, we can relate ψ and φ via Then the kinetic part of the Lagrangian becomes is the transformed capacitance matrix. Notice how the norms of the eigenvectors match with their corresponding eigenvalues in such a way that the non-zero entries are identical. The potential part of the Lagrangian becomes We see that φ CM does not enter at all in the Lagrangian. This mode is exactly like a center of mass-mode, or a constraint variable. It corresponds to a free particle that does not affect the rest of the system. As it does not enter at all, we simply ignore it, but without changing our notation (so for example now φ = (φ 1 , φ 2 , φ c ) T and C has lost its first row and column). We can now define momenta conjugate to the surviving coordinates With these we can write the kinetic Lagrangian in terms of the momenta where we have used that the inverse of C is symmetric, because C is itself symmetric. As C is diagonal, with all diagonal entries non-zero, it is trivially invertible. Incidentally C was not invertible before we removed φ CM . This is because the basis of node fluxes is overdetermined by exactly 1, i.e. there is one more coordinate than there are degrees of freedom. This is usually dealt with by naming one of the nodes the ground node, and setting the corresponding flux to zero. We have essentially done the same with φ CM , except it does not simply correspond to one of the nodes in the circuit, but in fact relates to all. As ground also has a more physical meaning, and it is difficult to avoid coupling to it, we will consider a small capacitive coupling to a ground node later on in the more realistic instance of this circuit. As mentioned in the main text this coupling introduces a kind of asymmetry in our system, which we must compensate for. Finally, with the Lagrangian fully transformed and written in terms of momenta, we can write up the Hamiltonian as where we have introduced the charging energy E C = 1 8C (in units = 2e = 1). Let us introduce H kin = L and H pot = −L pot , as we do not need to refer to the Lagrangian from this point on. If we consider the φ i to be positionlike coordinates, our system can be considered as three particles interacting through the cosine potentials. In that case we see that the diagonal elements of the capacitance matrix correspond to masses. This can be a helpful way to consider the system, and we will sometimes use language appropriate to it.

Recasting and truncating
We now move on to recast our Hamiltonian in terms of harmonic oscillators. The Hamiltonian of a harmonic oscillator is where φ is a position-like coordinate (in our case a flux) and p is its conjugate momentum. We promote φ and p to operators and give them the canonical commutation relation [φ, p] = i. We then define r = (K p /K φ ) 1/4 and introduce creation/annihilation operators a † /a related to φ and p via A brief calculation turns the Hamiltonian into the familiar form Let us do this for our Hamiltonian and all three flux coordinates. We promote all φ i and p i to operators and introduce for each of them creation/annihilation operators a † i /a i . As the Hamiltonian does not explicitly contain any φ 2 i -terms, we look at the second order terms from the expansion of the cosines in order to define r i . The second order terms alone are Notice, how the interaction terms cancel neatly. We will later see that all but the desired XXC-coupling and the ZZ-couplings cancel. With this we can define r i , which turn out to be identical in this very neat and symmetric case With this we can relate φ i and p i to their respective creation/annihilation operators as in the case of a single harmonic oscillator. Terms of the form a 2 and (a † ) 2 will cancel as in the case of a harmonic oscillator, but the fact that we are not actually working with harmonic oscillators is revealed by the presence of terms like a 4 and (a † ) 4 , which induce "rotations" between different energy levels, without any actual interaction taking place. This reveals that the harmonic oscillator basis is not the natural one for this system, but these rotations will be highly suppressed by a rotating wave-type argument. As we wish to be in the transmon regime, where the fluxes vary only slightly around the minimum of the potential, we will need to make r i small. This physically corresponds to having heavy particles in strong harmonic potentials, cf. the definition of the r i . In that case the particles move in such low depths of the potential that they can not perceive the cosine-behaviour of the potential, but simply see a harmonic potential with a slight anharmonicity. This is why the Fock basis works well for the system. More exactly the qubit transition frequencies will be much larger than the coupling strengths and so the number of excitations in the Fock basis will be approximately preserved, again by a rotating wave-type argument. The r i are exactly the parameters that define the sizes of the φ i , and so an expansion in φ i will in fact be an expansion in r i . The transmon regime is exactly when the r i are so small that an expansion to fourth order is justified. We will in this paper generally include all orders, but it is helpful to consider the r i as small to make the physics clearer. A usual value of the r i is < 1/6, and we will generally work only with even powers of them. Let us write the Hamiltonian in terms of the creation/annihilation operators We will now truncate this to the two lowest levels of each mode and represent the result via Pauli matrices. When we simulate the dynamics of the system we will use the full Hamiltonian above and include more than two levels for non-linear oscillator (normally four), in order to check whether we stay in the qubit space, etc. However, deriving the analytical Hamiltonian for two levels, will make it clearer what we can expect from our system and what we might be able to make it do, even though we will never actually use the two level Hamiltonian for any simulation. Before we reduce to two levels, let us introduce the idea of even and odd interactions. These are interactions that come from an even or odd number of creation or annihilation operators. Hence, XX-couplings are odd while ZZ-couplings are even, and when we include more than two levels there will be interactions that exchange 2, 3 or more excitations. We can see that the c-mode only has even interactions, because its only interaction terms come from a cosine which is an even function. The truncation is done with the usual method of calculating the matrix representation of an operator A by The truncation to the two lowest levels is justified by the anharmonicities (although as mentioned in the main text, these do turn out to be comparable in size to the coupling strengths), which we derive in the next section. The kinetic terms of the Hamiltonian are easiest to calculate and yield where matrix identities attached to numbers and other matrices are suppressed, i.e. any number added to a matrix is implicitly that number times the identities for each of the two level-systems I 1 ⊗ I 2 ⊗ I c , and for example σ z 1 is actually σ z 1 ⊗ I 2 ⊗ I c . To calculate the contributions from the cosines, we first derive the result quoted in the main text Eq. (4). We will write the cosines as complex exponential, and therefore need matrix elements of the form n|e ik(a † +a) |m (33) for some real constant k, and n, m = 0, 1. The operator is in fact a displacement operator, D(ξ) = e ξa † −ξ * a , with ξ = ik. We know of displacement operators that and we know its commutation relations, one of which is With this and the definition |n = 1 √ n! (a † ) n |0 , we can calculate the needed matrix elements. Incidentally, these identities are sufficient to calculate exactly the matrix elements of any product of sines and cosines of linear combinations of flux coordinates, with any number of energy levels. Let us also note that since the operator a † + a is symmetric, the operator D(ik) = e ik(a † +a) is also symmetric, i.e. n|e ik(a † +a) |m = m|e ik(a † +a) |n . We find With these we can write the operator e ik(a † +a) truncated to two levels in terms of Pauli matrices and the identity Let us now rewrite the potential part of our Hamiltonian in terms of exponentials and use the above identity to truncate it to two levels. First we do the following partial calculations, where we use the fact that the different a i commute to factorize the exponentials cos k(a † 1 + a 1 ) ± k(a † 2 + a 2 ) = 1 2 e ik(a † 1 +a1) e ±ik(a † 2 +a2) + e −ik(a † 1 +a1) e ∓ik(a † 2 +a2) cos k(a † c + a c ) = We can then rewrite the potential terms Defining the following constants we get the Hamiltonian as quoted in the main text. We see that there are no XX-type couplings other than the desired XXC-coupling.
The XXC-and ZZC-couplings could of course both be written as two independent couplings, a two-particle and a three-particle interaction. Writing the XX-type couplings as a single XXC-coupling, however, makes it clear that they disappear completely when the c-qubit is in the excited state. The XX-couplings are controlled by the state of the c-qubit. For the ZZC-coupling it is a bit more arbitrary, as there is an additional ZZ-coupling between qubits 1 and 2, i.e. even if we have the c-qubit in its ground state there is still a residual ZZ-coupling. The gate operation time, T g , defined as the time it takes for |0 1 1 2 1 c to evolve to |1 1 0 2 1 c and vice versa, can be easily calculated from considering a two level system. If we have two energy levels detuned by ∆ interacting with a coupling strength J, described by the Hamiltonian and initialise in one of the states, then it will evolve into the other in a time with a maximal amplitude Hence, as |0 1 1 2 1 c and |1 1 0 2 1 c can be considered two resonant levels (∆ = 0) that interact via J = J x 12c , we would expect the gate times to be given by the above formula, and the transfer to happen with an amplitude of e i π 2 , i.e. complete transfer of population with a phase of π 2 . Of course for a realistic system noise will decrease the amplitude and give it a different phase. Furthermore, the non-zero XX-coupling we find in the more realistic version of our circuit and higher order interactions induced by the strong couplings to higher levels will change the gate time and transfer amplitude.

Anharmonicities
We can calculate the anharmonicities exactly using the same method as above. The anharmonicity is defined as As there are multiple qubits in play here, there will be some contributions to the anharmonicities that come from interactions and are as such dependent on the state of the qubits. This is similar to how ZZ-couplings shift the energy levels depending on the states of the qubits. These contributions we ignore in the above definition, i.e. for this calculation we will only consider terms in the Hamiltonian that are not interactions. The contributions from interactions will always be smaller, as they are at least a factor r 2 smaller, and the non-interaction contribution are themselves at least of order r 4 . Remember that r is the small parameter and usually one expands to only fourth order, thus interaction contributions to the anharmonicities are usually ignored as well. In addition to the matrix elements found before, we will need To make sure we do not include contributions that come from ZZ-like couplings, we shall write the diagonal part of e ik(a † +a) as a 3x3 matrix in terms of the identity and the following two matrices The first is continuation of σ z and the second is the matrix that will carry the contribution to the anharmonicity (both the identity and Z do not contribute to the anharmonicity). With 2|e ik(a † +a) |2 calculated we can see that the diagonal part of e ik(a † +a) can be written as 1 − k 2 2 + k 2 2 Z + k 4 2 A e −k 2 /2 , which is consistent with the two level form of e ik(a † +a) we found previously. The idea is now that when calculating for example a 1 , and faced with product of exponentials, like e ik(a † 1 +a1) e ik(a † 2 +a2) , we shall take the second exponential to be 1 − k 2 2 times the identity, rather than for example just the identity alone. This factor is a result of the basis of operators we use (particularly the use of σ z ). For a different basis the anharmonicity would be differently distributed between interaction and non-interaction contributions. Let us finally note that we must necessarily have α 1 = α 2 because the Hamiltonian is symmetric under the exchange of the indices 1 and 2. Replacing exponentials with 1 − k 2 2 + k 2 2 Z + k 4 2 A e −k 2 /2 and expanding all products, we can simply sum the coefficients of the stand-alone A's for each mode. Doing this we get We can note here that Er 4 = E C 8 , revealing that the anharmonicities go as the charging energies, as usual with transmon qubits. Furthermore, we can see that the anharmonicities are very much comparable to the ZZ-and XXCcoupling strengths found above. Indeed, to fourth order in r, we have α 1 = α 2 = 3 4 α c and α c = J x qqc = 8 3 J z qq = 4J z qc . As mentioned in the main text, this means we cannot make use of all the available states in the qubit subspace, i.e. the space of just 0 or 1 excitation in each mode, as some of them couple strongly outside of the qubit subspace.

Gate
As the anharmonicities are not large enough to suppress interactions with higher levels we must restrict the available states further than just the qubit subspace. In particular the two states |1 1 1 2 0 c and |1 1 1 2 1 c with an excitation in both of qubits 1 and 2 couples non-negligibly to states out of the qubit subspace. In these states the excitations in qubits 1 and 2 can come together in just one of the modes, exciting it out of the qubits subspace. Similarly these two excitations can even go into the c-mode, via an XXZ-like interaction, because the c-mode is close to being resonant with the other two qubits. All other states, however, do not evolve out of the qubit subspace. In these states there are at most two excitations available and in that case one of them will be in the c-mode. As it is only coupled via even interactions this one excitation cannot leave the c-mode and so no mode can acquire both excitations.
This means that as long as we do not initialize the system with more than one excitation in total in qubits 1 and 2, we expect the state of the c-mode to be stationary. Hence, whether the XXC-coupling is turned on or off is completely determined by initialization and does not change during dynamics, barring noise effects. So working in the basis of |0 1 0 2 φ c , |1 1 0 2 φ c and |0 1 1 2 φ c with φ c = 0, 1, we may say that the system implements a gate represented by a diagonal matrix of phases when φ c = 0 and by when φ c = 1. The phases are a simple consequence of the fact that the different states have different energies, and so acquire different trivial dynamical phases. As mentioned in the main text, this gate could be used as a CNOT with |1 1 0 2 and |0 1 1 2 as the logical |0 and |1 states and the c-qubit as control. Naming one of the qubits 1 and 2 as input and the other as output, and always initializing with the output in the ground state, the gate can also implement a spin transistor that allows the transport of some superposition state from input to output conditioned on the state of the c-qubit.

External coupling
Although we are performing a change of coordinates we can still couple each of the new coordinates separately for external control. If we couple capacitively to each of the original node fluxes, the correspond terms in the Hamiltonian take the form for j = 1, 2, 3, 4, where K j are the capacitances of the couplings and D j are drives, i.e. essentially fluxes that we can control externally. Writing these capacitances in a diagonal matrix K D and the drives in a vector D ψ , we can write all the new terms as where we have ignored the terms that do not contain any node fluxes. After the coordinate transformation, ψ = V φ, we have The first term will induce undesired couplings between the modes if V T K D V is not diagonal. Calculations show that this happens only when K D has identical entries on its diagonal, i.e. all couplings to external control must be identical.
As for the second term, we would want it to be on the same form as before the coordinate transformation, such that we have a drive for each mode coupling to only that mode. If we introduce new drives D φ , and write D ψ = W D φ for some orthogonal matrix W, the second term in H D becomes We thus want W T K D V to be equal to some diagonal matrix C D with unit of capacitance, which it will be if we choose We have here used that V is orthogonal and both K D and C D are diagonal, which also means that K D is trivially invertible. Since the entries of K D ideally are identical, we can simply choose W = V (i.e. C D = K D ). Hence, identical couplings to the original coordinates with drives D ψ = V D φ , will result in identical couplings between the new coordinates φ and D φ . If the couplings are not identical we can still talk to each mode directly, but there will be induced XX-couplings between them. These will be small though, for small discrepancies.

Realistic circuit
Let us now go through all the same quantities but for the more realistic circuit, where we have included capacitive couplings to external control and readout and to ground. This circuit can be seen in Fig. 1 of the main text, where the loose ends pointing up collectively represent the external couplings to control, readout and ground. We assume that these are identical for each node, and contribute a capacitance K to each entry in the diagonal of the capacitance matrix. As mentioned this introduces an asymmetry between the modes, as we will go through in a moment, which we counter by adding a new branch to the circuit, connecting its two ends, consisting of an inductor with inductance L 0 and a parasitic capacitance C 0 . Furthermore to cancel interaction terms from this inductor we add an identical one to the middle block of the original circuit. To retain the same modes of the capacitance network, we must also add C 0 to the capacitance of the middle block. The capacitance matrix becomes This matrix has the same eigenvectors as before, i.e. the capacitance network has the same modes. If we change basis to φ, we get The center of mass-mode has now acquired a non-zero entry, but it still does not interact with any other mode, so we will still be ignoring it. More importantly the entries pertaining to φ 1 and φ 2 are no longer identical, even for C 0 = 0. This is the asymmetry that first and foremost creates a detuning between the modes and also makes the uncontrolled XX-coupling not cancel. Calculations show that even for C 0 = 0 and K = 1 fF qubits 1 and 2 are significantly detuned. This asymmetry, particularly the detuning, is countered by tuning L 0 . The kinetic part of the Hamiltonian can again be written in terms of the conjugate momenta where we have again removed φ CM -entries from the vectors and matrices. The potential part is = −E 2 cos √ 2φ 1 − √ 2φ 2 cos(2φ c ) + cos(2φ 1 + 2φ 2 ) where we have introduced the inductive energy E L0 = 1 2L0 . Notice how the inductor terms are not symmetric under the exchange of the indices 1 and 2. This is the asymmetry we will use to counter the asymmetry of the capacitance matrix. As before we can now extract the φ 2 i -terms and define the r i . Doing so we find These expressions indicate that the asymmetry between φ 1 and φ 2 may be contained in the sign of the factor 2 ± √ 2. Furthermore it seems justified that the asymmetry can be quantified by the ratios for capacitive and inductive contributions respectively. As we would like the asymmetries to cancel each other, we expect that these ratios have to be comparable in size.
We promote φ i and p i to operators with commutation relation [φ i , p i ] = i, and we introduce creation/annihilation operators a † i /a i for each mode, so that we may write We can then write the Hamiltonian in terms of creation/annihilation operators H = −2(E + 2(2 + √ 2)E L0 )r 2 1 (a † 1 − a 1 ) 2 − 2(E + 2(2 − √ 2)E L0 )r 2 2 (a † 2 − a 2 ) 2 − 2Er 2 c (a † c − a c ) 2 − E 2 cos r 1 (a † 1 + a 1 ) − r 2 (a † 2 + a 2 ) cos( √ 2r c (a † c + a c )) + cos( √ 2r 1 (a † 1 + a 1 ) + √ 2r 2 (a † 2 + a 2 )) + 4E L0 (2 + √ 2)r 2 1 (a † 1 + a 1 ) 2 + (2 − √ 2)r 2 2 (a † 2 + a 2 ) 2 Again, the asymmetry between φ 1 and φ 2 appears to be contained in the sign 2 ± √ 2. As before we can use Eq. (37) to truncate the cosines to two levels. Let us first note a generalization of Eq. (38) cos k 1 (a † 1 + a 1 ) + k 2 (a † 2 + a 2 ) and finally X ij for i, j = 0, 1, 2, 3 and i < j, whose a, b'th entries are (X ij ) ab = 0, for ab = ij, ji 1, for ab = ij, ji Together with the identity Z, A and B describe the contributions to the energy levels. In particular Z is similar to σ z , while A describes the anharmonicity and B describes a similar anharmonic energy shift beyond the regular anharmonicity only relevant for the third excited state and higher. In terms of the usual bosonic number operator N = a † a = diag(0, 1, 2, 3), we may say that Z corresponds to N , A corresponds to N 2 and B to N 3 . Alternatively, we may say that A corresponds to a † a † aa and B to a † a † a † aaa, which shows how A does not affect levels below the second excited one, while B does not matter for levels below the third excited. The X ij describe rotation or flipping between the i'th and j'th energy level, just like σ x describes rotation or flipping between the ground and excited state in the two-level case. In terms of these matrices, we can write e ik(a † +a) To reduce this to the two-level version one must merely exchange Z and X 01 with σ z and σ x respectively, and remove all other matrices except the identity.