Digital quantum simulation of spin models with circuit quantum electrodynamics

Systems of interacting quantum spins show a rich spectrum of quantum phases and display interesting many-body dynamics. Computing characteristics of even small systems on conventional computers poses significant challenges. A quantum simulator has the potential to outperform standard computers in calculating the evolution of complex quantum systems. Here, we perform a digital quantum simulation of the paradigmatic Heisenberg and Ising interacting spin models using a two transmon-qubit circuit quantum electrodynamics setup. We make use of the exchange interaction naturally present in the simulator to construct a digital decomposition of the model-specific evolution and extract its full dynamics. This approach is universal and efficient, employing only resources which are polynomial in the number of spins and indicates a path towards the controlled simulation of general spin dynamics in superconducting qubit platforms.

Quantum simulations using well controllable quantum systems to simulate the properties of another less tractable one [1,2] are expected to be able to predict the properties and dynamics of diverse systems in condensed matter [3,4], quantum chemistry [5] and high energy physics [6,7]. In particular, quantum simulations are expected to provide new insights into open problems such as modeling high-Tc superconductivity [8], thermalization [9] and non-equilibrium dynamics [10]. Up to now, several prototypical quantum simulations have been proposed and realized in trapped ions [11], cold atoms [12], and quantum photonics [13]. Examples include spin models [14][15][16], many-body physics [17], and relativistic quantum mechanics [18]. In the field of superconducting circuits quantum simulations are still in their infancy [19]. Topological properties [20,21] have been simulated recently, as have been fermionic models [22].
Quantum simulators are typically classified into two main categories, namely, analog and digital ones. Analog quantum simulators are designed to display intrinsic dynamics which are equivalent to those of the simulated system. While this approach is not universal it features control of the relevant Hamiltonian parameters better than in the system to be simulated. Instead, digital quantum simulators [2] can reproduce the dynamics of a quantum system via a universal digital decomposition of its Hamiltonian H = k H k into efficient elementary gates realizing H k . This approach is based on the Suzuki-Lie-Trotter expansion of the time evolution U (t) = e −iHt = lim n→∞ ( N k=1 e −iH k t/n ) n and was recently demonstrated experimentally in a trapped-ion digital quantum simulator [15].
Here we demonstrate digital quantum simulation of * ysalathe@phys.ethz.ch † Now at: IBM T. J. Watson Research Center, Yorktown Heights, NY 10598, United States spin systems [16] in an architecture known as circuit quantum electrodynamics (QED) [23].
Our experiments are carried out with two superconducting transmon qubits [24] coupled dispersively to a common mode of a coplanar waveguide resonator (see Appendix A for the device layout and setup diagram). We operate the circuit at 30 mK in a dilution refrigerator. The qubits Q1 and Q2 interact with a coplanar waveguide resonator with a fundamental resonance frequency at 7.14 GHz which serves both as a quantum bus [25] and for readout [26].
The natural two-qubit interaction is the XY exchange coupling [25] H xy 1,2 = J 2 (σ x 1 σ x 2 + σ y 1 σ y 2 ) mediated by virtual photons in a common cavity mode, which we also refer to as the XY interaction. Here, σ x,y i are the Pauli operators acting on qubit i and J denotes the effective qubit-qubit coupling strength [27]. The XY interaction is activated by tuning the transition frequency of qubit Q1 (5.44 GHz) into resonance with qubit Q2 (5.24 GHz) for a time τ using nanosecond time scale magnetic flux bias pulses [28] (see Appendix B). When the qubit transition frequencies are degenerate, the resonator-mediated coupling strength is spectroscopically determined to be J = −40.4 MHz. To make the presentation of the simulation results independent of the actual J, we express the interaction time τ for a given J in terms of the acquired quantum phase angle 2|J|τ . In our setup, the action of the XY gate (Fig. 1a) is characterized by full process tomography for a complete set of 16 initial twoqubit states and a series of 25 different interaction times τ finding process fidelities no lower than 89 % (see Appendix D).
In Fig. 2a,b we present non-stationary spin dynamics under the XY exchange interaction for a characteristic initial two-qubit state |↑ (|↑ + |↓ )/ √ 2 with spins pointing in perpendicular directions along +z and +x, respectively. During the XY interaction, the state of one spin is gradually swapped to the other spin and vice versa with (green) in the beginning and the final state is characterized using single-qubit basis rotations R tom 1,2 and joint two-qubit readout (yellow). (b) Digital quantum simulation of the twospin Heisenberg (XYZ) interaction for time τ . The first step after state-preparation is to apply the XY gate for a time τ (dashed box labeled as XY). In the second and third steps (dashed boxes with labels XZ and YZ), XZ and YZ gates are realized using single-qubit rotations R ±π/2 x,y (blue) by an angle ±π/2 about the x or y axis transforming the basis in which the XY gate acts. (c) Protocol to decompose and simulate Ising spin dynamics in a homogeneous transverse magnetic field. The circuit between the bold vertical bars with two dots is repeated n times, invoking each XY and phase gates for a time τ /n. See text for details. The actual pulse scheme is provided in Appendix C.
a phase angle of π/2. This corresponds to the iSWAP gate [29]. As a consequence, the measured Bloch vectors move along the YZ and XZ planes. For a quantum phase angle of 2|J|τ = π they point along the +y and +z directions respectively in good agreement with the ideal unitary time evolution indicated by dashed lines in Fig. 2a,b. We also find that the two-qubit entanglement characterized by the measured negativity [30] of 0.246 is close to the maximum expected value of 0.25 for this initial state at a quantum phase angle of π/2. As a consequence the Bloch vectors do not remain on the surface of the Bloch sphere but rather lie within the sphere. The anisotropic Heisenberg model describes spins interacting in three spatial dimensions where the sum is taken over pairs of neighbouring spins i and j. J x , J y and J z are the couplings of the spins along the x, y and z coordinates, respectively. Since it does not occur naturally in circuit QED we decompose the Heisenberg interaction into a sequence of XY and single-qubit gates, as shown in Fig. 1b. We combine three successive effective XY, XZ and YZ gates derived from the XY gate by basis transformations [16] to realize the isotropic Heisenberg model with J x = J y = J z = J versus interaction time τ . Since the XY, XZ and YZ operators commute for two spins the Trotter formula is exact after a single step.
To compare the Heisenberg (XYZ) interaction with the XY exchange interaction we have prepared the same initial state as presented in Fig. 2a,b. The isotropic Heisenberg interaction described by the scalar product between two vectorial spin 1/2 operators preserves the angle between the two spins. As a result, the initially perpendicular Bloch vectors of qubits Q1 and Q2 remain perpendicular during the interaction (Fig. 2c) and rotate clockwise along an elliptical path that spans a plane perpendicular to the diagonal at half angle between the two Bloch vectors (Fig. 2c).
In accordance with theory, the XYZ interaction leads to a full SWAP operation for a quantum phase angle of 2|J|τ = π/2 where the Bloch vectors point along the +x and +z directions. For the given initial state, we observed a maximum negativity of 0.210 close to the expected value of 0.25 for the Heisenberg interaction at a quantum phase angle of 2|J|τ = π/4. As for the XY interaction we have characterized the Heisenberg interaction with standard process tomography finding fidelities above 82 % for all quantum phase angles 2|J|τ .
Next, we consider the quantum simulation of the Ising model with a transverse homogeneous magnetic field where the magnetic field B pointing along the z axis is perpendicular to the interaction given by Jσ x i σ x j . Since the two-spin evolution ( Fig. 1c) is decomposed into twoqubit XY and single-qubit Z gates which do not commute, the transverse field Ising dynamics is only recovered using the Trotter expansion in the limit of a large number of steps n for an interaction time of τ /n in each step. To realize the Ising interaction term using the exchange interaction, the XY gate is applied twice for a time τ /n, once enclosed by a pair of π pulses on qubit Q1. This leads to a change of sign of the σ y 1 σ y 2 term which thus gets canceled when added to the bare XY gate. The external magnetic field part of the Hamiltonian is realized as single-qubit phase gates R φ z which rotate the Bloch vector about the z axis by an angle φ = Bτ /n per Trotter step. These gates are realized by detuning the respective qubit by an amount δ from its idle frequency corresponding to an effective B-field strength of B = 2πδ.
We experimentally simulate the non-stationary dynamics of two spins in this model for the initial state |↑ (|↑ − i|↓ )/ √ 2 which is well-suited to assess the simulation performance. In Fig. 3a expectation values for the digital simulation of the σ z 1,2 -components of the two spins are shown, as well as the two-point correlation function σ x 1 σ x 2 . The σ z 1,2 -components of the spins represented by the red and blue datasets in Fig. 3a, respectively, oscillate with a dominant frequency component of 2J due to the presence of the interaction term ∝ σ x 1 σ x 2 . Likewise, the XX correlation σ x 1 σ x 2 represented by the yellow dataset in Fig. 3a is non-stationary and oscillates at rate 2 √ B 2 + J 2 = 2 √ 10J ≈ 6.3J due to the presence of a magnetic field of strength B = 3J. The evolution of the measured final state shows agreement with a theoretical model (solid lines in Fig. 3a) which takes into account dissipation and decoherence with deviations being dominated by systematic gate errors (see Appendix E).
In Fig. 3b the fidelity of the simulated state is compared to the expected state at characteristic quantum phase angles both for the experimental realization (colored bars) and the ideal Trotter approximation (wire frames) after the nth step. In an ideal digital quantum simulator the theoretical fidelity (wire frame) converges for an increasing number of steps n (Fig. 3b). The experimental fidelity, however, reaches a maximum for a finite number of steps (Fig. 3b) after which it starts to decrease due to gate errors and decoherence [16]. As expected, the Trotter approximation converges faster for smaller quantum phase angles 2|J|τ . For 2|J|τ = π/4 the peak experimental fidelity (Fig. 3b) of 98.3 % is already observed for n = 1, whereas for 2|J|τ = 3π/2 the optimum of 80.7 % is observed for n = 5.
In future experiments, transmission line resonators may provide a means to design multi-qubit devices with non-local qubit-qubit couplings that directly reflect the lattice topology of spin systems such as frustrated magnets. Moreover, the incorporation of cavity modes as explicit degrees of freedom in the simulated models [31], following an analog-digital approach, and the integration of optimal control concepts, will be instrumental to scale the system to larger Hilbert-space dimensions. With this, the circuit QED architecture offers considerable potential for surpassing the limitations of classical simulations, which can be facilitated by using efficient digital decompositions of spin Hamiltonians as pursued in this work.
Acknowledgments The authors would like to thank Abdufarrukh Abdumalikov and Marek Pechal for helpful discussions. Furthermore we owe gratitude to Lars Steffen, Arkady Fedorov, Christopher Eichler, Mathias

