Tripartite quantum Rabi model with trapped Rydberg ions

We investigate a tripartite quantum Rabi model (TQRM) wherein a bosonic mode concurrently couples to two spin-$1/2$ particles through a spin-spin interaction, resulting in a spin-spin-boson coupling -- a departure from conventional quantum Rabi models featuring bipartite spin-boson couplings. The symmetries of the TQRM depend on the detuning parameter, representing the energy difference between the spin states. At zero detuning a parity symmetry renders the TQRM reducible to a quantum Rabi model. A subradiant to superradiant transition in the groundstate is predicted as the tripartite coupling strength increases. For non-zero detuning the total spin emerges as the sole conserved quantity in the TQRM. It is found that superradiance prevails in the groundstate as long as the tripartite coupling remains non-zero. We derive the Braak $\mathcal{G}$-function of the TQRM analytically, with which the eigenspectra are obtained. The TQRM can be realized in a viable trapped Rydberg ion quantum simulator, where the required tripartite couplings and single body interactions in the TQRM are naturally present. Our study opens opportunities to explore and create novel correlations and entanglement in the spin and motional degrees of freedoms with the TQRM.


I. INTRODUCTION
The quantum Rabi model (QRM) consists of a spin-1/2 particle coupled to a bosonic degree of freedom [1].Despite its simplicity the QRM exhibits rich physics [2], and finds applications in, e.g., benchmarking quantum computers [3], verifying the existence of supersymmetric quantum mechanics [4,5], investigating PT -symmetry breaking [6,7], and the generation of nonclassical states [8,9].An important feature of the QRM is that the validity of the rotating wave approximation breaks down as a result of the strong spin-boson coupling [10][11][12].This causes difficulties in obtaining the spectrum of the QRM analytically.Such a task was not accomplished until 2011 by Braak [13] through the G-function method, although the isolated doubly degenerate energies that represent the exceptional spectrum were analytically determined decades earlier [14].After this achievement, many studies have focused on understanding the mathematical structure and improving the numerical stability and performance of the G-function [15][16][17][18][19][20].When the coupling is strong the groundstate of the QRM exhibits a superradiant-like phase transition, where the population of the boson becomes sizable when increasing the coupling above a critical value [21][22][23][24].Motivated by the rich physics, the QRM has been experimentally realized with various physical systems [20], including trapped ions [25,26], cavity and circuit QED [27][28][29][30][31][32][33], cavity optomechanics [34], and quantum dots [35,36].Furthermore, superradiant transition has been experimentally demonstrated with a single trapped ion [26].
Modified QRMs have been extensively explored within the literature [20].The first type of extension concerns the modification of the bosonic operators.For example two photon QRMs, where the spin couples to the bosonic mode via a two photon process, have attracted a lot of attention [19,[37][38][39][40][41][42].The symmetry extends from a Z 2 group in the QRM to a Z 4 group in the two photon QRM [38].Despite the change of symmetry the spectra of the two photon QRM can still be found analytically [43,44].Another way to extend the QRM is to include more spins in the model, such as two qubit QRMs (TQQRMs) [45][46][47][48][49][50][51][52][53][54][55][56].Such extension leads to the Dicke model when number of spins is large [57].However, all of these existing studies have focused on bipartite couplings which couple individual spins with bosonic modes.
In this article we study an exotic tripartite quantum Rabi model (TQRM) where two identical spins couple to a monochromatic bosonic mode simultaneously through an Ising spin-spin interaction.To the best of our knowledge, such extension has not been studied before.Three-and multibody interactions have attracted a broad interest in the study of nuclear, atomic and many-body physics [58].Recent experimental and theoretical studies have realized multi-body interactions with trapped ions [59], neutral atoms [60], and superconducting circuits [61,62].A key feature of the TQRM is that the constituents of the TQRM are a bosonic mode and two identical spins.We show that the TQRM can be realized with two trapped ions in their highly excited Rydberg electronic states [63,64].In a linear Paul trap long-range dipole-dipole interactions between Rydberg ions couple to the breathing mode of a two-ion crystal [65], leading to the novel tripartite coupling between the two spins and bosonic mode (crystal phonons).The TQRM is achieved through modulating the Rydberg excitation and phonon-ion coupling with external laser fields.
The symmetries of the TQRM are controlled by the Rydberg excitation laser.At vanishing detuning the TQRM has a well defined parity symmetry that relies on the total population of the spins and phonon.When the detuning is finite, the parity symmetry disappears.We derive analytically the G-function that defines the eigenspectrum of the TQRM, and show that these analytical results indeed match those of numerical calculations.We then investigate the groundstate properties, focusing on the strong coupling regime.At zero detuning the TQRM reduces to an effective QRM, which differs from the QRM insofar that it describes a collective interaction between the spins and the bosonic mode.In the superradiant phase the spatial distribution of the phonon splits into two peaks symmetrically, and its state is captured by a classical mixture of two coherent states |±⟩ ( > 0).At finite detuning the superradiance gradually dominates the groundstate.It is found that the phonon spatial distribution has a single peak, which is shifted from the origin in the strong coupling regime.A careful examination shows that the phonon state is approximately described by coherent state |−⟩.
This article is organized as follows.In Section II we present the trapped ion system used as the quantum simulator and derive the TQRM Hamiltonian.In Sec.III we discuss the symmetries present in the model and how the single spin detuning affects the symmetry.The G-function is derived analytically in Sec.IV.In Sec.V we discuss superradiance and the phase space distribution function of phonons in the groundstate of the TQRM.We conclude in Sec.VI.
In our setting two ions are trapped in a linear Paul trap, as depicted in Fig. 1(a).Due to the harmonic trapping potential (axial trapping frequency ) and Coulomb repulsion the two ions form a Wigner crystal along the trap axis (-axis) at low temperatures, where coordinates of the -th ion are   ,   , and   (  = 1, 2).The ions vibrate around their equilibrium positions where N  is the net charge,  is the mass of the ion, and  0 the permittivity of free space.As we will show later, the electronic states (spin) will couple to the axial vibration, while the radial motions will not be involved.Along the trap axis we obtain a center-of-mass (c.m.) and breathing modes whose frequencies are  c.m. =  and  = √ 3, respectively [65].A Rydberg excitation laser couples the low energy metastable state |↓⟩ = |4 5/2 ⟩ and a Rydberg state |↑⟩ (e.g.| 1/2 ⟩ with  to be the principal quantum number) coherently [63].At the same time an off-resonant standing wave laser induces a state-dependent Stark shift   [84,85].The Hamiltonian of The dipolar interaction is induced through using the microwave dressed Rydberg states [74].In this example, we have considered 88

√
2 and  b = 1/ √ 2, respectively.For typical trap parameters, the length is in the order of a few tens of nanometers, while the distance between the two ions is several micrometers [64,65].This allows us to expand  d ( 2 − 1 ) around their equilibrium positions.To the linear order of   , we obtain where  d ( 0 ) and  ′ d ( 0 ) are the potential and its slope at the equilibrium distance  0 =  (0) 2 −  (0) 1 .  is the small deviation from the equilibrium position  (0)  .By modulating the Rydberg states with external microwave fields we can tune  d ( 0 ) = 0, i.e. the dipolar interaction vanishes at the equilibrium distance as shown in Refs.[73,74].See Fig. 1(c) and (d) for an example.As a result, the dipolar interaction couples directly to the breathing mode through gives the coupling strength.Note that the c.m. mode decouples with the electronic dynamics.We will exclusively focus on the coupled dynamics between the spin and breathing mode from now on.
With the above consideration one obtains ) (+ † ), with  = /4.Our aim is to achieve a spinspin-phonon coupled interaction ∝   1   2 ( +  † ).To achieve this, one can turn off the coupling between the individual spin and phonon through spin-dependent Stark shift [84] with The detail of the implementation can be found in the Appendix.The only remaining term is a collective coupling between the two spins and phonon, i.e. the breathing mode couples to the two spins simultaneously whose strength depends on the slope of the two-body dipolar interaction.This yields the Hamiltonian (1) Equation ( 1) represents the tripartite quantum Rabi model (TQRM) consisting of a monochromatic bosonic (phonon) mode of frequency  and two identical spins.The unique aspect of the TQRM compared to the QRM and its variants previously studied is that our model is characterized by a tripartite coupling, i.e. a spin-spin-boson coupling.The tripartite coupling is different from existing models.In conventional multi-spin QRMs the two-body interactions between spins are typically not present (see Refs. [45][46][47][48][49][50][51][52][53][54]), or the spin-spin interactions do not directly couple to the bosonic mode [55,56].

