Tuneable hopping and nonlinear cross-Kerr interactions in a high-coherence superconducting circuit

Analog quantum simulations offer rich opportunities for exploring complex quantum systems and phenomena through the use of specially engineered, well-controlled quantum systems. A critical element, increasing the scope and flexibility of such experimental platforms, is the ability to access and tune in situ different interaction regimes. Here, we present a superconducting circuit building block of two highly coherent transmons featuring in situ tuneable photon hopping and nonlinear cross-Kerr couplings. The interactions are mediated via a nonlinear coupler, consisting of a large capacitor in parallel with a tuneable superconducting quantum interference device (SQUID). We demonstrate the working principle by experimentally characterising the system in the single- and two-excitation manifolds, and derive a full theoretical model that accurately describes our measurements. Both qubits have high coherence properties, with typical relaxation times in the range of 15 to 40 microseconds at all bias points of the coupler. Our device could be used as a scalable building block in analog quantum simulators of extended Bose-Hubbard and Heisenberg XXZ models, and may also have applications in quantum computing such as realising fast two-qubit gates and perfect state transfer protocols.


INTRODUCTION
Analog quantum simulations, where engineered systems emulate the behaviour of other, less accessible quantum systems in a controllable and measurable way [1], show significant promise for improving our understanding of complex quantum phenomena without the need for a full fault-tolerant quantum computer [2][3][4][5]. In this paradigm, the versatility of the simulator is determined by the range of interaction types and complexity accessible to the emulating quantum system. Promising avenues for pushing beyond what can be simulated with a classical machine include the study of highly interacting manybody systems [6][7][8][9][10]. Superconducting circuit quantum electrodynamics (QED) is a very attractive platform for analog quantum simulation because of site-specific control and readout, and because of the flexible and engineerable system designs, which have led to the study of many interesting effects [11][12][13][14][15][16][17][18][19][20][21][22]. Adding new components to the circuit QED design toolbox such as novel types of interactions can dramatically increase the range of phenomena that can be simulated [23][24][25]. For example, for exploring exotic effects, like quantum phase transitions in systems of strongly correlated particles, it is important to be able to access and rapidly tune between different many-body interaction regimes.
In situ tuneable couplers have been successfully realised in a variety of circuit QED architectures [26][27][28][29][30][31][32][33][34][35], in particular using the more coherent transmon design [36][37][38][39]. In recent experiments, transmon arrays with tuneable exchange-type hopping interactions, have been employed to study many-body localisation phenomena of Bose-Hubbard and spin-1/2 XY models [19,20]. However, moving beyond linear couplings to incorporate addi-tional nonlinear interactions would enable the emulation of far more complex Hamiltonians. For example, nonlocal cross-Kerr interactions, present in extended Bose-Hubbard models [40,41], introduce much richer manybody phase diagrams, leading to intriguing phenomena such as crystalline and supersolid phases of light as the ratio of the hopping and cross-Kerr coupling strengths is varied [23,24]. In the qubit context, nonlinear cross-Kerr coupling, sometimes referred to as longitudinal coupling, is essential for engineering plaquette interactions in lattice gauge theories [25] and gives access to a large class of quantum-dimer and XYZ spin-model Hamiltonians.
Here, we demonstrate tuneable hopping and cross-Kerr interactions in a highly coherent two-transmon unit cell. Specifically, using a large capacitor in parallel with a tuneable nonlinear inductor as a coupling element, we are able to tune the ratio of the two coupling strengths, even suppressing hopping completely while maintaining a non-zero cross-Kerr coupling, giving access to different interaction regimes. We comprehensively characterise the energy landscape of this building block using different spectroscopic techniques. We show excellent agreement with a full theoretical model we have developed to describe the underlying circuit Hamiltonian including higher transmon excitation manifolds. Finally, we have thoroughly studied the qubit coherence as a function of the coupler bias, showing high relaxation times of 15 − 40 µs, and dephasing times reaching up to 40 µs at flux-insensitive operating points. Our work outlines a new recipe for building scalable analog quantum simulators of complex Hamiltonians using coupled transmon arrays, and our theoretical model provides an invaluable tool for designing and realising larger scale implementations.

