A trapped-ion simulator for spin-boson models with structured environments

We propose a method to simulate the dynamics of spin-boson models with small crystals of trapped ions where the electronic degree of freedom of one ion is used to encode the spin while the collective vibrational degrees of freedom are employed to form an effective harmonic environment. The key idea of our approach is that a single damped mode can be used to provide a harmonic environment with Lorentzian spectral density. More complex spectral functions can be tailored by combining several individually damped modes. The protocol is especially well-suited to simulate spin-boson models with structured environments. We propose to work with mixed-species crystals such that one species serves to encode the spin while the other species is used to cool the vibrational degrees of freedom to engineer the environment. The strength of the dissipation on the spin can be controlled by tuning the coupling between spin and vibrational degrees of freedom. In this way the dynamics of spin-boson models with macroscopic and non-Markovian environments can be simulated using only a few ions. We illustrate the approach by simulating an experiment with realistic parameters and show by computing quantitative measures that the dynamics is genuinely non-Markovian.


Introduction
The spin-boson model is an archetypical model of an open quantum system. It is applied in numerous contexts ranging from chemical reactions [1] and biological molecular aggregates [2] to solid state physics [3][4][5]. The model describes a single spin coupled to a dissipative environment formed by an infinite set of harmonic oscillators. It is well-known that the effect of thermal oscillator environments on a quantum system is fully described by a single scalar function, the spectral density (also spectral function) of the environment [4]. Although approximate analytic solutions have been found for some spectral densities [3,4], no closed analytic solution of the spin-boson model is known. Meanwhile, dynamics and thermodynamical properties of spinboson models have been investigated by a number of numerical approaches including techniques based on the numerical renormalization group [5], time-dependent density matrix renormalization group [6,7], path integral Monte Carlo [8], or the quasi-adiabatic propagator path integral approach [9]. Numerical simulations are especially needed for environments with spectral densities where the reorganization energy is of the order of the spectral width or for highly structured environments with long-lived vibrational modes. In these cases, strong system-environment correlations can lead to highly non-trivial dynamics and it is known that reduced effective models do not represent the dynamics faithfully [10,11]. Moreover, these types of spectral densities, which are of particular relevance for the excitonic and electronic dynamics in biomolecular systems [12], pose considerable challenges for the numerical methods. One example for this is the prediction of the results of nonlinear spectroscopy which is exceedingly hard on conventional computers already for small systems [13]. Therefore, an experimental simulator of spin-boson models with a high degree of control is desirable.
Trapped atomic ions provide a clean and highly controllable system where many dynamical quantities are directly accessible. They have proven to be a versatile platform for the simulation of a wide range of physical models including defect formation in classical phase transitions [14][15][16] as well as open and closed quantum systems [17][18][19][20][21][22]. The simulation of spin-boson models using trapped atomic ions has been proposed previously [23] requiring rather large crystals comprising 50-100 ions. These crystals feature a large number of vibrational modes which can be used to act as a mesoscopic environment for the spin. However, for such large crystals the level of control needed to simulate spin-boson models is experimentally very hard to achieve.
In this work, we develop a proposal to simulate the dynamics of spin-boson models with small crystals of trapped ions. Our procedure also relies on the vibrational degrees of freedom to act as the environment, but it makes use of the fact that a damped mode produces a Lorentzian spectral density [1]. This result was derived in [1] assuming that the damping is provided by an oscillator reservoir with Ohmic spectral density. Cooling of trapped ions, however, is usually described by a Lindblad equation. Here, we show that the Lorentzian spectral density can also be obtained for appropriate parameters if the damping is modeled by a Lindblad equation, extending the results of [24][25][26][27]. Combining several damped modes, arbitrary spectral densities can be constructed. The shape of the tailored spectral density is controlled by the couplings of the spin to the modes, the mode frequencies and the damping rates. Accordingly, both the shape of the spectral density and the strength of the dissipation can be tuned by the experimenter. As we show below, the effective Lorentzian spectral density can typically only be attributed to the damped oscillator in Lindblad description if the damping rate is considerably smaller than the mode frequency. Therefore, our protocol is especially well-suited for the simulation of spinboson models with environments that feature structured spectral densities. The spectral densities constructed with the protocol are continuous functions of frequency and can thus be identified with an environment made up of a macroscopic number of modes as it occurs in the condensed phase. Accordingly, with our method one can tailor environments with continuous and highly structured spectral densities using only a small number of oscillators to form the environment.
It is interesting to note that this approach is useful already for a small number of modes that are used to model the environment because the direct numerical simulation of a spin coupled to damped modes in Lindblad description becomes inefficient already for a few modes. The reduced overhead of our protocol brings the simulation of spin-boson models and the prediction of nonlinear spectroscopy of such systems to the realm of state-of-the-art trapped-ion setups.
The article is structured as follows. In section 2, we introduce the spin-boson model and the concept of the influence functional. Then, in section 3 we discuss the spectral densities generated by damped oscillators in different models of damping. Based on these results, we introduce our protocol for the simulation of spin-boson models and benchmark the procedure with a comparison to a numerically exact simulation of the full manybody dynamics in section 4. We then proceed to illustrate how the protocol can be implemented with currently available ion trap experiments. In section 5, we show that the protocol is robust to the most common experimental sources of noise. Finally, in section 6 we show that the dynamics in our simulations are truly non-Markovian by computing two quantitative measures of non-Markovianity. We close with a summary of our results and discuss future perspectives for simulations of spin-boson models with trapped ions in section 7.

