Quantum simulation of the quantum Rabi model in a trapped ion

The quantum Rabi model, involving a two-level system and a bosonic field mode, is arguably the simplest and most fundamental model describing quantum light-matter interactions. Historically, due to the restricted parameter regimes of natural light-matter processes, the richness of this model has been elusive in the lab. Here, we experimentally realize a quantum simulation of the quantum Rabi model in a single trapped ion, where the coupling strength between the simulated light mode and atom can be tuned at will. The versatility of the demonstrated quantum simulator enables us to experimentally explore the quantum Rabi model in detail, including a wide range of otherwise unaccessible phenomena, as those happening in the ultrastrong and deep strong coupling regimes. In this sense, we are able to adiabatically generate the ground state of the quantum Rabi model in the deep strong coupling regime, where we are able to detect the nontrivial entanglement between the bosonic field mode and the two-level system. Moreover, we observe the breakdown of the rotating-wave approximation when the coupling strength is increased, and the generation of phonon wave packets that bounce back and forth when the coupling reaches the deep strong coupling regime. Finally, we also measure the energy spectrum of the quantum Rabi model in the ultrastrong coupling regime.


INTRODUCTION
The interaction between light and matter is one of the most fundamental and ubiquitous physical processes. The semi-classical Rabi model was proposed in 1936 to analyze the effect of a varying, weak magnetic field on an oriented atom possessing nuclear spin [1]. It describes the dipolar interaction between a classical monochromatic field and a two-level system, successfully explaining the challenging experimental data in Ref. [2]. When the field is promoted to a quantum description, resulting in the simplest fully-quantum model of light-matter interaction, it is called the quantum Rabi model (QRM). Typically, the coupling strength in a light-matter system is much lower than the field frequency. In this scenario, the QRM can be simplified to the Jaynes-Cummings model (JCM) under the rotating-wave approximation (RWA) [3]. The JCM is an analytically solvable model that has been studied in cavity quantum electrodynamics(CQED) [4][5][6], atomic physics [7], quantum dots [8], circuit quantum electrodynamics (cQED) [9,10] and trapped ions [11][12][13][14], among other quantum platforms. Recent experimental achievements have shown the accessibility to the ultrastrong coupling (USC) regime [15,16] or even to the deep strong coupling (DSC) regime [17,18], where the coupling strength is comparable to or larger than the mode frequency. In these strong-coupling regimes, the RWA breaks down, rendering the JCM as a restricted description of the system, and requiring the use of the full QRM to correctly describe the emerging physical phenomena [19]. It is noteworthy to mention that, in such regimes, exotic dynamical properties of light-matter in-teraction [20] and potential applications to quantum information technologies [21] have been predicted and proposed.
The Hamiltonian associated with the QRM can be expressed as ( = 1) whereâ † (â) is the creation(annihilation) operator related to a bosonic mode with frequency ω m ,σ z,x are the Pauli operators acting on a two-level system with energy splitting ω 0 , and g is the coupling strength. Three major coupling regimes are defined depending on the ratio between the coupling strength g and the field mode frequency ω m , namely, the Jaynes-Cummings regime, with g ωm 0.1 and where the RWA is valid, the USC regime with 0.1 g ωm and where the RWA is not a valid approximation, and the DSC regime with 1 g ωm . Regardless of the apparent simplicity of its succint Hamiltonian form, an analytical solution of the QRM has only recently been found [22]. Nowadays, experimental efforts have been made to reach the DSC regime [17,18], which allows for the study of exotic physics such as the bouncing back and forth of the photon number wave packet along the parity chains, and the entangled nature of the ground state of the system [20]. The phenomenon of photon number wave packets bouncing back and forth has been observed in a classical simulator of a photonic waveguide system [23] and in an analogue and a digital quantum simulation with cQED systems [16,24]. However, the study of the ground state in DSC regime is still an open challenge.
In this work, we report the analog quantum simulation of the quantum Rabi model with a single trapped ion for all relevant coupling regimes. Among other results, we generate and observe the ground state of the QRM in the DSC regime in a trapped-ion quantum simulator for the first time. We demonstrate the full controllability and tunability of the QRM in a single trapped-ion system as proposed in Ref. [25], which enable us to generate the exotic ground state in the DSC regime by the adiabatic transfer from the simple ground state of the JCM. Moreover, we apply the capability of the ground state preparation to experimentally measure the energy spectrum of the QRM Hamiltonian (1).