Implementing nonlinear couplings
The working principle of the coupler is similar to that of a band-stop LC filter, relying on the fact that its impedance is infinite on resonance, as currents through the capacitor and the inductor interfere destructively. We implement a nonlinear analog of this circuit by using a nonlinear superconducting quantum interference device (SQUID) as the inductor, realising tuneable cross-Kerr and nearest-neighbour hopping interactions (Fig. 1a). The Josephson nonlinearity of the SQUID gives rise to tuneable higher-order nonlocal terms [40], including cross-Kerr interactions, which are equivalent to longitudinalσ zσz coupling in the qubit subspace. By contrast, the linear single-excitation hopping between the two sites is mediated by both the capacitor C c at a constant rate J C , and by the inductor at a tuneable rate −J L . Because of interference between these two processes, the hopping strength tunes in a different way from the cross-Kerr coupling, making different many-body interaction regimes accessible. In particular, at the point where the hopping rates cancel (J C = J L ), we can access a purely nonlinear regime with zero linear interaction.
The nonlinear coupler is implemented in a circuit QED device of two superconducting transmon qubits, the basic building block required for future lattice implementations (Fig. 1a). The optical micrograph along with a circuit diagram of the device are shown in Figs. 1b and 1c, respectively. Each transmon, consisting of two superconducting islands connected by an interdigitated capacitor C and a tuneable SQUID inductance, resonates at a plasma frequency ω ( √ 8E J E C − E C )/ , with Josephson energy E J and charging energy E C = e 2 2 C , where C is the effective transmon capacitance (see Supplementary Eq. (S12)). Transmon frequencies can be independently tuned using on-chip flux lines (Φ 1,2 ), and spectroscopy is performed through dispersively coupled readout resonators (R 1,2 ) measured via a common microwave feedline (see Supplementary Fig. S5 for the full measurement setup and Fig. S6 for qubit spectroscopy vs Φ 1,2 ). Microwave drives are applied via either the resonators or dedicated drivelines. The coupler capacitance C c and flux-tuneable SQUID (Φ 3 ) connect galvanically to the two transmons.
The physics of the two-transmon building block is, to a good approximation, well described by an extended Bose-Hubbard model, which is a Heisenberg XXZ spin model in the qubit regime. To achieve this, we need to operate the qubits detuned from coupler resonances, so that coupler excitations do not participate directly in the system dynamics. Under this condition, the system can be described by a simplified two-transmon Hamiltonian: whereâ ( †) ,b ( †) are bosonic annihilation (creation) operators for each transmon in the uncoupled basis, with on-site nonlinearity U = EC 2 . The interaction between the two sites can be described by hopping of single excitations at a rate The capacitive coupling J C is fixed and determined by the ratio of the coupling capacitor C c to an effective capacitance C eff , which depends on the circuit network (see Supplementary Eq. (S27)). The interaction strengths J L and V are determined by tuning the Josephson energy of the coupling SQUID, E c J = E c (max) J cos (πΦ 3 /Φ 0 ). Importantly, the cross-Kerr coupling V is different from the diagonal coupling that can be observed in linearly coupled transmon architectures, where the self-Kerr nonlinearity of each transmon leads to an effective cross-Kerr coupling between the normal modes of the system. Such effective diagonal coupling scales with the hopping interaction strength and vanishes at J = 0 [46], making it impossible to tune the ratio J/V independently. In our design, however, the cross-Kerr interaction results directly from the nonlinearity of the coupling junction and tunes to zero at a different coupler bias from the linear hopping interaction, giving access to different interaction regimes. As Φ 3 is tuned towards the filtering condition (J L ∼ J C ), the linear hopping term is suppressed more rapidly than the nonlinear cross-Kerr coupling, allowing access to the J V regime. By contrast, when Φ 3 ∼ 0.5 Φ 0 , the cross-Kerr coupling is suppressed (V = J L = 0) and the dynamics are dominated by singlephoton hopping at a rate J C .
In the full treatment of the quantum circuit, the nonlinear inductance also gives rise to correlated hopping and two-photon tunnelling correction terms, which play a role in the higher-excitation manifold. These terms may also lead to interesting physical phenomena, which we return to in the discussion section. We derive the full quantum model in the Supplementary Information [42], along with a classical normal mode analysis providing supporting intuition for the full system behaviour.