Spin-boson model
The spin-boson model describes a two-level system (spin 1/2) in a dissipative environment which is modeled by an infinite set of non-interacting harmonic oscillators. Denoting the energy splitting between the spin states ñ | and ñ | by ò and the coupling between them by  D, the Hamiltonian of the global system reads [3] a a n n denotes the raising (lowering) operator of environmental mode n and ω n its frequency while the real λ n describe the couplings of the spin to the environmental oscillators. The spectral density which determines the influence of the oscillator environment on the spin [3,4] reads with δ the Dirac δ-function. For a macroscopic environment, one assumes that the frequencies are so closely spaced that J(ω) becomes a continuous function of ω.
One is generally interested in finding the reduced dynamics of the spin for an environment with a certain spectral density. The path integral formalism [28] provides us with an exact expression for the propagator of the spin state where the effects of the environment are already included. For factorizing initial conditions r r r = Ä b s 0 with some spin state ρ s and the environmental modes in a thermal state ρ β at inverse temperature b = -( ) k T B 1 , the propagator for the spin reads [29] Here, the path integral ò Dq q q f 0 runs over all spin state trajectories connecting q(0)=q 0 and q(t)=q f , S 0 [q] is the action of the free spin evolution and F [q, q′] is the Feynman-Vernon influence functional [29]. The influence functional contains the effect of the environment on the spin dynamics. For an oscillator environment and the considered coupling it can be written as [3] * ò ò X a a n n n n . Note that the correlation function inequation (5) is evaluated with respect to the free evolution of the environmental oscillators. Alternatively, L(t) can be expressed in terms of the spectral density J(ω): Hence, we see that the influence of the environment on the spin is equivalently given either by the coordinate correlation function of the environment or by its spectral density.

Environments of damped harmonic oscillators
The key idea underlying our proposal for the simulation of spin-boson models is the result that a harmonic oscillator damped by an oscillator bath with Ohmic spectral density yields an effective environment with Lorentzian spectral density [1]. In this section, we will analyze this effective spectral density and compare it to that of a damped harmonic oscillator in Lindblad description. In particular, we show that the two models of damping yield the same effective spectral density for appropriate parameters. To this end, we analyze the correlation functions of the two reservoirs and show when the correlation functions and their spectral representations coincide. The observation that different environments that produce the same influence functional have the same effect on the dynamics of a reduced system [29] completes the argument that also a damped oscillator in Lindblad description can act as an effective harmonic environment. Let us start by considering an environment that consists of a single harmonic oscillator of free oscillation frequency Ω which is damped by an oscillator reservoir with Ohmic spectral function Here, K is a constant and ω c a high-frequency cutoff. For brevity, we call this damped harmonic oscillator 'Ohmic oscillator' in the following. For a strictly Ohmic environment, i.e. w  ¥ c , that causes damping at rate κ on the coordinate of the oscillator, the effective spectral density generated by the damped oscillator on the spin is Here, w k = Wm 2 2 is the reduced frequency of the damped oscillator and l the spin-oscillator coupling as inequation (1). Note that we restrict our considerations to the underdamped regime k < W.
In trapped-ion experiments, the motion of the ions is usually expressed in terms of a set of normal modes, each of which is a harmonic oscillator. Cooling of the modes is commonly described by a Lindblad equation [30,31]. Therefore, it is not immediately clear if we can obtain an effective spectral density as for the Ohmic oscillator,equation (7). We will now show that this is possible and that we obtain the same spectral function for appropriate parameters. It is easier to start by considering the reservoir correlation functions in the time domain in order to establish the correspondence between the effective environments for the two models of damping.

Reservoir correlation functions
The reservoir correlation function L(t) inequation (5) may be written explicitly in terms of the environmental coordinate correlation functions. To see this, note that we can write . Using the above identity, we may write the reservoir correlation function inequation (5) as where we have used that the oscillators are independent. We note again that for the coupling between spin and environment in the Hamiltonian inequation (1) the coordinate correlation function of the free environmental oscillators determines the influence on the spin. In the following, we consider only a single oscillator and therefore we omit the index n. Since the coordinate correlation function á , also L(t) is a complex-valued function with real and imaginary parts L′(t) and ( ) L t . For the oscillator damped by a strictly Ohmic bath, the coordinate correlation function and thus L(t) can be calculated analytically [3,32]. For a bath at inverse temperature β, we obtain Note that we again assumed the underdamped regime κ<Ω here. The reservoir correlation function for the Ohmic oscillator can also be written in the form ofequation (6) where J eff (ω) is given inequation (7). In Lindblad description, a damped harmonic oscillator coupled to a thermal reservoir at inverse temperature β evolves according to where a, a † are the ladder operators of the damped mode and the frequency ω m is taken to include possible renormalizations due to the damping. The dissipator reads [33] D r k r r k r r Note that the above form of the Lindblad equation is also used to describe laser cooling of trapped ions [30,31]. In the following, we call the damped oscillator where the damping is described by the Lindblad equation (16) 'Lindblad oscillator'. Employing the quantum regression theorem, we can compute the coordinate correlation of the Lindblad oscillator. Inserting the resulting expressions intoequation (9), we obtain the reservoir correlation function of the Lindblad oscillator The real part reads t L 2 m m and the imaginary part is given by Comparing equations (14) and (20), we find that the imaginary parts of the reservoir correlation functions L(t) and L L (t) are exactly equal. Note that here we have tacitly assumed that the damping rates κ and the reduced oscillator frequencies ω m are the same in both cases. However, despite this assumption ¢ ( ) L t L has a different functional form than L′(t) inequation (12).
we can also write ( ) L t L in terms of J eff (ω) ofequation (7) by taking the imaginary part ofequation (6). This is, however, not the case for ¢ L L (t). Writing ¢ L L (t) as inequation (6), we obtain  Hence, L L (t) in general cannot be written in terms of a single spectral density but takes the form Despite the differences, it is possible to obtain very good agreement also between the real parts L′(t) and ¢ L L (t) and their frequency space representations equations (7) and (21). References [32,34] estimate that the quantum regression theorem can only yield quantitatively correct predictions for the two-time correlation functions of the damped harmonic oscillator if k w  m and bk  1. Indeed, under these assumptions we find very good agreement between both L′(t) and L L ′(t) and their frequency space representations w ( ) J eff and ¢ J eff (ω). We can understand the two conditions as follows. The Lindblad equation (16) with the dissipator inequation (17) is a good description for the damped oscillator for weak coupling between the oscillator and its environment which is reflected by the condition k w  m . Furthermore, at very low temperatures the decay of L′(t) is dictated by the Matsubara frequencies in L 2 (t). This decay cannot be reproduced by ¢ L L (t) which only features a single decay rate. For 1, 23 1 the smallest Matsubara frequency is much larger than the decay rate κ and we can neglect L 2 (t) if we are interested in not too short time scales [3]. In this case, L′(t)≈L 1 (t) which can be correctly reproduced by the regression theorem result ¢ L L (t) in the considered parameter regime, as we will see shortly. Recalling that β=1/(k B T), we find that the condition inequation (23) puts a lower bound on the temperature where the identification of L(t) and L L (t) is possible for a fixed cooling rate. However, also too high temperatures lead to deviations such that there is an intermediate temperature range where the best agreement is achieved (see appendix A for a more detailed discussion).
Thus, we have established a regime where the coordinate correlation function of the Lindblad oscillator approximately coincides with that of the Ohmic oscillator. Summarizing, we require for the functions L L (t) and L(t) to coincide. In this regime, the Lindblad oscillator produces the same reservoir correlation function and thus influence functional as the Ohmic oscillator. According to the equivalence theorem in [35], in this regime the Lindblad oscillator acts as an effective macroscopic reservoir with Lorentzian spectral density as given inequation (7) above. For ion-trap experiments, one usually considers the mean occupation numbern of the bosonic modes rather than their temperature and therefore it is desirable to cast condition(23) in a form where it depends onn. Assuming a thermal state for a bosonic mode, we can associate the temperature In order to make our considerations more quantitative and to illustrate that the match of the reservoir correlation functions is indeed very good for parameters of interest, we make a numerical comparison of the functions L(t) and L L (t) in the regime k w n  , m 1 . Since the imaginary parts of the two functions are equal, we focus on the real parts L′(t) and L L ′(t). In figure 1 we plot L′(t)/λ 2 including the first 10 4 Matsubara frequencies together with L L (t)/λ 2 for the parameters specified in table 1 below. These parameters are realistic for an ion trap experiment. In part (a) of the figure we compare L′(t) and L L ′(t) on short and in part (b) on intermediate time scales. One can appreciate excellent agreement between the two functions.

