Spin-boson model as a simulator of non-Markovian multiphoton Jaynes-Cummings models

The paradigmatic spin-boson model considers a spin degree of freedom interacting with an environment typically constituted by a continuum of bosonic modes. This ubiquitous model is of relevance in a number of physical systems where, in general, one has neither control over the bosonic modes, nor the ability to tune distinct interaction mechanisms. Despite this apparent lack of control, we present a suitable transformation that approximately maps the spin-boson dynamics into that of a tunable multiphoton Jaynes-Cummings model undergoing dissipation. Interestingly, the latter model describes the coherent interaction between a spin and a single bosonic mode via the simultaneous exchange of n bosons per spin excitation. Resorting to the so-called reaction coordinate method, we identify a relevant collective bosonic mode in the environment which is then used to generate multiphoton interactions following the proposed theoretical framework. Moreover, we show that spin-boson models featuring structured environments can lead into non-Markovian multiphoton Jaynes-Cummings dynamics. We discuss the validity of the proposed method depending on the parameters and analyze its performance, which is supported by numerical simulations. In this manner, the spin-boson model serves as a good analog quantum simulator for the inspection and realization of multiphoton Jaynes-Cummings models as well as the interplay of non-Markovian effects, and thus, as a simulator of light-matter systems with tunable interaction mechanisms.