Tuneable single-photon hopping
We demonstrate our ability to tune the linear singlephoton hopping interaction between the two transmons, by measuring qubit-qubit avoided crossings, in Fig. 2. The top panels show example crossings, as qubit 1 is tuned on resonance with the other at ∼ 6.6 GHz, in three different coupling scenarios, J L > J C in Fig. 2a, Fig. 2b, and the J C -dominated regime in Fig. 2c. The measurements in Figs. 2a, c are performed via readout resonator R 2 , while in the zero coupling case (Fig. 2b) we measure via R 1 . The range of typical coupling strengths that can be achieved with this device is illustrated in Fig. 2d, where we plot |J|/2π vs calibrated coupler bias Φ 3 . We note that we have measured larger coupling strengths, up to 140 MHz, when operating at different qubit frequencies ∼ 5.  Table I). The normal-mode splitting gets suppressed at the crossover between inductively and capacitively dominated coupler regimes (Φ3/Φ0 0.3). Cross-talk effects between the different flux channels have been calibrated out experimentally.
Note that a higher transition of a lower-frequency sloshing mode of the circuit, crossing with the qubits around this point, limits the observed on-off ratio in this device to ∼ 10 (see Supplementary Fig. S3). This low-frequency mode, hybridising with the qubits, is associated with currents flowing only through the coupling inductor [42], and could be avoided with slightly different design parameters (see Supplementary Fig. S4). Additional avoided crossing measurements in the regime where the linear coupling gets suppressed and reverses sign are plotted in Supplementary Fig. S9, with spectroscopy on both qubits.
For a more complete characterisation of the tuneable hopping interaction, we fit the experimentally measured coupled-qubit spectrum with our theoretical model of the quantum circuit. More specifically, in Fig. 2e, we fix the qubits on resonance and plot the normal-mode splitting between the dressed states |+ = |01 +|10 . As J becomes comparable to the transmon anharmonicity, the |02 , |20 and |11 levels mix with each other, resulting in an effective repulsion of |11 . On the other hand, a qubit-qubit interaction with negative cross-Kerr coupling results in lowering the energy of |11 . (b) Combined two-tone spectroscopy data (black dots) showing The red curve is theory prediction assuming only hopping interaction between two transmons (V = 0), while the blue one shows simulation results obtained by taking into account also the higher-order nonlinear contributions, V 4 (a + a † ) 2 (b + b † ) 2 , which include the dominant cross-Kerr term. The parameters used are listed in Table I.
, as we tune the calibrated coupler flux bias. The blue and red curves are theory fits of the single-excitation qubit manifold to the quantum circuit Hamiltonian (Eq. (S11) in the Supplementary), showing excellent agreement with the experimentally obtained spectrum. The parameters obtained from this fit, which neglects higher-order couplings of each transmon to the low-frequency sloshing mode, are listed in Table I. Note that the antisymmetric mode frequency ω − /2π is unaffected by coupler tuning, which reflects the fact that this mode is only associated with charge oscillations across the qubit junctions (see Supplementary Fig. S2).

Tuneable nonlinear cross-Kerr coupling
As already discussed, a key feature of our implementation which differentiates it from previous tuneable couplers, is the nonlinear cross-Kerr interaction which can be tuned into different coupling regimes relative to the hopping strength and does not zero when J does. This cross-Kerr coupling, which in different contexts is referred to as σ z σ z , longitudinal [43], or dispersive [44], does not involve excitation hopping processes. Its presence, how-ever, does influence the dynamics as the occupation at one site can alter the energy level spectrum of a neighbouring site, in a process analogous to photon scattering [23].
In a coupled two-qubit system, the effect of cross-Kerr interaction can be seen as a shift of the energy level of the |11 state and can be determined spectroscopically from the transition energies relative to the ground state ω 11 − ω − − ω + . For weakly anharmonic systems, such as the transmon, this picture becomes more complicated in the presence of linear hopping (Fig. 3a). A three-state analysis at the two-excitation manifold reveals that |11 , |02 and |20 also couple to each other, resulting in an effective upwards repulsion of the |11 state [45], which scales as ∼ J 2 /E C [46]. Because the direction of this effect competes with the negative cross-Kerr shift, when the effects are similar in size, it can hinder the observation of cross-Kerr coupling in an individual spectroscopy measurement. To separate the two effects, it is therefore necessary to measure the shift for different coupling levels.
We experimentally demonstrate the presence of cross-Kerr interaction between the two transmons, by measuring transitions in the two-excitation manifold of the coupled system. More specifically, we extract the frequency shift of |11 at different couplings from a series of two-tone spectroscopy measurements (see Supplementary  Fig. S7), focusing on the inductively dominated regime 0 Φ 3 /Φ 0 0.25. In order to distinguish between the negative cross-Kerr shift and the positive shift from linear coupling J, we plot ω11−ω−−ω+ J as a function of J, in Fig. 3b. The red curve is theoretical prediction assuming only hopping interaction (V = 0) between the two transmons. The blue curve shows numerical results after diagonalising the transmon-transmon Hamiltonian with the full nonlinear coupling terms V 4 (a + a † ) 2 (b + b † ) 2 , which includes the dominant cross-Kerr interaction. We use the parameters listed in the second column of Table I, which differ slightly from the fitted parameters of Fig. 2e to accommodate the effects of extra higher-order terms (see later in Fig. 5). At J = 0 (Φ 3 ∼ 0.3 Φ 0 ) the cross-Kerr coupling |V |/2π is around 4 MHz, and it reaches a maximum of 10 MHz at Φ 3 = 0. We were unable to explore the region J/2π < 20 MHz in this device, due to a higher transition of the lower frequency sloshing mode hybridising with the qubits (see Supplementary Fig. S3).