Frequency space
In the previous section, we have seen that in the parameter regime specified inequation (24) the Lindblad description of the damped harmonic oscillator reproduces the reservoir correlation function L(t) of the Ohmic oscillator. In fact, in almost all cases environments are characterized by their spectral density rather than their correlation function L(t), which is obtained from the spectral density J(ω) throughequation (6).
The frequency space representation for the real part of the correlation function of the Lindblad oscillator is given inequation (22). Comparing the functions J eff (ω) and ¢ J eff (ω) from equations (7) and (21), we find that in general and that we hence cannot write L L (t) as a function of a single spectral density as inequation (6), in general. Yet, from our considerations in the previous section we expect that for parameters satisfyingequation (24) we have In this case, the reservoir correlation function can also be written in the form ofequation (6) as for a macroscopic oscillator environment.
In figure 2 we compare the left and right hand sides ofequation (26) for the parameters specified in table 1 for which we found excellent agreement between the correlation functions L L (t) and L(t) (see figure 1). Panel (a) shows J eff (ω) (solid line) and ¢ J eff (ω) (triangles) for small frequencies and part (b) shows the behavior around the resonance ω m /2π=100 kHz. Both parts of the figure show that we obtain very good agreement in frequency space, too. Part (c) of the figure shows the relative error  . Part (a) shows the behavior for small frequencies while part (b) depicts the two functions around the resonance ω m /2π=100 kHz. In part (c) we show the relative error ò J fromequation (28) over the relevant frequency range covered by the spectral density. Table 1. Parameters for the simulation of a spin-boson model with Lorentzian spectral density with trapped ions.
which is remarkably small over the whole range ω/2π=0-150 kHz. Note that the increase in the relative error for higher frequencies is because the spectral density J eff (ω) goes to zero more rapidly than w ¢ ( ) J eff . However, since both contributions are small, the effect of this difference should be negligible as long as the frequency of the spin coupled to this effective environment does not lie in this range.
In summary, we confirm the result of the previous section: for appropriate choices of mode frequency, cooling rate and temperature, the damped oscillator evolving according to the Lindblad equation can be attributed the effective spectral density J eff (ω),equation (7), of a macroscopic oscillator environment. Note that the treatment is not perturbative in the spin-motion coupling λ, so that this equivalence is valid for arbitrary values of λ as long as the Lindblad equation holds.

Trapped-ion simulations of spin-boson models
In this section, we introduce our protocol for the simulation of spin-boson models. We illustrate the procedure for a trapped-ion experiment but it can also be adapted to other experimental platforms.

