Analog quantum simulation of the Rabi model in the ultra-strong coupling regime

The quantum Rabi model describes the fundamental mechanism of light-matter interaction. It consists of a two-level atom or qubit coupled to a quantized harmonic mode via a transversal interaction. In the weak coupling regime, it reduces to the well-known Jaynes–Cummings model by applying a rotating wave approximation. The rotating wave approximation breaks down in the ultra-strong coupling regime, where the effective coupling strength g is comparable to the energy ω of the bosonic mode, and remarkable features in the system dynamics are revealed. Here we demonstrate an analog quantum simulation of an effective quantum Rabi model in the ultra-strong coupling regime, achieving a relative coupling ratio of g/ω ~ 0.6. The quantum hardware of the simulator is a superconducting circuit embedded in a cQED setup. We observe fast and periodic quantum state collapses and revivals of the initial qubit state, being the most distinct signature of the synthesized model.

Terms of the form e X Y e −X are calculated using the power series expansion of the exponential function, also called Hadamard lemma. For X = iω 1 tb †b , Y =b † , e i 2 ω1tσz andσ z clearly commute, while e i 2 ω1tσzσ x e − i 2 ω1tσz =σ + e iω1t +σ − e −iω1t (10) using e diag(a,b) = diag(e a , e b ). Terms in the last line in Supplementary Equation (6) are omitted in the following within the RWA, valid for η 1 /2ω 1 1. This is a good approximation as η 1 is bound in the experiment by the qubit anharmonicity |α| 350 MHz.
The η 1 term in the transformed Hamiltonian Supplementary Equation (6) is the most significant term, which justifies to move to its interaction picture. Under the constraint η 1 ≡ ω 1 − ω 2 , the Hamiltonian in the interaction picture becomes sin 2 η 1 t cos η 1 t + i sin 2η 1 t cos η 1 t − i sin 2η 1 t − sin 2 η 1 t Performing a time averaging in the spirit of a RWA casts Supplementary Equation (11) into the desired form of the quantum Rabi HamiltonianĤ with ω eff ≡ ω − ω 1 and noting that η 1 η 2 . In the experiment, the applied Rabi drives parasitically couple to the bosonic mode since the elements are degenerate during the simulation experiment and spatially close by in the circuit. This effect is accounted for by an additional term in the laboratory frame Hamiltonian, Supplementary Equation (3). η r denotes the effective amplitude of the parasitic drive and we take into account the effect of only the dominant Rabi drive frequency ω 1 . Performing the transformation U yields Terms rotating with e ±2iω1t ore omitted, resulting in a time independent drive term that is added to the Hamiltonian (12) in the rotating frame. The effective Hamiltonian including the parasitic drive term readŝ By applying the unitary displacement transformation we can cast the Hamiltonian (15) in the original form of the quantum Rabi Hamiltonian, including a qubit tunneling term ∝σ x ,D †Ĥ D = η 2 2σ

Supplementary Note 2. SAMPLE FABRICATION
Sample fabrication was carried out in one single electron beam lithography step. The Josephson junctions are formed by shadow angle evaporation and a Dolan bridge technique with electrode film thicknesses of 30 nm and 50 nm, respectively, resulting in a Al film thickness of 80 nm across the entire chip. The area of the Josephson junctions is designed to be 100 nm × 220 nm, resulting in a critical current I c = 45 nA for a single junction. Al is evaporated at a chamber background pressure of about 3 × 10 −8 mbar. We applied an Al metalization on the backside of the double-side polished intrinsic Si substrate as a ground reference for the microstrip elements on chip. Electric field coupling in the substrate due to the backside metalization accounts for roughly half of the qubit capacitance [4].

