A robust scheme for the implementation of the quantum Rabi model in trapped ions

We show that the technique known as concatenated continuous dynamical decoupling (CCD) can be applied to a trapped-ion setup for a robust implementation of the quantum Rabi model in a variety of parameter regimes. These include the case where the Dirac equation emerges, and the limit in which a quantum phase transition takes place. We discuss the applicability of the CCD scheme in terms of the fidelity between different initial states evolving under an ideal quantum Rabi model and their corresponding trapped-ion realization, and demonstrate the effectiveness of noise suppression of our method.


Introduction
Quantum coherence is an essential prerequisite to observe and exploit the intriguing phenomena in the quantum realm [1]. Indeed, technologies relying on those quantum properties are expected to surpass their classical counterparts in efficiency and performance. This new generation of quantum technologies encompasses a large diversity of possible applications which inlcude quantum simulation [2], quantum metrology [3], quantum communication [4] and quantum sensing [5], all of them requiring the preservation of quantum coherence for their correct functioning. In this respect, the loss of quantum coherence, or simply decoherence, is a crucial limitation as it occurs due to the unavoidable interaction of the quantum system with an uncontrolled environment as well as to the presence of experimental imperfections. Hence, the longtime maintenance of the quantum coherence of an evolving system is highly desired although its realization constitutes a formidable task.
During the past decades considerable efforts have been invested in the development of theoretical schemes to circumvent, as much as possible, the effect of the noise in the system with the goal of prolonging coherence times. Among them we find techniques such as decoherence-free subspaces [6], quantum error correction [7], or dynamical decoupling [8]. These are methods designed to handle specific noise scenarios, and present different benefits concerning noise supression. In particular, dynamical decoupling constitutes a promising tool to handle non-Markovian noise, and it is the central object of study in this article. In its continuous wave configuration, the effect of dynamical decoupling corresponds to the creation of a dressed basis with an energy gap such that, under certain circumstances that will be later developed, the effect of noise is suppressed. In addition, this technique allows for a concatenated configuration known as concatenated continuous decoupling (CCD) [9] that consists in applying concurrently different driving fields to eliminate further sources of noise, including those from imperfect driving fields themselves. Standard dynamical decoupling has been theoretically proposed in its continuous [10,11,12,13] and pulsed [8,14,15,16] configurations. Furthermore, these techniques have already been used in both radio frequency and Penning traps in [17,18] (continuous case) and in [19,20,21,22] (pulsed case) as a method to suppress noises on the registers and to drive robust single-and two-qubit gates. Furthermore, dynamical decoupling has been used to explore different models involving spin-spin interactions [23]. On the other hand, the CCD scheme has experimentally demonstrated its feasibility to preserve the coherence of an isolated nitrogen-vacancy center in diamond [9]. However, the convenience and possible benefits of the CCD method in an ion trap platform for quantum simulation purposes has not been proven yet.
In the present article we show how to apply the CCD scheme in a trapped-ion setting for a robust implementation of the paradigmatic quantum Rabi model that describes the interaction between a two-level system and one bosonic field mode. Despite of its apparent simplicity, this model exhibits a rich variety of physics, ranging from the relativistic Dirac equation [24,25,26,27] to critical phenomena as it can undergo a second-order quantum phase transition [28,29]. We demonstrate that, within the CCD scheme, high fidelities can be achieved and maintained during long evolution times in an ion trap setup in the presence of different noise sources and realistic conditions. While an experimental verification of such scheme in an ion trap is still required, the present theoretical results are promising and open the door to the study of robust and noise-resilient trapped-ion quantum simulations.
We exemplify and support by means of detailed numerics the applicability of the CCD scheme realizing the quantum Rabi model in three different parameter regimes. First, the case where the energy splitting of the two-level system matches the motional frequency and the rotating-wave approximation can be applied. In this situation the Jaynes-Cummings model [30] emerges and we can observe Rabi oscillations. Second, the realization of the Dirac equation [24,25,26,27] whose main hallmark is the Zitterbewegung, and finally, the extreme parameter regime [31] required to witness critical dynamics as a consequence of the emergence of a second-order quantum phase transition in the limit of strong coupling [28,29]. Additionally, we discuss possible drawbacks in the CCD scheme and identify particular situations where the method does not lead to an improved performance.
The present article is organized as follows. In Sec. 2 we introduce the Orstein-Uhlenbeck stochastic process [32,33], which we will use to model fluctuations in the trapped-ion setting as well as of the externally applied control fields. In Sec. 3 the CCD scheme is presented and explained. Furthermore, we show how CCD adapts to trapped-ion Hamiltonians giving rise to a noise protected quantum Rabi model in Sec. 4, while specific examples and their numerical simulations are shown in Sec. 5. Finally, we summarize the main conclusions in Sec. 6.