Simulation protocol
The Hamiltonian of the spin-boson model we want to simulate is given inequation (1). The influence of the environment on the spin dynamics is determined by the spectral density inequation (2) which we assume to be a smooth function of ω here.
In the previous section, we have seen that a Lindblad oscillator yields an effective environment with spectral density J eff (ω) fromequation (7) if the parameters satisfy the conditions inequation (24). Let us now assume that a spin is coupled to N independent damped harmonic oscillators in Lindblad description that satisfy the constraints inequation (24). Then, a spectral density J eff,n (ω) given byequation (7) with the corresponding λ n , κ n , ω n can be attributed to oscillator n, n=1, K, N.
The combined influence functional of different statistically and dynamically independent environments is given by the product of the individual influence functionals [29]. Using this property for the environment composed of the N independent Lindblad oscillators yields the reservoir correlation function Here, β n is the temperature of the reservoir associated with oscillator n. If all reservoirs have the same temperature, i.e. β n =β for n=1, K, N, their spectral densities add up and one can construct effective spectral densities Hence, the N independent Lindblad oscillators yield an effective environment with the spectral density J(ω) ofequation (30). If the number of available oscillators is not restricted, any spectral density can be decomposed as inequation (30). In the limiting case of infinitely many oscillators and vanishing damping,equation (30) yieldsequation (2). Of course, in practice the number of oscillators N is finite. Yet, this still allows us to create a large variety of spectral densities.
If we want to approximate a certain target spectral density J T (ω) with N oscillators, we can find the values for the λ n , κ n , ω n , n=1, K, N that reproduce the desired spectral density by minimizing the functional [36] Note that the optimization is subject to the conditions inequation (24) if we use Lindblad oscillators to form the environment. Therefore, if we use Lindblad oscillators and want to employ as few oscillators as possible, the procedure works best for structured environments.
Summarizing the above idea, in order to simulate a spin-boson model described by the Hamiltonian inequation (1) with a spectral density J(ω), we have to engineer a physical system in such a way that it evolves according to where the conditions inequation (24) are satisfied for each mode and the spectral densities associated with the damped modes fulfillequation (30). The HamiltonianH sb is the same as that inequation (1) but the index n now only runs over the set of damped modes. In order to confirm the above statement, we simulated the dynamics of s á ñ ( ) t z and s á ñ ( ) t x for the full spinboson Hamiltonian inequation (1) with spectral density J eff (ω) fromequation (7) using the numerically exact TEDOPA algorithm [6] and compared them with those given byequation (32) withH sb for a single mode. Details regarding the TEDOPA simulation technique and its implementation are given in appendix B. We considered an initial product state where ρ β is a thermal state with b =´-5.91 10 s 6 which corresponds to w = ( ) n 0.025 m for the Lindblad oscillator. Furthermore, we take ò=0, ω m /2π=100 kHz and κ/2π=1.25 kHz. These parameters are summarized in table 1. We chose a spin-mode coupling λ/2π=100 kHz and computed the evolution for spin energies Δ/2π=−50 and −100 kHz. The results are displayed in figure 3. For both values of Δ we obtain very good agreement which shows that the analogy to the macroscopic environment also holds when we probe the effective spectral density generated by the Lindblad oscillator away from the resonance and with a nonperturbative coupling. Note that one simulation for Δ/2π=−50 kHz took 15 days using 16 cores on a computing cluster which once more indicates the value of a trapped-ion simulator, especially for structured environments.

Ion trap implementation
Let us now proceed to illustrate how the ideas discussed above can be implemented in an ion-trap experiment. We consider N singly charged atomic ions with masses m j confined in a linear Paul trap with effective harmonic trapping potential. We assume trapping conditions such that laser cooled ions form a linear Coulomb crystal along z with equilibrium positions = ( ) z r 0, 0, j j T 0 0 . The motional degrees of freedom can then be described in terms of N uncoupled normal modes in each spatial direction [37,38]  where ω n, α is the frequency of mode n in spatial direction αä{x, y, z} with ladder operators a a † a a , n n , , . For simplicity we focus on the case of a spin coupled to a single damped mode which corresponds to a spinboson model with Lorentzian spectral density as inequation (7). This system already exhibits an interesting phenomenology and has been studied with a variety of numerical and analytical approaches, see e.g. [39][40][41][42]. For this purpose, we only need N=2 ions: one ion is used to encode the spin while the other ion provides sympathetic cooling of the shared modes of motion. In order to avoid that the cooling lasers couple to the spin transition, we choose to work with mixed-species ion crystals. Alternatively, one could rely on single site addressing. The internal levels of the spin ion are described by the Hamiltonian  For concreteness we consider a crystal composed of 24 Mg + and 25 Mg + . 25 Mg + has a nuclear spin and we can use the states = = ñº ñ of the 2 S 1/2 electronic hyperfine ground-state manifold to encode the spin. The spin can be driven either by a microwave or in a two-photon stimulated Raman configuration while the desired coupling of the spin to the motional degrees of freedom in the σ z basis is provided by a 'walking standing wave'. In this configuration the spin states are off-resonantly coupled to the P manifold by two laser beams near 280 nm whose beat note is tuned close to one of the motional mode frequencies [43]. The interaction of the spin ion with the applied fields is described by (see appendices C and D) where Ω d is the Rabi frequency of the applied microwave or stimulated Raman field and ω d ≈ω 0 its frequency. Ω odf , k L , ω L , f L are the effective laser Rabi frequency, wave vector, frequency and phase, respectively and r 2 denotes the position of the spin ion. Directing k L along z the laser only couples to the motion along this axis. A two-ion crystal features two axial modes, an in-phase and an out-of-phase mode of motion with frequencies ω 1,z ≡ ω 1 and ω 2,z ≡ ω 2 . The two modes are well separated in frequency such that choosing the laser frequency ω L ≈ω 2 the spin only couples to the ouf-of-phase mode. In an interaction picture rotating with the microwave and motional frequencies and under the rotating wave approximation, the system's Hamiltonian reads (see appendix D) where δ=ω 0 −ω d is the detuning of the field driving the spin transition and d w w w = - m 2 L 2 the detuning of the laser from the motional mode. The spin-motion coupling is given by 22 L . Note that the laser phase can be chosen such that λ is real.M 22 is the out-of-phase mode amplitude at the spin ion in mass weighted coordinates and m 2 its mass. Identifying   d = , Ω d =−Δ and δ m =ω m , we obtain the desired Hamiltonian, namely the spin-boson Hamiltonian ofequation (1) for a single mode. Adding the cooling on the second ion, the full system evolves according toequation (32). This is the desired time evolution for a simulation of the spin-boson model with a Lorentzian spectral density as inequation (7).
We simulate the dynamics of the system for experimentally realistic parameters. We consider an axial potential where a single 24 Mg + ion has a center-of-mass frequency ω com /2π=2.54 MHz. This potential leads to an out-of-phase mode frequency ω 2 /2π=4.36 MHz and η 2 ≈0.15 for the mixed crystal where we assumed that the lasers inducing the spin-dependent force are at right angles. Furthermore, we assume that EIT cooling [31] is applied to the 24 Mg + ion which has already been used to sympathetically cool mixed-species ion crystals [44]. We assume a cooling rate 2κ/2π=2.5 kHz and a steady-state population = n 0.025 of the mode which is realistic in light of the results in [44]. Note that one has to make sure that the conditions inequation (24) hold for the effective mode frequency ω m =δ m , which is the detuning of the spin-motion coupling and thus much smaller than the physical mode frequency. We choose the field driving the spin to be resonant, i.e. ò=0, and a detuning ω m /2π=100 kHz of the spin-motion coupling. Accordingly, we recover the parameters of table 1 and the correspondence holds. Note that experimentally a finite bias ò can easily be included by introducing a detuning to the field driving the spin transition. In the simulations, we truncate the motional Hilbert space at n max =15 excitations which makes truncation errors negligible.
In figure 4 we show the dynamics of s á ñ ( ) t z underequation (32) where = H H sb sb1 fromequation (37) with the parameters of table 1 for an initial state as inequation (33). We vary the spin-motion coupling λ/2π=10-200 kHz. In panel (a) we show the dynamics for Δ/2π=3 kHz. In this case the spin samples the low frequencies of the spectral density inequation (7). For small ω the spectral density shows Ohmic behavior J eff (ω)∼ω. We observe a transition from damped to overdamped oscillations with increasing spin-mode coupling λ. This behavior is expected for an Ohmic spectral density at finite temperatures [3,4]. Note, however, that our spectral density J eff (ω), even if Ohmic for small frequencies, does not yield the same correlation function as a strict Ohmic environment. Therefore we can only expect qualitatively similar dynamics [39,40]. In panel (b) we show s á ñ ( ) t z for Δ/2π=100 kHz such that the spin is resonant with the mode. The remaining parameters and the initial condition are the same as in part (a). In this regime, the spin dynamics shows a very complex behavior which one would intuitively call non-Markovian. In order to verify that the observed dynamics is truly non-Markovian, we computed two quantitative measures of non-Markovianity. The results are presented in section 6 below. For now, we remark that the two measures witness non-Markovianity for both the resonant and the Ohmic case for l ¹ 0.

