Connecting nth order generalised quantum Rabi models: Emergence of nonlinear spin-boson coupling via spin rotations

We establish an approximate equivalence between a generalised quantum Rabi model and its nth order counterparts where spin-boson interactions are nonlinear as they comprise a simultaneous exchange of $n$ bosonic excitations. Although there exists no unitary transformation between these models, we demonstrate their equivalence to a good approximation in a wide range of parameters. This shows that nonlinear spin-boson couplings, i.e. nth order quantum Rabi models, are accessible to quantum systems with only linear coupling between boson and spin modes by simply adding spin rotations and after an appropriate transformation. Furthermore, our result prompts novel approximate analytical solutions to the dynamics of the quantum Rabi model in the ultrastrong coupling regime improving previous approaches.

Introduction.-Thequantum Rabi model (QRM) lies not only at the heart of our understanding of light-matter interaction [1], but is also of considerable importance in diverse fields of research [2].Certainly, the Rabi model was primarily proposed to describe a nuclear spin interacting with classical radiation [3,4], whose first fully quantized version only appeared two decades later [5].The QRM contemplates a scenario which is of great generality as it encompasses two of the most basic, yet essential, ingredients in quantum physics, namely, a two-level system and a bosonic mode.Indeed, this model emerges in disparate settings, ranging from ion traps [6,7] to circuit or cavity QED [8,9], quantum optomechanical systems [10], color-centers in membranes [11], and cold atoms [12].
Even though the QRM has been exhaustively investigated during the last decades, a number of recent findings has brought it again into the research spotlight.Specifically, its integrability [13], the existence of a distinctive behavior beyond the Jaynes-Cummings doublets in the deep strong coupling regime [14], and the emergence of a quantum phase transition in a suitable limit [15][16][17][18], demonstrate that the physics of the QRM is far from being completely understood and exploited.Closely related to the QRM, we find the nth order QRM (nQRM) which differs from the QRM in that the nQRM comprises n-boson exchange interaction terms.This non-linear generalization of the QRM has recently attracted attention, mainly in its second-order form (2QRM) as it shows striking phenomena such as spectral collapse [19][20][21] and due to its relevance in preparing non-classical states of light in quantum optics [22,23].In addition, because of their different spectra, there is no unitary map between the QRM and the nQRM.
In this Letter, we demonstrate the existence of an approximate equivalence among a family of Hamiltonians comprising n-boson interaction terms, where the standard QRM or the two-photon QRM appear as special cases.As we will show, the dynamics of a nQRM is retrieved after an appropriate transformation of a generalized driven QRM (gQRM).As a proof of concept, we show how the dynamics of a 2QRM and a 3QRM can be captured without having access to the required twoand three-photon interaction terms.That is, a quantum system that contains a linear spin-boson coupling but lacks of non-linear interactions is sufficient for the simulation of models where non-linear terms are crucial.Indeed, because creating n-boson interactions is considered challenging in many quantum platforms, our method might open new avenues for their inspection.Moreover, we present a potentially scalable physical platform [24], a microwave-driven trapped ion setting [25][26][27][28], in which nQRMs are unattainable without resorting to our equivalence, which highlights its applicability.Finally, we specialize the general theory to analyze the standard QRM and find that our approach provides approximate analytical solutions that surpass in accuracy the standard Bloch-Siegert approximation of the QRM in the ultrastrong coupling regime [29][30][31].
Theory.-Our starting point is the Hamiltonian whose first two terms correspond to a bosonic mode of frequency ν and a two-level system with a frequency splitting ω, described by the usual annihilation (creation) operator a (a † ) and spin-1 2 Pauli matrices σ = (σ x , σ y , σ z ), respectively.Both subsystems interact through a set of coupling terms with amplitude Ω/2, considered here equal ∀j for simplicity, and α j being a phase that could be time dependent.We will show that Hamiltonian H s is central for our theory, as sketched in Fig. 1, which establishes an approximate map between QRM dynamics with those of nQRM.To this end, we perform a unitary transformation on is the displacement operator.Note that this transformation has been used in the context of trapped ions to derive the eigenstates of a system that comprises a laser interacting with a trapped ion, and for a fast implementation of the QRM [32,33].Now, choosing time dependent phases, α j = (ω + δ j )t, and moving to a rotating frame with respect to H T,0 = −(ω + δ 1 )σ x /2, the resulting Hamiltonian, denoted by H gQRM , reads (for more details see [34]) with p = i(a † − a) and U T,0 = e −itH T ,0 .The previous Hamiltonian corresponds to a more general version of the QRM, the gQRM, where the last term on Eq. ( 2) can be viewed as a classical driving acting on the system.In particular, we note that H gQRM adopts the form of a standard QRM in the case of having δ j = 0 ∀j.On the other hand, the starting Hamiltonian H s in Eq. ( 1) can be brought into the form of a H nQRM by properly choosing α j and in a suitable interaction picture.More specifically, by defining H s = H s,0 +H s,1 with H s,0 = (ν−ν)a † a+(ω− ω)σ z /2 and considering two interaction terms (i.e.j = 1, 2) such that δ 1,2 = ∓nν − ω ± nν (note that, as previously explained, α j and δ j are related through α j = (ω + δ j )t) with ω > 0 and ν > 0, we find that H I s,1 = e itHs,0 H s,1 e −itHs,0 approximately leads to with φ n = nπ/2 and g n = η n Ω/(2 n!).The validity of Eq. ( 3) is ensured when Ω ν, |ω + nν| nν together with |η| (a + a † ) 2  1 to safely perform a rotating wave approximation (RWA) [34].It is worth noting that the simulated nQRM can be brought into strong or ultrastrong coupling regimes as the parameters ω and ν can be tuned to frequencies comparable to g n .
In this manner, having access to H gQRM that includes only a linear spin-boson interaction, enables the exploration of a nQRM with non-linear spin-boson coupling (n > 1), whose physics is fundamentally different.For example, the most exotic hallmarks of the n = 2 counterpart of QRM, the so-called two-photon QRM (2QRM), reside in that the spectrum becomes a continuum for g 2 = ν/2 regardless of ω, and for g 2 > ν/2 the Hamiltonian is not longer lower bounded [19,20,35,36].The gQRM lacks of these features, and it is therefore not obvious that the physics of H nQRM can be accessed from H gQRM .Moreover, the H gQRM allows to simulate more exotic scenarios like a combined nQRM and mQRM, just considering enough interaction terms with appropriate values for δ j and α j [34].
Based on the previous transformations one can find the following expression among operators that establishes a relation between the gQRM and nQRM dynamics [34]: Here, U gQRM and U nQRM are the propagators of the gQRM and nQRM respectively, Γ(t) = U † s,0 T † (iη/2)U T,0 with U s,0 = e −itHs,0 , and the approximate character of Eq. ( 4) is only a consequence of the RWA performed to achieve H nQRM from H s .Note that the time dependence of the time-evolution propagators is not explicitly written.Hence, an initial state |ψ nQRM (0) after an evolution time t under H nQRM can be approximated as |ψ nQRM (t) ≈ Γ(t) |ψ gQRM (t) with the initial state |ψ gQRM (0) = T (iη/2) |ψ nQRM (0) .Remarkably, while the speed of the gQRM evolution is lead by the largest frequency ν, see Eq. ( 2), in the simulated nQRM the involved parameters are much smaller as it can be seen in Eq. ( 3) together with the previously commented conditions Ω ν, |ω + nν| nν, and g n = η n Ω/(2 n!).As a consequence, a long evolution time of gQRM is required to effectively reconstruct the dynamics of nQRM.Finally, our theory is completed with a mapping for the observables.As it can be derived from Eq. ( 4), the expectation value of an observable O nQRM , i.e. an observable of the nQRM, corresponds to evaluate O gQRM = Γ † (t)O nQRM Γ(t) in the gQRM.Because Γ(t) involves bosonic displacement and spin rotations, O gQRM may be in general intricate.Yet, for two relevant observables in nQRM, σ z and a † a, the mapping leads to simple operators, namely, σ z transforms into −σ x and a † a into a † a − η/2pσ x + η 2 /4 [34].Interestingly, it still possible to obtain good approximations for other observables by truncating bosonic operators.In this manner one can obtain the full qubit dynamics of nQRM by simply having access to σ x,y,z in gQRM.Indeed, e −η 2 /2 [σ z,y cos((ω + δ 1 )t) ∓ σ y,z sin((ω + δ 1 )t)] turns to be a good approximation of σ x,y in the gQRM frame [34].
To numerically confirm our approximate equivalence, in Fig. 2 we show the results of the simulated dynamics of a 2QRM and a 3QRM using a gQRM for a certain set of parameters and initial states |ψ(0 In addition, in Fig. 2(c) we show that the targeted σ x of a nQRM is retrieved by means of the previously mentioned bosonic truncation of σ x in the gQRM frame, i.e., e −η 2 /2 [σ z cos((ω + δ 1 )t) − σ y sin((ω + δ 1 )t)].Furthermore, in order to quantify the agreement among these models and the validity of the previous theory, we compute the fidelity between the ideal quantum state of the nQRM and the approximated state evolved in the gQRM and properly transformed with Γ(t), that is, F g,n (t) = ψ gQRM (t)| Γ † (t) |ψ nQRM (t) .The computed fidelities of the considered cases are well above 0.99, showing the good agreement among these two models.
Application for microwave driven ions.-Theproofof-concept of our method can be illustrated in a microwave driven ions platform.Note that the developed theory may be relevant in other quantum systems as circuit QED [9].The physics of a single microwave-driven trapped ion in presence of a magnetic field gradient is described by where ω is the qubit energy splitting with a value that depends on the ion species.For example, for the 171 Yb + ion, we have ω ≈ 12.4 GHz [37] plus a factor γB z with γ ≈ 1.4 MHz/G that depends on the applied static magnetic field B z .The coupling parameter ∆ determines the rate of the spin-boson coupling, while the last term corresponds to the action of microwave radiation on the system [38].In this setup the spin-boson coupling is restricted to be linear, see Eq. ( 5), and therefore, our theory appears as an alternative to introduce higher-order boson couplings in the dynamics.In order to take Eq. ( 5) into the form of Eq. ( 2), and subsequently (via the mapping T ) into the general expression in Eq. ( 1), we define ω = δ 1 + ω and move to a rotating frame with respect to the term ω 2 σ z .Considering two drivings such that ϕ 1,2 = π, ω 1 = ω and ω 2 = ω − (δ 2 − δ 1 ) and after eliminating terms that rotate at frequencies on the order of GHz, we find which equals H gQRM after a basis change on both the qubit and boson degrees of freedom, that is, e −i π 4 σy e −i π 2 a † a H I MW e i π 2 a † a e i π 4 σy = H gQRM , where H gQRM is given in Eq. ( 2) with η = 2∆/ν.Hence, it is possible to use a microwave-driven ion to simulate models with non-linear spin-boson couplings (see [34] for more details concerning the implementation in this setup).
Approximate analytical solution for the QRM.-Finding a solution to the QRM has been subject of a longstanding debate, which still attracts considerable attention [13,35,39,40].Based on our theory, we obtain a simple expression for the time-evolution propagator and expectation values of the QRM.In particular, the general expression given in Eq. ( 2) adopts the form of a standard QRM with a unique driving and δ 1 = 0, which, applying our method, approximately corresponds to we obtain now U † s,0 H s,1 U s,0 ≈ H aux instead of H nQRM , and where fast oscillating terms have been neglected performing a RWA, requiring again |η| (a + a † ) 2  1, and only considering resonant terms up to η 2 [34].As a consequence, the following analysis does not apply to the deep-strong coupling regime [14], found here when η ≥ 2. Therefore, the propagator for the QRM can be approximated as which is expected to hold even in the ultrastrong coupling regime of the QRM, although restricted to the condition Ω ν (see later for a numerical confirmation of regimes for Ω and ν in which our approach provides with insightful results).Because H aux has a simple form, the time evolution can be analytically solved, with an initial state |ψ aux (0) = T † (iη/2) |ψ QRM (0) .Indeed, x and E ± n = ±Ω/2(1−η 2 (n+1/2)).Now, employing the map between the two models, Eq. ( 8), we obtain the relation between observables.For example, a † a in the QRM translates to a † a + η 2 /4 + η/2(xσ z sin νt − pσ z cos νt) in H aux [34].In addition, we show that our method improves the typical Bloch-Siegert approximation of the QRM [30,31], e −S H QRM e S ≈ H BS , with H BS = (ν + gΛσ z )a † a+(Ω/2+ ξ)σ z − g(ia † σ − − iaσ + ), where the anti-Hermitian operator is given by S = iΛ(a † σ + + aσ − ) − ξ(a 2 − (a † ) 2 ), with parameters Λ = g/(ν + Ω), ξ = g2 /(2ν(ν + Ω)) and g = ην/2.In order to quantify the enhancement for different values of the coupling constant, we compute the overlap between time-evolved states for both approaches and the QRM, denoted by F QRM,aux and F QRM,BS .As shown in Fig. 3, the approximate solution reproduces correctly the time evolution of the QRM as the cou-pling enters in the non-perturbative ultrastrong regime, g/ν = 0.2 [31] with a fidelity F QRM,aux > 0.99, while the Bloch-Siegert approximation H BS fails as its fidelity drops to F QRM,BS ≈ 0.7.For smaller couplings both approaches lead to similar high fidelities (see Fig. 3(a)).
Conclusions.-Wepresented an approximate equivalence among a family of Hamiltonians, which includes the well-known QRM as well as its higher order counterparts such as the two-photon QRM.In particular, the generalized QRM allows to retrieve the nQRM dynamics with very high fidelities.This theoretical framework shows that these nQRMs can be accessed even in the absence of the required non-linear spin-boson exchange terms, as illustrated with a microwave-driven trapped ion.Moreover, we derived an approximate solution to the dynamics of the QRM even in the ultrastrong coupling regime.
This Supplemental Material: Equivalence Among Generalized nth Order Quantum Rabi Models

I. Transformation between Hs and HgQRM, and HnQRM
In this part we show the details of the unitary transformation between H s , given in the Eq. ( 1) of the main text, and the generalized quantum Rabi model (gQRM), H gQRM in Eq. ( 2).Later, we will discuss the approximate transformation that leads to the H nQRM .Let us start with the transformed Hamiltonian H T = T (iη/2)H s T † (iη/2), which adopts the following form In a straightforward manner, this procedure could be extended to a more general case, that is, when Ω j depends on j.However, for simplicity, we constrain ourselves here to the case in which Ω j ≡ Ω ∀j.As the last term in Eq. ( S2) is just a constant shift we will neglect it for the next developments.Then, we move to an interaction picture with respect to H T,0 = −(ω + δ 1 )/2 σ x .Defining H T = H T,0 + H T,1 , the resulting Hamiltonian H I T,1 reads with t = t − t 0 and U T,0 (t, t 0 ) the time-evolution operator associated to H T,0 .This readily leads to the Eq. ( 2) of the main text once the time-dependent phases are introduced, α j = ωt + δ j t (considering t 0 = 0) and after eliminating terms that rotate at a speed ∝ ω.
On the other hand, H s given in Eq. ( 1) can lead to the desired nQRM.For that, we move to an interaction picture with respect to H s,0 = (ν − ν)a † a + (ω − ω)/2σ z with H s = H s,0 + H s,1 .Then, the interacting part of H s can be written as with a(t) = ae −i(ν−ν)t and a † (t) = a † e i(ν−ν)t and U s,0 (t, t 0 ) the time-evolution operator associated to H s,0 .Then, expanding the exponential, considering that Ω ν and |η| (a + a † ) 2  1, and that δ 1,2 = ∓nν − ω ± nν with |ω + nν| nν, one can perform a rotating wave approximation just keeping those terms resonant with σ + a n and σ + a n .In general, with g n = η n Ω/(2 n!) and φ n = nπ/2.Hence, it is possible to achieve a H nQRM , although the corresponding attained coupling g n becomes smaller for increasing n, as it is proportional to η n /n!.In particular, for n = 2, H I s,1 can be approximated as Note that, while the Hamiltonians H s and H gQRM are related through a unitary transformation, the achievement of a n-photon QRM, H nQRM , from H s requires of certain relations between parameters, such as Ω ν, |ω + nν| nν and |η| (a + a † ) 2  1 to safely perform the rotating wave approximation.Because the frequencies ω and ν can be tuned, different combinations of Ω, ν, ν and η in H gQRM can lead to the same simulated nQRM.Therefore, it is important to recognize the main contribution that deteriorates the established equivalence as it may be overcome by properly tuning these free parameters.Moreover, besides the chosen parameters, the simulation of the dynamics of the 2QRM when a large number of bosonic excitations is involved is expected to breakdown as |η| (a + a † ) 2  1 is not longer satisfied.This is indeed the case for the 2QRM right at the spectral collapse as a † a blows up, however, the onset of the dynamics can be still reproduced.In Fig. S1 we show a comparison of the simulated 2QRM with different parameters and for g 2 /ν = 1/2 (at the spectral collapse) in (a) and (c), and g 2 /ν = 0.125 in (b) and (d), as shown in Fig. 2 of the main text.Note that decreasing ν the real evolution time becomes longer, although and at the same time, it leads to smaller values of η (and/or Ω/ν).For example, in Fig. S1(d) we observe that the fidelity when Ω/ν = 0.05 and ν/ν = 2.5 × 10 −4 slightly improves that of Ω/ν = 0.1 and ν/ν = 5 × 10 −4 , which together with the worst shown case (Ω/ν = 0.2), share the same value of η, namely, η = 0.05.This indicates that the main spurious contribution stems from the zeroth order in η, i.e., Ω/2σ + e ±2(ν−ν)t , as higher orders become smaller, ∝ η n Ω/(2 n!), while rotating at approximately equal frequency, ≈ ν.In addition, we explicitly show that the presented results are not affected by the Fock-space truncation; in particular, for the spectral collapse, the results do not change when doubling the number of Fock states, from N max = 100 to 200 (see Fig. S1(a) and (c)).

II. Combined nth and mth order QRM models from HgQRM
As stated in the main text, the approximately equivalence is not restricted to H nQRM and H gQRM .Indeed, more complex Hamiltonians than nQRM can be attained by a suitable H gQRM .Here we show how a Hamiltonian that comprises interaction terms of that of a nQRM and a mQRM can be accessed from H gQRM .This Hamiltonian is denoted here by H n,m and reads We first show how to achieve H n,m from H s (see Eq. ( 1)), following the same procedure as shown for H nQRM .For that, four interaction terms are now needed in Eq. ( 1), with α j = ωt + δ j t and δ j = ∓nν − ω ± nν for j = 1, 2 and δ j = ∓mν − ω ± mν for j = 3, 4.Then, assuming that Ω ν, |ω + nν| nν, |ω + mν| mν as well as |η| a + a † 1, H I s,1 approximately corresponds to H n,m , where we have dropped out non-resonant terms performing a RWA, i. e., terms rotating at frequencies ∼ ν have been neglected.Note that H I s,1 denotes H s,1 in the interaction picture with respect to H s,0 with H s,0 = (ν − ν)a † a + (ω − ω)σ z /2 and H s = H s,0 + H s,1 , as explained in the main text and in the previous section.In addition, the attained phases are φ k = kπ/2, while the couplings g k = η k Ω/(2 k!).It is worth mentioning that in this case one would gain tunability in the couplings by considering different frequencies Ω for each j; for example, one could achieve similar couplings g n ∼ g m with n = m.
On the other hand, following the procedure explained in the main text, we can bring H s into H gQRM , Eq. ( 2), by applying a unitary transformation.For the particular case considered here, i.e., four drivings with α j = ωt + δ j t and δ j = ∓nν − ω ± nν for j = 1, 2 and δ j = ∓mν − ω ± mν for j = 3, 4, H gQRM adopts the following form Finally, the time-evolution propagators of both models are related as given in Eq. ( 4), that is, where Γ(t) is defined as in the main text, Γ(t) = U † s,0 T † (iη/2)U T,0 .Therefore, the map between initial states and observables is identical as explained in the main text for the approximate equivalence between a nQRM and a gQRM.Therefore, the observables of H n,m transform in the same manner to the frame of gQRM.gQRM is not zero due to the approximate character of the equivalence.

III. Equivalence of observables and states
Here we show the derivation of the Eq. ( 4) of the main text.Having established the transformations that connect H gQRM with H s , and H nQRM with H s we can relate them in terms of the time-evolution operators, where U I refers to the time-evolution of H in an interaction picture and we have dropped the explicit time dependence for the sake of readability.Then, combining the Eqs.(S11), (S12) and (S13), we arrive to which is the Eq. ( 4) of the main text, with the relation between initial states |ψ gQRM (0) = T (iη/2) |ψ nQRM (0) .Finally, from Eq. (S15) it is straightforward to obtain the observable that must be measured in the gQRM frame in order to retrieve O nQRM of the nQRM, i.e., and thus, for O nQRM = σ z and a † a the transformation leads to while for other observables, like σ x and σ y , a more intricate expression is attained, be tuned to ≈ 9.25 kHz which is achievable with a magnetic field gradient smaller than 150 T m [S2].Note however that an evolution time of 20 ms is a rather long time to preserve the coherence of both qubit and bosonic mode from the inevitable presence of environmental noise sources, even when techniques to cope with noise are applied (see [S3] where a coherence time of ∼ 10 ms is measured).In this regard, we stress that although the results presented in Fig. 2 will be deteriorated by loss of coherence, a faithful simulation of nQRMs with H MW is still feasible at shorter times.Furthermore, depending on the specific platform, the parameters may be optimized to avoid spurious decoherence processes or be combined with techniques to extend quantum coherence as dynamical decoupling [S4].Finally, it is worth emphasizing that we have considered a microwave driven ion just to illustrate a direct application of our theory, as a proof of concept.

V. Relation between HQRM and Haux
The Eq. ( 2) of the main text reduces to a QRM by simply considering a driving with δ 1 = 0, that is, Making use of the derived approximate equivalence, the previous Hamiltonian can be approximately mapped into a simple Hamiltonian, H aux = Ω 2 1 − η 2 (a † a + 1/2) .This is accomplished by moving H s = H s,0 + H s,1 to an interaction picture with respect to H s,0 = νa † a + ωσ/2.Note that now ω = ν = 0 and recall that α 1 = (ω + δ 1 )t = ωt.Therefore, σ + e iη(a(t)+a † (t)) + H.c. , where a(t) = ae −iνt and a † (t) = a † e iνt .Requiring now |η| (a + a † ) 2 1 we expand the exponential, and assuming Ω ν, we can safely perform a RWA neglecting off-resonant terms (which rotate at frequencies larger or equal than ν) and keeping only the resonant terms up to η 2 .The following higher-order resonant term appears with η 4 .Hence, we obtain the relation given in the main text, In a straightforward manner as done for nQRM and gQRM, we can obtain the relation between the propagators, U QRM , U aux and U s,0 , U QRM ≈ U † T,0 T (iη/2)U s,0 U aux T † (iη/2), (S32) Note that this relation also follows from Eq. ( 4) considering now that gQRM reduces to QRM and H aux replaces H nQRM , and with H s,0 = νa † a + ωσ/2 since ω = ν = 0.

VI. Expectation values in the approximate QRM
As indicated in the main text and from Eq. ( 8), U QRM ≈ U † T,0 T (iη/2)U s,0 U aux T † (iη/2), (S33) which relates the time evolution of a QRM with the simple evolution in the Hamiltonian H aux , one can obtain the map between observables and initial states.In particular, we have U T,0 = e −it(−ω/2σx) and U s,0 = e −it(ω/2σz+νa † a) (S34) where we have set already ω = ν = 0 in H s,0 , i.e., H s,0 = νa † a + ωσ z /2.In addition, T (β) denotes the unitary transformation given in the main text, that is, T
work was supported by the ERC Synergy grant BioQ, the EU STREP project EQUAM and the CRC TRR21.The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 40/467-1 FUGG.J. C. acknowledges Universität Ulm for a Forschungsbonus.H. M.-C.thanks the Alexander von Humboldt Foundation for support.J. C. and R. P. contributed equally to this work.