Baur, Jonas
Mlynek who contributed to our experimental setup. We would also like to thank Tim Menke and Andreas Landig for contributions to the calibration software used in the present experiment.

Appendix A: Chip architecture and measurement setup
The present experiment was performed using two superconducting transmon [24] qubits Q1 and Q2 and one coplanar waveguide resonator R1 on a microchip (Fig. A1). The resonator R1 has a fundamental resonance frequency of ν r = 7.14 GHz. From spectroscopic measurements we have determined the maximum transition frequencies ν max = {5.55, 5.53} GHz and charging energies E C /h ≈ {260, 260} MHz of the qubits Q1 and Q2, respectively, where h is the Planck constant. The qubits Q1 and Q2 are coupled to resonator R1 with coupling strengths g/2π ≈ {120, 120} MHz. For this experiment the qubit transition frequencies in their idle state were offset to ν = {5.440, 5.240} GHz by applying a constant magnetic flux threading their SQUID loops with miniature superconducting coils mounted underneath the chip. At these idle frequencies, the measured energy relaxation and coherence times were T 1 = {7.1, 6.7} µs and T 2 = {5.4, 4.9} µs, respectively. The transition frequencies of the qubits Q3 and Q4 were tuned to 4.5 GHz and 6.1 GHz such that they do not interact with Q1 and Q2 during the experiment.
A schematic diagram of the measurement setup is shown in Fig. A2a. To realize two-qubit XY gates and single-qubit phase gates (Z), controlled voltage pulses generated by an arbitrary waveform generator (AWG) are used to tune the flux threading the SQUID loop of each qubit individually using flux bias lines [28]. The single-qubit microwave pulses (X,Y) are generated using sideband modulation of an up-conversion in-phase quadrature (IQ) mixer (Fig. A2b) driven by a local oscillator (LO) and modulated by an arbitrary waveform generator (AWG). The same up-conversion LO is used for the microwave pulses on both qubits to minimize the phase error introduced by phase drifts of microwave generators. We have used a quantum-limited parametric amplifier (PA) to amplify readout pulses at the output of R1 (Fig. A2c). Here the Josephson junction based amplifier in form of a Josephson parametric dimer (JPD) [32] is pumped by a strong pump drive through a directional coupler (D). To cancel the pump leakage, a phase (φ) and amplitude (A) controlled microwave cancelation tone is coupled to the other port of the directional coupler (D). Three circulators (C1-3) were used to isolate the sample from the pump tone. A circulator (C4) at base temperature followed by a cav- ity band-pass filter (BP) and another circulator (C5) at the still stage were used to isolate the sample and JPD from higher-temperature noise. The transmitted signal is further amplified by a high electron mobility transistor (HEMT) at the 4.2 K stage and a chain of ultra-low-noise (ULN) and low-noise (LN) amplifiers at room temperature as shown in Fig. A2d. The amplified readout pulse is down-converted to an intermediate frequency (IF) of 25 MHz using an IQ mixer (Fig. A2e) and digitally processed by field-programmable gate array (FPGA) logic for real-time data analysis.