TRAPPED-ION SYSTEM
In our experiment, a radio-frequency Paul trap is used to spatially confine an 171 Yb + ion that is then cooled down to its motional ground state by standard sideband cooling techniques [12] after Doppler cooling. In the lowenergy regime, the motion of the ion can be well approximated to that of a harmonic oscillator, and two energy levels of the hyperfine manifold of its electronic ground state can be used as a qubit. In particular, we encode the two-level system in the levels |F = 1, m F = 0 ≡ |↑ and |F = 0, m F = 0 ≡ |↓ of the S 1/2 hyperfine manifold, which have a transition frequency ω HF = (2π) 12.642812 GHz. We employ a radial vibrational mode of frequency ω X = (2π) 2.498 MHz as the bosonic degree of freedom of our simulator. The uncoupled Hamiltonian describing such a system is given byĤ 0 = ωHF 2σ z + ω Xâ †â . When a pair of counterpropagating Raman laser beams is driven onto the ion, the general ion-laser interaction is described byĤ Here, Ω l is the Rabi coupling strength proportional to the product of both laser field amplitudes, ∆k l is the net wave vector component of the Raman laser beams on the direction of the motion of the ion,x = x 0 (â+â † ) is the position operator of the ion, with ground-state wave-packet width x 0 = /2M ω X , where M is the mass of the 171 Yb + ion, and ω l and φ l are the differences of frequencies and phases of the Raman laser beams, respectively [14].
When both red and blue sideband interactions are simultaneously applied with equal strength, such that Ω = Ω r = Ω b , one can write the total Hamiltonian in the interaction picture aŝ Indeed, Eq. (4) corresponds to the Hamiltonian of the QRM in Eq. (1) considered in the interaction picture with respect to the uncoupled . Therefore, if we undo the interaction picture transformation, we havê where the parameters of the simulated QRM can be associated with the experimental ones as ω 0 = δ b +δr 2 , ω m = δ b −δr 2 and g = ηΩ 2 . Thus, such an experimental setup serves as a quantum simulator of the QRM, where one can simulate a wide range of coupling regimes by suitably tuning the laser intensities and detunings to match the desired ratio g ωm . It is important to point out that the observables of interest {â †â ,σ z , |n n|}, commute with all the adopted interaction-picture transformations, which are always with respect to a Hamiltonian of the form αâ †â + βσ z , such that their expectation values will remain unaltered in the laboratory reference frame [25].

COUPLING REGIMES AND BREAKDOWN OF THE RWA
For the experiment, we fix the coupling strength to g = (2π)12.5 kHz, and the detuning of the red sideband to δ r = 0, leaving δ b as a tunable parameter. In this manner we will be simulating a resonant QRM where the ratio g ωm will be determined by the selected detuning δ b . We experimentally explore three paradigmatic coupling regimes, namely the Jaynes-Cummings, the USC and the DSC regimes, accordingly selecting the values of the detuning for the blue sideband as δ b = 2ω m = (2π){625, 83.4, 41.6} kHz, which correspond to the ratios g/ω m = {0.04, 0.6, 1.2}, respectively. The experiment is carried out as follows. First, we perform standard Doppler and sideband cooling, which prepares the system in the state |↓, n = 0 [27], and then, we transfer the system to the initial state |↑, n = 0 by applying a carrier π pulse. After that, we turn on the red-sideband and blue-sideband transitions, with suitably chosen intensities and detunings, to implement the QRM Hamiltonian in the desired regime. We observe the dynamics of the QRM by measuring the average excitations of the spin σ +σ− and the phononic degrees of freedom â †â at specific evolution times t (see Methods).
In Fig. 1(a) and (d), the measurements for the simulation of the Jaynes-Cummings regime are plotted. Rabi oscillations, with a complete collapse and posterior revival of the excitation probability of the two-level system are clearly observed. In the same manner, the average number of phonons in the bosonic mode oscillates between 0 and 1, consistent with the notion that the wavefunction of the system should live in the space spanned by the corresponding JCM doublet {|0, ↑ , |1, ↓ }, as expected for such a regime. Figures 1(b) and (e) show the evolution of the same initial state in the USC regime for the coupling ratio g/ω m = 0.6. In this case, collapses and revivals of the excitation probability are not complete and the average number of phonons exceeds 1, indicating that the dynamics does not anymore happen exclusively in the JCM doublet. This departure from the JCM physics is associated with the breakdown of the RWA due to the large coupling ratio. In the DSC regime, plotted in Figs. 1(c) and (f) for the coupling ratio g/ω m = 1.2, the effects of the RWA breakdown are even clearer, where not even oscillations can be identified and where the average number of phonons grows above 6 for the plotted example. We also show in Fig. 1(g) the evolution of the total excitation number N = |↑ ↑| + n , which is a conserved quantity when the RWA is valid, but has a dynamical behavior as soon as the RWA breaks down.