Stochastic fluctuations: Orstein-Uhlenbeck process
A quantum system looses its quantum coherence due to an uncontrolled interaction with the environment. Such interaction introduces a stochastic noise or fluctuation in the system that we will model as an Orstein-Uhlenbeck (OU) stochastic process [32,33,34]. This effective description successfully reproduces the exponential decay of the quantum coherence due to dephasing noise as measured by Ramsey interferometry [35], as well as the behavior of a quantum system under fluctuations on the intensity of the applied radiation [9]. Moreover, as we will see later on, it also allows to vary the width of the spectral density, which quantifies the amount of power per unit of frequency. In this manner the OU process can describe different noise scenarios, and thus, it has been extensively used in the literature [10,11,36,37].
An OU process is characterized by two parameters, namely, τ and c, relaxation or correlation time and diffusion constant, respectively. While the former fixes the time in which the noise is correlated, the latter is proportional to the noise amplitude. A stochastic variable X(t) that obeys an OU process has an exact update formula [34], for an arbitrary value of ∆t. The term N(t) stands for a temporally uncorrelated normally distributed random variable, i.e., N(t) = 0 and N(t)N(t ′ ) = δ(t − t ′ ), where the overline denotes the stochastic average. The OU process is Gaussian, and hence, fully determined by its first and second moments, where σ 2 [X] denotes the variance of X, and thus, σ[X] its standard deviation. The power spectrum or spectral density, S X (f ), characterizes the nature of the noise, since it measures the amount of power per unit of frequency of X(t) at a frequency f . The stochastic variable X(t) can be written in Fourier series as X(t) = n P n e 2πifnt for t ∈ [0, T ] where P n are the corresponding Fourier coefficients at frequency f n . Then, the spectral density can be defined in the T → ∞ limit, as shown in [33], as The spectral density will be of importance in the next section, Sec. 3, for the understanding of the noise decoupling efficiency of the CCD method. Indeed, for the particular case of an OU process, S X (f ) can be analytically calculated giving rise to [33] S X (f ) = cτ 2 1 + 4π 2 τ 2 f 2 .
Therefore, the relaxation time τ sets a boundary in the frequency domain between white noise, i.e. S X (f ) ∝ f 0 , and Brownian or red noise, i.e. S X (f ) ∝ f −2 . This crossover frequency f cr can be estimated as S X (f cr )/S X (0) = 1/2, that is, f cr = 1/(2πτ ). In Fig. 1 we show a typical trajectory of an OU process for a fluctuating variable δ m (t) and its Fourier transform. Note that S X (f n ) ∝ |P n | 2 .
Here we are interested in magnetic-field fluctuations or simply dephasing noise, which can be written as H = δ m (t)/2 σ z where δ m (t) follows Eq. (1). The coherence time of the system depends then on the properties of δ m (t). For example, consider an initial state |↑ x at t = 0, i.e. σ x |↑ x = + |↑ x , evolving under H = δ m (t)/2 σ z , then it is easy to prove that where ϕ(t) = t 0 ds δ m (s) is the time integral of the stochastic variable δ m (t) and ϕ 2 (t) its autocorrelation function that can be written as [34] The coherence time T 2 is defined as the time instant at which σ x (T 2 ) = e −1 . Hence, from Eq. (6) and (5) it follows that that is, for a given τ and a coherence time T 2 , the diffusion constant can be determined. Nevertheless, depending on whether the noise is fast, i.e. with short memory, meaning τ ≪ T 2 , or slow, i.e. with long memory, which corresponds to τ T 2 , the coherence decays differently. Indeed, exponential decay is achieved when τ ≪ T 2 which is the typical scenario in ion traps [35]. In this case Eq. (7) acquires a simpler form: T 2 ≈ 2/cτ 2 . In contrast, for slow noise a Gaussian decay is observed. In Fig. 2 we plot σ x (t) as a function of the evolution time t for an initial state |↑ x evolved under fast and slow noise, considering T 2 = 3 ms, τ = 50 µs and τ = 5 ms, and c obtained according to Eq. (7). We can observe how the numerical stochastic average σ x (t) agrees with the exact expression in Eq. (5).