III. SYMMETRIES OF THE TQRM
We begin by examining symmetries in the TQRM and how the single body detuning affects the symmetries.It is convenient to rotate the Hamiltonian (1) around the    -axes by /2, yielding In the rotated basis the spin up and spin down states will be denoted with |⇑⟩ and |⇓⟩, respectively.After the rotation both the detuning and the tripartite coupling term of  R can induce transitions between spin states.In the following, we first analyze symmetries in the TQRM when the Rydberg excitation laser is resonant ( = 0), where transition between states is solely driven by the tripartite coupling.Next we consider finite detuning ( ≠ 0), where both terms will be taken into account.

A. Resonant TQRM
When  = 0, the spin states of the two ions are only coupled through the tripartite interaction, i.e.
Thus TQRM is equivalent to the conventional QRM in this sector.However, due to the collective nature the effective level separation becomes 4Ω, and the phonon mode couples electronic states of the two spins.In this regime  s exhibits a parity symmetry with parity operator Π p = exp( E ) [13], with the total excitation  E =  †  +   .Its eigenvalues are 1 (−1) when  E is an even (odd) integer.It can be shown that Π p  s Π † p =  s .Such parity, as depicted in Fig. 2(a), requires that when spin states change, phonon state |⟩ changes to | ± 1⟩.