I. INTRODUCTION
The rapid technological progress we have experienced during the last decades has made possible previously inconceivable experiments at the quantum regime, boosting their degree of precision, isolation and control to unprecedented limits [1]. Nowadays, quantum systems can be inspected in a very controllable manner in a number of distinct setups. This experimental breakthrough has therefore stimulated the emergence of research areas such as quantum information and computation and quantum simulation, where the exploitation of quantum effects will allow us to surpass both the capabilities of their classical counterparts in the near future [2]. In particular, quantum simulation considers a scenario in which a wellcontrolled quantum system serves as a simulator of other inaccessible systems [3][4][5]. In this manner, interesting quantum dynamics (i.e. the target dynamics) may be explored using, for example, optical lattices [6] or trapped ions [7]. The target dynamics can be obtained either by decomposing the timeevolution propagator in a set of simple quantum operations (digital quantum simulation) or by finding a map that brings the Hamiltonian into the desired form of the model to be simulated (analog quantum simulation) [5]. In this article we will consider the latter method, by using as a quantum simulator the paradigmatic spin-boson model [8,9].
The spin-boson model describes a spin immersed in an environment formed by a large, typically infinite, number of bosonic modes, in contrast to the quantum Rabi or Jaynes-Cummings models where the interaction comprises a single bosonic mode [10][11][12][13]. The spin-boson model encompasses very rich physics depending on how the spin couples with the * r.puebla@qub.ac.uk distinct bosonic modes. Hence, while it is a minimal model to scrutinize quantum effects of dissipation, it has application in a broad range of systems [8,9], ranging from defects in solid state platforms to quantum emitters in biological systems [14]. Moreover, much efforts inspecting the spin-boson model have dealt with its critical behavior, that is, with the emergence of a quantum phase transition between a delocalized and a localized phase of the spin degree of freedom as one increases the spin-environment coupling [8]. The simulation of the spin-boson model in the strong coupling regime is however computationally very demanding, as acknowledged in [15][16][17][18], since the spin and the bosonic modes become entangled forming a truly quantum many-body system. In some situations one can still resort to analytical methods which may simplify the problem considerably. Among these methods one finds the so-called reaction coordinate mapping [19][20][21][22][23][24][25], which can be viewed as a first step of the more general semiinfinite chain mapping of the environmental degrees of freedom [26,27]. The reaction coordinate is defined as a collective mode of the original environment oscillators. In this manner, one can bring the spin-boson model into the form of a generalized quantum Rabi model [10,11,13] whose bosonic mode undergoes dissipation as it interacts with the residual environment. In particular cases, upon rearranging the original environmental degrees of freedom, the dissipation acquires a Markovian character, and hence simplifying considerably the complexity of the problem (see for example [21]). It is also worth mentioning other attempts to effectively capture quantum dynamics with complex system-environment interactions, as for example the recent work relying on pseudo-modes [28] which builds on the proved equivalence for the dynamics of the system in both frames [29]. The quantum Rabi model (QRM), as well as its simplified version known as Jaynes-Cummings model (JCM) [12], arXiv:1904.07037v1 [quant-ph] 15 Apr 2019 play a central role in the description of light-matter interacting systems and in quantum information science [2,13]. In these models, the interaction mechanism between the spin and bosonic degrees of freedom has a linear form, namely, the spin gets excited or deexcited by absorbing or emitting one bosonic excitation. While this interaction is ubiquitous in quantum physics and with application in various experimental platforms [30], other forms of a spin-boson exchange mechanism beyond this simple case are also of interest. On the one hand, interactions beyond the linear fashion are of relevance for several applications in quantum computation and simulation (e.g. Kerr effect [31]). Furthermore, these exchange mechanisms may unveil interesting phenomena in light-matter systems [32,33] as well as in their multiple spin counterparts [34]. One possible generalization of the QRM or JCM consists in considering a spin-multiphoton interaction, where the spin exchange n excitations simultaneously with the bosonic mode. Such a generalization is often regarded as n-photon QRM or JCM, (nQRM or nJCM), and it has recently attracted attention mainly in its n = 2 form [32,33,[35][36][37][38], although models with n > 2 have been also analyzed [39]. From an experimental point of view however, such multiphoton terms are typically hard to attain. Thus, its realization may benefit from quantum simulation protocols, allowing for enough tunability and control over multiphoton interactions terms, as proposed using optical trapped ions [32,35] or superconducting qubits [37]. These latter schemes realize effective multiphoton exchange terms by exploiting the nonlinear fashion in which the spin and bosonic degrees of freedom couple. It is however still possible to realize such multiphoton models even when the setup comprises solely a linear, i.e. standard, interaction mechanism and thus, it is not suited for a direct simulation of these models, as shown in [40].
In this article we follow the theoretical framework developed in [40,41], combining the ideas of the reactioncoordinate mapping [19][20][21][22][23][24][25] to show that the paradigmatic spin-boson model, featuring a continuum of bosonic modes, can serve as an analog quantum simulator for the realization of dissipative multiphoton Jaynes-Cummings models. In particular, by considering a full spin-boson model we naturally extend the theoretical framework beyond the standard local master equation description of dissipation effects in the simulator, as considered in [41]. Furthermore, we show that the simulated multiphoton Jaynes-Cummings models may acquire non-Markovian behaviors when the spin-boson model features a structured environment, and thus, highlighting the suitability of the proposed theoretical framework to explore aspects of non-Markovianity in distinct light-matter interacting systems.
The article is organized as follows. In Sec. II we introduce the spin-boson model, while in Sec. III we explain how to map the spin-boson model into a different Hamiltonian comprising the desired spin-multiphoton interaction terms, and discuss how the dissipative effects must be transformed in to the aimed model. For that, we first introduce the reaction coordinate mapping in Sec. III A, while in Sec. III B we explain how to extend the theoretical framework to incorporate further bosonic modes in the realization of the desired multipho- σ + a n +H.c. σ + a †n +H.c.
(a) Spin-boson model in the customary star configuration, where each of the circles corresponds to a harmonic oscillator of the environment with frequency ω k interacting with the spin through σxf k (c k + c † k ), before the reaction coordinate mapping. In (b) we show an underdamped spin-boson spectral density JSB(ω), peaked at ω0 (cf. Eq. (8)). Upon the reaction coordinate mapping, a collective degree of freedom is included into the system, which in turn interacts with the residual environment, as sketched in (c) (see main text for further details). For an underdamped JSB(ω), JRC(ω) adopts a Markovian form, as depicted in (b). Such interaction with a collective coordinate can be exploited to realize Hamiltonians containing multiphoton interaction terms, as indicated in (c) and explained in detail in Sec. III. For structured environments one can still rearrange the original environment using more collective coordinates into the augmented system S , where each of them interacts now with its own residual environment, as sketched in (d) (see III B for further details). ton model. After having provided the theoretical derivation of how to perform the analog quantum simulation, we present examples and numerical results for the simulation of different multiphoton Jaynes-Cummings models in Sec. IV. Finally, we summarize the main conclusions of this article in Sec. V.

II. THE SPIN-BOSON MODEL
The spin-boson model describes a two-level system interacting with a large, typically infinite, number of bosonic modes, which constitute the environment. This model has been acknowledged as a paradigm for the inspection of quantum dissipation and quantum-to-classical transition [8,9]. As many physical systems can be well approximated as a twolevel system for sufficient low temperature, the spin-boson model has become a cornerstone in the description of quantum effects in diverse physical realizations, ranging from quantum-based setups [8,9] to biological complexes [14]. In addition, this model has played a key role in the development of the theory of open quantum systems [42], providing a suit-able test-bed to benchmark distinct approximations and tools aimed to efficiently deal with the large number of environment degrees of freedom. Moreover, the relevance of the spinboson model also encompasses the context of critical systems, as it features a quantum phase transition between spin localized and delocalized phases (see Refs. [43,44] and the references therein). Hence, the spin-boson model exhibits rich physics and it is of fundamental relevance in many different areas of research.
The Hamiltonian of the spin-boson model can be written as where each contribution reads as The first two terms represent the free-energy Hamiltonians of the spin and environment, while the last describes the interaction between them. Here we consider that the frequency splitting of the spin is given by ∆ 0 , while 0 accounts for the bias between the eigenstates of the two-level system |± , and with σ = (σ x , σ y , σ z ) the usual spin-1 2 Pauli matrices (see Fig. 1(a)). Hence, σ x |± = ± |± , σ z |e = |e and σ z |g = − |g . The interaction with the environment is dictated by H S−E , where the kth mode with energy ω k is coupled to the spin with a strength f k . These bosonic modes fulfill the usual commutation relation [c k , c † k ] = δ k,k . Remarkably, the system-environment interaction can be completely characterized in terms of the spectral density, J SB (ω) = k f 2 k δ(ω − ω k ). In addition, we comment that one could consider the application of n d drivings onto the spin. As discussed in [40,41], under certain conditions that we will explain in the following section, applying spin drivings enables the realization of different multiphoton Jaynes-Cummings interaction terms. In this manner, while a multiphoton Jaynes-Cummings model can be attained without the need of any additional driving, n d = 0, the realization of a multiphoton quantum Rabi model requires the application of at least one driving n d = 1. In general, the free-energy Hamiltonian of the spin under n d drivings with amplitude j and detuning ∆ j with respect to the spin frequency splitting ∆ 0 reads as Clearly, setting j>0 = 0 (or ∆ j = ∆ 0 ) we recover the form of the standard drivingless H S given in Eq. (2). For the sake of simplicity, in this article we will focus in cases with n d = 0, i.e., aiming to realize multiphoton Jaynes-Cummings models. However, we stress that the procedure explained in the following can be applied in a straightforward manner when n d > 0.

III. ANALOG SIMULATION OF MULTIPHOTON SPIN-BOSON INTERACTIONS
The task now consists in bringing the spin-boson Hamiltonian H SB into the form of a n-photon model, i.e, into a model containing interaction terms of the form σ ± a n and σ ± (a † ) n . For that, one could perform the approximate mapping used in [40,41] directly onto H SB . This would require the selection of a particular bosonic mode out of the environment with frequency ω q to now play the role of a in the interaction with the spin (c q → a), while treating the rest of c k =q as a residual environment. Here, however, we resort to a more sophisticated procedure, based on the so-called reaction coordinate (RC) mapping [19][20][21][22][23][24][25], which consists in rearranging the environment degrees of freedom, such that a small number of collective coordinates can be included in the Hamiltonian part, which in turn interact with the residual environment. In certain cases, the open-quantum system description of the augmented system is considerably simplified with respect to the original system plus environment. Clearly, if the spin-boson model involves just a discrete number of modes, the reaction-coordinate procedure then trivially retrieves the original discrete environment.

A. Reaction coordinate mapping
In the following we summarize how to make use of the RC mapping for a spin-boson model, which has been studied previously in different works [21,22], while referring to the Appendix A and Refs. [19][20][21][22][23][24][25] for further details of the calculations and of the RC mapping.
We shall start by defining a collective mode or reaction coordinate, described by the annhilation and creation operators a and a † , such that while the residual environmental degrees of freedom transform into b k and b † k , requiring that the latter appear in a normal form in the Hamiltonian. In this manner, the original spin-boson Hamiltonian adopts the form of where the former is given by and the other two terms are The reaction coordinate map is completed upon the identification of the parameters λ, Ω and g k or thus J RC (ω) = k g 2 k δ(ω − ω k ). For certain cases, such mapping allows for an exact relation between the original and transformed parameters [25].
Indeed, considering an underdamped spin-boson spectral density in the initial spin-boson model, one can show that the resulting spectral density for the residual environment interacting with the reaction coordinate reads as provided Λ/ω 1, and where the parameters are related according to γ = Γ/(2πω 0 ), Ω = ω 0 and λ = παω 0 /2 (see Appendix A or [19][20][21]25] for further details of this derivation). Here, the frequency ω 0 in J SB (ω) denotes the position at which the spectral density features a maximum, while Γ and α give account for its width and strength, respectively. For J RC (ω), the coupling strength is given by γ. In this manner, by augmenting the system incorporating a collective mode, the original spin-boson model with J SB (ω) is transformed into a spin plus reaction coordinate, which now in turn interacts with a Markovian environment, where the standard Born-Markov approximations can be performed [42]. Indeed, the master equation governing the dynamics of the augmented system, spin plus reaction coordinate, reads as (see Appendix A for the details of the calculation which closely follows [21]) with x = a + a † , while the quantities χ and Θ define the rates affecting the reaction coordinate. They are define as where Having obtained the reaction coordinate Hamiltonian, we undertake the transformation of H S+RC , and thus, of the Eq. (10), to achieve a model that comprises spin-multiphoton interaction terms. For that purpose, we will introduce two auxiliary Hamiltonians H a and H b which will arise in the intermediate steps by moving into suitable interaction picture and transforming them accordingly. The first step consists indeed in moving to a rotating frame in which while Eq. (10) transforms intȯ is the time-evolution operator of a Hamiltonian H x . Then, we perform a further transformation using the unitary operator T (α), defined as T where the Hamiltonian H b can be written as Hence, the dissipator acting on ρ b has the same form as in Eq. (14) but with transformed operators, namely T † xT , T †χ T and T †Θ T , where T ≡ T (−λ/Ω). Finally, by moving to an interaction picture with respect to H b,0 = (ν−ν)a † a−ωσ z /2, and expanding the exponential in Eq. (16) (the latter requires that |2λ/Ω| (a + a † ) 2 1 for truncating the exponential to a finite number of terms) we arrive to a Hamiltonian containing multiphoton interaction terms. The latter condition is commonly known as Lamb-Dicke regime. In addition, we consider the driving frequencies to be ∆ j = ±n j (ν − Ω) −ω with |Ω−ν| j /2 so that one can safely perform a rotatingwave approximation keeping only those terms that are resonant, i.e., time independent (see Appendix B for further details of the calculation). Note that, as H b is similar to the Hamiltonian describing an optical trapped ion under the action of lasers driving vibrational sidebands [45], the procedure to obtain Jaynes-Cummings or quantum Rabi models is analogous to those cases [32,46,47]. In this manner, we can approx- where H n contains the aimed multiphoton interactions, Note that the sets r and b encompass the terms with amplitude j driving red-and blue-sidebands, that is, those terms in Eq. (5) with frequency ∆ j∈r = +n j (ν − Ω) −ω and ∆ j∈b = −n j (ν − Ω) −ω. Each of these drivings will contribute with a multiphoton interaction, either σ + a nj +H.c. for j ∈ r or σ − a nj + H.c. for j ∈ b, which produce transitions between the states |m |g ↔ |m ∓ n j |e . In order to show how the dissipative part transforms, it is advisable to introduce the time-dependent unitary operator Then, one can see that, definingχ = ΦχΦ † ,Θ = ΦΘΦ † and x = Φ(a + a † )Φ † , the resulting master equation for ρ n (t) iṡ where the state ρ n (t) of the multiphoton model is related to the original spin-boson upon the reaction coordinate mapping, ρ S+RC (t), through a unitary transformation From the previous expression it follows that the purity of the total state ρ S+RC and that of ρ n are approximately equal. Moreover, the reduced spin state in the different frameworks are related according to and Tr RC [·] denote the trace over the environment degrees of freedom, and reaction coordinate, respectively. In this manner, having access to the spin degree of freedom one can have access to the dissipative spin dynamics dictated by the master equation (19) under a multiphoton Hamiltonian H n , given in Eq. (17), whose parameters can be tuned. In addition, we remark that the initial state at t 0 = 0 in the multiphoton frame is related to that of the spinboson model as ρ n (0) = T † ρ S+RC (0)T . At this stage, a few comments regarding the validity of Eq. (20) are in order. While the steps performed from H S+RC to H b are exact, H n is attained in an approximate manner. The good functioning of the simulation depends on how these approximations are met. That is, Eq. (20) holds within the Lamb-Dicke regime |2λ/Ω| (a + a † ) 2 1 and for parameters satisfying |Ω −ν| j /2 ∀j so that one can perform a rotating-wave approximation. Note that, as the parameters λ and Ω are directly related to the original spin-boson spectral density, these conditions set constraints onto the accessible parameters, as well as on the temperature of the environment. Furthermore, in order to observe coherent multiphoton dynamics the noise rates in Eq. (19) must be small compared to parameters involved in H n . For the considered shape of J SB (ω) this translates into Γ ν,g n whereg n = 0 (2λ) n /(2Ω n n!) for a n d = 0 and ∆ 0 = ±n(ν − Ω) −ω (cf. Eq. (17). Finally, we comment that previous scheme can be carried out beyond the Lamb-Dicke regime [41]. Admittedly, when the Lamb-Dicke approximation does not hold, the Hamiltonian H n is not longer a good approximation to the dynamics. In this case the Hamiltonian H n must be replaced by a suitable nonlinear Jaynes-Cummings or quantum Rabi model, whose coupling constants crucially depend on the Fock-state occupation number in a nonlinear fashion [48][49][50][51]. These nonlinear, yet multiphoton Hamiltonians appear then as a good approximation to H b , and thus to H SB whenever |2λ/Ω| (a + a † ) 2 1 is not fulfilled, as recently shown in [41]. In this article however we will constrain ourselves to parameters within the Lamb-Dicke regime.

B. Structured environments
As aforementioned, the simulation of multiphoton spinboson interactions is not restricted to a determined form of J SB (ω). Here we show the derivation of the procedure to obtain an effective multiphoton Hamiltonian when the initial spin-boson model features a more complicated interaction with the environment. For simplicity, we consider that J SB (ω) can be split in two parts, J SB (ω) = J SB,1 (ω) + J SB,2 (ω), although its generalization to more is straightforward. The first contribution, J SB,1 (ω), is considered here to be suitable for the realization of multiphoton interactions as described in III A. In addition, we will work under the assumption that the environment degrees of freedom corresponding to J SB,2 (ω) can be treated and simplified using again a collective or reaction coordinate, as sketched in Fig. 1(c).
As discussed previously, we identify a collective coordinate for each of the contributions to the spectral density J SB (ω). In this manner, we augment the system to include both reaction coordinates, denoted here by S = S + RC 1 + RC 2 . Hence, its Hamiltonian is given by where H S,d is the original spin Hamiltonian which may contain spin rotations, introduced in Eq. (5), while the subscripts denote the corresponding reaction coordinate. The parameters λ i and Ω i are determined by the spectral density J SB,i (ω).
The dynamics of the augmented system S is governed by the following master equation: where x i = a i + a † i for i = 1, 2, and χ i and Θ i are defined in analogy to Eqs. (11) and (12).
In order to find a suitable transformation to realize multiphoton interaction terms from H S we proceed in a similar manner as for a single reaction coordinate. That is, we first move to a rotating frame where H S ≡ H I a,1 , with H a = H a,0 + H a,1 and H a,0 = −∆ 0 /2σ x . Therefore, the transformed Hamiltonian reads as The next step is to perform the transformation using the unitary operator T (α). As aforementioned, we consider that the first reaction coordinate is suitable for the quantum simulation of multiphoton interaction terms, due to the form of its spectral density. This argument enables to choose α ≡ −λ 1 /Ω 1 , hence H b ≡ T † (−λ 1 /Ω 1 )H a T (−λ 1 /Ω 1 ). This transformation acts trivially on the second reaction coordinate, but it does affect the coupling between the latter and the spin. Finally, if we move to an interaction picture with respect to H b,0 = (Ω 1 −ν 1 )a † 1 a 1 −ωσ z /2 we obtain the Hamiltonian where we have considered ∆ j = ±n j (ν − Ω 1 ) −ω, and assumed the Lamb-Dicke regime |λ 1 /Ω 1 | (a 1 + a † 1 ) 2 1, and |Ω 1 −ν| j /2 to perform a rotating-wave approximation. Note that, while the multiphoton terms are identical to those of H n in Eq. (17), the second reaction coordinate interacts with the spin degree of freedom. Indeed, depending on the parameters of H n,2 , the effect of such interaction may effectively into non-Markovian effects for the reduced state of the spin and first reaction coordinate, ρ n = Tr 2 [ρ n,2 ]. The final master equation governing the dynamics of ρ n,2 iṡ where the operators involved are defined as in the case involving a single reaction coordinate [cf. Eq. (19)]. It is worth stressing that the relation between the states given in Eq. (20) still holds. From the previous derivation one can observe that the extension to more collective coordinates is straightforward.

IV. EXAMPLES AND NUMERICAL SIMULATIONS
In this section we provide examples of the previously explained general theoretical framework to investigate the performance of the quantum simulation of different multiphoton Hamiltonians H n , as well as to discuss the limitation in the parameter regime for their realization. In particular, in Sec. IV A we first consider the case in which the original spinboson model interacts just with a discrete number of modes, which can be viewed as a limit of vanishing spectral broadening Γ → 0. This scenario will allow us to examine the validity of the required approximations without the effect of dissipation. Then, in Sec. IV B, we will consider Γ = 0, where the reaction-coordinate mapping appears as a key step to realize a desired multiphoton Jaynes-Cummings model.
In all the cases, we assess the performance of the realization of the targeted multiphoton Jaynes-Cummings models by means of the fidelity F (t) between two states, In particular, we will analyze to what extent is the relation given in Eq. (20) satisfied. In other words, we will compare the aimed state of a multiphoton Jaynes-Cummings model ρ n (t) with the one retrieved using the analog simulator, Φρ S+RC (t)Φ † . We remark that when two reaction coordinates are included, the state ρ n (t) obeys the master equation given in Eq. (25), whose Hamiltonian is H n,2 , Eq. (24), while ρ S+RC (t) must be replaced by ρ S , as explained in III B. In addition, we will show that the theoretical framework allows us to realize non-Markovian multiphoton Jaynes-Cummings models. Among the different measures for non-Markovianity [52], we resort to the one based on the trace distance [53], defined as where |A| = √ A † A. Then, non-Markovian evolutions can be characterized as those for which D(ρ 1 (t), ρ 2 (t)) increases during certain time intervals, that is, for those for which the time-derivative of the trace distance for a pair of states ρ 1,2 , is σ(t, ρ 1,2 ) > 0. In general, one has to maximize over all possible pair of states ρ 1,2 in order to find a suitable non-Markovian measure [53]. For our purpose however, it will be sufficient with showing that σ(t, ρ 1,2 ) > 0 for a certain pair of states. In this manner, we offer a proof-of-principle that non-Markovian multiphoton models can be realized.

A. Dissipationless multiphoton Jaynes-Cummings models
We start considering the simplest case, namely, when the spin-boson model simply involves the interaction with a discrete number of modes. This corresponds to either consider Γ → 0 in the underdamped spectral density J SB (ω), or equivalently, assuming that dissipation effects are sufficiently small so that they can be discarded. Note that for a single bosonic mode with Γ = 0, the spin-boson model adopts the form of a generalized quantum Rabi model, which is indeed H S+RC as given in Eq. (7). Recall that in this particular case H SB ≡ H S+RC as there are no further modes in the system. In particular, we set n d = 0 in Eq. (5) as we aim to realize a single multiphoton Jaynes-Cummings interaction. The Hamiltonian for a nJCM can be written in general as H nJCM =ω 2 σ z +νa † a +g n σ + a n + σ − (a † ) n . (29) At resonant condition,ω = nν, the coupling constantg n fixes the time required to transfer population from the state |e |0 to |g |n , denoted as τ n = π/(2g n √ n!). Both are related to the spin-boson parameters as (cf. Eq. (17)) Clearly, as 2λ/Ω must be small to lie within the Lamb-Dicke regime, the couplingg n decreases considerably for increasing n, requiring longer evolution times under the spin-boson Hamiltonian to observe a significant effect, that is, an evolution time of the order of τ n . In Fig. 2 we show the results for the realization of 2JCM and 3JCM models using a spin-boson model interacting with a single bosonic mode. In order to observe the paradigmatic Rabi oscillations between the states |e |0 and |g |n , we choose ρ S+RC (0) = |− −| ⊗ ρ th RC as an initial state for the spin-boson model, where ρ th RC is a thermal state at temperature β −1 for the reaction coordinate mode, containing n th = (e βΩ − 1) −1 bosons. Recall that, as we consider here a single spectral density with Γ = 0, the reaction coordinate mode is simply the unique mode which interacts with the spin degree of freedom. In this manner, the initial state for the simulated multiphoton models reads as ρ nJCM (0) = T † ρ S+RC (0)T , which approximately amounts to ρ nJCM (0) ≈ |e e| ⊗ |0 0| for sufficiently low temperature and small 2λ/Ω. The chosen parameters for the simulation of the 2JCM, plotted in Fig. 2(a) and (b) are πα = 0 = 0.02ω 0 , recalling that Ω = ω 0 it results in 2λ/Ω = 0.2. Choosingν = 10 −3 Ω andω = 2ν, the coupling in 2JCM amounts tog 2 = 0.2ν. The initial reactioncoordinate thermal state, ρ th RC , contains n th = 10 −3 bosons. In Fig. 2(b) we show how the quantum simulation of the 2JCM model deteriorates for increasing number of bosons, as a large n th will eventually break down the Lamb-Dicke regime.
For the 3JCM we choose again πα = 0.02ω 0 , which leads in 2λ/Ω = 0.2. Then, we select the aimed coupling strength The considered initial state reads as ρ S (0) = |− −|⊗ρ th RC 1 ⊗ρ th RC 2 , with β very large such that ρ th ≈ |0 0|. In (b) we plot the time-derivative of the trace distance σ(t) after tracing out the second bosonic mode, and considering the initial states |e and |g for the spin in H S , while both reaction coordinates find themselves in their vacuum. Clearly, σ(t) > 0 during certain intervals, revealing the non-Markovianity introduced due to the interaction with the second mode. Panel (c) shows the evolution of purity for the state upon tracing the second mode, Tr[ρ 2 S+RC 1 (t)], and for the reduced state of the spin, Tr[ρ 2 S (t)], for the same case shown in (a). In panel (d) we compare the infidelity 1 − F (t) between the ideal state and the simulated one using H S for the three different initial states employed here. We refer to Sec. IV A for further details regarding the parameters and states considered in the simulation. of the multiphoton interaction to beg 3 = 0.1ν withω = 3ν, while we vary 0 /ω 0 . The temperature is set to βΩ ≈ 100 so that ρ th RC ≈ |0 0|. As in the previous case, the dynamics are well retrieved, see Fig. 2(c), where we have set 0 /ω 0 = 2 · 10 −3 . Note however that, as a consequence of the rotating-wave approximation performed to achieve a resonant third order (see Appendix B and cf. Eq. (17)), and due to the longer times required to simulate a 3JCM compared to the 2JCM, the condition |Ω −ν| 0 must be better satisfied. Indeed, for 0 /ω 0 = 10 −2 we already see a clear departure from the targeted dynamics, as indicated by a large infidelity 1 − F (t) 10 −1 , as shown in Fig. 2(d).
In the following we consider a spin interacting with two bosonic modes, again with Γ 1,2 = 0. As explained in III B, we perform the map onto the first bosonic mode to attain a multiphoton interaction. Upon suitable transformations and approximations, the spin-boson model will take the form of a multiphoton Jaynes-Cummings model H nJCM,2 , where the subscript 2 indicates the presence of a second reaction coordi- . In (c) we compare the different behavior as Γ/ν varies for the von Neumann entropy of the reduced spin state, SvN(ρS(t)). The values of Γ/ν are indicated close to each curve. Finally, the state infidelity 1 − F (t) between the targeted ρ2JCM and its approximate simulation, ΦρS+RC(t)Φ † , is plotted in panel (d) for different Γ/ν. See main text for further details on the parameters employed for the simulation. nate in the system. The Hamiltonian H nJCM,2 reads as H nJCM,2 =ω 2 σ z +νa † 1 a 1 + Ω 2 a † 2 a 2 +g n σ + a n 1 + σ − (a † 1 ) n − λ 2 σ z (a 2 + a † 2 ). (32) In this manner, the spin exchanges n quanta with the first bosonic mode as in H nJCM , while the last term effectively shifts the spin frequency depending on the state of the second mode. The reduced state for the spin and first bosonic mode is given then by ρ nJCM (t) = Tr 2 [ρ nJCM,2 (t)]. Indeed, due to the interaction with the second bosonic mode, the multiphoton Jaynes-Cummings model may exhibit non-Markovian features. For that, we consider the spin-boson Hamiltonian H S given in Eq. (21), which then approximately realizes H nJCM,2 . In particular, we select ∆ 0 = −2Ω 1 , so that the simulated model involves two-photon interaction terms, i.e., a 2JCM. The results are plotted in Fig. 3, while the parameters are πα i = 0.02Ω i such that 2λ i /Ω i = 0.2 for i = 1, 2, 0 /Ω 1 = 10 −2 . The coupling strength in H 2JCM,2 is given byg 2 = 0.2ν withν = Ω 2 . As in the single-mode case, Rabi oscillations will be clearly visible selecting ρ S (0) = |− −| ⊗ ρ th RC1 ⊗ ρ th RC2 . After its transformation, this state corresponds approximately to an initial spin state |e in the nJCM frame. In the same manner, in order to analyze the emergence non-Markovian behavior we consider the initial states |g g| and |e e| for the spin in H S . This implies initial spin states |± in the nJCM frame, which for pure dephasing noise it has been shown to be the pair of states maximizing σ(t) [53]. The results plotted in Fig. 3 have been performed considering a sufficiently low temperature such that ρ th RC1,2 ≈ |0 0|. We then compute the trace distance D(ρ 1 , ρ 2 ) using the states upon tracing out the second mode, ρ 2JCM (t) = Tr 2 [ρ 2JCM,2 ]. As shown in Fig. 3(b) the time-derivative of the trace distance, σ(t), becomes positive during certain intervals -a clear indication of the non-Markovian behavior of the simulated multiphoton Jaynes-Cummings model. In addition, we also calculate the non-trivial evolution of the purity for the states ρ S+RC1 (t) and ρ S (t) = Tr RC1 [ρ S+RC1 (t)], which is shown in Fig. 3(c). According to our theoretical framework, their purity is approximately equal to that of ρ 2JCM (t) and the reduced spin state upon tracing both bosonic degree of freedom in the 2JCM, Tr[ρ 2JCM (t)], respectively. Finally, the infidelity 1 − F (t) between the targeted state ρ 2JCM,2 (t) and its reconstructed one Φρ S+RC1+RC2 (t)Φ † in Fig. 3(d).

B. Dissipative multiphoton Jaynes-Cummings models
We now consider a more realistic scenario in which the spin-boson model interacts with an environment whose spectral density has an underdamped shape, i.e. J SB (ω) has the form of Eq. (8) with Γ = 0. In this manner, we extend the theoretical framework beyond the standard local master equation description [41]. As explained in Sec. III A, this situation can be mapped using a reaction coordinate, which now in turn interacts with a Markovian residual environment. The evolution of the state of the augmented system, spin and reaction coordinate, evolves according to the master equation given in (10). Indeed, the effect of spectral broadening, Γ = 0, introduces dissipation into the simulated multiphoton Jaynes-Cummings model, whose state now obeys the master equation (19). We remark that the performance of the simulated dissipative model is not altered when the effect of dissipation is taken into account correctly. Nevertheless, whenever Γ ν, dissipation dominates the dynamics, and the paradigmatic Rabi oscillations will eventually fade away. In Fig. 4 we show the results of numerical simulations aimed to retrieve a 2JCM with different Γ/ν values, and for different quantities. As for Fig 2, we used πα = 0 = 0.02ω 0 , so that 2λ/Ω = 0.2. We chose againν = 10 −3 Ω andω = 2ν, and therefore the coupling in 2JCM amounts tog 2 = 0.2ν, while the temperature is such that ρ th RC contains n th = 10 −3 bosons. The spin is initialized in the |− state, so that ρ S+RC (0) = |− −|⊗ρ th RC . In particular, the value Γ/ν = 2 · 10 −1 considered in Fig. 4(a) already produces a significant departure from the Rabi oscillation between the states |e |0 and |g |2 in the dissipationless 2JCM (cf. Fig. 2(a) for Γ = 0). Note that the results plotted in Fig. 4(a) corresponds to a critically damped 2JCM since Γ =g 2 . As plotted in Fig. 4(b), the effect of the dissipation is clearly visible in the evolution of the purity for both the total state (spin plus bosonic mode) and the reduced spin state, namely, Tr[ρ 2 S+RC (t)] and Tr[ρ 2 S (t)]. As in pre-vious cases, the purity of these states is directly related to those of the simulated model as a consequence of the relation ρ 2JCM (t) ≈ Φρ S+RC (t)Φ † . Furthermore, Rabi oscillations or population revivals appear in the evolution of von Neumann entropy, S vN (ρ) = −ρ log 2 ρ for the reduced spin state. This is plotted in Fig. 4(c) for different Γ/ν values. Finally, we note that the performance of the quantum simulation is independent of the dissipation as demonstrated by the good fidelities attained in these cases [cf. Fig. 4(d)], allowing for the simulation of different parameter regimes in a nJCM.

V. CONCLUSIONS
We have proposed a theoretical scheme to realize multiphoton Jaynes-Cummings models using the paradigmatic spin-boson model, which contains a continuum of bosonic modes, as an analog quantum simulator. While the spinboson model lacks naturally of these multiphoton interaction terms, we make use of a suitable transformation that approximately maps the spin-boson model into a dissipative multiphoton Jaynes-Cummings model, whose parameters can be tuned. In order to bring the spin-boson model, typically interacting with an infinite number of bosonic modes, into the form of the aimed multiphoton model, we first rearrange the environment degrees of freedom using the so-called reactioncoordinate method [19][20][21][22][23][24][25]. This method allows us to include a set of collective bosonic modes into the coherent description of the problem, which then in turn interact with the residual environment. For certain types of interactions between the spin and the environment, characterized by the spectral density, the reaction coordinate mapping emerges as a powerful tool to reduce the complexity of the problem. In particular, for an underdamped spectral density, the reaction coordinate takes a simple form as it interacts with the residual environment in Markovian fashion. The resulting Hamiltonian is then used to generate multiphoton interaction terms, following the theory explained in [40,41], while the dissipation effects must be transformed accordingly. Furthermore, we extend the scheme to spin-boson models with structured environments. In these cases, the original spin-boson Hamiltonian can be mapped onto the one of a spin interacting with more reaction coordinates. In this manner, we show how to extend the theoretical framework to account for these additional modes. In particular, due to presence of two or more reaction coordinates, the attained multiphoton Jaynes-Cummings model can exhibit non-Markovian features. We perform numerical simulations starting from the spin plus reaction-coordinate Hamiltonians and aiming to realize different multiphoton Jaynes-Cummings. We first perform simulations considering one reaction coordinate without dissipation to better illustrate the performance of the required approximations to achieve two-and three-photon Jaynes-Cummings models. We then demonstrate that non-Markovian multiphoton Jaynes-Cummings can be indeed attained when a second reaction coordinate is included, as unveiled by the standard trace distance measure [53]. Finally, we provide numerical simulations investigating the interplay between spectral broadening, dissipation and the decoherence in the targeted multiphoton models. the so-called Leggett prescription, one get One can reproduce the same calculation also after performing the RC mapping and express J RC (ω) in terms of the corresponding kernelL 0 (z). However, since at this stage we are just rearranging the environment in a more convenient way by using a suitable normal mode transformation, the integral kernel must be the same before and after the mapping, hence one can useL 0 (z) instead ofL SB (z) in Eq. (A2). By considering the Ohmic spectral density J RC (ω) = γωe −ω/Λ , one obtains J SB (ω) = 4γΩ 2 λ 2 ω (Ω 2 − ω 2 ) 2 + (2πγΩω) 2 . (A3) It is easy to see that one exactly recovers the underdamped spectral density given by Eq. (8) by simply requiring that γ = Γ/(2πω 0 ), Ω = ω 0 and λ = παω 0 /2. Furthermore, one also needs to solve the dynamics, i.e. writing down the corresponding master equation for the mapped system, system plus reaction coordinate. The guiding idea is to treat exactly the coupling between the spin and the RC, while the interaction between the latter and the residual environment is treated perturbatively up to the second order. This enables to rely on the standard Born-Markov approximation, provided that either the coupling between the augmented system and the residual environment is weak or the residual environment correlation time is short compared to the relevant timescale of the system. Within this approximation, one can work out a master equation that, in the Schrödinger picture, reads aṡ where ρ ≡ ρ S+RC , A = a + a † and the residual environment is assumed to be in a thermal state, i.e. ρ E = e −βH E /Tr E {e −βH E }.
In order to obtain an expression for the interaction picture operators, one can proceed by truncating the space of the augmented system up to n basis states and numerically diagonalising the Hamiltonian H S+RC . To this end, let |φ n be an eigenstate of H S+RC , i.e. H S+RC |φ j = ϕ j |φ j , therefore the operator A can be expanded as A = jk A jk |φ j φ k |, while in the interaction picture one has where A jk = φ j | A |φ k and ξ jk = ϕ j − ϕ k . Finally, by plugging Eq. (A4) into Eq. (A4), and assuming the imaginary parts to be negligible, one gets the final form of the master equation given by Eq. (10).