Concatenated Continuous Decoupling (CCD)
In this section we explain the technique known as dynamical decoupling in a concatenated scheme (CCD) [9] that corresponds to the addition of several continuous decoupling fields. Note that the use of continuous fields, not pulsed, will be maintained throughout the article. Consider a situation where the Hamiltonian is H = ω 0 (t)/2 σ z where ω 0 (t) = ω 0 + δ m (t) with δ m (t) the stochastic fluctuation of ω 0 , which strongly affects the quantum coherence of the system. Then, in order to eliminate its effects a continuous driving field with Rabi frequency Ω is introduced. This situation is described by the Hamiltonian 1/e T 2 =3ms Figure 2. Decoherence due to a fluctuating σ z term (dephasing noise), which follows a OU process. Two different noises have been considered: fast or short memory noise (solid-blue) τ = 50 µs and slow noise (solid-red) τ = 5 ms. The coherence time T 2 = 3 ms is fulfilled in both cases when the diffusion constant c is calculated from Eq. (7). An initial state |↑ x is prepared, and then the expectation value σ x (t) is calculated as an average over 1000 stochastic trajectories, and presented together with their corresponding sample variance (dashed lines). The solid black lines show the exact σ x (t) according to Eq. (5). As expected, for slow noise a Gaussian behavior is observed e −(t/T 2 ) 2 , while the case with fast noise decays exponentially e −t/T 2 .
In an interaction picture w.r.t. ω 0 /2σ z we have thus, selecting ω = ω 0 and invoking the rotating-wave approximation (RWA), the previous Hamiltonian (in the case Ω ≪ ω 0 ) reads The first term on the r.h.s of the above equation produces no transition in the basis {|↑ x , |↓ x } as long as the fluctuating term, δ m (t), has vanishing Fourier coefficients, |P n | ≪ 1, in the vicinity of frequencies f n ≈ Ω. In other words, to protect the system against the noise, the Rabi frequency Ω must lie in the region in which the noise spectrum is negligible. In this manner, transitions in the dressed basis {|↑ x , |↓ x } as a consequence of the stochastic term δ m (t)/2 σ z have an energy penalty and can be neglected. We will denote this first step as the first layer of protection, since only one additional driving has been introduced. From a more rigorous point of view, that noise elimination is achieved after the application of a RWA on each of the noise components as a consequence of the presence of the term Ω 2 σ x . In addition, and because the RWA presents a slightly different behavior depending on the initial state of the system, the proposed method inherits its dependence. Note however the existence of certain states for which its evolution under noise and the Hamiltonian just gives rise to a global phase.
For such dark states, introducing a first layer deteriorates the coherent evolution since, in the rotated basis, noisy terms are able to produce transitions. In 5.1 we will comment more about this scenario and show an example. Now one should also consider that the Rabi frequency Ω is not completely stable and represents another source of fluctuations, that is, Ω ≡ Ω[1+δ Ω (t)] with δ Ω (t) another stochastic fluctuation with a small amplitude. However, the CCD scheme offers the possibility to further protect the system against δ Ω (t) with a second layer by introducing one additional driving to cancel δ Ω (t) [9].
In Fig. 3 we sketch the main idea behind the effectiveness of dynamical decoupling to cancel interfering stochastic processes. In Fig. 3 (b) the evolution of the coherences as a function of the evolution time is plotted for three different drivings. The success depends on the properties of the noise (a): when the Rabi frequency of the driving does not exceed the crossover frequency of the noise (Ω 1 < f cr ) no protection is achieved. On the contrary, as the Rabi frequency gets larger, Ω 2,3 f cr , the quantum coherence is preserved during longer times since transitions due to the original noise occur with a smaller probability in the new dressed basis. This shows the crucial interplay between noise properties and driving frequencies in a dynamical decoupling scheme. Then, one can apply the same criteria to cancel further fluctuations of additional drivings fields in the CCD scheme. Note that the same techniques can be applied to other noise models that present a similar behavior, i.e. models exhibiting a spectral density that vanishes for asymptotically large frequencies.