Qubit coherence
Maintaining high coherence for all interaction strengths is an essential requirement for future implementations based on our building block device. In Fig. 4, we investigate the individual qubit properties as a function of the coupler bias, with the transmons far detuned from each other by ∼ 1. top and bottom sweetspots (Φ 1 = 0, Φ 2 /Φ 0 = 0.5). In Figs. 4a, b, we demonstrate high relaxation times T 1 (15 − 40 µs) over the entire coupling range. We also report a systematic study of dephasing times in our device, obtained from spin-echo measurements (Figs. 4c, d). Full theoretical model of the higher excitation manifold of the device. Measurement of the coupled system eigenspectrum (as Fig. 2e) at higher powers. Dots show theoretical calculations using the full quantum circuit Hamiltonian including all next-to-leading order transmon-transmon coupling terms and the sloshing mode contributions, for the parameters listed in Table I. Simulations are performed using a Hilbert space dimension of N = 15 3 .

DISCUSSION
Our work demonstrates a key building block for circuit QED devices capable of exploring a rich vein of manybody physics in extended Hubbard models. In the context of driven nonlinear arrays, a chemical potential term, µ = ω q − ω d , could be straightforwardly implemented by coherently driving the transmons through the drive lines at a frequency ω d , which would enable the study of rich many-body phase diagrams with all J, V, µ tuneable [23,24]. It may also be possible to implement topological pumping of interacting photons, by modulating the frequency of each transmon, to study bosonic transport of Fock states in a nonlinear array configuration [48]. In realisations where higher-excitation manifolds might be explored, additional higher-order terms arising from the junction nonlinearity should be considered. For example, in our implementation, the nonlinearity of the medium leads to correlated hopping terms, such that a photon can hop between sites, on the condition that another photon is present. Additional contributions at higher excitation manifolds involve photon-pair tunnelling processes which might lead to exotic phenomena such as fractional Bloch oscillations [49]. These contributions are explicitly derived in the Supplementary Information [42], following a full quantum mechanical treatment of our circuit. In Fig. 5, we plot the coupled system eigenspectrum up to the two-excitation manifold, based on our full theoretical model including all next-to-leading order terms, which is found to be in excellent agreement with our data obtained at high powers. Our circuit can also be used to study many-body effects in spin models. When the transmon anharmonicity is much larger than the coupling strength (E C J), a truncation to the qubit subspace is justified, and the transmon-transmon interaction is described by a Heisenberg XXZ Hamiltonian 2J(σ xσx +σ yσy ) + Vσ zσz .
The coupling strengths available in this device are J/2π ∼ 8 − 140 MHz, V /2π ∼ 0 − 10 MHz, with orders of magnitude lower qubit decay rates (3 − 15 kHz). In a slightly different design, with a larger coupling capacitor C c , it would also be possible to further explore the J V regime (see Supplementary Fig. S4). One could then simulate an Isingσ zσz interaction Hamiltonian (J = 0), with antiferromagnetic couplings of ∼ 10 MHz. Additionally, time modulated magnetic fluxes threaded through the coupler SQUID can enable a large set of spin-spin interactions (e.g. pureσ xσx orσ yσy ) [50], therefore, giving access to emulating the dynamics of almost any spin model and exploring their phase diagrams. Connecting the coupler to four transmons on a lattice could enable simulating models with topological order such as the famous toric code [50]. Our circuit could also be employed to engineer plaquette terms in lattice gauge theories or ring-exchange interactions in dimer models, where a longitudinal coupling much larger than the hopping term is required in order to emulate effective fields on the lattice [25]. Moreover, a similar architecture, featuringσ zσz coupling between transmons has been proposed theoretically for the realisation of a microwave single-photon transistor [51].
In order to scale this circuit to larger lattice sizes, future experiments could take an approach where each transmon is connected to couplers via the same superconducting island, with the other island used for drive control and readout. Using a two-island transmon design has the advantages of reducing or eliminating the number of possible current loops involving current flow across qubit junctions, as well as allowing spurious flux crosstalk to be eliminated by linear compensation techniques (see Methods). Our coherence measurements (Fig. 4) suggest that this coupler design can be realised without significantly limiting qubit coherence times, showing promise for scaling up to larger lattice sizes.
In conclusion, the implemented circuit increases the range of available interactions and phenomena that can be explored with circuit QED analog quantum simulators. We have demonstrated hopping and cross-Kerr interactions with in situ tuneability between two transmon qubits in a flexible and scalable superconducting circuit. The observed decay rates are orders of magnitude lower than the coupling strengths, making this a viable platform for analog quantum simulation experiments. Moreover, our full theoretical model of the quantum circuit is in excellent agreement with the measurements, providing a powerful tool for designing future larger scale implementations.

