Experimentally simulating the dynamics of quantum light and matter at deep-strong coupling

The quantum Rabi model describing the fundamental interaction between light and matter is a cornerstone of quantum physics. It predicts exotic phenomena like quantum phase transitions and ground-state entanglement in ultrastrong and deep-strong coupling regimes, where coupling strengths are comparable to or larger than subsystem energies. Demonstrating dynamics remains an outstanding challenge, the few experiments reaching these regimes being limited to spectroscopy. Here, we employ a circuit quantum electrodynamics chip with moderate coupling between a resonator and transmon qubit to realise accurate digital quantum simulation of deep-strong coupling dynamics. We advance the state of the art in solid-state digital quantum simulation by using up to 90 second-order Trotter steps and probing both subsystems in a combined Hilbert space dimension of ∼80, demonstrating characteristic Schrödinger-cat-like entanglement and large photon build-up. Our approach will enable exploration of extreme coupling regimes and quantum phase transitions, and demonstrates a clear first step towards larger complexities such as in the Dicke model.

Digital quantum simulations [1] promise a quantum advantage without a universal, fault-tolerant quantum computer, with applications in fields such as quantum chemistry [2,3] and condensed-matter physics [4][5][6][7]. In a digital quantum simulator, sequences of discrete interaction components synthesise the evolution of an artificial Hamiltonian, allowing access to more exotic dynamics than the simulator can realise naturally. Systems involving ultrastrong light-matter interactions raise significant challenges for both theoretical analysis [8][9][10][11][12][13] and experimental study [14], making them ripe candidates for exploration via quantum simulation.
The standard QRM [8] describes the dynamics of a single two-level atom (energy ω R q ) coupled to a single quantum harmonic field mode (energy ω R r ) by field-dipole interaction: where a (a † ) and σ − (σ + ) are annihilation (creation) operators for field mode and atom, respectively, σ j are the Pauli spin operators, and g R is the coupling strength.
Under small coupling (g R ω R q , ω R r ), this reduces to the Jaynes-Cummings (JC) model via the rotating-wave approximation: H JC = − ω q 2 σ z + ω r a † a + g aσ + + a † σ − , which contains only the excitation-number-conserving interaction terms, aσ + and a † σ − , and has an exact solution. In the USC regime (g R ∼ ω R q , ω R r ), however, the excitation-nonconserving terms aσ − and a † σ + cannot be neglected and only total parity [σ z n (−1) n |n n|] is conserved [29]. Without the strong symmetry of number conservation, the combination of an infinitedimensional oscillator with an explicitly quantum (twolevel) atomic system makes the full QRM difficult to solve [10]. Demonstrating the ground-state entanglement and large ground-state photon numbers that can arise in the QRM is an open challenge in USC research.
Ultrastrong-coupling dynamics can also produce nontrivial quantum states and significant build-up of photon numbers [29]. Many characteristic dynamical features of USC can already be seen in the degenerate-qubit limit ω R q = 0. Here, the interaction-picture Hamiltonian is a coherent drive on the oscillator mode, with an amplitude ±g R e iω R r t conditioned on the σ x basis state of the atom. The conditional coupling ±g R coherently displaces the field, but in a continuously rotating direction given by e iω R r t , creating two diametrically opposite cir-cular trajectories in phase space [see Fig. 3(a)]. Relating the diameter and circumference of these trajectories, πα max =αT R , with field displacement rateα = g R and period T R = 2π/ω R r , gives a maximum amplitude α max = 2r set by the relative coupling ratio r ≡ g R /ω R r . Figure 1(a) illustrates the atomic and photonic parity dynamics for characteristic USC regimes, starting in an eigenstate of the uncoupled system, |1 q ⊗ |0 r . Because this is a superposition of σ x eigenstates, evolution gives rise to an atom-field entangled state (Bell-cat state) [30], |+, +α q,r −|−, −α q,r . For r 1, the two trajectories remain virtually indistinguishable, giving evolution closely approximating simple JC dynamics with an atom-field detuning equal to ω R r [31]. As r increases, the curves start distorting from the sinusoidal JC exchange oscillations (USC regime), until dUSC is reached, where the parities exhibit a characteristic Gaussian-shaped "collapse", followed by flat plateaus and periodic revivals at multiples of T R . The cross-over between these dynamical regimes is related to the maximum distinguishability of the two coherent states of the field. When the paths separate completely, the qubit appears to be in a mixed state, with parity 0.5.
In our circuit QED simulator, the Rabi atom and field mode are simulated by a transmon qubit (Q R ) [32] and a coplanar waveguide resonator (R R ) with energies ω q and ω r , respectively [device shown in Fig. 1(b)]. Because the transmon is only weakly anharmonic (ω 0-1 q −ω 1-2 q ω 0-1 q ), directly increasing the qubitresonator coupling g leads to a breakdown in its qubit behaviour at small r, and full circuit quantization shows that dUSC cannot be reached for any circuit parameters [33]. Instead, building on the proposal in Ref. 23, we perform a digital simulation of the QRM for arbitrarily large r using a coupling in the manifestly non-USC regime (r < 10 −3 ). The full Rabi Hamiltonian can be decomposed into two JC-like interactions [23]: where H AJC = σ x H JC σ x contains only counter-rotating interaction terms, and the effective Rabi parameters g R = g, ω R r = 2∆ r and ω R q = ∆ q ≡ ∆ JC q − ∆ AJC q are not related to the natural circuit frequencies, but defined relative to a nearby rotating frame (∆ = ω − ω RF ), and can be arbitrarily small. Using the standard method of Trotterization [1], Rabi dynamics can therefore be simulated deep into the USC regime by decreasing ∆ r and ∆ q . Figure 1(c) illustrates the second-order Trotter step used here (see Methods and [31]). An asymmetric transmon with two flux-insensitive "sweet" spots [31,32] is driven and measured at its lower sweet spot (5.452 GHz) far below the resonator (6.381 GHz), and short JC interaction gates are applied by fast frequency-tuning flux pulses [34]. Experimentally, a rotating frame is usually defined by the frequency of a drive tone. Here, the choice of rotating frame specifies the required rotation axis of the π pulses which create the AJC interaction. By appropriately updating the pulse phases, which are controlled with high precision, we can therefore arbitrarily select the rotating frame detuning from the resonator, even though these pulses are applied far from both resonator and rotating frame (see Methods).
Numerical modelling of the digital Rabi protocol highlighted several challenges for device design and fabrication [31]. Most significantly, due to practical fluxpulsing bandwidths which limit the shortest achievable Trotter step, it is challenging to digitise fast compared with the dynamics. Reaching acceptably low Trotter error in interesting regimes of r therefore required small qubit-resonator coupling (here, g/2π = 1.95 MHz). This also placed constraints on other device parameters, including coherence (for long simulation times), flux-tuning precision and qubit-resonator frequency targetting (due to a very narrow resonance). An extra qubit Q W was strongly and dispersively coupled to R R to probe the intraresonator quantum state via its photon-dependent frequency shift (−1.26 MHz per photon) using pulse sequences based on Ramsey interferometry. We used Q W to implement a range of photon measurements: average photon number with a controllable dynamic range, average photon parity and, combined with coherent field displacements through an external input coupler, direct Wigner tomography of the resonator [9,31,35,36]. Qubits were driven and measured through dedicated read-out resonators.
We first experimentally simulate the QRM for the degenerate-qubit case over a wide range of r covering the full USC regime, from r ∼ 0.3 to r → ∞ (Fig. 1). We use 60 Trotter steps to simulate 1.2 µs of dynamics (gt = 4.68π) and measure either qubit or photon parity after each step. (Simulations start in the state |1 q ⊗ |0 r for all results in the main text.) A simplified pulse sequence is illustrated in Fig. 1(d). The qubit and photon parity dynamics [ Fig. 1(e, g)] show very similar qualitative behaviour, consistent with parity conservation. The revival periods T R are in excellent agreement with the predictions of USC Rabi dynamics, and strikingly different from those predicted for a pure Jaynes-Cummings interaction with the equivalent qubit-resonator detuning T JC = 2π/ 4g 2 + ∆ 2 q-r [31]. The experimental simulations also show the Gaussian-shaped parity collapse set by the simulated g R , which is a key signature of dUSC dynamics. From fits to the initial points of the qubit data, we calculate an average g R ≈ 2π × 1.79 MHz, slightly lower than the expected g R = g ≈ 2π × 1.95 MHz determined from independent spectroscopy and vacuum Rabi oscillations. This is consistent with a small residual fluxpulse distortion and provides the best estimate for the simulated g R achieved in these experiments.
From the observation of parity revivals, combined with the simulated g R , we can estimate the range of r reached Digital quantum simulation of quantum Rabi model parity dynamics in the degenerate-qubit case. (a) Parity dynamics of the ideal quantum Rabi model for qubit (green) and resonator (red) in coupling regimes: r = g R /ω R q = 0.1 (dotted), 0.5 (dashed) and 1.0 (solid). (b) Two-transmon, three-resonator circuit QED chip (detailed description in supplement [31]). (c) Sequence schematic for second-order Trotterisation. The rotating frame defining the simulated resonator frequency (ωr) is controlled via the QR bit-flip pulse phases. (d) Example simplified experimental pulse sequence for 5 Trotter steps followed by a photon parity measurement. (e-h) Measured qubit and photon parity dynamics for up to 60 Trotter steps, with the extreme dUSC regime in the centre decreasing to weaker USC near the edges. The data show clear Gaussian-shaped collapses for all r, along with the characteristic plateaus of the ultrastrong coupling regime. Qubit revivals are observed up to r ∼ 0.8, while photon parity shows clear revivals up to r ∼ 1.8. Slices are plotted in (f, h) for evenly spaced ω R r /g R between the red and blue dashed lines in (e, g), respectively. For r 1.5, some deviation from the expected revival time in the photon parity results from a small residual Kerr-type nonlinearity in the resonator [31] and is correlated with significant photon populations. Arrows in (f) and (h) show expected revival times for each slice.  [31]. Photon decay becomes increasingly critical in the USC regime, because even a single decay destroys the qubit-resonator entanglement, and losing a photon becomes increasingly likely for larger photon numbers. The qubit parity revivals rely on entanglement being maintained. This is supported by measurements of reduced qubit entropy, which show that the qubit state collapses to the mixed state, before displaying a revival in purity [31]. The resonator parity dynamics, however, are more robust to decay and provide a more direct measure of the dUSC dynamics. Photon parity collapses and revivals prove the field un-dergoes large-amplitude excursions through phase space even during a single cycle of the resonator period. The difference between qubit and photon parity dynamics is a signature of break-down in parity conservation, caused by resonator decay. We next directly explore the build-up of large photon populations (Fig. 2), another feature of USC dynamics that contrasts strikingly with the excitation-conserving dynamics expected under weak coupling. Using a Ramsey pulse sequence with small separation τ , the excitation probability in Q W becomes a measure of average photon number in the resonator [31]. The dynamic range and sensitivity of this photon meter are controlled via τ [Figs 2(a, b)]. Measured with a linear range of ∼ 0-8 photons [ Fig. 2(c)], the resonator displays the complementary build-up of photons which causes the collapse of qubit and photon parity, clearly demonstrating the violation of number conservation expected for the QRM. As with photon parity, clear oscillations can be seen out to r ∼ 1.8 [ Fig. 2(e)]. The large central feature appears to deviate from the expected trend, but is in fact due to photon number exceeding the dynamic range of the photon meter. To explore this region further, we extended the linear range to ∼ 0-20 photons using a photon meter with a non-centred refocussing pulse [ Fig. 2(d)] and simulated up to 90 Trotter steps (gt = 7.0π), allowing photon oscillations beyond 1.5 µs to be observed. This range operated at the limits of approximately uniform driving given the bandwidth of the 12 ns (4σ) Q W pulses. At r 2, the photon dynamics in Figs 2(c) and (e) are clearly skewed, causing the observed oscillations to deviate from the expected revival period T R (also observable in the photon parity [ Fig. 1(g)]. This results from a residual Kerr nonlinearity in R R inherited from the dispersively coupled ancilla qubit [38].
Exploring the resonator oscillations more quantitatively, the maximum photon number in each vertical (constant-r) slice [ Fig. 2(f)] compares well with the expected ideal behaviour. The discrepancy between the two curves in the overlapping region results from bandwidth limitations in the high-dynamic-range (HDR) photon meter and the limits in linearity of the number-toprobability mapping for Q W . Because of the sinusoidal conversion, the calibrated value at either end of the range compressed slightly towards the centre from the real photon number. The measurement saturates at the highest r even for the HDR meter, suggesting that we observe more than 30 photons (average) building up in the resonator for the strongest USC regions. Given the Poissonian statistics expected for coherent states, this accesses a resonator subspace of dimension ∼ 40.
Combining the parity measurement with coherent displacements from an external drive allows observation of resonator phase-space dynamics using direct Wigner tomography [9,35]. Figure 3(a) shows unconditional maximum-likelihood tomograms (ignoring the state of Q R ; see Methods) measured after each Trotter step with r ∼ 0.9 (full movie available [31]), with the full trajectory obtained from two-dimensional double-Gaussian fits of the raw data. The resonator state displays the clear signatures of USC dynamics, first separating into two distinct Gaussian (coherent-state) peaks which follow opposite circular trajectories before re-coalescing at the origin. The peaks do not return perfectly to the origin because of photon decay, in agreement with a numerical simulation at g R /2π = 1.79 which includes T 1,r = 3.5 µs (green curves).
By capturing the complete resonator quantum state, the Wigner function also enabled the demonstration of coherence in dUSC dynamics, by contrast with photon parity and number measurements, which were largely insensitive to coherence. Observing this requires correlating the resonator and qubit states, because the coherence is stored in entanglement. We did this in two ways. First, we measured the Wigner function after 10 Trot-ter steps for r ∼ 0.9 with Q R initialised in states |0 , |1 , |+ and |− [Figs 3(b-e)]. This showed that the resonator and qubit were correlated, consistent with the expected Bell-cat entanglement. Second, we ran the simulation for r ∼ 0.9 and 2.1 (8 Trotter steps) with the qubit prepared in the excited state, conditioning the Q W measurement on the state of Q R in the σ z basis (Fig. 4). For the expected Bell-cat state, an outcome of 0 (1) for Q R leaves the resonator in an odd (even) Schrödinger cat state (|α ∓ |−α ). Numerical modelling shows that only in the USC regime is negativity in the Wigner function observed for both Q R measurement outcomes. The negative regions observed in all the Wigner functions demonstrate nonclassicality for all resonator cat states, which arises from coherence in the underlying Bell-cat entanglement. Reduced visibility is again caused primarily by photon decay, but also by single-shot readout infidelity (here, ∼ 85-90%) and experimental drift over the long measurements. These different measurements provide clear evidence of qubit-resonator entanglement arising from coherent dUSC dynamics.
Finally, by detuning the qubit frequency during the JC half of the Trotter steps [ Fig. 1(c)], we also experimentally simulate dynamics for the nondegenerate-qubit case of the QRM for effective qubit frequencies g R /ω R q ∼ 4, 2 and 1 (Fig. 5). These regimes access the full complexity of QRM dynamics. Comparison with numerical modelling of the ideal QRM (no decoherence) shows that the experimental simulations capture many features of the ideal dynamics, even up to r 1. However, the main deviation from the degenerate-qubit case occurs primarily in the regime ω R r ω R q [39], making the measurements here even more susceptible to photon decay. Numerical modelling of the digital QRM simulation including the measured T 1,r confirms that simulation fidelity is primarily limited by decay.
Demonstrating stabilisation by decreasing step sizes will be an important part of validating the behaviour of future complex digital simulators achieving quantum advantage [40]. In the supplementary material, we showed that using second-order Trotterisation and decreasing the Trotter step size both significantly improved performance [31]. This indicates that the simulation is not limited by an error-per-gate noise floor as in previous circuit QED simulations [7], and enables us to linearly increase the number of Trotter steps for increasing simulated time, rather than keeping the number fixed [3,6,7]. This is a crucial step towards the quadratic scaling needed for universal quantum simulation [1]. In combination, these achievements advance the state of the art in solid-state digital quantum simulation, bringing circuit QED simulators to a level previously attained only in trapped-ion systems [5].
A QRM simulator has direct advantages over natural USC systems. Although USC can lead to ground-state entanglement and significant ground-state photon pop- (a) Selected frames from a "movie" (measured over ∼ 40 hours) showing the phase-space evolution of the resonator reduced state for r ∼ 0.9, with the final panel showing the full trajectories determined from 2D double-Gaussian fits to the raw data (the full movie is provided in the supplemental material [31]). Plotted tomograms are maximum-likelihood reconstructions of direct Wigner tomography measured data with a systematic phase correction (see Methods). When the effective drive on the intracavity field created by the Rabi interaction has a strength comparable to the resonator's natural frequency (i.e., g R ∼ ω R r ), this drive is able to create a significant displacement of the cavity field before the phase-space rotation caused by ω R r brings the field back towards the origin. This effect is observed clearly here in the creation of two well-resolved, rotating peaks and subsequent re-coalescence which are characteristic signatures of dUSC dynamics. Deviation from the ideal circular trajectories (orange curves) arises from photon decay. The measured trajectory shows excellent agreement with a numerical Trotter simulation at g R /2π = 1.79 MHz which includes resonator T1,r = 3.5 µs (green curves). From the fits, we calculate an estimated Wigner function width σ = 0.526 ± 0.003, instead of the predicted 0.5, indicating a displacement calibration error of ∼ 5% [31]. Background noise arises from phase instability of microwave sources and frequency stability of the Wigner qubit over the long measurement.
Nonclassical Schrödinger cat states of the Rabi resonator from conditioned dUSC-driven entanglement. The plots show resonator Wigner functions from maximum-likelihood state reconstructions for two different ultrastrong coupling strengths with g R /ω R r ∼ 0.9 (top, 10 Trotter steps) and g R /ω R r ∼ 2.1 (bottom, 8 Trotter steps), conditioned on measuring QR in |0 (left) and |1 (right). The regions of negativity and visibility of several fringes between the well-resolved coherent state peaks are clear signatures of nonclassicality in the Rabi field mode and demonstrates the coherence and entanglement of the underlying qubit-resonator state. Combined with the qubit conditioning shown in Fig. 3, observing clear cat states for both outcomes of the QR measurement is a clear signature of coherent USC dynamics.
ulations, these potentially interesting ground states are not readily accessible in natural USC systems [14,33,41] without the ability to rapidly (nonadiabatically) tune or switch off the ultrastrong coupling. In systems where the coupling reaches many gigahertz, tuning system parameters on this timescale represents a significant technical challenge [16,17]. In our simulator, however, cavity photons are always real (not virtual), detectable and usable, and it is straightforward to nonadiabatically tune system parameters to implement quantum quenches [42]. This makes a circuit QED chip with natural JC interactions an ideal platform to explore the preparation of interesting ground states in future experiments. The challenge is that the simulator decay processes differ from those in a natural USC system and do not move the system towards the USC ground state [11]. This highlights the need to improve T 1,r so that photon decay does not limit the dynamics. It should be possible to improve T 1,r ten-fold using novel processing methods [43].
Finally, the phase technique we have developed to define a rotating frame via single-qubit pulses introduces a precise and flexible paradigm for engineering artificial Hamiltonians which can be applied across architectures such as trapped ions and cold atoms [5,26,27]. In combination with the number of Trotter steps demonstrated, the technique will allow accurate simulation of the time-dependent Hamiltonians [5,7,44] required to perform adiabatic preparation of USC ground states. It is therefore ideally suited for exploring novel quantum phase transitions relying on extreme coupling regimes recently identified for the QRM [26,45,46]. Further, by extending to small-scale Dicke-model systems [23,25], it will avoid the problem of additional nonlinear evolution terms [25] which have been suggested to prevent the onset of a long-predicted superradiant phase transition in a range of physical systems [12,13,33,47].

Phase-controlled Trotterisation of the quantum Rabi model
In the digital QRM simulation proposed in Ref. [23], the effective parameters of the simulated Rabi Hamiltonian are g R = g, ω R r = 2∆ r and ω R q = ∆ JC q −∆ AJC q , where ∆ r = ω r − ω RF and ∆ q = ω q − ω RF are defined relative to a rotating frame. This rotating frame is essential to reaching the USC regime with weakly anharmonic transmon qubits, by allowing us to tune the simulated ω R r and ω R q . Typically, the frequency of a rotating frame is set by a physical generator or drive signal that defines a rotation or a measurement basis. In the digital simulation, the rotating frame is still abstract, since no drive is used to induce an interaction. Here, we describe a method we have developed for controlling the frequency of the rotating frame which is simple, high-resolution and flexible.
In the experiment, the frequency of the rotating frame is defined by the rotation axes of the bit-flip pulses (set by the pulse phase) that convert every second JC interaction into an effective AJC interaction. The qubit, however, is driven at its bottom sweet spot, around 1 GHz below the resonator. The drive generator phase therefore changes rapidly by comparison with the target rotating frame and it is necessary to reset the generator back in phase with the target rotating frame each time a pulse is applied.
We now derive the relation between the bit-flip pulse phases and the rotating-frame frequency. The symmetric, second-order Trotter step for the digital QRM simulation is: where U JC (τ ) = exp(−iH JC τ / ) and an arbitrary AJC step is defined by the phases used to set the rotation axes φ 1,2 of the bit flips R φ (π). Writing the JC Hamiltonian in the The cases implemented are g R /ω R q ∼ 4 (top), ∼ 2 (middle) and ∼ 1 (bottom), with the plots showing measured qubit dynamics (left), numerically simulated dynamics of a Trotterised QRM with the measured T1,r ∼ 3.5 µs included (centre), and ideal Rabi dynamics (right). The results illustrate that the nondegenerate-qubit dynamics do not deviate significantly from the degenerate-qubit case in the regime where ω R r ω R q . The measured dynamics exhibit many qualitative features in good agreement with the ideal QRM and show excellent agreement with the numerical Trotter simulation with decay, indicating that the fidelity of the measured results to the ideal case is limited primarily by resonator decay.
rotating frame of the resonator, and using the identity Next, noting that ∆φ = πω R r τ 1 if τ 1/ω R r , and providing the Trotter condition = gτ 1 is fulfilled, we can combine exponentials using Trotter approximations to give first: and then the full Trotter step So far, we have considered arbitrary φ 1 and φ 2 . In the experiment, however, we keep ∆φ constant for all sequential pairs of bit flips. Specifically, for the nth Trotter step, the two phases are φ 1 = φ 0 + (2n−2) ∆φ and φ 2 = φ 0 +(2n−1) ∆φ, where the choice of φ 0 has no effect on the dynamics. Setting φ 0 = 3∆φ/2 gives φ Σ = 4n∆φ, and the nth Trotter step can be rewritten in terms of a frequency ω 0 = 2∆φ/τ and a time t n = nτ : which corresponds to an effective Hamiltonian: Until this point, the calculation has been carried out with both qubit and resonator in a frame rotating with the resonator. We now transform H eff into a rotating frame where both qubit and resonator are rotating at frequency Thus, in the new rotating frame, the final effective Hamiltonian implemented by the Trotterisation is: This completes the mapping of the phase-controlled Trotterisation into the form of a simulated Rabi Hamiltonian and we can now identify the effective simulated parameters g R = g, ω R q = ∆ JC q-r and ω R r = −ω 0 = −2∆φ/τ . It is worth emphasising here that, by controlling the phase difference between successive bit-flip pulses on the qubit, we are able to define the rotating frame frequency ω 0 , and hence the effective resonator frequency in the simulated Hamiltonian.

Trotter step
For a second-order Trotter step with simulated time τ , the Trotter step consists of 3 flux pulses (τ /2, τ and τ /2) and 2 single-qubit rotations with buffers separating the different gates. Adjacent τ /2 flux pulses from neighbouring Trotter steps are implemented as a single flux pulse of length τ . Each flux pulse was followed by a 5 ns phase-compensation flux pulse [31]. For most of the data presented in this work, the simulated τ = 20 ns. The qubit drive pulses on Q R were 16 ns total duration (4σ) and the pulses buffers were 10 ns. The total Trotter step for τ = 20 ns was therefore τ step = 122 ns. In addition to the drive-pulse phase advance required to define ω R r , another linear phase advance ∆φ = (ω drive q − ω r )τ step /2 is required to compensate the rapid rotation of the qubit drive with respect to the resonator frequency.

Qubit control
Qubit rotations were implemented using DRAG pulses [48,49], with a Gaussian envelope in the X quadrature and a derivative-of-Gaussian envelope in the Y quadrature. The 4σ pulse durations were 16 ns for Q R and 12 ns for Q W . The performance of the Trotter sequences, which contained up to 180 bit-flip pulses, was very sensitive to details of the Q R pulse calibrations. In particular, the drive amplitude was calibrated using a sequence of 50 π-pulse pairs preceding a single π/2 pulse. All parameters were typically calibrated just before launching a long measurement. The drive amplitude was intermittently recalibrated during the scans. Because only 2 or 3 pulses were applied to Q W for the photon measurements, it was optimised using the AllXY sequence [50] of 21 combinations of two σ x and σ y rotations (either π/2 or π). The frequency of Q W was regularly calibrated during photon measurements using Ramsey sequences.

Wigner tomography reconstructions
Tomograms shown in Figs 3 and 4 are maximum likelihood reconstructions [51,52] of the resonator quantum state from direct Wigner tomography measurements [9]. The Wigner function at a phase-space position α is: where ρ r is the resonator density matrix, Π = n (−1) n |n n| is the photon parity operator and D (α) is the coherent displacement operator. For each measured α, we calculated M α = D (α) ΠD † (α) using an operator dimension much larger than the largest |α| 2 in the measured phase space, to avoid edge effects when calculating D (α). The M α were then truncated to a maximum photon number sufficient to capture all of the reconstructed state, but small enough to allow fast reconstructions and ensure an informationally complete set of operators (n max = 12 and 8 for tomograms in Figs 3 and 4, respectively). The maximum likelihood reconstruction was carried out using convex optimisation [53,54]. In Fig. 3, a systematic phase correction was applied to the density matrices to correct for a miscalibration of the resonator drive phase used in the coherent displacement. Finally, the reconstructed density matrix was then used to calculate the plotted Wigner functions.

SUPPLEMENTARY INFORMATION
This supplement provides experimental details and additional data supporting the claims in the main text.

Experimental setup
The sample and low-temperature microwave components were mounted inside magnetic and infrared radiation shielding consisting of two layers of cryogenic mu metal around a layer of aluminium, with an internal layer of copper foil coated in a mixture of silicon carbide and Stycast (2850 FT) [1]. Microwave coaxial cables are connected to the PCB-mounted chip via non-magnetic SMP connectors (Rosenberger).
The qubit drive and read-out tones are sent through two dedicated feedlines which are connected via a short coaxial cable off-chip. The input line for the qubit drives is filtered at the mixing chamber with 30 dB cold attenuation, a small home-built inline eccosorb filter and a 10 GHz low-pass filter (K&L 6L250-10000/T20000-0/0). [The resonator input line filter is 8 GHz low-pass (K&L 6L250-8000/T18000-0/0).] The output line passes through two 3-12 GHz isolators (Pamtech CWJ1019K) and a circulator (Quinstar CTH0408KCS) mounted above the mixing chamber on the way to a 4-8 GHz cryogenic HEMT amplifier (Low-Noise Factory LNF-LNC4 8A), two room-temperature amplifiers (Miteq AFS3-04000800-10-ULN, then AFS3-00101200-35-ULN-R), RF demodulation (Marki 0618LXP IQ mixer) and amplification, and finally digitised in a data acquisition card (AlazarTech ATS9870). The flux-bias lines are filtered at the mixing chamber with 1.35 GHz low-pass filters (Minicircuits VLFX-1350) followed by home-built eccosorb filters. All input lines are thermalised with 20 dB attenuators mounted at the 4 K plate. The microwave input lines and output line are connected to the fridge through a DC block.
Qubit and resonator drive pulses are created via singlesideband modulation with IQ mixers and generated by two arbitrary waveform generators (AWGs; Tektronix AWG5014). We use a 3-7 GHz IQ mixer (Marki 0307MXP) for the resonator and two custom-built 4-8.5 GHz IQ mixers (QuTech F1c: DC-3.5 GHz IF bandwidth) for the qubit drives. The qubit drive pulses were amplified by a high-power (35 dB) microwave amplifier (Minicircuits ZV-3W-183) before passing through a 5.5 GHz low-pass filter (Minicircuits LFCN 5500+) to minimise amplifier noise at the readout resonator frequencies.
Most microwave units receive a 10 MHz reference from a microwave generator (Agilent E8257D) via a homebuilt distribution unit. However, the generators used for driving Q R and R R (R&S SGS100A) synchronised directly via a 1 GHz reference. This was critical to achieving the phase stability required to measure R R Wigner functions during measurement runs lasting up to 40 hours. The frequencies for these two generators were also always set to a multiple of the trigger repetition rate (5 kHz), to ensure a stable phase relationship. For phasesensitive measurements, a 500 MHz scope (Rigol DS4034) monitored the relative trigger timing between the master and slave AWGs to select consistent delay configurations between the AWG outputs.
Home-built low-noise current sources mounted in a TU Delft IVVI-DAC2 rack provided precision DC bias currents for flux tuning of the qubit frequencies. The DC bias for Q R was combined with the amplified output of one channel of the master AWG (the same as used for generating Q R drive pulses) using a reactive bias tee (Minicircuits ZFBT-6GW+). The flux pulses from the AWG were amplified using a home-built 2 V/V flux-pulse amplifier.

Device fabrication
The device was fabricated using a method similar to that of Ref. 2, but with several specific improvements: 1. The transmon design includes a rounded spacing between the shunt capacitor plates [ Fig. S2(a)] to avoid the regions of high electric field which can increase sensitivity to interface two-level fluctuators [3].
2. The flux-bias line was centred between the transmon capacitor plates to symmetrise the capacitive coupling with the goal of decoupling the qubits from possible decay-inducing effects of voltage noise fluctuations on the flux-bias lines.
3. As in our previous work [2], the transmon qubits were patterned with niobium titanium nitride (NbTiN) capacitor plates to further reduce susceptibility to noise from two-level fluctuators. Prior to evaporation of the aluminium (Al) junction layers, a short hydrogen-fluoride (HF) dip removed surface oxides to facilitate a good contact between the evaporated Al and NbTiN thin film. To avoid contact problems caused by unwanted etching into the silicon substrate during patterning of the NbTiN, we: 1) optimised the reactive-ion etch (RIE) recipe and duration to minimise the substrate etch and eliminate underetch (under the NbTiN); and 2) introduced a narrow bay in the NbTiN fingers at the contact point to create a softer etch for more reliable contact [ Fig. S2(c)].
4. The junction development process and doubleangle evaporation parameters were optimised to improve the reliability of the very small junction sizes needed for the asymmetric qubit [ Fig. S2(b)]. FIG. S1. Experimental schematic showing the connectivity of microwave electronics and components in and outside the dilution refrigerator. The sample mounted below the mixing chamber typically remained at around 30 mK. Qubit and resonator drive lines and flux-bias lines were thermalised and attenuated at the 4-K and 30-mK stages and were low-pass filtered before arriving at the sample. The qubits and resonator drive pulses were generated by AWGs and IQ mixers. Home-built low-noise current sources provided DC bias currents for qubit frequency tuning, which were combined with fast frequency-tuning bias pulses using reactive bias tees. AWG markers provided the gating for pulse-modulated measurement pulses.
Device operating parameters and qubit performance Figure S3(a) shows the frequencies for the two qubits and three resonators on the device as a function of the applied qubit flux in units of the flux quantum Φ 0 = h/2e, along with the operating points for both qubits during the quantum simulation experiments. Measured device parameters are summarised in Tab. S1. Qubit T 1 , T 2,echo and T * 2 decay times are shown as a function of qubit frequency in Fig. S3(b,c,d).
At the operating point, the Rabi qubit Q R was designed to sit below the resonator R R and be pulsed up into resonance with it to avoid continually crossing the resonator with the Q R 's 1-2 transition during the  [4]. T * 2 s reported here were measured by fitting a decaying double sinusoid to a long, beating Ramsey signal and represents the underlying coherence of the qubits. At the operating point for QW far from the sweet spot, no beating was observed in the Ramsey measurements.
long flux-pulse sequence. Because of significant protocol times and two operating points, an asymmetric qubit design with two flux-insensitive "sweet" spots was used for Q R [5], with drive pulses applied at its bottom sweet spot. The first-order flux insensitivity at this point also mitigated some of the impact of rapid, long-range fluxpulsing on the qubit pulse tuning. The maximum and minimum frequencies for Q R in the final cooldown were 6.670 GHz and 5.451 GHz, respectively.
The asymmetric design also minimised the stringent challenge of targetting the qubit frequency to resonator closely on the scale of the very small coupling frequency. Ideally, the resonator would have been closer to the qubit top sweet spot to maximise phase coherence also during the interaction pulses. However, with the asymmetric design, the reduced flux gradient relaxes this constraint. With an asymmetry parameter of α = (E J,max −E J,min )/(E J,max +E J,min ) ∼ 0.68, the Ramsey time T * 2 for Q R did not typically drop below a few microseconds, even at the positions with steepest flux gradient.
The asymmetry of Q R was smaller than targetted, with the result that the bottom sweet spot was also lower in frequency than intended. The ancilla qubit Q W (a standard symmetric transmon) was therefore operated around 650 MHz below its own maximum-frequency sweet spot of 5.653 GHz. At this operating point, its T * 2 was typically 1.5 µs. Because we were able to drive Q W and achieve good photon-sensitive operation at this lower position, we chose not to rapidly tune its frequency up to the sweet spot to perform the photon meter measurements.
To identify the flux operating point that positioned Q R precisely at the bottom sweet spot, we applied the following procedure. We first decoupled the applied DC qubit fluxes, applying the appropriate linear correction to compensate for flux cross-talk. Then, after positioning Q W roughly at its selected operating point, we applied a simple excitation swapping sequence for Q R with R R with fixed swap time (near a full swap) and varying amplitudes of positively and negatively directed pulses. Finally, we varied the applied flux on Q R and identified the operating point as the symmetric flux point where the qubit hit the resonance for positive and negative pulses of equal amplitude. We were able to identify this point to 1 part in 5000. Because the precise choice of operating frequency for Q W was not critical, any slight shift in frequency due to residual DC cross-talk remaining after the flux decoupling measurements was unimportant.