Trapped-ion Hamiltonian and CCD
Consider a trapped-ion with internal electronic structure described by ω I /2 σ z and νa † a representing the motional mode energy with ν the trap frequency. The interaction created by a laser irradiation is captured in the term Ω j /2σ x e i(kx−ω j t−φ j ) + H.c. . Hence, under the influence of applied radiation the trapped-ion Hamiltonian reads [38] where k j is the wave vector of each laser field, ω j its frequency, φ j an initial phase, Ω j the Rabi frequency of the jth laser, andx the ion position operator. Before starting with further developments, let us introduce some typical values of the parameters in the previous equation according to the state-of-the-art in experiments with 40 Ca + [25,27]. Here, the axial trap frequency is ν = 2π × 1.36 MHz, ω I is on the optical regime at 729 nm, i.e. ω I = 2π × 4 · 10 14 Hz, and the Rabi frequency is typically on the order of several kHz [25,27]. Additionally, we should consider the coherence time of the internal levels of the ions as the main limiting factor that affects to the quality of the experiments with 40 Ca + [25,27]. As we already commented this is caused by magnetic-field fluctuations which give rise to a coherence time T 2 ≈ 3 ms, see [27]. We will consider this value throughout the present article. Note however that, by using a cryogenic setup [39], a longer coherence time of T 2 ≈ 18 ms has already been achieved. Additionally, laser-intensity fluctuations are present in any realistic ion trap experiment, Figure 3. Schematic representation of the CCD scheme. In (a) the normalized power spectrum of the noise is plotted. Depending on the Rabi frequency Ω i of the additional pulse, different evolution of the coherences is observed (b). As sketched in (c), the original basis suffers dephasing. Then, if the introduced Ω i is small compared to the characteristic frequency of the noise, there is essentially no protection, while Ω i f cr coherence times are enhanced significantly as the noise term δ(t)σ z is not enough to produce transitions in the new dressed basis. Noise parameters are τ = 50 µs and T 2 = 3 ms, while the Rabi frequencies Ω 1 = 2π × 0.5 kHz, Ω 2 = 2π × 5 kHz and Ω 3 = 2π × 50 kHz, and ω 0 ≫ Ω 3 such that a RWA can be safely applied.
while its frequency ω j and phase φ j can be very accurate. Although these magnetic and intensity fluctuations are the main limiting factor for the coherence time of the system, there are still another sources of noise which will be not considered here as they will produce significant effects only on time scales significantly longer than T 2 = 3 ms. In this respect, phonon dephasing has been measured with an incidence of few Hz [40]. This provides a limit of the time scale across which the dynamics can be observed, which is, approximately, two orders of magnitude larger than the one we could consider if the magnetic noise is not eliminated. Concerning the heating rate it can be estimated that, on average, one phonon is gathered in ∼ 100 ms [40], or in ∼ 500 ms for a cryogenic setup [39]. Furthermore, the lifetime of the qubit for the D 5/2 state of 40 Ca + is ∼ 1s [40].
Regarding the trapped-ion Hamiltonian, in the interaction picture w.r.t.
where we have already performed the optical RWA, i.e., we neglect the terms that rotate at frequency ω I + ω j (counter rotating terms). Since ω j will be chosen such that ω j ≈ ω I and because Ω j ≪ ω I + ω j , this approximation can be safely carried out. We denote ∆ j = ω I − ω j , thus, choosing ∆ = 0, ν or −ν one arrives to a carrier, red sideband or blue sideband interaction, respectively, when the system is adjusted to lie within the Lamb-Dicke regime (η j (a + a † ) 2 ≪ 1). Here, the Lamb-Dicke parameter η j is , m the mass of the ion and = 1 throughout the whole article; thus,x = x 0 a + a † . Finally, we would like to remark that all the numerical simulations of trapped-ion Hamiltonians presented in this article have been performed after the optical RWA and without further assumptions.

CCD for a single trapped-ion setup
We discuss now how to employ a CCD scheme in a single trapped-ion setup. In [41] it is demonstrated that, by using two traveling waves to excite the red-and blue sideband transitions, and by setting properly the parameters Ω 1,2 , φ 1,2 and ω 1,2 , the Rabi model can be simulated in a variety of parameter regimes which includes the Dirac equation as a particular case. However, the presence of different noise sources could significantly deteriorate its realization. Therefore, a noise-resilient implementation is desired to enhance coherence control and fidelity. For that reason, in the following we apply a CCD scheme to a single trapped-ion setup. We use the first layer (4.1.1) to tackle the dephasing noise as it is the main limiting factor for the coherence time of the system, while the second layer is introduced to handle laser-intensity fluctuations (4.1.2).