A. Microwave setup
The schematic diagram of the measurement and microwave setup is depicted in Supplementary Figure 1. We readout the qubit state by observing a dispersive shift of the readout resonator which is acquired via a 400 ns long readout pulse. The resonator shift is extracted from the microwave reflection signal at a single-ended 50 Ω matched transmission line that capacitively couples to the readout resonator. We frequency-convert the readout pulse by heterodyne single sideband mixing to eliminate parasitic population of the readout device during pulse-off time.
Qubit excitation and Rabi driving are performed by heterodyne mixing with respective IQ frequencies for different drive and excitation frequencies, using one single microwave source and one IQ mixer, see Supplementary Figure 1(a). This allows for the required phase control on the phases ϕ 1 , ϕ 2 of the Rabi drives. In particular, we fix the idling time between initial excitation pulse and the onset of the Rabi drive, such that the acquired phase during that time is constant, and apply the Rabi drives with a constant relative phase ϕ i with respect to the phase used for the excitation pulse. We chose an LO frequency located 20 MHz or 65 MHz above the qubit control frequency, located at about ω/2π + 95 MHz. In the experiment with the second Rabi drive added, we generate the drives initially in separate IQ mixers that share a common LO input. We suppress phase errors by employing identical coaxial cables for the high-frequency lines prior to the combination of the drive pulses. The microwave pulses forXŶ control of the qubit are applied via the same transmission line used for readout.
The qubit transition frequency is adjusted by a dc current applied to the on-chip flux coil. High frequency noise is filtered at the 4 K stage with RCR type π-filters at about 25 kHz and on the base plate via an RC-element enclosed in copper powder [5]. Fast flux pulses for fastẐ pulsing of the qubit are sent through a separate microwave line and combined with the offset current by means of a bias tee located at the base plate.

B. Compensation for finite bias tee time constant
In order to combine dc and ac flux signals that are applied to the flux bias line without breaking the 50 Ω impedance matching, we make use of a bias tee in the experiment. Due to a finite time constant τ of the bias tee, a continuous compensation for decharging effects is required in order to produce the desired pulse sequence in the flux bias line on chip. Supplementary Figure 2(a) shows the schematic circuit diagram present during flux pulse generation. V ex denotes the amplitude of the voltage pulse, R is the line resistance in front of the bias tee and C is the relevant capacitance of the bias tee. From Kirchhoff's law we can write with I the current that is admitted to the flux bias line. The condition of constant current follows from requesting which yields Since I is not a function of time according to the initial claim, Supplementary Equation (22) can be integrated yielding This shows that the required correction of the externally applied voltage V ex is linear in time t with a slope proportional to 1/τ . See the pulse sequence applied (blue) and the resulting pattern (red) in Supplementary Figure 2(b). During pulse-off time, V ex has to be kept constant and no further correction is required. The described compensation cannot be performed for longer than approximately 1 µs within one continuous pulse sequence, since the output stage of our pulse generator eventually saturates. With additional amplification of the generated signal, the current is increasing linearly during compensation time, implying an experimental limitation as higher currents can ultimately heat the cryostat. The utilized pulse sequence therefore scales nicely for perspective longer simulation times ∆t since a compensation during ∆t is not required.
We calibrated the time constant of the bias tee used in the experiment to be τ = 0.7 µs.

C. Summary of the device parameters
Supplementary Table 1 lists the relevant parameters of the simulation device that were found experimentally and used in the master equation simulations, together with drive parameters typically used for the quantum simulation.
Values for the geometric coupling strength g are summarized in Supplementary Table 2. Values from the spectroscopy experiment and the vacuum Rabi experiment are in reasonable agreement. The geometric coupling appears to be effectively enhanced during the simulation experiment, which could be due to a parameter drift during successive cool-downs. In order to achieve a better agreement between experiment and master equation simulation, we use a slightly increased value for the geometric coupling strength in classical simulations. We believe that the effective modification of this parameter is a consequence of the discussed parasitic driving of the bosonic mode in the experiment. We checked that this has no qualitative influence on the results, in particular the non-conservation of the total excitation number remains evident.

A. Calibration of the qubit basis
Prior to the experiments presented in the main text, we calibrated the qubit basis spanned by eigenstates |g , |e . This is of relevance for the data analysis detailed in Supplementary Note 5 B and to demonstrate the base line shifts observed in Fig. 3 in the main text. The calibration is performed by Rabi spectroscopy including the fast z pulsing scheme as adopted in successive quantum simulation experiments. The applied pulse sequence and measured Rabi oscillations are depicted in Supplementary Figure 3. The length of the applied excitation pulse is given on the horizontal axis in Supplementary Figure 3(b). Fundamental qubit states (black lines) can be mapped to dispersive shifts of the readout resonator with the bosonic mode in its groundstate and off-resonant.