IV. ENERGY SPECTRUM
Eigenspectra of many variants of QRMs have been obtained analytically [20], typically through the Braak G-function.The G-functions are derived either through a transformation into the Segal-Bargmann space of complex analytic functions [87], or through the Bogoliubov operator approach (BOA) introduced by Chen et al. [44].We will employ the BOA as it generates algebraically simpler G-functions than the Segal-Bargmann approach [44,48,51].When  = 0 our model reduces to an effective QRM, whose analytical spectrum is known [13].In the following we will focus on the spectrum of Hamiltonian (1) in the triplet sector for  ≠ 0.
i.e., |0  ⟩ is a coherent state, where the sign is − (+) when  =  ( = ).The relation between the coherent states |0  ⟩ [88] and |0⟩ is important in a later step of our derivation Using the Bogoliubov operators  and  † , and in the triplet manifold, the Hamiltonian  can be written as where , where   ,   , and   are the expansion coefficients (proportional to probability amplitudes), respectively.Substituting this expansion into the Schrödinger equation for eigenenergy .Multiplying by the state ⟨  | from the left hand side, we derive coupled equations of the coefficients, From Eq. ( 5) and ( 7) one finds expressions of   and   with respect to the common coefficient   .Substituting these into Eq.( 6) yields a recursion relation for   , where the coefficient C  is defined as , with  =  −  +  2 /.Using the initial coefficients  0 = 1 and  1 = C 0 , the coefficients   can be evaluated iteratively when parameters {, , , Ω, } are specified.
If  is a nondegenerate eigenenergy of  the states |Ψ  ⟩ and |Ψ  ⟩ must differ only by a complex coefficient , i.e. |Ψ  ⟩ =  |Ψ  ⟩.Left-multiplying by the vacuum bra-state ⟨0| on both sides, we obtain where we have used relation [44].Cross-multiplying to eliminate the arbitrary constant  an analytical expression for the Braak Gfunction is found, The roots of the G-function are used to evaluate eigenenergies of Hamiltonian (1).
Examples of the function G() are plotted in Fig. 3(a) and (c).Unlike the TQQRM [48,51], there are no regular pole structures in our model.This results exclusively from the detuning term that breaks the symmetry in the QRM [13,91].Roots of the G() function are evaluated numerically.In Fig. 3(b) eigenvalues of  as a function / are shown.The eigenenergy obtained from the G() function agrees with the numerical diagonalization of the Hamiltonian.Anticrossings of the spectra are found at finite .However the groundstate energy clearly separates from the first excited state.The energy gap between the ground and first excited state remains finite when / > 1. Increasing  the gap increases too, as can be seen in Fig. 3(b) and (d).