Appendix B: Implementation of the XY gate
The interaction between two qubits with degenerate transition frequencies dispersively coupled to the same CPW resonator is described by the exchange coupling [33] J(σ + 1 σ − 2 + σ − 1 σ + 2 ) which can also be written in terms of Pauli operators as J 2 (σ x 1 σ x 2 + σ y 1 σ y 2 ). We activate this interaction by tuning the transition frequency of qubit Q1 into resonance with qubit Q2 with a flux pulse (Fig. A3) for an interaction time τ which we varied from 0 to 60 ns. At the frequency of qubit Q2, we obtain a coupling strength J = −40.4 MHz from a fit to the spectroscopically measured avoided crossing. To compensate overshoots of the flux pulse due to the limited bandwidth of the flux line channel, we use an inverted linear filter based on room-temperature response measurements of the flux line channel and in-situ Ramsey measurements of the residual detuning of qubit Q1 in the time interval from 0 to 2 µs after the flux pulse.
Since the outcome of the XY gate depends strongly on the relative phase of the two-qubit input state, we have used the same LO signal for the upconversion of the single-qubit pulses acting on both qubits Q1 and Q2 (green lines in Fig. A2a). Then the initial relative phase between the qubits is defined solely by the pulse sequence generated by the AWG and the cable lengths. In addition, we choose the shape of the flux pulse that realizes the XY gate such that the dynamic phase acquired by qubit Q1 during the idle time and the rising edge of the flux pulse cancels any unwanted relative phase offset of the initial state. We satisfy this condition by tuning the frequency of Q1 to an intermediate level (buffer) for a fixed time of 16 ns before and after the XY gate (Fig. A3). A suitable buffer level is found by performing Ramseytype experiments with a single XY gate while sweeping the buffer amplitudes. This calibration procedure is carried out for each interaction length of the XY gate. The second buffer at the falling edge of the flux pulse is used to ascertain that the relative phase between the qubits after tuning qubit Q1 back to its original position is the same as the initial relative phase. interaction. To decouple this undesired effect we have used a standard refocusing technique [36] implemented by two consecutive π pulses on qubit Q2 (magenta boxes in Fig. A4). In the end of each pulse sequence we perform dispersive joint two-qubit state-tomography [37] by single-qubit basis transformations followed by a pulsed microwave transmission measurement through resonator R1.