DSC REGIME AND PHONON WAVEPACKETS
We focus now on the DSC regime and explore two scenarios, namely the case of the degenerate QRM with ω 0 = 0 and the non-degenerate case with ω 0 = 0. For this experiment, we fixed the coupling strength to g = (2π)12.5 kHz, and vary δ r and δ b , while keeping always a ratio g/ω m = 1.25. For the degenerate case, we use detunings δ r = δ b = (2π)10 kHz, while for the nondegenerate case we use δ r = 0 and δ b = (2π)20 kHz, which corresponds to ω 0 = 0.8g. For the initial state, we choose the ground state of the JC model |↓, 0 , which should have no dynamical properties when the coupling strength g/ω m is small enough. In Fig. 2, we show that the situation is different when considering the DSC regime. First of all, in panels (a) and (b), we show the evolution of the population of Fock state |0 , after tracing out the spin degree of freedom. For the spin degenerate case, we can clearly observe that this population collapses to zero and that it is stabilised at zero except for every one period of the mode frequency, that is to say, except for times k 2π ωm , k being an integer, where a full revival of the population is detected. On the other hand,  the non-degenerate case shows a degradation of these revivals for long times, as it was predicted in Ref. [20]. Additionally, we sample several points during one period T = 2π ωm = 100µs of the evolution of the degenerate case and measure its phonon distribution, as shown in panels (c-h). We obtain the phonon number distributions by fitting the spin-excitation evolution under the bluesideband transition with the function shown in the Methods section. At time zero, the population is concentrated on Fock state |0 , and as time elapses higher Fock states are populated. The evolution resembles a wavepacket that travels along a chain of Fock states up to a maximum determined by ∼ 4(g/ω m ) 2 and then comes back to the initial states at one period of the mode frequency. This phenomenon was theoretically predicted in Ref. [20] as characteristic of the DSC regime, and it is referred to as the bouncing back and forth of phonon-number wave packets. We note that the ω m = 0 case has already been studied in the literature for its application as a geometric phase gate, where the effects over the spin degrees of freedom are the main interest [28,29].