First layer
In order to achieve the Rabi model within the CCD scheme, we apply an extra laser, denoted by the subscript a, with the objective to introduce a term Ω a cos(ω I t)σ x into the dynamics. This is accomplished by setting ω a = ω I (resonant with the frequency splitting of the ion), φ a = 0 and a Rabi frequency Ω a ≪ ω I . Then, the trapped-ion Hamiltonian in a rotating frame w.r.t. H 0 = ω I /2σ z + νa † a and after the optical RWA reads where δ m (t) follows an OU process and is responsible of the dephasing noise, ∆ j = ω I −ω j is the detuning and η j the Lamb-Dicke of the jth laser. Note that the additional laser a has zero detuning, ∆ a = 0 which ensures a carrier interaction (i.e. a σ x proportional term) within the Lamb-Dicke regime, and when other terms, i.e. the ones with a linear dependence in the Lamb-Dicke parameter, can be averaged out because of the condition Ω a η ≪ ν. Hence, only the first term of the following expansion is considered, In this way the additional continuous driving a, provides a dressed spin-basis, {|↑ x , |↓ x }, in which the system is protected against the magnetic-field fluctuation or dephasing noise, δ m (t)/2σ z , as long as Ω a fulfills the criteria given in Sec. 3. Then, the magnetic-field fluctuation can be eliminated and the Hamiltonian (13) is Furthermore, choosing properly the detunings and phases, ∆ j and φ j , a tunable Rabi model can be obtained from the previous effective Hamiltonian. This can be accomplished by setting two lasers j = 1, 2 with ∆ 1 = ν − ξ and ∆ 2 = −ν + ξ (detuned red and blue sideband), for which only the terms at first order in η (η 1,2 = η) of the expansion in Eq. (14) survive, provided by ξ ≪ ν and Ω j ≪ ν; that is, we are applying the vibrational RWA. Finally, the Rabi model is achieved when an interaction term is orthogonal to the free energy term of the two-level system, which in this case is σ x . Therefore, it suffices to set the phases as φ 1 = φ 2 = 0 and the Rabi frequencies Ω 1,2 = Ω, The previous Hamiltonian corresponds to a Rabi model in a rotating frame w.r.t. ξa † a, i.e.
We remark that the previous effective Hamiltonian is only valid under both optical and vibrational RWA, within the Lamb-Dicke regime and when Ω a is such that the noise δ m (t) has vanishing small component at that frequency. Under the same approximations, the Dirac equation can be obtained. The corresponding Hamiltonian of the (1 + 1) Dirac equation [24,26] reads H D = c Dp σ x + m D c 2 σ z , where c D is the speed of light, m D the mass of the 1 2 -spin particle, andp the momentum operator. To realize such a Hamiltonian from Eq. (15), we select ∆ 1 = ν, ∆ 2 = −ν (red and blue sideband), φ 1 = 3π/2, φ 2 = π/2 considering η 1,2 = η and Ω 1,2 = Ω (together with ∆ a = 0 and φ a = 0). Then, Eq. (15) reads wherep = i(a † − a)/2. This is equivalent to the Dirac equation with the following parameters c D = ηΩ and m D = Ω a /(2η 2 Ω 2 ).

Second layer
Once the main source of noise, magnetic field fluctuations, is overcome by means of the first layer, the following step consists in facing laser-intensity fluctuations which can still spoil quantum coherence. The intensity of a jth laser is now modeled as Ω j (t) = Ω j 1 + δ Ω j (t) , where Ω j is the desired Rabi frequency and δ Ω j (t) describes a small stochastic fluctuation. Such fluctuation will be present for all the lasers used in the setup. That is, the laser intensities are not completely stable, but fluctuate around its mean value Ω j . We characterize these fluctuations as an OU process with τ Ω = 1 ms following [42], and an amplitude of 0.1% (p = 0.001) of the laser intensity Ω j . Thus, one can characterize this as σ[δ Ω ] = p, which leads to c Ω = 2p 2 /τ Ω . Note that the laser-amplitude noise is chosen to be slow, compared to δ m (t). This fact can be seen as a technological requirement as otherwise the noise might not be easily handled within the CCD scheme as we will discuss later on. In this way, once δ m (t)/2 σ z is overcome, the main fluctuation in Eq. (15) appears in the free energy term of the two-level system (i.e. as dephasing noise). Note that the rest of the Rabi frequencies, Ω j , are multiplied by a Lamb-Dicke parameter which reduces the influence of the errors introduced into the system by their fluctuating character. Therefore, we can proceed as for the first layer to deal with the term Ω a δ Ωa (t)/2 σ x . To eliminate its contribution an additional continuous driving, denoted by the subscript b, is introduced, but with a time-dependent Rabi frequency Ω b 2 cos(Ω a t). The Hamiltonian describing this situation in a rotating frame w.r.t.
where we have already fixed ∆ b = 0. By simplicity, we only write down explicitly the fluctuation δ m (t) and δ Ωa (t), although all the functions δ Ω j (t) have been taken into account in our numerical simulations, see next Section. As we need an orthogonal carrier with respect to σ x for Ω b , we select φ b = π/2 which leads to Ω b cos(Ω a t)σ y . Now we move to a rotating frame w.r.t. Ω a /2σ x obtaining The spin raising and lowering operators have contributions of σ x and σ y , i.e. σ ± = 1 2 (σ x ± iσ y ), which in a rotating frame with respect to Ω a /2 σ x makes σ y to rotate at frequencies ±Ω a while it does not affect σ x . We then invoke the RWA to average out those rotating terms. Note that this is valid under the assumption Ω b ≪ Ω a . The free energy term of the effective two-level system is given now by σ y , and hence, the new dressed spin-basis is |↑ y , |↓ y . In this basis the fluctuating term Ω a δ Ωa (t)/2 σ x can be depreciated following the same arguments given in Sec. 3, as well as δ m (t). Hence, the Hamiltonian can be approximated by We can summarize the operating regime on the second layer as Ω b ≪ Ω a ≪ ω I . Additionally, Ω a has to be large enough to ensure decoupling with respect to δ m (t), this condition is Ω a 1/(2πτ m ) or, in different words, Ω a has to be larger than the crossover frequency, see Sec. 2. At the same time, and following the same arguments, Ω b needs to handle the fluctuation Ω a δ Ωa (t)/2 σ x , and hence, Ω b 1/(2πτ Ω ) which implies the relation τ Ωa ≫ τ m . Yet, both the intensity of the noise and the RWA (Ω b ≪ Ω a ) play a decisive role to successfully apply a second layer of protection in the CCD scheme.
We note that now we may use only one traveling wave to produce the Rabi-like interaction. Setting ∆ 1 = +ν − ξ, φ 1 = 3π/2 , we arrive to after using the vibrational RWA. The previous equation is equivalent to the Rabi model in a rotating frame w.r.t. ξa † a, As in the case of the first layer, the Dirac equation can be realized in a straightforward manner. Choosing Ω 1 = Ω, η 1 = η, ∆ 1 = ν and φ 1 = π the Eq. (21) reduces to which is equivalent to the Dirac Hamiltonian with c D = ηΩ/2 and m D = 2Ω b /(η 2 Ω 2 ). Note that the effective Hamiltonians given in Eqs. (22) and (24) are valid under a number of approximations, as for the first layer. Additionally, we now require Ω b ≪ Ω a due to a RWA, but at the same time Ω b must be still large enough to decouple with respect to the noisy term Ω a δ Ωa (t)σ x .