Chip fabrication
The capacitive network of superconducting islands and ground plane, together with readout and control lines are defined on a thin NbTiN film [52] on top of a high resistivity Si substrate. The film is patterned using e-beam lithography on ARN7700 resist and etched with SF6/O2 plasma reactive-ion etching. Josephson junctions are then fabricated on each SQUID using Al-AlOx-Al shadow evaporation following e-beam lithography patterning on a PMGI/PMMA lift-off mask and HF dip to remove surface oxides on the NbTiN contact pads. Finally, after an e-beam patterning and reflowing PMGI step, followed by a second e-beam patterning of a MAA/PMMA resist stack, we evaporate Al airbridges which are used as crossovers above all lines in order to ensure a uniform ground plane.

On-chip flux cross-talk calibration
Due to the compact geometry of our device, on-chip cross-talk between all flux channels is quite significant and extra care is required in order to decouple them. This requirement is vital for independent control, especially for larger scale implementations where such effects could become a major experimental challenge. We employ a systematic calibration procedure (see Supplementary Fig. S8) which is enabled by the fact that the frequency of the lower circuit sloshing mode (3.2 GHz at flux insensitive point) is directly associated with tuning the coupling strength via Φ 3 . There is therefore one circuit degree of freedom corresponding to each bias channel. We track spectroscopically the frequency of each degree of freedom (e.g. qubit 1) around its flux insensitive point (top sweetspot) and determine the flux offset as we vary the other two channels (2 and 3). Repeating this for all three degrees of freedom and flux channels we were able to measure and calibrate all cross-talk effects, making all flux bias channels Φ 1,2,3 orthogonal. Note that the calibration method employed here allows us to dis-tinguish between the on-chip flux cross-talk effects from the intrinsic qubit frequency shifts that are expected by varying the coupling inductance [42]. The latter are deliberately not calibrated out in order to be able to fit the measurement data in Fig. 2e and Fig. 5 with the full circuit theory Hamiltonian, however we could straightforwardly compensate for them if required.

Device parameters
Fitting of the Full nonlinear Parameter single-excitation circuit model manifold (Fig. 2e)  The device parameters are presented in Table I. In the first column, we list the circuit parameters obtained by fitting the resonantly coupled transmon-transmon spectrum in the single-excitation manifold (Fig. 2e) with a simplified circuit model that neglects any higher-order couplings to the sloshing mode. The parameters in the second column are used in our full numerical model that includes all next-to-leading order terms in the circuit Hamiltonian (Supplementary Eq. (S11)), to describe the obtained data at higher excitation manifolds (Fig. 5).

DATA AVAILABILITY STATEMENT
The experimental data and numerical code are available by the corresponding author upon reasonable request.
The authors declare no conflict of interest.

I. DERIVATION OF FULL CIRCUIT HAMILTONIAN
Here, we analytically derive the full Hamiltonian that describes the two-transmon coupled system including all circuit degrees of freedom. We begin with a Lagrangian description of the circuit in the harmonic limit and proceed with the Hamiltonian formulation and canonical quantisation, following the methodology described in Ref. [S1, S2].