Impact of experimental sources of noise
In this section, we present an estimate of the impact of typical experimental sources of noise on the quality of the simulations of the spin-boson model dynamics. We consider two different types of noise. First, we consider the impact of a generic dephasing noise on the experimental results. Although it can be suppressed very well in many experiments, dephasing is ubiquitous in trapped-ion experiments. The second source of noise that we consider is related to the concrete implementation of the protocol that we propose. The σ z spin-motion coupling via optical fields utilizes coupling of the spin states via a decaying state. In general, this type of coupling leads to decoherence due to residual off-resonant excitation of the upper level (see appendix C).

Dephasing noise
We consider the implementation of the spin-boson Hamiltonian as presented in the previous section and assume that the spin levels are additionally subject to dephasing noise. Including noise effects on the spin, the state ρ of the spin ion and the relevant damped mode evolves according to where H sb1 is the spin-boson Hamiltonian for a single mode ofequation (37). The dissipator now consists of two parts: D kn , describes the damping of the mode and is given inequation (17) figure 5, we compare the noise free dynamics of the spin, given byequation (32) above, with those given byequation (38) including the dephasing noise inequation (39). The parameters are again those of table 1 and we use the initial condition ofequation (33). We vary the spin-motion coupling λ/2π=10-200 kHz. In parts (a) and (b) of figure 5 the spin energy is Δ/2π=3 kHz. Part (a) shows the dynamics for * = T 1 2 ms and part (b) shows the dynamics for * = T 10 2 ms. These coherence times have already by far been surpassed with magnetic field sensitive trapped-ion qubits, see e.g. [45], where coherence times of 300 ms have been observed. The figure shows that for * = T 10 2 ms there is already no appreciable effect in the dynamics for the Ohmic case. For the resonant case Δ/2π=100 kHz there is no visible effect already for * = T 1 2 ms as can be appreciated in part (c) of the figure. This is due to the much shorter time scale in this case. In the light of these results, it seems fair to neglect dephasing noise.

Decoherence due to σ z spin-motion coupling
The use of optical fields to obtain the σ z spin-motion coupling introduces additional decoherence on the spin through off-resonant scattering of photons from the applied beams. We analyze the effects of this type of decoherence on the quality of the spin-boson model simulations in this section. Motivated by the results of the previous section, we neglect dephasing noise in our analysis. The dissipative effects of the considered spin- fromequation (37) for the simulation parameters specified in table 1 and the initial condition ofequation (33). This setting corresponds to a spin-boson model with a Lorentzian spectral density as inequation (7). We vary the spin-motion coupling λ. The values of λ/2π in kHz are given in the bars on top of the plots. In part (a) the spin energy Δ/2π=3 kHz is much smaller than the mode frequency ω m /2π=100 kHz, so that the environment is approximately Ohmic. In part (b) the mode is resonant with the spin (Δ/2π=100 kHz). motion coupling only act on the spin. Accordingly, the system of spin and mode again evolves according toequation (38) but with a different dissipator D s than that ofequation (39).
In order to find expressions for the Lindblad operators that describe the dissipative effects due to the σ z spinmotion coupling, we use a simplified three-level model for the internal structure of 25 Mg + . The model and the calculations to obtain the effective Lindblad operators are summarized in appendix C. With the effective Lindblad operators we find, the spin dissipator now reads Here, the Lindblad operators L jk are given by å å where the index l runs over the applied laser beams. We consider that two laser beams are applied and Ω l,s is the Rabi frequency of laser l coupling spin state s to the upper level. w D W G  , , l s R , 0 denotes the detuning of the beams from the excited state and the Γ s are the decay rates from the upper level to spin state s.
The Lindblad operators inequation (41) describe Rayleigh scattering where the spin state is not altered upon a scattering event but can introduce dephasing. Those inequation (42) describe Raman scattering where the spin state is changed upon a scattering event. If we assume that the modulus of the Rabi frequencies of the two lasers providing the spin-dependent force is approximately equal W » W | | l s , 0 , we can estimate the effective scattering rate Γ eff ≈Γ Ω L /Δ R , where W = W D ( ) 2 L 0 2 R is the approximate effective laser Rabi frequency. Hence, decoherence can be largely suppressed if we choose Δ R large enough.
For the spontaneous emission rate and the Lamb-Dicke parameter we take the parameters of 25 Mg + Γ/2π= 41.4 MHz and η≈0.15, which we also used in the previous section. Furthermore, we assume an equal branching ratio for the decay of the excited level to the spin states and assume that all laser Rabi frequencies have the same modulus W = W | | l s , 0 . In figure 6 we present the dynamics resulting fromequation (38) without the spin dissipator D s , as presented in figure 4, and compare them to the dynamics incorporating the effects of D s from equation (40). Again, we use the parameters of table 1 with the initial condition inequation (33) and vary the spin-motion coupling λ/2π=10-200 kHz. Parts (a) and (b) of figure 6 compare the dynamics for the quasi Ohmic case Δ/2π=3 kHz and laser detunings Δ R /2π=100 GHz and Δ R /2π=1 THz, respectively. The solid lines represent the unperturbed dynamics while the symbols incorporate the effects of the additional decoherence on the spin. For a detuning of 100 GHz, the spin dynamics is noticeably perturbed already for weak spin-motion couplings. Yet, the right qualitative behavior of the dynamics is preserved. For the large detuning, there is only an appreciable effect for the strongest spin-motion coupling. Finally, part (c) of the figure shows the dynamics for the resonant case Δ/2π=100 kHz and a laser detuning Δ R /2π=100 GHz. In this case, there is no appreciable effect on the spin dynamics, again due to the much shorter time scale.
The results of figure 6 show that, as we expected, errors due to the σ z spin-motion coupling can be suppressed to a large extent if the laser detuning can be chosen large enough. In order to avoid this source of error, one could also rotate the spin basis and provide spin-motion coupling in a different basis, for instance by a Mølmer-Sørensen interaction [46]. Another way to circumvent the considered type of noise would be to use the spin-motion coupling induced in the near field of microwave currents [47], where the considered type of spinmotion coupling is also available but spontaneous emission is negligible.

