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

TL;DR

This paper demonstrates how the spin-boson model can be approximately transformed into a tunable multiphoton Jaynes-Cummings model, enabling simulation of complex light-matter interactions and non-Markovian effects using existing quantum systems.

Contribution

The authors introduce a method to map the spin-boson model onto a multiphoton Jaynes-Cummings model, facilitating the simulation of non-Markovian multiphoton interactions in a controlled manner.

Findings

The transformation effectively simulates multiphoton interactions.

Structured environments induce non-Markovian multiphoton dynamics.

Numerical simulations support the validity of the method.

Abstract

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.

Equations

$$H_{\mathrm{SB}}=H_{\mathrm{S}}+H_{\mathrm{E}}+H_{\mathrm{S}-\mathrm{E}}$$

$$H_{\mathrm{S}}$$

$$H_{\mathrm{E}}$$

$$H_{\mathrm{S}-\mathrm{E}}$$

$$H_{\mathrm{S},\mathrm{d}}=\frac{{\Delta }_0}{2}{\sigma }_x$$

$$+\sin ({\Delta }_j-{\Delta }_0)t{\sigma }_y].$$

$$\lambda (a+a^{†})=\sum _k{f}_k(c_k+c_k^{†}),$$

$$H_{\mathrm{S}+\mathrm{RC}}=\frac{{\Delta }_0}{2}{\sigma }_x+\Omega a^{†}a+\lambda {\sigma }_x(a+a^{†})$$

$$+\sum _{j=0}^{n_d}\frac{{ϵ}_j}{2}\left[\cos ({\Delta }_j-{\Delta }_0)t{\sigma }_z+\sin ({\Delta }_j-{\Delta }_0)t{\sigma }_y\right],$$

$$J_{\mathrm{SB}}(\omega )=\frac{\alpha \Gamma {\omega }_0^2\omega }{({\omega }_0^2-{\omega }^2{)}^2+{\Gamma }^2{\omega }^2},$$

$$J_{\mathrm{RC}}(\omega )=\gamma \omega e^{-\omega /\Lambda }$$

$${\dot{\rho }}_{\mathrm{S}+\mathrm{RC}}(t)=$$

$$+\left[x,\left\{\Theta ,{\rho }_{\mathrm{S}+\mathrm{RC}}(t)\right\}\right].$$

$$\chi$$

$$\Theta$$

$$H_a=\Omega a^{†}a$$

$$+\sum _{j=0}^{n_d}\frac{{ϵ}_j}{2}\left[\cos {\Delta }_{j}t{\sigma }_z+\sin {\Delta }_{j}t{\sigma }_y\right].$$

$${\dot{\rho }}_a(t)=-i\left[H_a,{\rho }_a(t)\right]-[x,[\hat{\chi },{\rho }_a(t)]]+[x,\{\hat{\Theta },{\rho }_a(t)\}].$$

$${\dot{\rho }}_b=-i\left[H_b,{\rho }_b\right]$$

$$+\left[T^{†}xT,\left\{T^{†}\hat{\Theta }T,{\rho }_b(t)\right\}\right],$$

$$H_b=\Omega a^{†}a+\sum _{j=0}^{n_d}\frac{{ϵ}_j}{2}\left[{\sigma }^{+}{e}^{2\lambda (a-a^{†})/\Omega }{e}^{-i{\Delta }_{j}t}+H.c.\right].$$

$$H_{\mathrm{n}}=\frac{\tilde{\omega }}{2}{\sigma }_z$$

$$+\sum _{j\in b}\frac{{ϵ}_j(2\lambda {)}^{n_j}}{2{\Omega }^{n_j}{n}_j!}\left[{\sigma }^{+}(-a^{†}{)}^{n_j}+H.c.\right].$$

$$\Phi =U_{b,0}^{†}{T}^{†}{U}_{a,0}.$$

$${\dot{\rho }}_{\mathrm{n}}(t)=-i[H_{\mathrm{n}},{\rho }_{\mathrm{n}}(t)]-\left[\tilde{x},\left[\tilde{\chi },{\rho }_{\mathrm{n}}(t)\right]\right]+[\tilde{x},\{\tilde{\Theta },{\rho }_{\mathrm{n}}(t)\}]$$

$${\rho }_{\mathrm{n}}(t)\approx \Phi {\rho }_{\mathrm{S}+\mathrm{RC}}(t){\Phi }^{†}.$$

$$H_{{\mathrm{S}}^{\prime }}=H_{S,d}$$

$$+{\Omega }_2{a}_2^{†}{a}_2+{\lambda }_2{\sigma }_x(a_2+a_2^{†}),$$

$${\dot{\rho }}_{{\mathrm{S}}^{\prime }}(t)=$$

$$-\left[x_1,\left[{\chi }_1,{\rho }_{{\mathrm{S}}^{\prime }}(t)\right]\right]-\left[x_2,\left[{\chi }_2,{\rho }_{{\mathrm{S}}^{\prime }}(t)\right]\right]$$

$$+\left[x_1,\left\{{\Theta }_1,{\rho }_{{\mathrm{S}}^{\prime }}(t)\right\}\right]+\left[x_2,\left\{{\Theta }_2,{\rho }_{{\mathrm{S}}^{\prime }}(t)\right\}\right],$$

$$H_a=$$

$$+$$

$$H_{\mathrm{n},2}$$

$$+\sum _{j\in r}\frac{{ϵ}_j}{2n_j!}{\left(\frac{2{\lambda }_1}{{\Omega }_1}\right)}^{n_j}\left[{\sigma }^{+}{a}_1^{n_j}+H.c.\right]$$

$$+\sum _{j\in b}\frac{{ϵ}_j}{2n_j!}{\left(\frac{2{\lambda }_1}{{\Omega }_1}\right)}^{n_j}\left[{\sigma }^{+}(-a_1^{†}{)}^{n_j}+H.c.\right],$$

$${\dot{\rho }}_{\mathrm{n},2}(t)$$

$$-\left[{\tilde{x}}_1,\left[{\tilde{\chi }}_1,{\rho }_{\mathrm{n},2}(t)\right]\right]-\left[{\tilde{x}}_2,\left[{\tilde{\chi }}_2,{\rho }_{\mathrm{n},2}(t)\right]\right]$$

Full text

I Introduction