ADIABATIC GROUND-STATE PREPARATION
As mentioned in the previous section, the ground state of the QRM in the Jaynes-Cummings regime (g ω m ) is given by the state |↓, 0 , while the ground state of the QRM in the DSC regime is a nontrivial state where spin and field are entangled, and which to the best of our knowledge has never been implemented in a physical quantum platform.
In our experiment, we generate the ground state of the QRM in the DSC regime by starting in the ground state of the Jaynes-Cummings regime, the state |↓, 0 , and adiabatically increasing the coupling ratio g ωm towards the DSC regime. To achieve this, one could choose to either increase g = ηΩ/2, which can be done by rising the laser intensity, or decrease ω m = (δ r − δ b )/2. Because it is experimentally more feasible to manipulate the detuning of the Raman lasers than their power, we choose the latter in our experiment. For that, we fix coupling strength g to be (2π)12.5 kHz and δ r = 0, leaving δ b as the only tunable parameter, which we manipulate with an exponential time dependence of the form δ b (t) = (δ Max − δ Tar )e − t τ + δ Tar . Here, we set δ Max = (2π)0.2 MHz, while δ Tar is determined by the ratio g ωm we want to reach, and τ = tTar 10 , t Tar = 300µs being the total duration of the adiabatic process. The adiabaticity of our scheme is guaranteed by the numerical computation of the fidelity between the instantaneous ground state of the Hamiltonian and the adiabatically evolved state, which is shown in Fig. 3a.
In panel (b) of Fig. 3, we show the spin evolution dur- In panel (a), we show the adiabatic scheme for the preparation of the ground state of the QRM at the DSC ratios g ωm = 1.2, 1.5, and 2.0, as starting from the initial JC ratio g ωm = 0.125. The fidelities between the instantaneous ground state and the numerically evolved state quantify the adiabaticity of our process. In panel (b), we show the evolution of the excitation probability and the purity of spin state during the adiabatic ground state preparation (0 ≤ t ≤ tTar) and during the reverse process (tTar ≤ t ≤ tRev). The plot corresponds to the preparation of the ground state at the ratio g ωm = 1.2. The red dashed line and the green solid line are obtained by direct diagonalization of the QRM Hamiltonian and numerical simulation of the adiabatic process, respectively, including heating and dephasing of the motional mode, which are expected experimental imperfections. The blue circles with error bars correspond to the experimental results. The purple line represents the numerically computed purity of the spin, defined as Tr(ρ 2 spin ), where ρspin is the reduced density matrix of the spin after tracing out the motional degree of freedom. The orange squares are the corresponding experimental results, computed from the spin-tomography. Panels (d) to (f ) show the phonon-number distributions correlated with |↓ (lower panel) and |↑ (upper panel), which are obtained by fitting the standard blue-sideband signals after the spin-projective measurement for g ωm = 1.2, 1.5, and 2.0, respectively. Finally, panel (f ) shows the purity of the spin and estimated lower bounds for the purity of the whole system, for each tested coupling ratio. The purity of the spin is obtained from the measured reduced density matrix, which was done by spin-tomography. The lower bounds of the purities of the whole system are estimated by measuring the probability of being in |↓, 0 after the reverse process at tRev.
ing the adiabatic process for the time interval (0-t Tar ). The plot corresponds to the case g/ω m = 1.2, with the cases for other ratios showing similar behavior. At time t Tar , the system is expected to be in the ground state of the QRM for the selected coupling regime. In panels (d) to (f), we plot the outcome of the phonon distributions correlated with the spin performed at this time for coupling ratios g/ω m = 1.2, 1.5 and 2.0, respectively (see Methods). To verify the quantum coherence maintained within the preparation of the ground state, we reverse the adiabatic process in an attempt to recover the initial ground state |↓, 0 . In panel (b), we can observe how the spin returns to state |↓ . As a complementary proof, we plot the purity of the spin state. To this aim, we trace out the phononic degrees of freedom and measure the density matrix ρ associated with the spin degree of freedom [30], from which we calculate the purity, defined as Tr(ρ 2 ), during the whole process. The degradation of the purity during the preparation of the ground state of the QRM in the DSC regime confirms that the adiabatically prepared ground state is indeed an entangled state, and the subsequent revival of the purity when the adiabatic process is inverted proves that we are able to recover the initial state and therefore that the whole process preserves quantum coherence.
On the other hand, the expectation values of the parity operator Π = σ z e −iπâ †â , for the states in panels (d) to (f) of Fig. 3, are respectively 0.74(0.08), 0.70(0.08) and 0.52(0.13), showing that the ground states in the DSC regime belong mostly to the same parity chain due to the Z 2 symmetry of the QRM [20]. As the coupling ratio increases, the deviation from the ideal parity value +1 becomes more prominent due to imperfections of the adiabatic process and the motional heating arising from the occupation of larger Fock states.
By measuring the probability of recovering the initial state |↓, 0 after the ground state preparation and inverse process at time t Rev , we estimate a lower bound of the purity of the prepared ground state. The revival probabilities are 0.89(0.024), 0.87(0.027) and 0.79(0.03) for the three ratios g/ω m = 1.2, 1.5 and 2.0, which give the lower bounds 0.79(0.042), 0.75(0.047) and 0.62(0.047), respectively. As shown in Fig. 3c, the reduced-spin purities, taking values 0.545(0.006), 0.514(0.003) and 0.505(0.002), are significantly smaller than the lower bounds of the total system, which prove the existence of entanglement within the prepared ground state at t Tar (See Methods).