Quantification of the degree of non-Markovianity of the dynamics
In this section, we investigate the non-Markovian character of the dynamics presented in figure 4. There are several different ways to define non-Markovian dynamics. Here, we compute two quantitative measures of non-Markovianity: N RHP and N BLP as presented in [48,49], respectively.
First we analyze the non-Markovian character of the dynamics presented in figure 4 according to the measure N RHP . To this end, let us consider an open quantum system of finite dimension d whose time evolution is described by a completely positive and trace preserving dynamical map E t t , 0 . For an initial state ρ(t 0 ), the system's state at a later time  t t 0 is given by According to [48], the dynamical map describes a Markovian evolution if and only if the map E t t where Here,   ... 1 denotes the trace norm and yñ = å ñ = | |n n , is a maximally entangled state of the open system with an ancillary system of the same size. E   y y Ä ñá is the so-called Choi matrix and is positive if and , is completely positive [50]. Note that is completely positive. Thus, for a Markovian dynamics g(t)=0 for all times and N RHP evaluates to zero.
We evaluated N RHP numerically for the spin-boson system consisting of a spin coupled to a damped mode described byequation (32) with the Hamiltonian inequation (37). We considered the parameters that we used to produce figure 4, which are given in table 1, and the initial state inequation (33). The numerical computation of N RHP requires the evaluation of a discrete version ofequation (45). Since the evaluation of the measure is numerically demanding for the considered case, we restricted the time intervals that we inspected to T=0.01/Δ for the 'Ohmic' case (Δ/2π=3 kHz) and to T=0.1/Δ for the resonant case (Δ/2π=100 kHz). In both cases, we divided the time intervals in N=10 4 steps to approximateequation (45). More details regarding the numerical evaluation of N RHP are given in appendix E. The results of these computations are shown in figure 7(a). Figure 6. Impact of the dissipative effects of the σ z spin-motion coupling using optical fields on the spin-boson model simulation. In the figure the dynamics without noise on the spin, given byequation (38) without  s , are represented by solid lines and are the same as that in figure 4. The symbols depict the dynamics incorporating noise as given byequation (40). Part (a) shows the dynamics for the quasi Ohmic case Δ/2π=3 kHz and a Raman beam detuning Δ R p /2 =100 GHz. Part (b) shows the same as part (a) for Δ R /2π=1 THz. In part (c) we show the dynamics for the resonant case Δ/2π=100 kHz and Δ R /2π=100 GHz. The values of λ/ 2π in kHz are given in the bars on top of the plots. The remaining simulation parameters are those used for figure 4 and are reported in table 1 and the text.
Let us now turn to the measure of non-Markovianity N BLP [49]. The computation of N BLP is somewhat easier than that of N RHP . N BLP was originally proposed as a measure of non-Markovianity based on the monotonicity of the trace distance under completely positive and trace preserving evolutions and is given by and (· ·) D , is the trace distance. The integral extends over those subintervals of I where σ (t)>0. Thus, N BLP detects non-Markovianity of a dynamical map E t t , 0 if the trace distance between two initial states ρ 1 and ρ 2 increases in the course of the dynamics induced by E t t , 0 . A nonzero value of N BLP can be associated with a backflow of information from the environment to the system [49]. It is known that optimal state pairs ρ 1 , ρ 2 that saturate the maximum inequation (46) are orthogonal and lie on the boundary of state space [51]. However, since we only want to witness non-Markovian dynamics we do not need to perform the maximization inequation (46). We can provide a useful lower bound on N BLP by computing the measure for the eigenstates  ñ | , ñ | x and ñ | y of the Pauli matrices σ z , σ x and σ y , respectively, as initial states. For the numerical computation of N BLP we considered the whole interval [0, 20/Δ] for both values of Δ. We considered N=10 4 equally spaced points t i in that interval and computed the time evolution for the spin starting in each of the eigenstates of the Pauli matrices. We then computed the discrete version of N BLP for the pairs of eigenstates belonging to the same Pauli matrix. Here, the sum runs over those i where the term in brackets is larger than zero and . We note that due to the finite number of 'measurements' there will be a small deviation from the true value of N RHP [52]. The results for the initial state pairs that led to the largest values of N BLP are shown in part (b) of figure 7. The values for the 'Ohmic' case are obtained for the initial spin states r = ñá ( ) | | 0 x s and in the resonant case for the initial spin states r = ñá ñá ( ) | | | | 0 , s . In both cases the measure is non-zero for all couplings λ/2π>0 . An evaluation of N RHP requires process tomography and is therefore experimentally time-consuming already for a single spin. Hence, it might be easier to experimentally detect non-Markovian dynamics using N BLP which only requires state tomography. We remark that N BLP witnesses non-Markovianity in all regions where N RHP does. The somewhat discontinuous behavior of N BLP for the resonant case is due to the finite time interval we are sampling. Also note that N BLP , as we consider it here, is not normalized. Accordingly, since the time interval we are considering in the 'Ohmic' case is much longer than that for the resonant case, we cannot compare the degree of non-Markovianity of the two cases for N BLP .