B. Qubit basis shift
The qubit state is expected to assume an incoherent steady state on the equator of the Bloch sphere for long simulation times ∆t due to the finite coherence in our circuit. However, we notice a shift in the steady state qubit population dependent on the initially prepared qubit state. This is apparent in Supplementary Figure 3(c), showing measurements for the qubit prepared in |e (blue) and |g (green) and in Fig. 3(f)-(i) in the main text. We attribute this effect to a change in the effective qubit basis, which is not captured by the master equation simulations performed. We additionally conjecture that the basis shifts are caused by an effective tilt of the qubit Bloch sphere as an artifact of the frequency tuning in experiment prior to applying the Rabi drives. We isolate the effect as an initialization issue since the basis shift cancels out in good approximation when averaging the mutually antiparallel simulation sequences in Fig. 3(f)-(i).

Supplementary Note 5. CLASSICAL MASTER EQUATION SIMULATION
Numerical simulations are based on a master equation solver provided by the QuTiP package [6,7] for python. The time evolution of a given initial state or density matrix is calculated by solving the von Neumann equation associated with the given system Hamiltonian in the absence of dissipation. For including losses to an imperfect environment of the system Hamiltonian, the time evolution of the density matrix is calculated via the Lindblad master equation.
Classical simulations in the main text include a finite lifetime of the qubit of about 5 µs, a dephasing time of about 0.5 µs and a bosonic mode inverse lifetime of κ = 3.9 × 10 6 s −1 . The Fock space of the bosonic mode is truncated at a photon number of 25, since higher excitation numbers were found to not play a significant role. The transmon qubit is treated as a three-level system with the experimentally found anharmonicity α/2π = −350 MHz. We use a transversal qubit coupling operator with the coupling matrix elements g ij found by evaluating the Cooper pair number operator in the charge basis [8,9].

A. Verification of the simulation scheme
Since a rotating frame is not an inertial system, the laws of physics may be drastically altered when describing a physical system in a rotating frame. In the framework of analog quantum simulation, this offers a rich toolbox which allows to access intriguing effective parameter regimes that are hard or impossible to access in the laboratory frame. While the qubit dynamics in the rotating frame is well reproduced as demonstrated in Fig. 3 in the main text, it is not a priori clear that the bosonic mode population evolution of the quantum Rabi model is well reflected by the laboratory frame simulation. To verify the simulation scheme for the ideal system we compare the bosonic mode population in the driven laboratory frame, and the expected evolution of the ideal quantum Rabi model via classical master equation simulations. Here, we neglect dissipation and parasitic driving of the bosonic mode. Supplementary Figure  4(a) demonstrates good agreement when comparing the time evolutions of the ideal (solid line) and the constructed (dashed line) Hamiltonians. The violation of excitation number conservation in the quantum Rabi model manifests in excitation numbers of the bosonic mode of larger than one for small ω eff . We find that the maximum photon excitation in the bosonic mode roughly equals one for choosing simulation conditions where ω eff ≈ g ∼ 2π × 5 MHz.
For ω eff /2π = 5 MHz we demonstrate that the photon population in the bosonic mode is independent of the applied drive amplitude η 1 , reflecting the fact that it does not appear in the synthesized Hamiltonian, Supplementary Equation (12). Simulations for various η 1 are depicted in Supplementary Figure 4(b).

B. Qubit population retrieval from measured dispersive shifts of the readout resonator
In measuring the dispersive shift of the readout resonator during the quantum simulation experiment, we observe a bulged and shifted equatorial baseline following the expected population evolution of the bosonic mode as obtained from the master equation simulation. We attribute this to a photon exchange coupling f between the bosonic mode and the readout resonator, potentially mediated by the qubit. By inheriting nonlinearity from the qubit it gives rise to a dispersive shift on the readout resonator dependent on the photon number in the bosonic mode. The complete Hamiltonian including the readout resonator of resonance frequency ω r with creation (annihilation) operatorâ † (â) and with the RWA applied takes the form in the subspace spanned by the bosonic mode and the readout resonator. We isolate this additional dispersive shift ∝ f 2 by adding up measured data for the qubit prepared either in |g or |e , as described in the Methods of the main text. Comparison with measured data suggests f ∼ MHz.