A. Analytical description in the harmonic limit
An intuitive picture about the circuit and the relevant degrees of freedom can be drawn in the harmonic limit, where all Josephson elements are approximated by linear inductors, L i = φ0 FIG. S1. Harmonic limit of the circuit diagram shown in Fig. 1, where each Josephson element is approximated by a linear inductor.
inductive energies of the linearised system are and where the node flux φ i (expressed in units of the reduced flux quantum φ 0 = /2e) is related to the potential at node i asφ i = V i . The system Lagrangian, therefore, is where is the capacitance matrix, and is the inverse of the inductance matrix of the circuit, expressed in the node flux basis φ T= [φ 1 , φ 2 , φ 3 , φ 4 ].


  + + -- Normal modes of the resonantly coupled transmon-transmon system in the harmonic limit. (a) Normal-mode frequencies as a function of the coupling inductance. The normal-mode splitting (blue and red curves) is suppressed at the point where the filter frequency (dashed curve) is resonant with both transmons. The sloshing mode frequency (green curve) vanishes at the limit of no inductive coupling (Lc → ∞). (b) -(d) Schematic representation of the circuit charge oscillations associated with each normal mode. Note that the antisymmetric mode is associated with charge oscillations across the transmon inductors only and therefore its frequency ω−/2π is independent of the coupling inductance Lc as shown in (a).
Solving the characteristic/secular equation in the resonantly coupled case (L 1 = L 2 ) yields the following normal modes (up to normalisation factors) where the coefficients a, b, c, d depend on the geometry of the circuit. The first two normal modes correspond to symmetric and antisymmetric combinations of the coupled transmons as schematically depicted in Figs. S2b, c. The third term corresponds to a sloshing mode associated with charge oscillations across the coupling element (Fig. S2d). The last term describes a zero-frequency rigid mode, corresponding to all capacitors charging up simultaneously. The normal-mode frequencies of the linearised circuit are plotted in Fig. S2a as a function of the coupling inductance L c , for our device parameters. Note that the normal-mode splitting and, subsequently, the coupling between the two transmons is suppressed at the point where they are both on resonance with the filter frequency 1/ √ L c C c (dashed curve in Fig. S2a).

B. Hamiltonian description of the nonlinear circuit
The inductive energy of the nonlinear circuit (Fig.1c in the main text) is where, in the second step, we have expressed it using the transmon and sloshing mode variables ψ A= φ 1 − φ 2 , ψ B= φ 4 − φ 3 and ψ S= ). The first two terms describe the Josephson energy of each transmon, while the third one describes the inductive coupling between them, as well as the inductive energy of the sloshing mode. Here, ψ A,B correspond to the uncoupled transmon modes in the limit C c = 0, L c → ∞. Defining ψ R= 1 2 (φ 1 + φ 2 + φ 3 + φ 4 ) and performing a change of basis from φ T to ψ T = [ψ R , ψ S , ψ B , ψ A ], the capacitance matrix transforms as Therefore, the system Lagrangian, expressed in the mode variable basis ψ, is The conjugate momenta Q i = ∂L ∂ψi , describing the charges associated with each mode, can be determined by inverting the capacitance matrix, since Q = [C ]ψ. Performing a Legendre transformation H = iψ i Q i − L and promoting all variables to operators, satisfying [ψ i ,Q i ] = i , we obtain the full circuit Hamiltonian: where , (S12) and is the determinant of the capacitance matrix.

Inductive coupling contributions
The linear contributions to the inductive coupling, arising from the first-order expansion of E c where the first term corresponds to a direct dipole-dipole interaction between the transmons due to the coupling inductor, and the second results in an indirect transmon-transmon coupling via the off-resonant sloshing mode. The last terms describe a linear contribution to the inductive energy of each mode.
where the first direct transmon-transmon coupling term describes correlated hopping processes including small corrections to the linear coupling, and the second term contains the cross-Kerr contribution. In our desing, the sloshing mode is at much lower frequency than the transmons, therefore the third term is not directly relevant, however, the fourth term results in qubit coupling to the 0 − 3 transition of the sloshing mode, which was observed in spectroscopy near the zero coupling region. As we show in Fig. S3, this transition always crosses around the zero coupling region, even when operating at different qubit frequencies. This is an undesired effect which arises due to the fact that the two-island transmon design introduces four nodes in our circuit. A possible solution to this, would be to carefully engineer where these higher transitions occur and make sure that they do not cross with the qubits at the critical operation points. For example, as we show in Fig. S4, a zero linear coupling regime could be reached for a slightly different circuit design with a larger coupling capacitor C c . Notice that, all the self-inductive terms to infinite order O(ψ 2 A , ψ 2 B , ψ 2 S , ψ 4 A , ψ 4 B , ψ 4 S , ...) result in the following contribution to the Hamiltonian The first term describes the Josephson energy of the sloshing mode, while the rest two terms account for a correction to the transmon Josephson energies as E c J is tuned. This change accounts for the measured transmon frequency shifts vs coupler bias Φ 3 , shown in Figs. 4e, f of the main text.