Conclusions and outlook
In summary, our work provides a route towards the physical simulation of spin-boson models with continuous spectral densities using damped oscillators in Lindblad description. Due to the constraints that have to be satisfied such that we can attribute a continuous Lorentzian spectral density to the damped oscillator in Lindblad description, the protocol we have developed is most promising for the simulation of structured environments. For these environments, our protocol has the potential to achieve a significant reduction of the technical requirements for the implementation of this paradigmatic model for decoherence and dissipation employing trapped ions. The joint effect of different damped modes allows one to tailor a large variety of spectral densities with rich non-Markovian features. We showed that it is possible to carry out simulations of non-trivial dynamics making use of just one motional mode, and illustrated the practicality of our approach by simulating an experiment with realistic parameters.
In order to tailor more complex spectral densities than in this simulated proof-of-principle experiment, one would need to couple the spin to two or more damped modes with the appropriate couplings and cooling rates that match the effective spectral density to the desired one. In case several modes are used, it could be advantageous to use the transverse modes of motion. Due to the smaller bandwidth of the transverse phonon frequencies it is easier to couple to and cool several modes at the same time. It should be borne in mind that the cooling rates should be considerably smaller than the spacing between modes. Only then the damping of each mode can be described by a dissipator as inequation (17). In order to fill possibly unwanted gaps in the effective spectral density, one could then use the modes of the second transverse direction of motion and place the effective frequencies of these modes between those of the first direction.
Let us finally note that the model can be extended in two ways. More complex spectral densities can be obtained by including more modes by either adding more coolant ions or coupling the spin to the modes of more than one spatial direction. More spins can be included by adding more spin ions. Spin-spin interactions are nowadays routinely implemented such that models of interacting spins coupled to a dissipative environment can be realized. Then, trapped ions could be used as a testbed for the dynamics of exciton transport in complex spectral densities as it occurs in photosynthetic pigment protein complexes. Especially higher order spectral responses, e.g. 2D electronic spectroscopy, of these systems are exceedingly hard to compute numerically even for only a few electronic sites coupled to an environment with structured spectral density [13]. In this way, trapped ions could contribute to the understanding of the physical mechanisms underlying photosynthesis.