V. SUBRADIANT AND SUPERRADIANT PHASES IN THE GROUNDSTATE
In this section we investigate properties of the groundstate of the TQRM.Through mean-field calculations and numerical diagonalization of the full Hamiltonian we will show that a subradiant to superradiant transition can be induced by increasing  in the resonant TQRM.Phase space densities of the phonon state are symmetric but highly delocalized in the superradiant regime.The phonon state can be described by a classical mixture of coherent states |±⟩.In the detuned TQRM we will show that the phonon state is a displaced oscillator.The subradiance disappears, as the mean phonon number is non-zero as long as  > 0 in the detuned TQRM.

A. Delocalized superradiant phase of the resonant TQRM
The groundstate exhibits a subradiant to superradiant phase transition when  is larger than a critical value.We first determine the critical value within the mean-field approach.This is done by replacing operators  ( † ) with their mean value  ( * ) in Hamiltonian Eq. (1).By assuming that  is real we diagonalize the Hamiltonian, and the lowest energy obtained as a function of .In the groundstate E G / = 0, which allows the definition of a critical coupling  c = √ Ω.When  ≤  c we find  = 0.The groundstate is subradiant with zero phonon occupation.When  >  c , two non-zero solutions,  = ± 0 with , are obtained indicating superradiance with a finite phonon occupation.As shown in Fig. 4(a), E G is an even function of .Two minima are found when  >   , consistent with the above analysis.
At the mean-field level this superradiant transition is characteristic of a spontaneous symmetry breaking occurring related to a second order transition [21,92].Here the energy (Fig. 4(b)) and its first derivative (Fig. 4(c)) are continuous when varying .However the second derivative of the groundstate energy shows a discontinuity at the critical coupling, as seen in Fig. 4(d).Both the mean-field and the exact diagonal- ization calculation agree well in the subradiant and superradiant phases.Around the critical point  =  c , the full quantum model shows a crossover, which is different from the sharp transition in the mean-field calculation.
We now investigate the phonon states in the subradiant and superradiant phase.We are particularly interested in the phonon properties in the position representation.The position and momentum are given through phonon operators  =  b ( † + ) and  = i  b ( † − ), respectively.We first numerically evaluate the reduced density matrix of the phonon  b = tr[ t ], by tracing the spin degree of freedom from the total density matrix  t = | G ⟩⟨ G |.This allows us to obtain spatial density distribution  b () = ⟨| b |⟩ in the position representation.As shown in Figs.5(a) the spatial density has a Gaussian profile centered at  = 0 in the subradiant phase.When  ≳  c the density deviates from the Gaussian.It first becomes non-Gaussian around the critical value, and then splits into two separated Gaussian shapes.
The Wigner quasiprobability distribution  (, ) [93], on the other hand, exhibits distinctive features in both phases.As shown in Figs.5(b) and (c) the Wigner distribution is stretched along the -axis when  ∼  c .Two separate peaks are observed when  >  c , depicted in Fig. 5(d)-(e).These peaks are centered around  ≈ ±2 b  0 , similar to the spatial density shown in Fig. 5(a).The Wigner distribution is stretched horizontally (i.e.along the  axis) in the superradiant state.In certain regions, the Wigner function becomes negative, indi-  r [94,95].When  <  c , the fidelity is high when  r = |0⟩⟨0|, as shown in Fig. 6(a).This is a direct manifestation of the subradiance in this regime.When  ∼  c , however, the fidelity becomes low in general when projecting  b to |± 0 ⟩ or their superpositions | ± ⟩ =  ± (|+ 0 ⟩ ± |− 0 ⟩), where ) is the normalization constant.This is consistent with our numerical calculation shown in Fig. 5(c), where the Wigner function exhibits negative values.In this regime, the phonon state is a non-Gaussian state, and hence can not be accurately described by coherent states.
On the other hand, the purity of the phonon density matrix decreases rapidly with increasing  when  ≫  c , as shown in Fig. 6(b).It turns out that the  b becomes a classical mixture When  ≫  c , the corresponding fidelity  approaches unity, as shown in Fig. 6(c).Using  c , we can evaluate the uncertainty of  and  analytically, (Δ) 2 = (4 2 +1)/2 2 and (Δ) As  is a real number in the mean-field calculation, this leads to (Δ) 2 =  2 /2.The uncertainty (Δ) 2 is therefore same with that of the coherent state |± 0 ⟩.In other words the phonon state is strongly stretched along the -axis, while the uncertainty along -axis remains a constant.Numerical diagonalization of the full Hamiltonian shows that (Δ) 2 only deviates from  2 /2 around  ∼  c , as can be seen in Fig. 6(d).We first investigate the groundstate of the detuned TQRM with the mean-field approach.In Fig. 7(a) the parameter  as a function of  is shown.When increasing  > 0 and restricting  > 0 we find  is always negative and its magnitude || increases monotonically.This is different from the resonant TQRM where both positive and negative branch of  are found.Another important difference is that the sharp change of  at the subradiant-superradiant transition disappears when  > 0. Instead a smooth crossover emerges around  ∼  c .We then calculate the mean phonon population, shown in Fig. 7(b).It is found that the crossover at  ∼  c persists.The numerical diagonalization and mean-field results deviate apparently when  is small.This trend changes when  is large, where the mean-field result agrees with the numerical calculation and are largely independent of .Another feature is that when  is large the phonon number is non zero as long as  > 0, indicating that subradiance does not exist any more.In other words the subradiant to superradiant transition is removed by finite .
The negativity of  can be understood through the meanfield analysis.When  ≫  and  ≫ Ω the groundstate energy functional is well approximated by When plotting E G as a function of  [see Figs.7(c)-(d)] it is apparent that the minima locate at  < 0. The value of  can be determined through solving E G / = 0.One finds groundstate energy  G ≈ − 2 / − 2 when  1 = −/.This explains qualitatively why  becomes negative in the groundstate.
The c.m. of the phonon density is shifted from  = 0 (when  = 0) towards  < 0 (when  ≠ 0).Our numerical simulations show that the center of the spatial density [Fig.8(a)] and  In this regime the reduced density matrix of the phonon is well represented by the coherent state | 1 ⟩.The Uhlmann-Jozsa fidelity between  b and coherent state | 1 ⟩ is shown in Fig. 9(a).The corresponding fidelity increases when increasing .As a result the phonon state population is well approximated by  2  1 , which is non-negligible even when  <  c .This means that subradiance becomes impossible when  ≠ 0, which is consistent with the numerical data shown in Fig. 7(b).The fidelity becomes relatively low around  ∼  c .Importantly the fidelity already approaches unity when / c ≈ 2 for all  shown in the figure.Due to the high fidelity we find that the purity of the phonon density matrix approaches unity, as shown in Fig. 9(b).As a result  b approaches to a pure state when  ≠ 0 and  >  c .This is in sharp contrast to the resonant TQRM, where  b is a mixed state in the superradiant regime.