Calibration of the flux distortions
Implementing the digital Trotterisation of the Rabi model proposed in Ref. [6] required tuning the qubit frequency with a long series of square interaction pulses. To achieve this, it was necessary to compensate for the filtering effects of electronics and microwave components in the line (Fig. S1) [7]. One of the particular challenges of an experiment using a long train (up to 10 µs) of very short pulses (10-20 ns) is that the system is sensitive to both short-and long-time pulse distortions. These effects included the intrinsic bandwidth of the AWG and the flux-pulse amplifier, the high-pass characteristics of the bias tee, a range of low-pass effects including the Minicircuits and eccosorb filters and filtering from the skin QR is coupled to RR near its shorted end in order to achieve the required small coupling g. (b, c) Josephson junctions are contacted to the NbTiN SQUID loop fingers using small bays to achieve better contact. In (b), it is possible to see the large asymmetry in junction size, with a zoom on the small junction in (c). effect of the coaxial cabling, pulse bounces at impedance mismatches, as well as more intangible effects such as transient decays in step responses. Subject to the system operating in a linear regime (e.g., the AWG operating in a comfortable amplitude range), this could be achieved by applying predistortions to the target fluxing sequence. Figure S4 illustrates the calibration process used in this experiment. Rather than building a single, comprehensive model for all flux distortions, we took a divideand-conquer approach, applying a series of corrections to compensate individual effects. For processes outside the fridge, we calculated the required compensations by directly measuring the system step response using a fast oscilloscope (R&S RTO1024, 10 Gs/s sampling rate and 2 GHz bandwidth). We applied predistortion corrections sequentially, at each step correcting the longest-time behaviour and zooming in to shorter time scales once the longer-time response is successfully corrected. Once measuring through the fridge, we optimised on the shape of the two-dimensional flux-pulse resonance, the so-called "chevron". Again, we typically focussed initially on correcting the coarse features before zooming in to finer details.
The procedure we used to calculate the external cor-rections was: • Sample a measured step response at a period τ : • Invert H to find the transfer matrix of the so-called predistortion kernel and calculate the step response of the predistortion kernel as Hu[n], where u[n] is the discrete Heaviside function. This numerical matrix inversion step limits the length of the step response that can be treated in this way. The sampling period τ is chosen to ensure the sampled step response covers the region of interest.
• Fit the numerically inverted kernel step response using a simple functional form which can then be used to construct a high-resolution predistortion kernel (the impulse response calculated as above from a high-resolution step response). The downsampling of the step response reduces the fit function dependence on high-frequency effects. For each step, we varied the sampling period to check that the fit parameters were relatively robust to details of the sampling. Figure S4(a) shows the step response from the AWG measured after the home-built flux-pulse amplifier (see Fig. S1), with a zoom into the top of the step in (b). In this case, the longest-time response was actually an effectively linear ramp over the long step response. Here, we used a slightly modified procedure to the one above, fitting a linear function directly to the measured step response. Using Laplace transformations, it is possible to show that a step response with a linear ramp, (1 + αt) u(t), can be corrected using a predistortion kernel with an exponentially decaying step response exp(−αt) u(t). After this linear correction, we then implemented a series of three corrections with "exponential-approach" predistortion step responses of the form (1 + α exp(−t/τ ))u(t) with τ values between 5 µs and 500 ns (various amplitudes), determined using the above procedure. Figure S4(c) shows the corrected step function measured after applying the four initial corrections. The small but distinct sawtooth structure in the otherwise flat step response is due to the vertical resolution of the AWG.
After correcting for distortions from the AWG and flux-pulse amplifier, we measured the step response after the bias tee, at the fridge input. Figures S4(d, e) show the measured step response and sampled predistortional Note that, at the sweet spots, measured qubit T * 2 times here are limited by slow frequency-switching processes in the qubits such as quasiparticle tunnelling [4].
kernel step response calculated using the above procedure (with τ = 50 ns). The high-pass characteristics of a reactive bias tee's RF input naïvely predict a kernel step response with a full initial step followed by a continually increasing linear voltage ramp. From Fig. S4(e), however, it is clear that the kernel step response is not completely linear. We instead fit the step response to a quadratic form and proceed as above. The step response measured after compensating for the bias tee is shown in Fig. S4(f).
Inside the fridge, we calibrated the flux-pulse predistortions to optimize the shape of the flux chevron [ Figs S4(g-i)], which probes the excitation-swapping exchange interaction between qubit Q R and resonator R R as a function of flux-pulse amplitude and interaction time. When the qubit is exactly on resonance, the swapping interactions are expected to be slowest and strongest. As it moves off resonance, the oscillations speed up and reduce in amplitude. Interestingly, despite the good performance of the bias-tee correction when measured outside the fridge, the chevron measured with the same corrections [ Fig. S4(g)] showed a clear ramp in the start of the interaction signal (the lateral skew), consistent with an under-compensated bias tee. We do not understand the cause of this discrepancy, but corrected it empirically by adjusting the linear coefficient of the bias-tee correction. The chevron measured after optimising this correction (final linear coefficient corresponded to a time constant τ = 9.7 µs) showed the characteristic asymmetric signature of low-pass filtering from the skin effect [ Fig. S4(h)]. This was corrected by applying a kernel numerically calculated from a step response of the form (1 − erf(α 1GHz /21 √ t + 1)) u(t) [8], using α 1GHz = 1.7 dB. Finally, we implemented another series of exponential-approach kernels with values of τ between 1500 ns and 30 ns, to achieve the result in Fig. S4(i).