The rapid technological progress we have experienced during the last few decades has made possible previously inconceivable experiments at the quantum regime, boosting their degree of precision, isolation and control to unprecedented limits Dowling and Milburn (2003). Currently, 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 Nielsen and Chuang (2000). In particular, quantum simulation considers a scenario in which a well-controlled quantum system serves as a simulator of other inaccessible systems Feynman (1982); Johnson et al. (2014); Georgescu et al. (2014). In this manner, interesting quantum dynamics (i.e., the target dynamics) may be explored using, for example, optical lattices Bloch et al. (2012) or trapped ions Blatt and Roos (2012). The target dynamics can be obtained either by decomposing the time-evolution 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 (analogous to quantum simulation) Georgescu et al. (2014). In this article, we will consider the latter method, by using as a quantum simulator the paradigmatic spin-boson model Leggett et al. (1987); Weiss (2008).

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 Rabi (1936, 1937); Jaynes and Cummings (1963); Scully and Zubairy (1997). The spin-boson model encompasses very rich physics depending on how the spin couples with the distinct bosonic modes. Hence, while it is a minimal model to scrutinize the quantum effects of dissipation, it has application in a broad range of systems Weiss (2008); Leggett et al. (1987), ranging from defects in solid state platforms to quantum emitters in biological systems Huelga and Plenio (2013). Moreover, much effort inspecting the spin-boson model has dealt with its critical behaviour, 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 Leggett et al. (1987). The simulation of the spin-boson model (or of a generic open quantum system) in the strong coupling regime is however computationally very demanding, as acknowledged in Tanimura and Kubo (1989); Tanimura (1990); Prior et al. (2010); Dattani et al. (2012); Dattani (2013); Wilkins and Dattani (2015); Strathearn et al. (2018), 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 Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014, 2016); Strasberg et al. (2016, 2018); Nazir and Schaller (2018), which can be viewed as a first step of the more general semi-infinite chain mapping of the environmental degrees of freedom Chin et al. (2010); Woods et al. (2014). 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 Rabi (1936, 1937); Scully and Zubairy (1997) 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, hence simplifying considerably the complexity of the problem (see for example Iles-Smith et al. (2014)). It is also worth mentioning other attempts to capture quantum dynamics effectively with complex system-environment interactions, as for example the recent work relying on pseudo-modes Mascherpa et al. , which builds on the proven equivalence for the dynamics of the system in both frames Tamascelli et al. (2018).

The quantum Rabi model (QRM), as well as its simplified version known as the Jaynes-Cummings model (JCM) Jaynes and Cummings (1963) play a central role in the description of light-matter interacting systems and in quantum information science Scully and Zubairy (1997); Nielsen and Chuang (2000). 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 Braak et al. (2016), other forms of a spin-boson exchange mechanisms 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., the Kerr effect Lloyd and Braunstein (1999)). Furthermore, these exchange mechanisms may unveil interesting phenomena in light-matter systems Felicetti et al. (2015); Pedernales et al. (2018), as well as in their multiple spin counterparts Garbe et al. (2017). One possible generalization of the QRM or JCM consists of considering a spin-multiphoton interaction, where the spin exchanges $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 Felicetti et al. (2015); Puebla et al. (2017); Cui et al. (2017); Pedernales et al. (2018); Felicetti et al. (2018); Xie et al. (2019), although models with $n>2$ have been also analysed Lo et al. (1998). 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 interaction terms, as proposed using optical trapped ions Felicetti et al. (2015); Puebla et al. (2017) or superconducting qubits Felicetti et al. (2018). 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 Casanova et al. (2018).

In this article, we follow the theoretical framework developed in Casanova et al. (2018); Puebla et al. (2019), combining the ideas of the reaction-coordinate mapping Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014, 2016); Strasberg et al. (2016, 2018); Nazir and Schaller (2018) to show that the paradigmatic spin-boson model, featuring a continuum of bosonic modes, can serve as an analogue quantum simulator for the realization of different dissipative multiphoton Jaynes-Cummings models by tuning the frequency and bias parameter of the spin. In this manner, we demonstrate the emergence of a connection between the dynamics of these paradigmatic and fundamental quantum models, which was not previously unveiled. Moreover, as the spin-boson model is of considerable experimental significance, i.e., it describes the ubiquitous scenario of a two-level system interacting with an arbitrary environment, our method paves the way for the simulation of multiphoton Jaynes-Cummings models in distinct setups. 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 Puebla et al. (2019). Furthermore, we show that the simulated multiphoton Jaynes-Cummings models may acquire non-Markovian behaviours when the spin-boson model features a structured environment, 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 Section II, we introduce the spin-boson model, while in Section 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 into the aimed model. For that, we first introduce the reaction coordinate mapping in Section III.1, while in Section III.2, we explain how to extend the theoretical framework to incorporate further bosonic modes in the realization of the desired multiphoton model. After having provided the theoretical derivation of how to perform the analogue quantum simulation, we present examples and numerical results for the simulation of different multiphoton Jaynes-Cummings models in Section IV. Finally, we summarize the main conclusions of this article in Section 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 Weiss (2008); Leggett et al. (1987). As many physical systems can be well approximated as a two-level 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 Leggett et al. (1987); Weiss (2008) to biological complexes Huelga and Plenio (2013). In addition, this model has played a key role in the development of the theory of open quantum systems Breuer and Petruccione (2002), providing a suitable test-bed to benchmark distinct approximations and tools aimed to deal with the large number of environment degrees of freedom efficiently. Moreover, the relevance of the spin-boson model also encompasses the context of critical systems, as it features a quantum phase transition between spin localized and delocalized phases (see Refs. Vojta (2006); Hur (2010) 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:

[TABLE]

where each contribution reads as:

[TABLE]

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 $\Delta_{0}$, while $\epsilon_{0}$ accounts for the bias between the eigenstates of the two-level system $\left|\pm\right\rangle$ and with $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ the usual spin-$\frac{1}{2}$ Pauli matrices (see Figure 1(a)). Hence, $\sigma_{x}\left|\pm\right\rangle=\pm\left|\pm\right\rangle$, $\sigma_{z}\left|e\right\rangle=\left|e\right\rangle$ and $\sigma_{z}\left|g\right\rangle=-\left|g\right\rangle$. The interaction with the environment is dictated by $H_{\rm S-E}$, where the $k^{\text{th}}$ mode with energy $\omega_{k}$ is coupled to the spin with a strength $f_{k}$. These bosonic modes fulfil the usual commutation relation $[c_{k},c_{k^{\prime}}^{\dagger}]=\delta_{k,k^{\prime}}$. Remarkably, the system–environment interaction can be completely characterized in terms of the spectral density, $J_{\rm SB}(\omega)=\sum_{k}f_{k}^{2}\delta(\omega-\omega_{k})$, which here is assumed to be known. In anticipation of the developed theoretical framework that allows us to bring $H_{\rm SB}$ into the form of a multiphoton Jaynes–Cummings model, we comment that while the frequency splitting $\Delta_{0}$ tunes the multiphoton order of the interaction, the bias parameter $\epsilon_{0}$ will be proportional to the interacting strength of the simulated model (see Section III).

In addition, we comment that one could consider the application of $n_{d}$ drivings onto the spin. As discussed in Casanova et al. (2018); Puebla et al. (2019), under certain conditions that we will explain in the following section, applying spin drivings enables the simultaneous realization of different multiphoton Jaynes-Cummings interaction terms. In this manner, while a multiphoton Jaynes-Cummings model can be attained without the need for any driving, $n_{d}=0$, the realization of a multiphoton quantum Rabi model requires the application of at least one, i.e., $n_{d}=1$. In general, the free-energy Hamiltonian of the spin under $n_{d}$ drivings with amplitude $\epsilon_{j}$ and detuning $\Delta_{j}$ with respect to the spin frequency splitting $\Delta_{0}$ reads as:

[TABLE]

Clearly, setting $\epsilon_{j>0}=0$ (or $\Delta_{j}=\Delta_{0}$), we recover the form of the standard drivingless $H_{\rm S}$ given in Equation (2). For the sake of simplicity, in this article, we will focus on 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 Analogue simulation of multiphoton spin–boson interactions

The task now consists of bringing the spin-boson Hamiltonian $H_{\rm SB}$ into the form of a $n$-photon model, i.e, into a model containing interaction terms of the form $\sigma^{\pm}a^{n}$ and $\sigma^{\pm}(a^{\dagger})^{n}$. For that, one could perform the approximate mapping used in Casanova et al. (2018); Puebla et al. (2019) directly onto $H_{\rm SB}$. This would require the selection of a particular bosonic mode out of the environment with frequency $\omega_{q}$ to now play the role of $a$ in the interaction with the spin ($c_{q}\rightarrow a$), while treating the rest of $c_{k\neq q}$ as a residual environment. Here, however, we resort to a more sophisticated procedure, based on the so-called reaction coordinate (RC) mapping Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014, 2016); Strasberg et al. (2016, 2018); Nazir and Schaller (2018), which consists of 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.