VI. DISCUSSION AND CONCLUSION
In this work, we have shown that the TQRM can be realized with a pair of trapped Rydberg ions, in conjunction with laser induced spin-dependent force.In this setting the single-body terms in the Hamiltonian can be controlled by tuning the laser and Paul trap electric field.In typical experiments the detuning, Rabi frequency, and trap frequency can be varied flexibly from kHz to MHz.The challenging part is to control the tripartite coupling strength , e.g. to probe subradiant and superradiant phases.In the Rydberg ion system, / c =  ′ d ( 0 )/4 √ 2Ω, which depends on the trap frequency implicitly through the equilibrium positions.The dipolar interaction profile is realized by using the microwave dressing scheme discussed in Ref. [74].In this method, coherent coupling of multiple Rydberg states of different parities can shift the dipolar interaction globally, such that the interaction strength is zero at the twoion equilibrium distance, leading to strong potential gradients that couple the two ionic spins and the phonon mode.For example, we can use same parameters given in Ref. [74] to couple states |50⟩ and |50⟩ of 88 Sr + ions.Considering trap frequency  = 2 × 1 MHz, and Ω = 2 × 25 kHz, we obtain / c ≈ 0.13, which falls in the weak coupling (subradiant) regime.To access the strong coupling (superradiant) regime one can choose higher Rydberg states with  = 80, and tighter trap with  = 2 × 2.02 MHz, which leads to / c ≈ 1.04.Hence the tunable laser and trap parameters allows us to explore the physics of TQRM, as well as strong and collective spin-phonon coupled dynamics.
In conclusion, we have studied a novel TQRM where the phonon interacts with two spins simultaneously through the tripartite coupling.The Braak -function is derived analytically for the TQRM, which determines the regular eigenspectra of the model.We have analyzed groundstate properties of the TQRM by varying the detuning .In the resonant case ( = 0), the TQRM can be reduced to the conventional QRM.In the subradiant to superradiant transition region, the phonon becomes a non-Gaussian state whose Wigner function exhibits negative values.In the strong coupling (superradiant) regime, it is found that the reduced density matrix of the phonon is a classical mixture of two coherent states |± 0 ⟩ with equal probability.In the detuned TQRM ( ≠ 0), the subradiance disappears for finite  > 0. When  ≫   , we find that  b becomes a pure state | 1 ⟩.This work opens a new way to explore exotic physics through multi-partite couplings between spins and bosons.One can extend the setting to a longer ion chain, where multiple spins can couple to multiple phonon modes through the tripartite coupling [76].Such tripartite coupling permits to create, e.g.hyperentanglement in the spin and phonon degrees of freedom simultaneously, which finds quantum information processing applications [96][97][98][99][100][101][102].