Calibration of the photon meters
Using a photon meter based on a Ramsey sequence's sensitivity to qubit frequency and Q W 's dispersive frequency dependence on resonator photon number allows   detection of average photon number with controllable sensitivity and dynamic range. Suppose the resonator is in the state ψ = j α j |j . To implement the photon meter, we apply a Ramsey pair of π/2 pulses with pulse separation τ on Q W at a frequency Ω d W = Ω 0 W − d2χ, corresponding to the d th photon peak. Different photonnumber frequency components accrue different phases during the variable delay between pulses, given by θ j = (j−d) 2χτ . By driving first around σ x and then around σ y , the d th photon term ends up on the equator of the Bloch sphere. Measuring the excitation of Q W then gives a measurement probability Provided τ is chosen such that θ j is small for all photon components j present in the photon state, Increasing τ therefore increases the sensitivity of measured probability to average photon number, but decreases the accessible range of photon numbers for which the linearity condition sin θ j ≈ θ j holds. An accurate calibration of the photon meter also requires an accurate calibration of the single-photon dispersive frequency shift 2χ and Q W 's zero-photon frequency (which determines Ω d W ). Here, we describe a self-consistent calibration of our photon meters which does not rely on quantities derived from other measurements, such as spectroscopy, and relies primarily on knowing drive-pulse frequencies, probably the most accurate control parameter we have in the experiment. At each stage, we first calibrate Q W 's zero-photon frequency using a standard Ramsey sequence. With the performance of Q W at the operating point (dephasing time T * 2 ∼ 1.5 µs), we routinely achieved frequency accuracy better than 10 kHz.
To calibrate the single-photon dispersive shift [sequence shown in Fig. S5(a)], a calibrated SWAP pulse on Q R transfers an excitation into R R , before the resonator photon number is probed via Q W . The single-photon excitation in R R dispersively shifts the frequency of Q W by 2χ. Driving Q W at the calibrated zero-photon frequency around σ x and then σ y , the correct parity condition corresponds to the point where the curve crosses 0.5 excitation probability [ Fig. S5(d): 383 ns wait time]. This measurement is robust to both the relatively short resonator photon decay time T 1,r ∼ 3.5 µs and the short dephasing time of Q W at its operating point (T * 2 ∼ 1.5-1.8 µs at ∼ −650 MHz detuned from its top sweet spot), because these processes both reduce the visibility of the curve, but not the oscillation period, and therefore do not affect the value of the crossing point. The zero-photon frequency calibration is the main limitation, because that calibration limits the accuracy with which the crossing point represents the correct delay time between π/2 pulses.
The wait time identified above specifies the time between the end of the first pulse and the beginning of the second required to realise a photon parity measurement, but this does not account for the finite pulse duration. To calibrate the effective value of τ , we fix the pulse separation and sweep the frequency of the Q W drive generator this time without loading any photons into the resonator [ Fig. S5(b)]. For a pulse separation of 383 ns, the effective τ is ∼ 398 ns [ Fig. S5(e)]. Note that the difference here is not quite the same as the drive pulse width used in the experiment (4σ = 12 ns). This value of τ is related to the dispersive shift of R R on Q W in the usual way: τ = π/2χ, giving 2χ/2π = −1.26 MHz. Note that, when used directly as a parity meter, the read-out of Q W was calibrated using a parity pulse pair either with the usual phase on the second pulse, or a phase shifted by π radians. This accounted for the reduced parity visibility from the short T * 2 of Q W at its operating point and helped to track any fluctuations in the correct parity extremes as a result of drift in qubit frequency and T * 2 . Figures 2(a, b) show the pulse sequences for two different photon meters used in the experiment, one with the standard Ramsey sequence [calibrations in Figs S5(f, g)] and one an unbalanced "echo"-like sequence with an offcentre refocussing pulse (calibrations not shown). The mapping between average photon number and qubit excitation is approximately valid provided the phase advance/delay is less than 30 degrees, which corresponds to a qubit excitation of 0.25. We select the appropriate Ramsey pulse separation by driving the qubit at the frequency corresponding to the mid-point of the desired range (here, the 4-photon position), calculated from the dispersive shift and the calibrated zero-photon frequency, and choosing the separation which gives the target excitation probability of 0.25 [ Fig. S5(f)], here 4 ns. The effective τ was calibrated, as above, to be ∼ 19 ns. Moving to the smaller τ necessary for a higher photon number dynamic range requires frequency refocussing. Ultimately, the main limitation to the range achievable with such a photon meter is set by the bandwidth of the drive pulse.
We used a photon number meter calibrated using the above procedure to follow the excitation-swapping oscillations of a vacuum-Rabi exchange between Q R and R R , plotted as a function of the duration of the flux pulse on Q R [ Figure S5(h); sequence in Fig. S5(a)]. The drifting baseline results from pulsed flux cross-talk between Q R and Q W . To correct this, we repeated the same measurement without initially exciting Q R in order to avoid exciting photons in R R [Fig. S5(c)]. This curve was compensated by adjusting the drive phase of the second Ramsey pulse in the photon meter (on Q W ), leading to the compensated measurement in Fig. S5(j). To maximise the sensitivity of the cross-talk calibration, during the calibration, Q W can be driven at the zero-photon frequency, which then places the expected "null" measurement result on the equator of the Bloch sphere. A modified version of this procedure can be carried out for all flux-pulse sequences of interest. Note that cross-talk compensation was also necessary to ensure an accurate calibration of the parity condition in Fig. S5(d) above.