Numerical results
Here we present numerical simulations of the previous derived effective Hamiltonians. We compare the usefulness of CCD scheme in contrast to the bare realization, denoted here as zeroth layer (see for example [41] and Appendix A for a derivation), i.e., when no protection against noise is provided. We explore two physical regimes in the realized quantum Rabi model, namely, the paradigmatic resonant case to observe Rabi oscillations, and the limiting case where a quantum phase transition takes place [28,29]. Then, we present the case of the evolution of a Dirac particle. We emphasize that all the numerical simulations involving trapped-ion Hamiltonians have been carried out after the optical RWA without performing further approximations.
The bare realization or zeroth layer is accomplished by two lasers while the first layer involves and additional laser for protection purposes, Finally, the second layer adds a time-dependent Rabi frequency, The effective magnetic-field fluctuation is described by δ m (t), as shown in Sec. 2 and 3, with parameters τ m = 50 µs and T 2 = 3 ms. Note that distinct experimental setups may suffer different magnetic-field fluctuation, and thus τ m may differ. In this respect, depending on the correlation time τ m , our scheme can be adapted to suppress magneticfield fluctuations by setting properly the Rabi frequencies Ω j , as discussed in Sec. 3. However, for a too short noise correlation time, i.e. in the limit of Markovian noise τ m /T 2 → 0, the tunability of the simulated Rabi models using CCD scheme is reduced as the Rabi frequency must fulfill Ω a > 1/(2πτ m ) to ensure decoupling. We recall that the characteristic frequency from which the spectral density starts to decay as 1/f 2 corresponds to f cr = 1/(2πτ m ), and therefore Ω a > f cr , as explained in Sec. 3. In addition, the fluctuation of the jth laser's amplitude, denoted as δ Ω j (t), is parametrized with τ Ω = 1 ms and c Ω = 2p 2 /τ Ω as it describes a relative amplitude fluctuation, with p = 0.1%. We have considered an equal noise for the lasers with intensities Ω 1 and Ω 2 , i.e. δ Ω 1 (t) = δ Ω 2 (t), while the fluctuations of the rest are completely independent. However, we also performed simulations with uncorrelated noise between Ω 1 and Ω 2 and no significant differences have been observed. In all the simulations, the trap frequency has been chosen as ν = 2π × 1.36 MHz, the Lamb-Dicke parameter as η 1,2 = 0.06 and η a,b = 0.01 [25,27]. Table 1. Trapped-ion parameters to simulate the quantum Rabi model using CCD scheme.
Zeroth layer First layer Second layer