C. Effect of the qubit tunneling term in the effective Hamiltonian
For a realistic parasitic coupling η r ∼ 0.1η 1 of the dominant Rabi drive to the bosonic mode and ω eff /2π = 5 MHz, we obtain an effective qubit tunneling term We demonstrate that the position of the quantum revival corresponds to 2π/ω eff in quantum simulations with initial state prepared in |0 ⊗ |g , where both the qubit and the bosonic mode are in its groundstate. Classical simulations are done for κ ∼ 3.9 × 10 6 s −1 and 1/T 1 = 0.2 × 10 6 s −1 , 1/T 2 = 2.0 × 10 6 s −1 and in the absence of a parasitic drive of the bosonic mode.
We find a better agreement of experimental data with classical simulations for a slightly increased geometric coupling g/2π ∼ 5.5 MHz as compared to the measured value in Fig. 2 in the main text of g/2π = 4.3 MHz. This can be explained by the slight excess population due to the parasitic coupling of the Rabi drives to the bosonic mode and the effect may be additionally enhanced by a population of higher transmon levels. Increased decay and decoherence rates, accounting for a changed environmental spectral density in the rotating frame, did not lead to the correct ratio of revival and idling amplitudes. Simulations using g/2π = 4.3 MHz show results with a ∼ 30% decrease in the maximum population of the bosonic mode with no qualitative consequences. With the measured geometric coupling we approach an USC regime with g eff /ω eff ∼ 0.6 for the experiment depicted in Supplementary Figure 7(d). With the slightly increased value for g eff we reach a relative coupling ratio of 0.7. Measured values of g eff and the effective value are summarized in Supplementary We verify the simulation data for the full quantum Rabi model presented in the main text by comparison to measured data with intentionally departing from the required parameter constraints. While an increase in revival amplitude and an increased amplitude of the fast oscillations for larger ∆t is visible by comparing Supplementary  Figure 8(a), (b) presented in the main text, these signatures vanish for violating the parameter conditions required by the simulation scheme. We chose ϕ 1 ∼ ϕ 2 + π but η 1 = ω 1 − ω 2 in Supplementary Figure 8(c) and ω 2 = ω 1 − η 1 − 2π × 10 MHz while ϕ 1 = ϕ 2 in Supplementary Figure 8(d) in order to demonstrate that the desired experimental features vanish.
An additional limitation of the simulation quality is imposed by an uncertainty of the effective Rabi frequency η 1 , which is extracted from a rather broad peak in the Fourier transformed qubit evolution, see Supplementary Figure 5. This causes the main constraint of the simulation scheme η 1 = ω 1 − ω 2 to be poorly satisfied in particular at small simulation times.

Supplementary Note 8. SPECTROSCOPY OF THE AVOIDED CROSSING BETWEEN QUBIT AND BOSONIC MODE
The coupling between qubit and bosonic mode is pre-characterized in a two tone spectroscopy measurement using a vector network analyzer and a microwave source. The dispersive readout resonator shift is measured with a continuous microwave probe tone and an additional microwave drive tone is applied through the same transmission line to excite the qubit transition. The fit of the observed avoided crossing yields a minimum line separation of 2g/2π = 7.8 MHz.    The qubit is prepared in its groundstate |g and the bosonic mode is initially in the vacuum state |0 . The first revival appears at 2π/ω eff , respectively, and the different plots correspond to a varying ω eff . The blue line shows a classical master equation simulation of the ideal effective Hamiltonian in the rotating frame, while the red data points are the qubit population evolution in the effective qubit frame. The black line shows the classically simulated bosonic mode population in the rotating frame, which increases with decreasing ω eff . Note that the scale on the horizontal axis is not equal for each plot. The depicted qubit signal is extracted from measured data with the protocol described in the Methods in the main text and based on the classical simulations (black) in the absence of a parasitic drive of the bosonic mode. Error bars denote a statistical standard deviation as detailed in the Methods of the main text.  Fig. 4 in the main text with η1/2π = 58 MHz. The measurement satisfying the required parameter constraints (b) is compared to measurements with parameters intentionally violating the phase matching condition, ϕ1 = ϕ2 (c) and the condition η1 = ω1 − ω2 (d). When the required constraints are not satisfied, the expected signatures are suppressed. The dispersive shift of the readout resonator induced by the bosonic mode is subtracted based on classically simulated data. Error bars throughout the figure denote a statistical standard deviation as detailed in the Methods of the main text. Supplementary Figure 9. Avoided crossing between qubit and bosonic mode in spectroscopy The qubit transition frequency is tuned by a dc current applied to the flux coil. The dispersive shift of the readout resonator is proportional to the excitation number of the qubit and is depicted in colors.