SPECTRUM
The ground-state preparation can be extended to study the low-lying energy spectrum of the QRM by coherent spectroscopy [31]. In particular, we have measured the energy spectrum in the region g ωm ∈ [0, 1]. A Z 2 parity exists in the QRM model, which divides the Hilbert space in two, namely a subspace of parity +1 and other of parity −1 [20]. Here, we focus on the energy splittings between the ground state and the first three excited states of opposite parity to the ground state [20]. For that, we have used a relatively weak modulated field as a probe on top of the simulation of the QRM, with the system initially in the ground state of the corresponding regime. We sweep the frequency of the probe pulse until we detect a transition, and we associate the frequency of the probe to the energy difference of the transition. To generate transitions between states of opposite parity, we use the probe pulse of the form where g p ( g) is the strength of the modulation field, and ν p is swept to find the resonant frequencies. In the region g ωm = 0.1 to 0.3, g p /g is 0.02, while the pulse duration is 350 µs. For the ratios g ωm = 0.4 to 1.0, the ratio g p /g is 0.01, with a pulse duration of 450 µs. Population transfer is clearly seen when ν p is resonant with the energy splittings as shown in Fig. 4.

CONCLUSION
We have implemented the quantum simulation of all relevant coupling regimes of the QRM in a single trapped ion, obtaining direct evidence of the breakdown of the RWA. Historically, trapped ions have been always linked to the JCM physics, which has been enhanced here toward the more general QRM. In the DSC regime, we observe the phonon number wave-packets bounce back and forth as well as collapses and revivals of the initial state, confirming previous theoretical predictions. The adiabatic preparation of the ground state of the QRM in the DSC regime was produced for the first time in a quantum platform, and its reconstruction has enabled us to demonstrate the entanglement present in its ground state. As a direct application of this adiabatic method, we have been able to measure the energy splittings between states of different parity and recreate the energy spectrum of the QRM in the USC regime. In conclusion, our work presents a detailed experimental exploration of the QRM in a wide range of physical regimes. Our experimental methods can be directly extended to the study of the phase transition in the QRM [32][33][34] or to the simulation of the Dicke model [35][36][37] by considering the presence of more ions.

Calibration of the detuning of the blue and red sideband transitions
For this experiment, since we fix the coupling strength g, the key ingredient for the simulation is to precisely set the detuning of the lasers. For the ratio g ωm = 0.04 case, since we set δ r = 0, and the δ b is much larger than the corresponding coupling strength, we obtain the resonance frequency of red-sideband transition with the detuning of the blue-sideband transition fixed at δ b = (2π)625 kHz. For the USC/DSC regime, the coupling strength is comparable to the effective mode frequency, such that we need to carefully deal with the ac-Stark shift mainly introduced by the off-resonant excitation on the carrier transition. We measure the ac-Stark shift with a Ramsey experiment and calibrate the shift in the bichromatic pulse within 1 kHz accuracy. We further improve the frequency precision within 0.15 kHz range by giving the same detuning δ = (2π)10 kHz with different signs for both beams similarly to the scheme in Refs. [38,39].

Phonon number state population distribution
In Fig. 3, we obtain the phonon number distribution. This is performed by driving the resonant blue-sideband transition |↓, n ← |↑, n + 1 after the spin projective measurement and fitting the obtained spin population evolution with the formula [11,26,40] where p(n) is the phonon number state population, γ is the empirical decay constant, and t is the pulse duration of the blue sideband.
Verification of entanglement for the ground-state of the quantum Rabi model In the main text, we use Tr ρ 2 spin − P 2 Rev < 0 to verify the existence of the entanglement between the spin and the phonon degrees of freedom. This can be understood as follows. First, we introduce the purity-based entanglement witness [41,42] W for the target stateρ Tar at time It can be proved that W [ρ] ≥ 0 for arbitrary separable state. Thus W [ρ] < 0 serves as the sufficient condition for the inseparability ofρ. However, the purity of the whole system Tr ρ 2 Tar requires full information ofρ tot , which is quite demanding in our current experimental setup. Instead, we perform the disentangling operation, which is the time reversal of the adiabatic ground-state preparation in this particular case, and measure a single component P Rev ≡ T r[| ↓, 0 ↓, 0|, ρ Rev ] of the spectral decomposition of the final state at time t Rev , which corresponds to the probability of the initial state . It's straightforward that P 2 Rev ≤ Tr [ρ Rev ]. With reasonable assumption that the evolution performed in this experiment never increases purity, we can push the inequality to an earlier time. Then we have the following inequality Tr ρ 2 Tar ≥ Tr ρ 2 Rev ≥ P 2 Rev .
In other word, P Rev serves a lower bound for Tr ρ 2 Tar . Putting Eqs. (8) and (9)