Appendix A. Correlation functions of Ohmic and Lindblad oscillator
In this appendix, we show how one can see that the reservoir correlation functions L(t) and L L (t) for the Ohmic and Lindblad oscillators coincide when the conditions inequation (24) are satisfied. The real and imaginary parts of L(t) are specified in equations (13) and (14), while those of L L (t) are given in equations (19) and (20). Since the imaginary parts are equal, we focus on the real parts.
¢ L L (t) only features a single decay rate and therefore cannot reproduce the contribution L 2 (t) in L′(t) which incorporates the Matsubara frequencies. Furthermore, ¢ L L (t) cannot reproduce the sine component in ( ) L t 1 . Accordingly, we need to be able to neglect L 2 (t) and the sine contribution in L 1 (t) to identify ¢ ( ) L t L and L′(t). We start by considering L 2 (t). The Matsubara frequencies ν n determine the time scale on which L 2 (t) decays, the smallest decay rate being ν 1 . Accordingly, if the decay rate κ is much smaller than the smallest Matsubara frequency n 1 , L 2 (t) drops to zero much faster than L 1 (t) [32,34]. This is the regime where 1, A1 1 and we recoverequation (23) of the main text. In this regime, one expects that L 2 (t) contributes to L(t) only on very short time scales and is negligible if we are interested in not too short time scales [3]. This is the case in our considerations. If , we can neglect L 2 (t) completely. Assuming that we can neglect L 2 (t) on the time scales of interest, the decay of correlations is given by L 1 (t) which only features a single decay rate. Now, we need to find the regime where In the limit  b k  1, we can expand the sine and cosine terms in L 1 (t) in this small parameter. To first order we obtain if the reservoirs are at the same inverse temperature β. Accordingly, we assume that the reservoir in the Lindblad description and the Ohmic oscillator bath have the same inverse temperature β.
Let us now assume that the values of κ and ω m are fixed. We should note now that the condition inequation (A1), which characterizes the regime where L 2 (t) is negligible, favors higher temperatures. However, in order to suppress the sine component in L 1 (t), lower temperatures are more favorable. Accordingly, we estimate that the approximation is best in some intermediate temperature regime.
In order to illustrate the above statement, we compute the distance where we used the abbreviations   The nearest-neighbor geometry as well as the coefficients ω n and t n are directly related to the recurrence coefficients of the three-term recurrence relation defining the orthogonal polynomials p n (ω) [6]. This transformation from the spin-boson model to a one-dimensional geometry is depicted in figure B1.
In the second step, this emerging configuration is treated by the time evolving block decimation (TEBD) method. TEBD generates a high fidelity approximation of the time evolution of a one-dimensional system subject to a nearest-neighbor Hamiltonian with polynomially scaling computational resources. TEBD does so by dynamically restricting the exponentially large Hilbert space to its most relevant subspace thus rendering the computation feasible [55,56]. TEBD is essentially a combination of an MPS description for a one-dimensional quantum system and an algorithm that applies two-site gates that are necessary to implement a Suzuki-Trotter time evolution. Together with MPS operations such as the application of measurements this yields a powerful simulation framework. An extension to mixed states is possible by introducing a matrix product operator to describe the density matrix, in complete analogy to an MPS describing a state [55]. Such an extension is needed in our simulations in order to build the thermal state of the oscillator chain.
A last step is necessary to adjust this configuration further to suit numerical needs. The number of levels for the environment oscillators is restricted to a value d max to reduce the required computational resources. A suitable value for d max is related to the sites average occupation which, in turn, depends on the environment structure and temperature. In our simulations, we set d max =5: this value provides converged results for all examples provided. The Hilbert space dynamical reduction performed by TEBD is determined by the bond dimension. The optimal choice of this parameter depends on the amount of long range correlations in the system. For all the simulations used in this work, a bond dimension χ=200 provided converged results. At last, we observe that the mapping described above produces a semi-infinite chain that must be truncated in order to enable simulations. In order to avoid unphysical back-action on the system due to finite-size effects, i.e. reflections from the end of the chain, the chain has to be sufficiently long to completely give the appearance of a 'large' reservoir. These truncations can be rigorously certified by analytical bounds [53]. For the examples provided in the main text, chains of n=15 sites are enough to avoid boundary effects. In order to further optimize our simulations, we augmented our TEDOPA code with a reduced-rank randomized singular value decomposition (RRSVD) routine [57,58]. Singular value decomposition (SVD) is at the heart of the dimensionality reduction TEBD relies on. RRSVD is a randomized version of the SVD that provides an improved-scaling SVD, with the same accuracy as the standard state-of-the-art deterministic SVD routines.
In order to benchmark the quality of the effective model presented in the main text, we compared the dynamics of the full spin-boson model inequation (B1) with spectral density as inequation (7) of the main text with those of a spin coupled to a damped harmonic oscillator in Lindblad description inequation (32). The simulation parameters are found in table 1 and the results are presented in figure 3. We set the hard cutoff of the macroscopic environment inequation (B3) to ω max /2π=200 kHz. For both cases, one can appreciate very good agreement between the two dynamics.
Note that the simulation of one curve for the case Δ/2π=−50 kHz took 15 days with 16 cores on the bwForCluster JUSTUS such that simulations for the case Δ/2π=3 kHz as presented in figure 4 of the main text are out of reach. . In this case, the lasers are far offresonant for all transitions and the ground states are only weakly coupled to the excited state. We can then adiabatically eliminate the excited state from the dynamics and obtain an effective dynamics in the ground state manifold. Applying the formalism of [59] to our system, we obtain the effective Lindblad equation The effective Hamiltonian H eff has three contributions H eff =H g +H sr +H odf . The first part contains the shifted ground state levels where the Δò s are the ac-Stark shifts of the spin levels due to the applied laser beams The second part, H sr , describes two-photon stimulated Raman transitions between the spin states where a photon is absorbed from one laser beam followed by stimulated emission into the other beam Thus, we obtain three effects on the spin states. The first is an ac-Stark shift of the spin levels due to the laser fields. The differential ac-Stark shift between spin levels can usually be canceled in experiments by adjusting polarization and intensity of the lasers [60]. Hence, we ignore this contribution. Alternatively, it could be absorbed into ω 0 . Figure C1. The figure shows a three level Λ-system consisting of the ground states ñ | and ñ | which are separated in frequency by ω 0 and both feature a dipole-allowed transition to the decaying excited state ñ |e . The transitions are driven by two lasers with frequencies ω 1/2 and Ω l,s denotes the Rabi frequency of laser l on transition ñ  ñ | | s e. Δ R is roughly the detuning of the lasers from the excited state. Depending on the effective laser frequency ω L =ω 1 −ω 2 , different operations on the ground states can be implemented (see text). Spontaneous emission from the excited to the ground states happens at rates Γ s and is indicated by the curly lines.
If one chooses the frequency difference between lasers close to the transition frequency between the spin states ω 1 −ω 2 ≈ω 0 , the second part of the Hamiltonian is resonant and one can drive coherent two-photon stimulated Raman transitions between the spin states. In this case, we usually have Ω odf , Ω rw = ω 0 , the third contribution H odf is highly off-resonant and can be neglected in a rotating wave approximation.
Finally, there is the regime of the spin-dependent optical dipole forces where the beatnote between the two lasers matches one of the motional frequencies ω 1 −ω 2 ≈ω k . Usually w w  k 0 such that now the stimulated Raman processes in H sr are highly off-resonant and can be neglected in a rotating wave approximation. with the effective laser frequency ω L =ω 1 −ω 2 and phase f L =f 1 −f 2 . Furthermore, we have written the phases e k r i l explicitly again and introduced the effective laser wave vector k L =k 1 −k 2 . Note that we have omitted the first part of H odf inequation (C13). For our choice of laser frequency this term would couple to the motion but it can usually be canceled choosing the appropriate laser intensities, polarizations and detunings [60].
Let us turn to the dissipative part. The effective Lindblad operators are found to read: By keeping only the dominant contributions, i.e. those parts of the action of the Lindblad operators that are time-independent, and using d D 

Appendix D. Spin-boson Hamiltonian with trapped ions
In this section, we show how to obtain the spin-boson Hamiltonian inequation (37) of the main text in an ion trap experiment. In the main text we consider a 24 Mg + -25 Mg + crystal. 25 Mg + has electronic hyperfine ground states with total angular momentum F=2, 3 for the valence electron in the 2 S 1/2 state whose degeneracy can be lifted by a magnetic field.  0 . Their motion is then conveniently described in terms of normal modes [37,38]. For a crystal of N ions, we obtain N modes in each direction such that, taking into account the coupled harmonic motion, the system's Hamiltonian becomes . For the 'Ohmic' case (Δ/2π=3 kHz) we chose T=0.01/Δ and N=10 4 and for the resonant case (Δ/2π=100 kHz) T=0.1/Δ and N=10 4 . Note that taking a too small time step eventually leads to discontinuous behavior in N RHP .