Appendix D: Process tomography
We perform standard two-qubit process tomography [38,39] of the XY gate and of the simulated isotropic Heisenberg (XYZ) model for a varying interaction time τ . Fig. A5 shows the process χ matrices characterizing the XY gate for a quantum phase angle π/2 (Fig. A5a) and π (Fig. A5b) corresponding to a √ iSWAP gate [29,40] and an iSWAP gate [41,42] with process fidelities of 97.8 % and 95.3 %, respectively. Heisenberg interaction with a quantum phase angle π/2 leads to a SWAP gate (Fig. A6a) with a process fidelity of 86.1 %. While the SWAP gate belongs to the two-qubit Clifford group, there is no natural interaction in standard circuit QED architecture to directly implement the SWAP gate [43,44]. For a phase angle π, the Heisenberg interaction is an identity gate (Fig. A6b) with a process fidelity of 83.6 %.

Appendix E: Error contributions
The single-qubit gate fidelities measured by randomized benchmarking [45][46][47] amount to 99.7 %. The dominant contribution to the loss in fidelity originates from the two-qubit XY gates for which a process fidelity F p,XY = 95.7 % is obtained from process tomography averaging over all quantum phase angles. This indicates that the errors in the implementation of the XY gate limit the fidelity of the final state of the quantum simulation. To confirm this, we calculate the expected process fidelity for the Heisenberg and Ising protocol from the observed XY gate fidelity by assuming independent gate errors in all three steps. For the Heisenberg (XYZ) model simulation neglecting the small singlequbit gate errors, we expect a mean process fidelity F p,XYZ ≈ 1−3(1−F p,XY ) = 87.1 %, which is close to the observed value of 86.3 %. For the Ising model simulation we expect a process fidelity of F p,Ising ≈ 1−2n(1−F p,XY ). From the relation F s = (dF p + 1)/(d + 1) between state (F s ) and process fidelity (F p ), we obtain the expected mean state fidelities of {93. To estimate the dominant source of systematic errors, we consider a model which includes relaxation (T 1 ) and dephasing (T 2 ) and state-dependent phase errors described by an effectiveJ z σ z 1 σ z 2 term with interaction strengthJ z . In addition, we include an extra offset in the single-qubit phase gate acting on qubit Q2 from cross talk of the flux pulses acting on qubit Q1 in each Trotter step. By fitting the final state predicted by this model to the observed states, we estimate an unwanted interaction angleJ z τ z of approximately 2.3 • and a constant phase offset of 4.6 • .  Pulse sequences that are applied on qubit Q1 (red) and qubit Q2 (blue) to implement the Heisenberg (a) and Ising spin (b) models. The Gaussian-shaped DRAG microwave pulses are applied to the charge lines of the respective qubits to implement single-qubit rotations R φ x,y about the x or y axis of the Bloch vector by an angle φ. Each sequence starts with the preparation of an initial state (green boxes) and ends with microwave pulses for basis rotations to perform state-tomography (yellow boxes). The microwave pulses marked with magenta boxes are used for refocussing. The black vertical bars with the two dots in panel (b) indicate that the enclosed pulse sequence is repeated n times. The XY gates are realized by applying flux pulses to the flux line of qubit Q1 for a time τ /n. The phase gates R φ/n z are implemented by detuning the transition frequency of each qubit from their idle frequencies applying flux pulses for a time τ /n. The numbers stated below the pulses on qubit Q1 represent timescales in ns. A5. (a) Measured real and imaginary part of the XY process χ matrix (Re χ, Im χ), in the basis {I = identity, X = σx,Ỹ = −iσy, Z = σz}, describing the mapping from any initial state to the final state for a quantum phase angle of 2|J|τ = π/2. The dashed wire frames represent the theoretically optimal matrix elements and the colored bars represent measured positive (blue) and negative (red) matrix elements. The fidelity of the experimentally observed process with respect to the ideal process is indicated in the black boxes. (b) As in (a) for a phase angle π.

FIG.
FIG. A6. (a) Measured real and imaginary part of the Heisenberg (XYZ) process χ matrix (Re χ, Im χ), in the basis {I = identity, X = σx,Ỹ = −iσy, Z = σz}, describing the mapping from any initial state to the final state for a quantum phase angle of 2|J|τ = π/2. The dashed wire frames represent the theoretically optimal matrix elements and the colored bars represent measured positive (blue) and negative (red) matrix elements. The fidelity of the experimentally observed process with respect to the ideal process is indicated in the black boxes. (b) As in (a) for a phase angle π.