Calibration of Wigner tomography
We implement Wigner tomography using the direct method of Ref. 9. After the algorithm part of the pulse sequence [represented in Fig. S6(a)  resonator photon state before the usual parity readout pulses. The phase-sensitive resonator drive tone is created via single-sideband modulation in an IQ mixer. We calibrate the drive frequency and amplitude using the already calibrated photon meter (Figs S6(b, c), respectively). The drive amplitude is calibrated in the middle of the linear range, where we expect the best performance. Figure S6(c) illustrates the breakdown of the linear mapping between average photon number and Q W excitation probability both towards the edge of the linear regime and above the range, as the higher photon components wrap around in phase. In the digital QRM simulation, for phase-sensitive Wigner tomograms (e.g., Figs 3 and 4), it was critical to maintain phase stability between the drives on Q R and R R during the measurement. To achieve this, the two microwave generators were synchronised using a 1 GHz reference, with frequencies set as a multiple of the 5 kHz experimental repetition rate. Figure S6 shows one-and two-dimensional Wigner tomograms of a zero-photon (d, e) and one-photon (f, g) state (scaled in terms of photon parity). The maximum visibilities in Figs S6(f, g) do not reach the expected values, because these tomograms were measured without an accompanying full set of parity meter calibrations. However, the radial symmetry observed in these tomograms demonstrates the correct behaviour of the coherent resonator drive.
The curves in Figs S6(d, g) show fits to the data of a classical mixture of zero-photon and one-photon Wigner function cross-sections, with a free x-axis scaling parameter has been included in the fits. These fits demonstrate that the measured tomograms agree well with theoretical expectations, subject to an x-axis scaling error of ∼ 5%. That is, the fits indicate that the amplitude calibrations result in a small systematic overestimate in displacement by 5%. This also agrees with two-dimensional double Gaussian fits of individual frames of the unconditional Wigner movie in Fig. 3(a) of the main text, which give an average Gaussian widthσ = 0.526 ± 0.003, compared with the expected value of 0.5.

Analog vs Digital Jaynes-Cummings Dynamics
Simple modelling of the Trotterised version of the full Rabi model shows that high-quality simulations require both slow dynamics and short Trotter steps (i.e., fast flux pulsing). Such an experiment is sensitive to both short-time and long-time effects in the flux-pulse shaping. A simpler experiment which verifies the performance of this flux pulsing is to implement a digital simulation of the standard Jaynes-Cummings (JC) interaction underlying the standard excitation-swapping experiments demonstrated with single flux pulses (Fig. S4).
In the standard continuous-wave (single-pulse) version of a JC excitation-swapping interaction, resonance between the qubit and resonator frequencies gives rise to maximum visibility oscillations of the excitation moving between the two components. When detuned, the different phases accrued by the qubit and resonator during the interaction decrease the oscillation visibility, while increasing the oscillation frequency. This gives rise to the characteristic shape of the flux chevron. Significant care is required, however, to accurately reproduce the (analog) JC interaction with a digital pulse train. Figures S7(a, b) show analog and digital versions of the JC interaction (viewed through the qubit excitation) under otherwise identical conditions. The digital chevron shows a series of resonances which do not appear in analog measurements (not shown), and there is also no chevron visible at the natural resonance condition around 2.45 Vpp.
The new features relate to the extra "interaction off" times in the digital version. The regular spacing between neighbouring satellite resonances is around 50 MHz (after converting AWG amplitude to qubit frequency), which is the inverse pulse duration. During the interaction time, the qubit-resonator relative phase evolves as expected. However, in the "off" time between interaction pulses, the qubit accrues phase at a different rate, and will hence not have the required phase at the beginning of the next pulse for the interaction to pick up where it left off at the end of the previous pulse. Therefore, the necessary condition for observing a chevron feature at exactly the position of the natural resonance is that the qubit phase accrued (relative to the resonator) during the "off" time should be a multiple of 2π. The observation of multiple satellite resonances is a form of digital aliasing, where the interaction will build up constructively from pulse to pulse provided the relative phase accrued between qubit and resonator during the "on" time of the pulse again differs only by an integer multiple of 2π. However, this is an aliasing of the dynamics itself, not just an aliasing of the measurement, which could also occur in natural continuous-wave (CW) chevrons and would never lead to the observation of extra satellite peaks.
This pulsed interaction can also be viewed as a Trot-terised simulation of the CW interaction. While successive interaction pulses obviously commute with each other, they do not necessarily commute with the "off" pulses. The condition on qubit-resonator phase during the "off" pulse can be understood as the condition where the Trotter error vanishes, because the Hamiltonian term resulting from the qubit detuning coincides with the identity. The satellites arise because the phase contribution from the qubit detuning in the "on" pulse is identical if the frequency change matches a multiple of 2π phase.
To compensate for the phase error accrued in the qubit during the "off" pulses, we apply a 5 ns compensation flux pulse between interaction pulses. Using the fluxpulse amplitude which corresponds to the centre of the CW chevron, the amplitude of the compensation pulse was swept to identify the correct compensation point. In this way, very good agreement was achieved between the digital JC dynamics and the traditional analog version [Figs S7(d-i)]. The main differences are a slightly reduced visibility because of the increased experiment time, and a slighly lower effective coupling frequency (g/2π ∼ 1.8 MHz, instead of ∼ 1.95 MHz). The latter most likely arises from residual short-time pulse imperfections which do not contribute significantly to the long interactions in the analog form.

Trotter simulation with excited and ground initial states
In the degenerate-qubit case, when understood in terms of the cavity trajectories in phase space, it is clear that the structure of the dynamics of the full Rabi model with USC should not depend on whether the qubit starts in the ground or excited state. This contrasts with the JC interaction, where the |g, 0 state is decoupled from the rest of the system and the system will only undergo nontrivial dynamics if an initial excitation is loaded in the system. Indeed, in a natural USC system, if it were possible to turn the coupling on and off rapidly, it would be extremely interesting to watch an uncoupled-system ground state evolve into a state with excitations in the qubit and cavity. In this digital simulation, however, this is less satisfying, since the protocol in any case involves regularly injecting excitation into the system in the form of qubit flipping pulses. Most of the results reported here therefore take the more conservative position of initialising the system with an excitation, with the motivation that observing a difference between the simulated dynamics and what would be expected in a weakcoupling scenario could then only result from the simulated counter-rotating terms. Although there were some stability issues during the measurement with groundstate initialisation, there is nevertheless extremely good agreement between the two cases, for example with the timing of the revivals in both cases agreeing with the theoretical predictions. For this particular measurement of ground-state initialisation, qubit revivals are observed even out to r ≡ g R /ω R q ∼ 1.

Trotterisation performance vs Trotter order
As discussed already, initial modelling of a Trotterised Rabi simulation showed that unusually low qubitresonator coupling between Q R and R R was required to be able to achieve reasonable simulation fidelities given the hard bandwidth limitations of flux-based fast frequency tuning. This, however, required longer experimental times for the simulations, which in turn placed significant constraints on qubit and resonator coherence. Indeed, the shorter-than-anticipated resonator coherence time proved to be the biggest limitation. As a result, it was critical to use all available measures to minimize the Trotter error in our simulations, given the limits on the shortest achievable Trotter step sizes.
The accuracy of the Trotter approximation is set by the amount of non-commutativity between different components in the step [10]. While first-order Trotterisations [exp(A + B) ≈ exp(A) exp(B)] lead to Trotter errors that scale with single commutators (quadratically with simulation time), higher-order Trotterisations can be used to eliminate lower orders of Trotter error. For example, the symmetry of a second-order Trotterisation [exp(A + B) ≈ exp(A/2) exp(B) exp(A/2)] ensures that first-order error terms (related to single commutators) cancel, pushing the largest Trotter error terms out to third order in simulation time. For two-part Hamiltonians, however, second-order Trotterisation in practice only involves modifying the pulses in the first and last Trotter steps. All the results in the main text were obtained using a second-order Trotterisation. The plots in Fig. S9 illustrate that this was absolutely critical in order to extend the simulations deep into the ultrastrong coupling regime. The first-order and second-order Trotterisation agree reasonably well at r < 0.5, but behave fundamental differently at the higher values. The first-order simulation starts to show qualitatively different behaviour for relative coupling strengths r 0.5. In particular, only in the second-order case are the characteristic plateaus and revivals of the USC regime observable.
Trotterisation performance vs Trotter step size As illustrated in Fig. S9, the effects of Trotter error are most visible in the high r regimes, which is reasonable, considering that for low r, the Rabi model is well approximated by the JC model where the excitation-nonconserving terms (non-commuting with the excitation-conserving terms) do not play a significant role. This was also visible when studying the performance of the simulation as a function of the Trotter step size.
Measurements and numerical simulations show significant reduction in Trotter error as the number of Trotter steps over 1.2 µs increased from 24 to 60. The Trotter error shows up in two ways, namely the central features departing from the expected plateaus, and a tendency for the dynamical landscape to "break apart", even out into the lower coupling regimes. In the measured results and the simulation with decay, the fine details do not appear as strongly, but the effect appears to wash out the oscillation dynamics more rapidly. Only at the smallest step size are these effects absent from the measured results, and in the ideal simulations (without decoherence) there are even then central features which only disappear at a still smaller 10 ns step size. The measured results agree very closely with the numerical Trotter dynamics which include only the effect of photon decay, again highlighting that the primary limiting factor in our experiments was T 1,r . It is clear from these results that moving towards the smallest possible Trotter steps will be a key challenge for reaching quantum supremacy in complex quantum simulations.

Qubit entropy dynamics
In the Rabi model, as the resonator states separate, the qubit-resonator entanglement causes the reduced qubit state to collapse towards the maximally mixed state. entanglement is still present when the resonator states re-coalesce at the origin in phase space. While many possible uninteresting effects may cause an initial collapse in qubit purity, a revival in purity is a signature of entanglement with another system, in this case the resonator. After each Trotter step, a tomographically complete set of measurements on Q R was used to reconstruct its reduced state using maximum-likelihood tomography. We use the von Neumann entropy to characterise the purity of the reduced qubit state and observe revivals in qubit purity out to r > 0.8 [ Fig. S11(a)], consistent with the observed revivals in qubit parity. While the observed revivals shown in the slices [ Fig. S11(b)] appear smaller than the qubit parity revivals, in fact this is deceiving, resulting from the fact that purity (as with other entropy measures) is a quadratic function of the qubit population difference. The inset shows that the background noise of this signal is small and that the revivals are quite distinct. Moreover, plotting an appropriate square root of the entropy (not shown) shows that the revivals are consistent with the qubit parity case.