Capacitive coupling contributions and mode variable elimination
The two transmons couple directly via their charge degrees of freedom as 4Det[C ]Q AQB . Additionally they couple together indirectly via the off-resonant sloshing mode 1 The mode variable Q R has no inductive energy associated with it and only contributes to the analytical dynamics in (S11) by coupling with the transmons via 1 C ABRQ R (Q A +Q B ). Notably, this coupling would vanish if the two islands of each transmon had exactly the same capacitances to ground (C 1g = C 2g ). In our design there is a small asymmetry (C 1g = C 2g ), therefore, this contribution should be taken into account for an accurate theoretical description of the circuit. It is possible, however, to eliminate this variable, while taking into account the stray contributions to the coupling between the two transmons. To illustrate this, we rewrite these coupling terms as The first term is derived following a redefinition of the rigid modeQ R = f (Q R ,Q A ,Q B ). This zero frequency mode has no contribution to the equations of motion of all the other modes and can, therefore, be eliminated. The second term accounts for a renormalisation of the transmon charging energies, while the last term accounts for a small contribution to the direct capacitive coupling, which can be absorbed into a redefinition of J C . We have therefore eliminated one variable which reduces the computational power required for numerical modelling and makes the theoretical description of the circuit more elegant.

C. Circuit quantisation in the harmonic oscillator basis
The Hamiltonian describing the uncoupled system is given by The measurement setup including the wiring diagram and an optical micrograph of the full 2 mm × 7 mm chip is shown in Fig. S5. Our device was mounted on a printed circuit board (PCB) at the mixing chamber plate of a dilution refrigerator (fridge temperature ∼ 20 mK). The sample was shielded against external magnetic fields, using a pair of cryogenic mu-metal shields as well as aluminium layers, and infrared radiation, using an additional copper can coated with a mixture of silicon carbide and Stycast 2850 [S4]. Cryogenic attenuators and home-build eccosorb filters were used in all microwave lines, resulting in a total attenuation of ∼ 60 dB for each line. For the DC flux bias lines we use home-build eccosorb filters and commercial low-pass filters (VLFX) with 1.35 GHz frequency cutoff. The input line additionally has a 10 GHz low-pass filter, and the output signal passes through 3 circulators, connected in series, before being amplified using cryogenic (HEMT) and room-temperature amplifiers. The signal is then demodulated, amplified and registered with a data acquisition card.

 
FIG. S6. Frequency tunability of (a) qubit 1 and (b) qubit 2 vs dedicated flux bias (raw data including theory fits). All flux channels have been made orthogonal following the calibration procedure described in Fig. S8.

III. TWO-TONE SPECTROSCOPY
In order to determine the diagonal and cross-Kerr couplings ( Fig. 3 in main text), we need to access the twoexcitation manifold in the resonantly coupled system (Fig. S7a). However, directly exciting the even states with one source is difficult because of the transmon selection rules [S5]. This could, in principle, be achieved via single-tone spectroscopy at very high powers (as in Fig. 5), however this results in broadening and Stark shifting of the coupled qubits frequencies, causing some uncertainty in our measurements. We therefore perform two-tone spectroscopy, where we read out via one resonator and sweep the drive frequencies of two microwave tones (f S1 and f S2 ) applied through each driveline, as shown in Fig. S7b. The bright vertical lines, correspond to the dressed state transition frequencies ω + /2π, ω − /2π, which can be easily excited with each source. The diagonal 45 • lines correspond to twophoton transitions to |02 , |20 and |11 as different photons from each source combine together to excite these states. Each data point in Fig. 3b of the main text is extracted from a series of similar two-tone spectroscopy measurements.