Quantum Rabi model realization
Here we present the numerical simulations of the trapped-ion Hamiltonian realizing the quantum Rabi model to observe the paradigmatic Rabi oscillations. The simulated quantum Rabi model in the ith layer can be written as where σ i TLS and σ i ⊥ stand for the Pauli matrices of the free energy term of the two-level system and the orthogonal direction of the interaction, respectively. The parameters used to simulate this model using Eqs. (25), (26) and (27) are gathered in Table 1, as well as their relation with the effective frequencies given in Eq. (28),Ω i ,ω i andλ i . Note that Ω 1,2 = Ω and η 1,2 = η for zeroth and first layer. In order to achieve the same effective model, regardless of the layer, we will introduce dimensionless constants to define a target Hamiltonian. These are R ≡Ω i /ω i and g ≡ 2λ i /(ω i √ R). Hence, fixing R and g, H R,i /ω i represents the same effective quantum Rabi model. We set ω 0,1,2 =Ω 0,1,2 = 2π × 5 kHz to simulate a resonant case R = 1, and a dimensionless coupling constant g = 1/4. This implies that: (i) for H I 0 , i.e. for the bare realization, δ 2 = 2π × 10 kHz, δ 1 = 0 and Ω 1,2 = 2π × 20.83 kHz; (ii) for H I 1 (first layer) ω 1 = 2π × 5 kHz, Ω a = 2π × 5 kHz and Ω 1,2 = 2π × 20.83 kHz; (iii) for H I 2 (second layer) ω 2 = 2π × 5 kHz, Ω b = 2π × 5 kHz and Ω a = 40Ω b = 2π × 200 kHz, Ω 1 = 2π × 41.67 kHz.
We illustrate how CCD improves the realization of the Rabi model by means of the fidelity among the wavefunction of the ideal Rabi model, |ψ R,i (t) , and its noisy trapped-ion realization |ψ i (t) for the ith layer of protection, which reads We will also compare the oscillations of the population on the excited state of the qubit which is given by σ i TLS + 1 /2 in both cases, ideal and the trapped ion realization with   where σ TLS |↑ TLS = + |↑ TLS and σ ⊥ |↑ ⊥ = + |↑ ⊥ . To the contrary, there are specific situations in which CCD scheme could deteriorate the desired realization. In particular, if the considered initial state is parallel to both magnetic noise, δ m σ z and Hamiltonian (i.e. when we deal with the dark state), to apply CCD scheme is counterproductive since it changes a source of noise, that originally just gives rise to a global phase, to an orthogonal noise producing transitions and distorting the dynamics. This is the case for |ψ(0) = |0 |↓ TLS in the Rabi model when g ≪ 1 i.e. when the Jaynes-Cummings model arises. As we see in Fig. 5, for R = 1 and g = 1/4 the fidelity of the first layer is noticeably worse that an unprotected realization, while the second layer is just as good as the original. This reveals that CCD scheme does not necessarily lead to an improved realization; it depends on several factors which have to be taken into account beforehand.

Critical dynamics of the superradiant quantum phase transition in the Rabi model
In order to illustrate the versatility of the CCD scheme, we analyze the realization of a time-dependent Rabi Hamiltonian in the ultra-strong coupling regime. In this respect, it has been recently shown that the Rabi model (Eq. (28)) undergoes a quantum phase transition in the R = Ω/ω 0 → ∞ limit at the critical point g c = 2λ c / √ Ωω 0 = 1 despite of consisting only of a single two-level system and a single-mode bosonic field [28]. For finite R, critical behavior is revealed in the form of finite-frequency scaling functions, in an approach that is equivalent to finite-size scaling in traditional phase transitions [43,44]. As shown in [29], the presence of the quantum phase transition can be observed with a single trapped-ion that interacts with one of its vibrational modes. This can be achieved resorting to non-equilibrium universal scaling functions [45,29] in terms of the expectation value σ i TLS of Eq. (28), which can be measured with high-fidelity in a trapped-ion system [46,47]. To obtain such non-equilibrium universal scaling functions one can proceed as follows. Prepare an initial state |ψ(0) = |0 |↓ TLS at g = 0 for a fixed R, such that σ i TLS |↓ TLS = − |↓ TLS , and then quench continuously in a time τ Q the coupling constant g until g = g c = 1 is reached. Then, at g(τ Q ) = 1 for a frequency ratio R we calculate the quantity is the ground-state expectation value of σ i TLS at g = 1 and R. The non-equilibrium universal function is found as S(T ) = R µ σ i TLS R where T ≡ R −γ/(µ(1+ζ)) τ Q . The critical exponents are µ = 2/3, γ = 1 and ζ = 1/2 [28,29]. Note however that the driving time τ Q cannot be arbitrarily short since S(T ) is obtained assuming adiabatic dynamics away from the critical point. On the other hand, in an iontrap realization, the duration of the dynamics to reconstruct S(T ) is severely restricted due to the presence of various sources of noise [29].
Here, by applying the CCD scheme, we offer a way to overcome these noises, which facilitates the observation of universal scaling functions, and illustrate that the CCD scheme is valid in an extreme parameter regime and even when quench dynamics is considered. Note however that, due to the large desired value of R, the second layer is expected to fail as R ∝ Ω b but Ω b ≪ Ω a is required to fulfill the RWA. Hence, for this specific case the approximations leading to the quantum Rabi model will break down. The Fig. 6 shows the universal non-equilibrium function S(T ) as a function of the rescaled driving time T . The solid black line corresponds to the ideal quantum Rabi model, while the points to the trapped-ion realization using a first layer protection with R = 50 (circles) and R = 100 (squares) for 0.02 ≤ τ Q ≤ 8.6 in units of 2π/ω i . In the inset the results using zeroth and second layer are plotted. Observe the remarkable improvement compared to the zeroth layer, and the failure of the second layer as Ω b becomes comparable to Ω a . The simulation parameters areω 0,1 = 2π × 1 kHz, ω 2 = 2π × 400Hz, whileΩ i = Rω i . For the second layer Ω a is set to 2π × 200 kHz, and hence Ω a /Ω b = 10 and 5 for R = 50 and 100, respectively, which already provides evidence of the expected failure of the RWA. Additionally, the quench is attained by tuning linearly in time the laser intensities from 0 to Ω f . For the zeroth and first layer, Ω f results in 2π × 117.8 kHz and 2π × 166.7 kHz for R = 50 and R = 100, respectively. For the second layer Ω f amounts to 2π × 94.3 kHz and 2π × 133.3 kHz for R = 50 and R = 100, respectively.