III.1 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 Iles-Smith et al. (2014, 2016), while referring to Appendix A and References Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014, 2016); Strasberg et al. (2016, 2018); Nazir and Schaller (2018) 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 annihilation and creation operators $a$ and $a^{\dagger}$, such that:

[TABLE]

while the residual environmental degrees of freedom transform into $b_{k}$ and $b^{\dagger}_{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 $H_{\rm SB}=H_{\rm S+RC}+H_{\rm RC-E^{\prime}}+H_{\rm E^{\prime}}$, where the former is given by:

[TABLE]

and the other two terms are $H_{\rm RC-E^{\prime}}+H_{\rm E^{\prime}}=(a+a^{\dagger})\sum_{k}g_{k}(b_{k}+b_{k}^{\dagger})+(a+a^{\dagger})^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}}+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$. The reaction coordinate map is completed upon the identification of the parameters $\lambda$, $\Omega$, and $g_{k}$ or, thus, $J_{\rm RC}(\omega)=\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k})$. For certain cases, such mapping allow for an exact relation between the original and transformed parameters Nazir and Schaller (2018). Indeed, considering an underdamped spin-boson spectral density in the initial spin-boson model,

[TABLE]

one can show that the resulting spectral density for the residual environment interacting with the reaction coordinate reads as:

[TABLE]

provided $\Lambda/\omega\gg 1$ and where the parameters are related according to $\gamma=\Gamma/(2\pi\omega_{0})$, $\Omega=\omega_{0}$, and $\lambda=\sqrt{\pi\alpha\omega_{0}/2}$ (see Appendix A or Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014); Nazir and Schaller (2018) for further details of this derivation). Here, the frequency $\omega_{0}$ in $J_{\rm SB}(\omega)$ denotes the position at which the spectral density features a maximum, while $\Gamma$ and $\alpha$ account for its width and strength, respectively. For $J_{\rm RC}(\omega)$, the coupling strength is given by $\gamma$. In this manner, by augmenting the system incorporating a collective mode, the original spin–boson model with $J_{\rm SB}(\omega)$ 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 Breuer and Petruccione (2002). 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 Iles-Smith et al. (2014)):

[TABLE]

with $x=a+a^{\dagger}$, while the quantities $\chi$ and $\Theta$ define the rates affecting the reaction coordinate. They are defined as:

[TABLE]

where $x_{jk}=\left\langle\phi_{j}\right|x\left|\phi_{k}\right\rangle$, $H_{\rm S+RC}\left|\phi_{j}\right\rangle=\varphi_{j}\left|\phi_{j}\right\rangle$, and $\xi_{jk}=\varphi_{j}-\varphi_{k}$.

Having obtained the reaction coordinate Hamiltonian, we undertake the transformation of $H_{\rm S+RC}$, and thus, of Equation (III.1), 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 a suitable interaction picture and transforming them accordingly. The first step consists indeed of moving to a rotating frame in which $H_{\rm S+RC}\equiv H_{a,1}^{I}$ where $H_{a}=H_{a,0}+H_{a,1}$ with $H_{a,0}=-\Delta_{0}/2\sigma_{x}$. In this manner, we find:

[TABLE]

while Equation (III.1) transforms into:

[TABLE]

where $\hat{\chi}=U_{a,0}\chi U_{a,0}^{\dagger}$ and $\hat{\Theta}=U_{a,0}\Theta U_{a,0}^{\dagger}$, such that $U_{x}=\mathcal{T}e^{-i\int_{0}^{t}dsH_{x}(s)}$ is the time-evolution operator of a Hamiltonian $H_{x}$. Then, we perform a further transformation using the unitary operator $T(\alpha)$, defined as $T(\alpha)=1/\sqrt{2}\left[D^{\dagger}(\alpha)\left(\left|e\right\rangle\left\langle e\right|-\left|g\right\rangle\left\langle e\right|\right)+D(\alpha)\left(\left|g\right\rangle\left\langle g\right|+\left|e\right\rangle\left\langle g\right|\right)\right]$ with $D(\alpha)=e^{\alpha a^{\dagger}-\alpha^{*}a}$ the standard displacement operator. Hence, $H_{b}\equiv T^{\dagger}(-\lambda/\Omega)H_{a}T(-\lambda/\Omega)$ such that $\rho_{b}=T^{\dagger}\rho_{a}T$, which leads to (see Appendix B for further details):

[TABLE]

where the Hamiltonian $H_{b}$ can be written as:

[TABLE]

Hence, the dissipator acting on $\rho_{b}$ has the same form as in Equation (14), but with transformed operators, namely $T^{\dagger}xT$, $T^{\dagger}\hat{\chi}T$, and $T^{\dagger}\hat{\Theta}T$, where $T\equiv T(-\lambda/\Omega)$. Finally, by moving to an interaction picture with respect to $H_{b,0}=(\Omega-\tilde{\nu})a^{\dagger}a-\tilde{\omega}\sigma_{z}/2$ and expanding the exponential in Equation (16) (the latter requires that $|2\lambda/\Omega|\sqrt{\left<(a+a^{\dagger})^{2}\right>}\ll 1$ for truncating the exponential to a finite number of terms), we arrive at a Hamiltonian containing multiphoton interaction terms. The latter condition is commonly known as the Lamb-Dicke regime. In addition, we consider the driving frequencies to be $\Delta_{j}=\pm n_{j}(\tilde{\nu}-\Omega)-\tilde{\omega}$ with $|\Omega-\tilde{\nu}|\gg\epsilon_{j}/2$, so that one can safely perform a rotating-wave 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 Leibfried et al. (2003), the procedure to obtain Jaynes-Cummings or quantum Rabi models is analogous to those cases Pedernales et al. (2015); Felicetti et al. (2015); Lv et al. (2018). In this manner, we can approximate $H_{b,1}^{I}\equiv U_{b,0}^{\dagger}H_{b,1}U_{b,0}\approx H_{\rm n}$, where $H_{\rm n}$ contains the aimed at multiphoton interactions,