FIG. 1 .
FIG. 1.(a) Two ions are confined in a linear Paul trap with trap frequency .The ions are coherently laser excited to Rydberg states.In their Rydberg states the ions interact through dipolar interaction depending on their separation  0 =  (0) 2 −  (0) 1 .A standing wave laser field (blue arrows) induces a site dependent Stark shift.(b) The Rydberg excitation laser excites the ions from groundstate |↓⟩ to Rydberg state |↑⟩ with Rabi frequency Ω and detuning 2.(c) The dipolar interaction.At  0 = 2.7 µm the interaction  d ( 0 ) = 0 (marked by the dashed perpendicular line).(d) The slope  ′ d of the dipolar interaction. ′ d = −2 × 174.7 MHz/µm when  0 = 2.7 µm.The dipolar interaction is induced through using the microwave dressed Rydberg states[74].In this example, we have considered88 Sr + with principal quantum number  = 80, and  = 2 × 2.02 MHz.See text for details.

FIG. 4 .
FIG. 4. (a) Mean-field energy functional E G for / c = 0.0 (solid blue), 1.1 (dashed red), and 2.0 (dash-dotted purple).The groundstate energy is obtained at the minima of the curves.(b) The groundstate energy, (c) its first and (d) second derivatives, as well as (e) the population of the bosonic mode  †  =  2 with respect to the coupling.Solid (dashed) lines represent the mean-field (quantum) results.In the following calculations, we scale the Hamiltonian with respect to  and set Ω = .

FIG. 6 .𝜌 r 2 between
FIG. 6.(a) The Uhlmann-Jozsa fidelity  and (b) purity of the reduced density matrix  b .In Panel (a) fidelities  of  b and the pure states |± 0 ⟩ (solid blue), | + ⟩ (dashed red), and | − ⟩ (dash-dotted purple) are considered for  >  c .When  <  c , the projection to the vacuum state is examined.In Panel (c) the fidelity between  b and  c is shown.Panels (d) and (e) show the quadrature variances obtained by the mean-field (solid red) and full numerical calculations (dashed blue).

FIG. 7 .
FIG. 7. (a) Mean-field  and (b) phonon population ⟨ † ⟩ (|| 2 ) of the resonant TQRM.In Panel (b) the quantum and mean-field results are denoted by lines and lines with symbols, correspondingly.Energy functional E G for (c) / c = 0.95 and (d) / c = 1.4 are shown.In both the subradiant and superradiant regime, minima of the energy functional locate at negative .We have considered  > 0 in the calculation, where /Ω = 0.2 (solid blue), /Ω = 0.5 (dashed red), and /Ω = 1.0 (dash-dotted purple) in all the panels.
EP/W015641/1 and the University of Nottingham.IL acknowledges funding from the European Union's Horizon Europe research and innovation program under Grant Agreement No. 101046968 (BRISQ).This work was supported by the University of Nottingham and the University of Tübingen's funding as part of the Excellence Strategy of the German Federal and State Governments, in close collaboration with the University of Nottingham.This work is partially funded by the Going Global Partnerships Programme of the British Council (Contract No. IND/CONT/G/22-23/26).The data used to create the figures in this article can be accessed via an online repository.