Dirac equation realization in a trapped-ion setting
The parameters to realize the Dirac equation, H D,i /c D = rσ i TLS +pσ i ⊥ with r ≡ m D c D , using Eqs. (25), (26) and (27) are gathered in the Table 2.
In order to observe the paradigmatic Zitterbewegung [25] we calculate the expectation value of the position operatorx = (a + a † ) as a function of time for an initial state |ψ(0) , eigenstate of σ i ⊥ (in particular we consider |↑ ⊥ ). We then set a value m D and c D , or equivalently, r. Note that the presented scheme for first and second layer does not allow for a realization of the strict massless limit, r = 0, since r is proportional to Ω a or Ω b and Ω a,b = 0 does not provide a protected Hamiltonian against fluctuations, while in the zeroth layer, r is just proportional to the detuning δ. Nevertheless, for r > 0, CCD scheme still improves the simulated Dirac equation, as we illustrate in the following. We set r = 2, (i) δ = 2π × 5 kHz, (ii) Ω a = 2π × 5 kHz, (iii) Ω b = 2π × 5 kHz and Ω a = 2π ×200 kHz. This implies (i) for Eq. (25) Ω 1,2 = 2π ×20.8 kHz and ∆ 1,2 = ±ν +δ, (ii) for Eq. (26) Ω 1,2 = 2π × 20.8 kHz and (iii) for Eq. (27) Ω 1 = 2π × 41.7 kHz. In Fig. 7 we plot the fidelity F 0,1,2 (t) (a) and position expectation value x(t) (b) as a function of time. The fidelity corresponds to F i (t) = | ψ D,i (t)| ψ i (t) |, where |ψ i (t) and |ψ D,i (t) are the wave-function of the trapped-ion and ideal Dirac equation of the ith layer, respectively. Note that the final time corresponds to t = 3(2π/c D ) = 2.4 ms. The improvement is clearly shown in Fig. 7. The second layer works worse at longer times than the first one, which is mainly due to laser-amplitude fluctuations and breakdown of RWA (note that Ω a = 40Ω b ). Nevertheless, for shorter times, the simulation of Dirac equation in the second layer is considerably enhanced. Finally we want to comment that the access to motional variables is achieved by, for example, adding a second ion to the trap and computing the time derivative of the qubit expectation value [25,27], see Appendix B for more details. In principle, this protocol requires to prepare the ancillary ion in a certain quantum state that we will select as parallel to the magnetic noise δ m (t). Hence, during the realization of the dynamics, this ion is not affected by external fluctuations, while, for the reconstruction of the time derivatives, a fast evolution is required. In this manner the noise will have an small incidence in the reconstruction of x(t) .

Appendix B. Measurement of vibrational operators
After the system evolution within the CCD scheme we have that the final state is |ψ(t ′ ) . Then, we can use another ion which is initialized into the state |↑ , and therefore does not suffer from the action of a noisy term δ m (t)/2σ A z , where σ A i are the Pauli operators of the ancillary ion. Hence, it does not require CCD protection. After the final time t ′ , a short evolution of time t of the form U = e −iΩtσ A xx is applied to the state |ψ(t ′ ) |↑ . Then, it is easy to demonstrate that ∂ t σ A y t=0 = 2Ω ψ(t ′ )|x|ψ(t ′ ) . (B.1)