[TABLE]

Note that the sets $r$ and $b$ encompass the terms with amplitude $\epsilon_{j}$ driving red- and blue-sidebands, that is those terms in Equation (II) with frequency $\Delta_{j\in r}=+n_{j}(\tilde{\nu}-\Omega)-\tilde{\omega}$ and $\Delta_{j\in b}=-n_{j}(\tilde{\nu}-\Omega)-\tilde{\omega}$. Each of these drivings will contribute with a multiphoton interaction, either $\sigma^{+}a^{n_{j}}+{\rm H.c.}$ for $j\in r$ or $\sigma^{-}a^{n_{j}}+{\rm H.c.}$ for $j\in b$, which produce transitions between the states $\left|m\right\rangle\left|g\right\rangle\leftrightarrow\left|m\mp n_{j}\right\rangle\left|e\right\rangle$.We stress that for a time-independent spin-boson model, as given in Equations (1)–(4) (or equivalently with $n_{d}=0$ in $H_{S,d}$ as given in Equation (II), one obtains a single $n$-photon [anti]-Jaynes–Cummings interaction term, $\sigma^{+}a^{n}+{\rm H.c.}$ [$\sigma^{+}(-a^{\dagger})^{n}+{\rm H.c.}$], by choosing $\Delta_{0}=n(\tilde{\nu}-\Omega)-\tilde{\omega}$ [$\Delta_{0}=-n(\tilde{\nu}-\Omega)-\tilde{\omega}$] in the original spin–boson Hamiltonian $H_{\rm SB}$. Thus, one needs the knowledge of the relevant bosonic frequency $\Omega$ to simulate multiphoton interaction terms properly.

In order to show how the dissipative part transforms, it is advisable to introduce the time-dependent unitary operator:

[TABLE]

Then, one can see that, defining $\tilde{\chi}=\Phi\chi\Phi^{\dagger}$, $\tilde{\Theta}=\Phi\Theta\Phi^{\dagger}$ and $\tilde{x}=\Phi(a+a^{\dagger})\Phi^{\dagger}$, the resulting master equation for $\rho_{\rm n}(t)$ is:

[TABLE]

where the state $\rho_{\rm n}(t)$ of the multiphoton model is related to the original spin-boson upon the reaction coordinate mapping, $\rho_{\rm S+RC}(t)$, through a unitary transformation:

[TABLE]

From the previous expression, it follows that the purity of the total state $\rho_{\rm S+RC}$ and that of $\rho_{\rm n}$ are approximately equal. Moreover, the reduced spin state in the different frameworks are related according to ${\rm Tr}_{\rm B}[\rho_{\rm SB}(t)]={\rm Tr}_{\rm RC}[\rho_{\rm S+RC}(t)]\approx{\rm Tr}_{\rm RC}[\Phi^{\dagger}\rho_{\rm n}(t)\Phi]$, where ${\rm Tr}_{\rm B}[\cdot]$ and ${\rm Tr}_{\rm RC}[\cdot]$ 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_{\rm n}$, given in Equation (III.1), 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 spin–boson model as $\rho_{\rm n}(0)=T^{\dagger}\rho_{\rm S+RC}(0)T$.

At this stage, a few comments regarding the validity of Equation (20) are in order. While the steps performed from $H_{\rm S+RC}$ to $H_{b}$ are exact, $H_{\rm n}$ is attained in an approximate manner. The good functioning of the simulation depends on how these approximations are met. That is, Equation (20) holds within the Lamb-Dicke regime $|2\lambda/\Omega|\sqrt{\left<(a+a^{\dagger})^{2}\right>}\ll 1$ and for parameters satisfying $|\Omega-\tilde{\nu}|\gg\epsilon_{j}/2\ \forall j$, so that one can perform a rotating-wave approximation. As a consequence, this approximation also sets a constraint on the total duration for a good simulation (see Appendix B). Note that, as the parameters $\lambda$ and $\Omega$ 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 Equation (19) must be small compared to the parameters involved in $H_{\rm n}$. For the considered shape of $J_{\rm SB}(\omega)$, this translates into $\Gamma\ll\tilde{\nu},\tilde{g}_{n}$, where $\tilde{g}_{n}=\epsilon_{0}(2\lambda)^{n}/(2\Omega^{n}n!)$ for an $n_{d}=0$ and $\Delta_{0}=\pm n(\tilde{\nu}-\Omega)-\tilde{\omega}$ (cf. Equation (III.1).

Finally, we comment that the previous scheme can be carried out beyond the Lamb-Dicke regime Puebla et al. (2019). Admittedly, when the Lamb-Dicke approximation does not hold, the Hamiltonian $H_{\rm n}$ is no longer a good approximation to the dynamics. In this case, the Hamiltonian $H_{\rm 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 de Matos Filho and Vogel (1994, 1996); Vogel and Filho (1995); Cheng et al. (2018). These nonlinear, yet multiphoton Hamiltonians appear then as a good approximation to $H_{b}$, and thus to $H_{\rm SB}$ whenever $|2\lambda/\Omega|\sqrt{\left<(a+a^{\dagger})^{2}\right>}\ll 1$ is not fulfilled, as recently shown in Puebla et al. (2019). In this article, however, we will constrain ourselves to parameters within the Lamb-Dicke regime.

III.2 Structured environments

As previously mentioned, the simulation of multiphoton spin–boson interactions is not restricted to a determined form of $J_{\rm SB}(\omega)$. 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_{\rm SB}(\omega)$ can be split in two parts, $J_{\rm SB}(\omega)=J_{\rm SB,1}(\omega)+J_{\rm SB,2}(\omega)$, although its generalization to more is straightforward. The first contribution, $J_{\rm SB,1}(\omega)$, is considered here to be suitable for the realization of multiphoton interactions as described in III.1. In addition, we will work under the assumption that the environment degrees of freedom corresponding to $J_{\rm SB,2}(\omega)$ can be treated and simplified using again a collective or reaction coordinate, as sketched in Figure 1(c).

As discussed previously, we identify a collective coordinate for each of the contributions to the spectral density $J_{\rm SB}(\omega)$. In this manner, we augment the system to include both reaction coordinates, denoted here by ${\rm S^{\prime}}={\rm{S}}+{\rm{RC}}_{1}+{\rm{RC}}_{2}$. Hence, its Hamiltonian is given by:

[TABLE]

where $H_{\rm S,d}$ is the original spin Hamiltonian, which may contain spin rotations, introduced in Equation (II), while the subscripts denote the corresponding reaction coordinate. The parameters $\lambda_{i}$ and $\Omega_{i}$ are determined by the spectral density $J_{\rm SB,i}(\omega)$. The dynamics of the augmented system $\rm S^{\prime}$ is governed by the following master equation:

[TABLE]

where $x_{i}=a_{i}+a^{\dagger}_{i}$ for $i=1,2$, and $\chi_{i}$ and $\Theta_{i}$ are defined in analogy to Equations (11)-(12).

In order to find a suitable transformation to realize multiphoton interaction terms from $H_{\rm S^{\prime}}$, we proceed in a similar manner as for a single reaction coordinate. That is, we first move to a rotating frame where $H_{\rm S^{\prime}}\equiv H_{a,1}^{I}$, with $H_{a}=H_{a,0}+H_{a,1}$ and $H_{a,0}=-\Delta_{0}/2\sigma_{x}$. Therefore, the transformed Hamiltonian reads as:

[TABLE]

The next step is to perform the transformation using the unitary operator $T(\alpha)$. As previously mentioned, 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 one to choose $\alpha\equiv-\lambda_{1}/\Omega_{1}$, hence $H_{b}\equiv T^{\dagger}(-\lambda_{1}/\Omega_{1})H_{a}T(-\lambda_{1}/\Omega_{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}=(\Omega_{1}-\tilde{\nu}_{1})a^{\dagger}_{1}a_{1}-\tilde{\omega}\sigma_{z}/2$, we obtain the Hamiltonian $H_{\rm n,2}\approx H_{b,1}^{I}\equiv U_{b,0}^{\dagger}H_{b,1}U_{b,0}$,

[TABLE]

where we have considered $\Delta_{j}=\pm n_{j}(\tilde{\nu}-\Omega_{1})-\tilde{\omega}$ and assumed the Lamb-Dicke regime $|\lambda_{1}/\Omega_{1}|\sqrt{\langle(a_{1}+a^{\dagger}_{1})^{2}\rangle}\ll 1$, and $|\Omega_{1}-\tilde{\nu}|\gg\epsilon_{j}/2$ to perform a rotating-wave approximation. Note that, while the multiphoton terms are identical to those of $H_{\rm n}$ in Equation (III.1), the second reaction coordinate interacts with the spin degree of freedom. Indeed, depending on the parameters of $H_{\rm n,2}$, the effect of such an interaction may turn effectively into non-Markovian effects for the reduced state of the spin and first reaction coordinate, $\rho_{\rm n}={\rm Tr}_{2}[\rho_{\rm n,2}]$. The final master equation governing the dynamics of $\rho_{\rm n,2}$ is:

[TABLE]

where the operators involved are defined as in the case involving a single reaction coordinate (cf. Equation (19)). It is worth stressing that the relation between the states given in Equation (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_{\rm n}$, as well as to discuss the limitation in the parameter regime for their realization. In particular, in Section IV.1, we first consider the case in which the original spin-boson model interacts just with a discrete number of modes, which can be viewed as a limit of vanishing spectral broadening $\Gamma\rightarrow 0$. This scenario will allow us to examine the validity of the required approximations without the effect of dissipation. Then, in Section IV.2, we will consider $\Gamma\neq 0$, where the reaction-coordinate mapping appears as a key step to realize a desired multiphoton Jaynes-Cummings model. The dynamics of each model is obtained by a standard numerical integration (fourth-order Runge-Kutta) of the corresponding master equation, namely Equations (III.1) and (19) for the spin-boson and multiphoton Jaynes-Cummings model, respectively. Note that for a structured environment, the master equations are given in Equations (III.2) and (III.2).

In all 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,

[TABLE]

In particular, we will analyse to what extent is the relation given in Equation (20) satisfied. In other words, we will compare the aimed state of a multiphoton Jaynes-Cummings model $\rho_{\rm n}(t)$ with the one retrieved using the analogue simulator, $\Phi\rho_{\rm S+RC}(t)\Phi^{\dagger}$, that is $\rho_{1}(t)\rightarrow\rho_{\rm n}(t)$ and $\rho_{2}(t)\rightarrow\Phi\rho_{\rm S+RC}(t)\Phi^{\dagger}$ in Equation (26). We remark that when two reaction coordinates are included, the state $\rho_{\rm n}(t)$ obeys the master equation given in Equation (III.2), whose Hamiltonian is $H_{\rm n,2}$, Equation (III.2), while $\rho_{\rm S+RC}(t)$ must be replaced by $\rho_{\rm S^{\prime}}$, as explained in III.2.

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 de Vega and Alonso (2017), we resort to the one based on the trace distance Breuer et al. (2009), defined as:

[TABLE]

where $|A|=\sqrt{A^{\dagger}A}$. Then, non-Markovian evolutions can be characterized as those for which $\mathcal{D}(\rho_{x}(t),\rho_{y}(t))$ increases during certain time intervals, that is for those for which the time-derivative of the trace distance for a pair of states $\rho_{x,y}$,

[TABLE]

is $\sigma(t,\rho_{x,y})>0$. In general, one has to maximize over all possible pairs of states $\rho_{x,y}$ in order to find a suitable non-Markovian measure Breuer et al. (2009). For our purpose, however, it will be sufficient to show that $\sigma(t,\rho_{x,y})>0$ for a certain pair of states in a multiphoton Jaynes-Cummings model and that it can be retrieved using a spin-boson model. That is, we calculate $\sigma(t,\rho_{x,y})$ using two initial states $\rho_{x,y}$ in the multiphoton Jaynes-Cummings model and corroborate that $\sigma(t,\rho_{x,y})$ is obtained to a very good approximation when the states $\rho_{x,y}(t)$ are replaced by their simulated ones using the spin-boson model, namely $\rho_{x}(t)\rightarrow\Phi\rho_{\rm x,S+RC}(t)\Phi^{\dagger}$ and $\rho_{y}(t)\rightarrow\Phi\rho_{\rm y,S+RC}(t)\Phi^{\dagger}$. In this manner, we offer a proof-of-principle that non-Markovian multiphoton models can be realized.

IV.1 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 considering $\Gamma\rightarrow 0$ in the underdamped spectral density $J_{\rm SB}(\omega)$ or, equivalently, assuming that dissipation effects are sufficiently small so that they can be discarded. Note that for a single bosonic mode with $\Gamma=0$, the spin–boson model adopts the form of a generalized quantum Rabi model, which is indeed $H_{\rm S+RC}$, as given in Equation (III.1). Recall that in this particular case, $H_{\rm SB}\equiv H_{\rm S+RC}$, as there are no further modes in the system. In particular, we set $n_{d}=0$ in Equation (II) as we aim to realize a single multiphoton Jaynes-Cummings interaction. The Hamiltonian for a nJCM can be written in general as:

[TABLE]

At resonant condition, $\tilde{\omega}=n\tilde{\nu}$, the coupling constant $\tilde{g}_{n}$ fixes the time required to transfer the population from the state $\left|e\right\rangle\left|0\right\rangle$ to $\left|g\right\rangle\left|n\right\rangle$, denoted as $\tau_{n}=\pi/(2\tilde{g}_{n}\sqrt{n!})$. Both are related to the spin-boson parameters as (cf. Equation (III.1)):

[TABLE]

Clearly, as $2\lambda/\Omega$ must be small to lie within the Lamb-Dicke regime, the coupling $\tilde{g}_{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 $\tau_{n}$.

In Figure 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 $\left|e\right\rangle\left|0\right\rangle$ and $\left|g\right\rangle\left|n\right\rangle$, we choose $\rho_{\rm S+RC}(0)=\left|-\right\rangle\left\langle-\right|\otimes\rho_{\rm RC}^{\rm th}$ as an initial state for the spin-boson model, where $\rho_{\rm RC}^{\rm th}$ is a thermal state at temperature $\beta^{-1}$ for the reaction coordinate mode, containing $n^{\rm th}=(e^{\beta\Omega}-1)^{-1}$ bosons. Recall that, as we consider here a single spectral density with $\Gamma=0$, the reaction coordinate mode is simply the unique mode that interacts with the spin degree of freedom. In this manner, the initial state for the simulated multiphoton models reads as $\rho_{\rm nJCM}(0)=T^{\dagger}\rho_{\rm S+RC}(0)T$, which approximately amounts to $\rho_{\rm nJCM}(0)\approx\left|e\right\rangle\left\langle e\right|\otimes\left|0\right\rangle\left\langle 0\right|$ for sufficiently low temperature and small $2\lambda/\Omega$. The chosen parameters for the simulation of the 2JCM, plotted in Figure 2(a) and (b), are $\pi\alpha=\epsilon_{0}=0.02\omega_{0}$; recalling that $\Omega=\omega_{0}$, it results in $2\lambda/\Omega=0.2$. Choosing $\tilde{\nu}=10^{-3}\Omega$ and $\tilde{\omega}=2\tilde{\nu}$, the coupling in 2JCM amounts to $\tilde{g}_{2}=0.2\tilde{\nu}$. The initial reaction-coordinate thermal state, $\rho_{\rm RC}^{\rm th}$, contains $n^{\rm th}=10^{-3}$ bosons. In Figure 2(b), we show how the quantum simulation of the 2JCM model deteriorates for increasing number of bosons, as a large $n^{\rm th}$ will eventually break down the Lamb-Dicke regime.

For the 3JCM, we choose again $\pi\alpha=0.02\omega_{0}$, which leads to $2\lambda/\Omega=0.2$. Then, we select the aimed coupling strength of the multiphoton interaction to be $\tilde{g}_{3}=0.1\tilde{\nu}$ with $\tilde{\omega}=3\tilde{\nu}$, while we vary $\epsilon_{0}/\omega_{0}$. The temperature is set to $\beta\Omega\approx 100$ so that $\rho_{\rm RC}^{\rm th}\approx\left|0\right\rangle\left\langle 0\right|$. As in the previous case, the dynamics are well retrieved; see Figure 2(c), where we have set $\epsilon_{0}/\omega_{0}=2\cdot 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. Equation (III.1)) and due to the longer times required to simulate a 3JCM compared to the 2JCM, the condition $|\Omega-\tilde{\nu}|\gg\epsilon_{0}$ must be better satisfied. Indeed, for $\epsilon_{0}/\omega_{0}=10^{-2}$, we already see a clear departure from the targeted dynamics, as indicated by a large infidelity $1-F(t)\gtrsim 10^{-1}$, as shown in Figure 2(d).

In the following, we consider a spin interacting with two bosonic modes, again with $\Gamma_{1,2}=0$. As explained in III.2, 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_{\rm nJCM,2}$, where the subscript $2$ indicates the presence of a second reaction coordinate in the system. The Hamiltonian $H_{\rm nJCM,2}$ reads as:

[TABLE]

In this manner, the spin exchanges $n$ quanta with the first bosonic mode as in $H_{\rm 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 $\rho_{\rm nJCM}(t)={\rm Tr}_{2}[\rho_{\rm 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_{\rm S^{\prime}}$ given in Equation (III.2), which then approximately realizes $H_{\rm nJCM,2}$. In particular, we select $\Delta_{0}=-2\Omega_{1}$, so that the simulated model involves two-photon interaction terms, i.e., a 2JCM. The results are plotted in Figure 3, while the parameters are $\pi\alpha_{i}=0.02\Omega_{i}$ such that $2\lambda_{i}/\Omega_{i}=0.2$ for $i=1,2$, $\epsilon_{0}/\Omega_{1}=10^{-2}$. The coupling strength in $H_{\rm 2JCM,2}$ is given by $\tilde{g}_{2}=0.2\tilde{\nu}$ with $\tilde{\nu}=\Omega_{2}$. As in the single-mode case, Rabi oscillations will be clearly visible selecting $\rho_{\rm S^{\prime}}(0)=\left|-\right\rangle\left\langle-\right|\otimes\rho_{\rm RC_{1}}^{\rm th}\otimes\rho_{\rm RC_{2}}^{\rm th}$. After its transformation, this state corresponds approximately to an initial spin state $\left|e\right\rangle$ in the nJCM frame. In the same manner, in order to analyse the emergence of non-Markovian behaviour, we consider the initial states $\left|g\right\rangle\left\langle g\right|$ and $\left|e\right\rangle\left\langle e\right|$ for the spin in $H_{\rm S^{\prime}}$. This implies initial spin states $\left|\pm\right\rangle$ in the nJCM frame, which for pure dephasing noise, it has been shown to be the pair of states maximizing $\sigma(t)$ Breuer et al. (2009). The results plotted in Figure 3 have been performed considering a sufficiently low temperature such that $\rho^{\rm th}_{\rm RC_{1,2}}\approx\left|0\right\rangle\left\langle 0\right|$. We then compute the trace distance $\mathcal{D}(\rho_{x},\rho_{y})$ using the states $\rho_{x,y}(t)$ resulting in tracing out the second mode, ${\rm Tr}_{\rm 2}[\rho_{\rm 2JCM,2}(t)]$, for the two different initial states $\rho_{2JCM,2}(0)\approx\left|\pm\right\rangle\left\langle\pm\right|\otimes\rho_{\rm RC_{1}}^{\rm th}\otimes\rho_{\rm RC_{2}}^{\rm th}$. As shown in Figure 3(b), the time-derivative of the trace distance, $\sigma(t)$, becomes positive during certain intervals, a clear indication of the non-Markovian behaviour of the simulated multiphoton Jaynes-Cummings model. In addition, we also calculate the non-trivial evolution of the purity for the states $\rho_{\rm S+RC_{1}}(t)$ and $\rho_{\rm S}(t)={\rm Tr}_{\rm RC_{1}}[\rho_{\rm S+RC_{1}}(t)]$, which is shown in Figure 3(c). According to our theoretical framework, their purity is approximately equal to that of $\rho_{\rm 2JCM}(t)$ and the reduced spin state upon tracing both bosonic degree of freedom in the 2JCM, ${\rm Tr}[\rho_{\rm 2JCM}(t)]$, respectively. Finally, the infidelity $1-F(t)$ between the targeted state $\rho_{\rm 2JCM,2}(t)$ and its reconstructed one $\Phi\rho_{\rm S+RC_{1}+RC_{2}}(t)\Phi^{\dagger}$ in Figure 3(d).

IV.2 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_{\rm SB}(\omega)$ has the form of Equation (8) with $\Gamma\neq 0$. In this manner, we extend the theoretical framework beyond the standard local master equation description Puebla et al. (2019). As explained in Section III.1, 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 (III.1). Indeed, the effect of spectral broadening, $\Gamma\neq 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 $\Gamma\gg\tilde{\nu}$, dissipation dominates the dynamics, and the paradigmatic Rabi oscillations will eventually fade away. In Figure 4, we show the results of numerical simulations aimed to retrieve a 2JCM with different $\Gamma/\tilde{\nu}$ values and for different quantities. As for Figure 2, we used $\pi\alpha=\epsilon_{0}=0.02\omega_{0}$, so that $2\lambda/\Omega=0.2$. We chose again $\tilde{\nu}=10^{-3}\Omega$ and $\tilde{\omega}=2\tilde{\nu}$, and therefore, the coupling in 2JCM amounts to $\tilde{g}_{2}=0.2\tilde{\nu}$, while the temperature is such that $\rho_{\rm RC}^{\rm th}$ contains $n^{\rm th}=10^{-3}$ bosons. The spin is initialized in the $\left|-\right\rangle$ state, so that $\rho_{\rm S+RC}(0)=\left|-\right\rangle\left\langle-\right|\otimes\rho_{\rm RC}^{\rm th}$. In particular, the value $\Gamma/\tilde{\nu}=2\cdot 10^{-1}$ considered in Figure 4(a) already produces a significant departure from the Rabi oscillation between the states $\left|e\right\rangle\left|0\right\rangle$ and $\left|g\right\rangle\left|2\right\rangle$ in the dissipationless 2JCM (cf. Figure 2(a) for $\Gamma=0$). Note that the results plotted in Figure 4(a) correspond to a critically-damped 2JCM since $\Gamma=\tilde{g}_{2}$. As plotted in Figure 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 ${\rm Tr}[\rho_{\rm S+RC}^{2}(t)]$ and ${\rm Tr}[\rho_{\rm S}^{2}(t)]$. As in previous cases, the purity of these states is directly related to those of the simulated model as a consequence of the relation $\rho_{\rm 2JCM}(t)\approx\Phi\rho_{S+RC}(t)\Phi^{\dagger}$. Furthermore, Rabi oscillations or population revivals appear in the evolution of von Neumann entropy, $S_{\rm vN}(\rho)=-\rho\log_{2}\rho$ for the reduced spin state. In particular, for an initial state $\rho_{nJCM}(0)\approx\left|e\right\rangle\left\langle e\right|\otimes\left|0\right\rangle\left\langle 0\right|$ and due to the $n$-photon interaction with a bosonic degree of freedom, the spin state oscillates between a pure ($S_{\rm vN}=0$) and a maximally-mixed state ($S_{\rm vN}=1$) in a time $\tau_{n}/2$. This further corroborates that one can witness the multiphoton transitions of the aimed multiphoton Jaynes-Cummings model monitoring the spin even without access or control on the bosonic environment. This is plotted in Figure 4(c) for different $\Gamma/\tilde{\nu}$ 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. Figure 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 analogue quantum simulator. While the spin-boson model naturally lacks 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. Importantly, the parameters of the multiphoton model, as well as the order of the interaction can be controlled by tuning the frequency splitting and bias parameter of the spin in the original spin-boson model. 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 reaction-coordinate method Thoss et al. (2001); Martinazzo et al. (2011); Iles-Smith et al. (2014, 2016); Strasberg et al. (2016, 2018); Nazir and Schaller (2018). 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 a Markovian fashion. The resulting Hamiltonian is then used to generate multiphoton interaction terms, following the theory explained in Casanova et al. (2018); Puebla et al. (2019), 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 the 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 models. 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 models can be indeed attained when a second reaction coordinate is included, as unveiled by the standard trace distance measure Breuer et al. (2009). Finally, we provide numerical simulations investigating the interplay between spectral broadening, dissipation and the decoherence in the targeted multiphoton models.

Acknowledgements.

G.Z. is supported by the H2020-MSCA-COFUND-2016 project SPARK(Grant No. 754507). R.P. and M.P. acknowledge the support by the SFI-DfEInvestigator Programme (Grant 15/IA/2864). M.P. acknowledges the H2020 Collaborative Project TEQ(Grant Agreement 766900), the Leverhulme Trust Research Project Grant UltraQuTe (Grant No. RGP-2018-266) and the Royal Society Wolfson Fellowship (RSWF\R3\183013). J.C. acknowledges support by the Juan de la Cierva Grant IJCI-2016-29681. I.A. acknowledges support by Basque Government Ph.D. Grant No. PRE-2015-1-0394. We also acknowledge funding from Spanish MINECO/FEDER FIS2015-69983-P and Basque Government IT986-16. This material is also based on work supported by the U.S. Department of Energy, Office of Science, Office of Advance Scientific Computing Research (ASCR), Quantum Algorithm Teams (QAT) Program under Field Work Proposal Number ERKJ333. J.C. and E.S. acknowledge support from the projects QMiCS(820505) and OpenSuperQ(820363) of the EU Flagship on Quantum Technologies.

Appendix A Reaction coordinate mapping

In this Appendix, we provide the necessary steps for the reaction coordinate mapping, as well as for the derivation of the master equation given in Equation (III.1), following closely Iles-Smith et al. (2014). As outlined in Section III.1, given the Hamiltonian of the spin–boson system $H_{\rm SB}=\frac{\epsilon_{0}}{2}\sigma_{z}+\frac{\Delta_{0}}{2}\sigma_{x}+\sigma_{x}\sum_{k}f_{k}(c_{k}+c_{k}^{\dagger})+\sum_{k}\omega_{k}c_{k}^{\dagger}c_{k}$, one can achieve the RC mapping by defining a collective coordinate such that $\lambda(a+a^{\dagger})=\sum_{k}f_{k}(c_{k}+c^{\dagger}_{k})$, where $a$ and $a^{\dagger}$ are respectively the annihilation and creation operators of the RC. This transformation leads to a new Hamiltonian where the original system interacts with the residual environment only through the RC,

[TABLE]

where $H_{\rm S+RC}$ is given by Equation (III.1), while $H_{\rm RC-E^{\prime}}=(a+a^{\dagger})\sum_{k}g_{k}(b_{k}+b_{k}^{\dagger})+(a+a^{\dagger})^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}}$, $H_{\rm E^{\prime}}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$.

The crucial point of this procedure is to find an explicit relation between the spectral density of the original configuration, i.e. $J_{\rm SB}(\omega)=\sum_{k}f_{k}^{2}\delta(\omega-\omega_{k})$, and the analogue quantity of the transformed system $J_{\rm RC}(\omega)=\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k})$. In order to obtain this relation, one can rephrase the problem classically. Indeed, since the spectral density only depends on the interaction between the system and the environment, one can momentarily regard the spin as a continuous coordinate $q$ subject to a potential $V(q)$. After solving the corresponding Hamilton equations of motion in the Fourier space, one obtains an equation of the form $\hat{L}_{\rm SB}(z)\hat{q}(z)=-\hat{V}^{\prime}(z)$, where $\hat{L}_{\rm SB}(z)=-z^{2}\left(1+\int_{0}^{+\mathcal{1}}d\omega\frac{2J_{\rm SB}(\omega)}{\omega(\omega^{2}-z^{2})}\right)$. Therefore, using the so-called Leggett prescription, one gets:

[TABLE]

One can reproduce the same calculation also after performing the RC mapping and express $J_{\rm RC}(\omega)$ in terms of the corresponding kernel $\hat{L}_{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 use $\hat{L}_{0}(z)$ instead of $\hat{L}_{\rm SB}(z)$ in Equation (34). By considering the Ohmic spectral density $J_{\rm RC}(\omega)=\gamma\omega e^{-\omega/\Lambda}$, one obtains:

[TABLE]

It is easy to see that one exactly recovers the underdamped spectral density given by Equation (8) by simply requiring that $\gamma=\Gamma/(2\pi\omega_{0})$, $\Omega=\omega_{0}$, and $\lambda=\sqrt{\pi\alpha\omega_{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 us 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 time scale of the system. Within this approximation, one can work out a master equation that, in the Schrödinger picture, reads as:

[TABLE]

where $\rho\equiv\rho_{\rm S+RC}$, $A=a+a^{\dagger}$, and the residual environment is assumed to be in a thermal state, i.e., $\rho_{E^{\prime}}=e^{-\beta H_{\rm E^{\prime}}}/{\rm{Tr}}_{\rm E^{\prime}}\{e^{-\beta H_{\rm E^{\prime}}}\}$.

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_{\rm S+RC}$. To this end, let $\left|\phi_{n}\right\rangle$ be an eigenstate of $H_{\rm S+RC}$, i.e., $H_{\rm S+RC}\left|\phi_{j}\right\rangle=\varphi_{j}\left|\phi_{j}\right\rangle$; therefore, the operator $A$ can be expanded as $A=\sum_{jk}A_{jk}\left|\phi_{j}\right\rangle\left\langle\phi_{k}\right|$, while in the interaction picture, one has:

[TABLE]

where $A_{jk}=\left\langle\phi_{j}\right|A\left|\phi_{k}\right\rangle$, and $\xi_{jk}=\varphi_{j}-\varphi_{k}$. Finally, by plugging Equation (37) into Equation (36) and assuming the imaginary parts to be negligible, one gets the final form of the master equation given by Equation (III.1).

Appendix B Derivation of $H_{b}$ and $H_{\rm n}$

In this Appendix, we show how to obtain the Hamiltonians $H_{b}$ and $H_{\rm n}$, given in Equations (16) and (III.1), respectively. In particular, for $H_{b}$, the following expressions are needed:

[TABLE]

Thus, the resulting Hamiltonian $H_{b}=T^{\dagger}H_{a}T$, with $H_{a}=\Omega a^{\dagger}a+\lambda\sigma_{x}(a+a^{\dagger})+\sum_{j}\epsilon_{j}/2(\cos\Delta_{j}t\sigma_{z}+\sin\Delta_{j}t\sigma_{y})$, reads:

[TABLE]

where we have neglected a constant energy shift. Therefore, by selecting $\alpha=-\lambda/\Omega$, we obtain a simple Hamiltonian to pursue multiphoton interactions, namely:

[TABLE]

which is indeed Equation (16). Moving now to an interaction picture w.r.t. $H_{b,0}=(\Omega-\tilde{\nu})a^{\dagger}a-\tilde{\omega}\sigma_{z}/2$, we obtain:

[TABLE]

with $a(t)=ae^{-i(\Omega-\tilde{\nu})t}$. Requiring $|2\lambda/\Omega|\sqrt{\left<(a+a^{\dagger})^{2}\right>}\ll 1$ and selecting $\Delta_{j}=\Delta_{n}^{\pm}\equiv\pm n(\tilde{\nu}-\Omega)-\tilde{\omega}$, we resonantly drive multiphoton Jaynes–Cummings interaction terms, while the rest of the terms in the expansion of the exponential term are off-resonant and rotating with a large frequency compared to its amplitude, i.e., $n|\Omega-\tilde{\nu}|\gg\epsilon_{j}/2$ (for zeroth order) where $n$ is the selected order of the interaction $\sigma^{\pm}a^{n}$. In this manner, performing these two approximations, one obtains:

[TABLE]

where $\Delta_{j\in r}=n_{j}(\tilde{\nu}-\Omega)-\tilde{\omega}$ and $\Delta_{j\in b}=-n_{j}(\tilde{\nu}-\Omega)-\tilde{\omega}$, which corresponds to Equation (III.1). The largest error committed in the previous approximation stems from the zeroth order in the expansion of the exponential. These contributions are of the form $\epsilon_{j}/2(\sigma^{+}e^{in_{j}(\Omega-\tilde{\nu})t}+{\rm H.c.})$, which will produce a significant effect after a time $t\approx n_{j}(\Omega-\tilde{\nu})/\epsilon^{2}_{j}$. For a single $n$-photon interaction term, population transfer occurs in a characteristic time $\tau_{n}=\sqrt{n!}(\Omega/2\lambda)^{n}/\epsilon_{0}$ (see Section IV.1). Hence, we can provide a rough estimate for the duration of a correct simulation of the desired multiphoton Jaynes–Cummings model to be $t=k\tau_{n}$ with $k\approx(2\lambda/\Omega)^{n}n(\Omega-\tilde{\nu})/(\epsilon_{0}\sqrt{n!})$.