GHZ-like states in the Qubit-Qudit Rabi Model

We study a Rabi type Hamiltonian system in which a qubit and a d-level quantum system (qudit) are coupled through a common resonator. In the weak and strong coupling limits the spectrum is analysed through suitable perturbative schemes. The analysis show that the presence of the multilevels of the qudit effectively enhance the qubit-qudit interaction. The ground state of the strongly coupled system is a found of Greenberger-Horne-Zeilinger (GHZ) type. Therefore, despite the qubit-qudit strong coupling, the nature of the specific tripartite entanglement of the GHZ state suppress the bipartite entanglement. We analyze the system dynamics under quenching and adiabatic switching of the qubit-resonator and qudit-resonator couplings. In the quench case, we found that the non-adiabatic generations of photons in the resonator is enhanced by the number of levels in the qudit. The adiabatic control represents a possible route for preparation of GHZ states. Our analysis provides relevant information for future studies on coherent state transfer in qubit-qudit systems.


I. INTRODUCTION
The quantum Rabi model (QRM) describes the interaction between a two-level system and a single quantized harmonic oscillator mode. It is one of the most celebrated models in atomic physics for light-matter interaction [1]. In quantum technology, Rabi-like models are widely employed to describe the effective physics emerging in a variety of different contexts ranging from spintronics [2,3] to trapped ions [4], and from circuit quantum electrodynamics (cQED) [5] to atom-superconducting qubit hybrid schemes [6]. Despite its simple form, the Rabi model was solved exactly only recently [7,8]. The ground state of the quantum Rabi model consists of a non-classical highly entangled state of two-level system and bosonic mode [7,9]. In cQED, different regimes of interaction between the two-level system and the bosonic field can be explored. In particular, weak and strong coupling regimes are routinely exploited for read-out and coherent state transfer [10]. Recent research has demonstrated the possibility of reaching the ultrastrong and deep strong coupling regimes, too [11,12].
We shall see that the physics of our model is particularly interesting in the ultrastrong coupling regime [18,[25][26][27][28]. In particular, the ground state of the system turns out to be defining a Greenberger-Horne-Zeilinger (GHZ) entangled state. GHZ states present great significance among all types of multipartite entanglement [29]. These states exhibit maximal correlations between three or more quantum systems. GHZ states have been considered a key resource in fundamental physics since the early stages of quantum information. They have also been proven useful in various quantum technologies, including quantum error-correcting codes [30] and quantum metrology beyond the Heisenberg limit [31].
We point out that the generation of GHZ hybrid entanglement defines a challenging problem in quantum technology. In cQED, for example, GHZ hybrid entanglement has been achieved by state-dependent phase shift operations which involve complicated control and feedback sequences [23,32]. In this context, exploration of the ultrastrong coupling regime has been demonstrated beneficial for GHZ state preparation [33,34]. Indeed, our scheme guarantees a straightforward preparation of hybrid GHZ states, as such state could appear in the ground state of the QRM at large coupling strength.
The article is organized as follows: in Sec. II, we introduce a generalization of the quantum Rabi model to describe the qubit-qudit interaction through the resonator bosonic field. We demonstrate, through an analytical approach based on adiabatic approximation and perturbative expansion that hybrid GHZ states constitute the ground state solutions in the ultrastrong coupling limit. In Sec. III, we study the bipartite entanglement between the qubit and qudit mediated by the common resonator, quantified by negativity [35]. We show that the presence of the GHZ state induces an exponential suppression of the negativity for large values of the coupling strengths. Dynamical features are investigated in Sec. IV. In particular, we show the dynamics after quenching the coupling strength, and propose adiabatic state preparation schemes for the hybrid GHZ states. A short discussion of the main results of the manuscript is presented in Sec. V.

II. THE QUBIT-QUDIT RABI MODEL
We investigate the physical system schematically pictured in Fig. 1. The scheme features a qubit (two-level quantum system) and a qudit (d-level quantum system) individually coupled to a common quantum resonator described by bosonic degrees of freedom. The Hamiltonian reads (here, and in the rest of this manuscript, we work in units = 1) (1) Here,σ x,y,z are the Pauli matrices for the qubit with transition frequency Ω 1 ,Ĵ z,± d are the spin (d − 1)/2 operators with level spacing Ω 2 , andâ(â † ) is the annihilation (creation) operator for the bosonic field with frequency ω. The coupling strengths g 1,2 in Eq. (1) quantify the vacuum-Rabi splittings. Employing a jargon that is widely used in the literature, we will denote the qubit and the qudit as "artificial atoms". For d = 2, our model is equivalent (up to a sign convention) to the two-qubit quantum Rabi model [36][37][38][39]. In contrast to the single qubit Rabi model, this generalized model is not integrable for general parameter values [37,38].
The eigenvalues and eigenstates ofĤ can be readily obtained numerically. Here, we devise analytical approximation schemes both in the weak-coupling and in the ultrastrong coupling regimes.
The weak coupling limit (g 1 , g 2 ≪ ω), in the presence of strong qubit/qudit-resonator detuning (Ω 1 , Ω 2 ≪ ω), can be treated by means of a Schrieffer-Wolff transformation [5,40]. In particular, we apply the following unitary transformation to the Hamiltonian Eq. (1): FIG. 1. Model schematics. The system is composed of a qubit with level spacing Ω1, an harmonic oscillator (resonator) with characteristic frequency ω, and a d-level quantum system (qudit) with level spacing Ω2. The qubit and the qudit are coupled to the resonator through the coupling constants g1,2.
where we choose In the weak coupling limit considered here ǫ i , ξ i ≪ 1.
In particular, at the lowest order in the expansion, the effective HamiltonianĤ eff =VĤV † reads: with renormalized frequencies [41]: and effective coupling g eff = [g 1 (ǫ 2 + ξ 2 ) + g 2 (ǫ 1 + ξ 1 )]/2. In Eq. (5),Ĥ 0 = ωâ †â − Ω1 2σ z + Ω 2Ĵ z d is the uncoupled Hamiltonian and [..., ...] denotes the commutator. Within our approximation, the energy spectrum consists of different manifolds characterized by a fixed value of resonator photon number operatorN =â †â (the interactions between different manifolds can be neglected due to the large resonator frequency compared to other energy scales). For the qubit (d = 2), we have (Ĵ ± 2 ) 2 = 0, (Ĵ z 2 ) 2 = 1/4, and the spectrum of the Hamiltonian in Eq. (5) can be found by diagonalizing a 4 × 4 matrix consisting of two 2 × 2 blocks. In the lowest manifold (N = 0), the four eigenenergies are given by In the general qudit case (d ≥ 3), the eigenenergies can be obtained by computing the roots of the product of two In the ultrastrong coupling regime, the numerical results to be presented below are corroborated by an analytical approach combining the adiabatic approximation in the displaced oscillator basis [42] and degenerate perturbation theory. More precisely, we first obtain the exact spectrum of the reduced Hamiltoniañ neglecting the free Hamiltonian of the qubit and qudit in the limit where Ω 1 , Ω 2 ≪ ω, and g 1 , g 2 < ∼ ω. Then, these terms are restored within a perturbative expansion.
where δ i,j is the Kronecker delta [43]. Notably, the twofold degeneration of ground state is only resolved at dorder perturbation theory inĤ ′ , with a correction proportional to Ω 1 Ω d−1

2
(see Appendix C). The corresponding eigenstates read where C is the normalization factor, and the functions C σ,m (g 1 , g 2 ) ∝ e −g 2 1 /ω 2 , e −g 2 2 /ω 2 are exponentially suppressed with the coupling strengths. As discussed above, the analytical expression in Eqs. (11)-(12) are expected to hold in the regime where the free Hamiltonian terms of the atoms are treated as perturbations to the interacting system, i.e. Ω 1 , Ω 2 ≪ ω, g 1 , g 2 .
In Fig. 2, we display the low-energy spectrum of the Hamiltonian in Eq. (1) as a function of the coupling strength (with g 1 = g 2 ) for d = 2 (panels a,d), d = 3 (panels b,e), and d = 4 (panels c,f). For visualization purposes, we plot the energy as a function of dg 1 (with d number of levels in the qudit), since the ground state energy for g 1 = g 2 scales as E 0 /ω ∼ −(dg 1 /ω) 2 for large values of g 1 . The analytical expressions obtained in the lowcoupling limit (dotted-dashed) and through perturbation theory in the ultrastrong coupling regime (dashed) are compared with the solutions obtained through numerical diagonalization (solid). In the low coupling regime, the analytical expression of Eq. (9) gives a very accurate description of the spectrum for dg 1 < ∼ 0.4Ω (see Fig. 2a), and Fig. 2d)). For the general qudit case, the expressions obtained through RWA approximation reproduce the numerical results in a less satisfactory way. Still, they correctly reproduce the spectrum for dg 1 < ∼ 0.3Ω. For Ω 1 = 0.15ω, Ω 2 = 0.1ω (top panels), an excellent agreement between the strong coupling regime approximations and numerical solutions arises for dg 1 , dg 2 > ∼ 0.3ω. For higher Ω 1 and Ω 2 values (Ω 1 = 0.55ω, Ω 2 = 0.5ω, bottom panels), the adiabatic approximation breaks down below dg 1 , dg 2 ∼ 0.75ω (Fig. 2b).
With increasing coupling strength g 1 , g 2 , the higher order correction terms in Eq. (12) are suppressed exponentially, and the states |Ψ ± approach the GHZ-type states. Note that the states |0 ↑, |0 are coherent states with opposite displacement in the phase space, which are asymptotically orthogonal in the limit g 1 , g 2 → ∞. As a result, the ground state under such large coupling assumption can be approximated as: The validity of this approximation is investigated in Fig. 3, where the fidelity between the GHZ state and the ground state of the Hamiltonian obtained through numerical diagonalization is plotted as a function of the coupling strength g 1 = g 2 . In this manuscript, we define the fidelity between two pure states |φ , |ψ as F = | φ|ψ |.
In agreement with our perturbative calculation, the fidelity of the ground state with the GHZ state approaches the unit value in the limit g 1 ≫ Ω 1 , Ω 2 .

III. NEGATIVITY
In our scheme, it is interesting to investigate the bipartite entanglement between the qubit and qudit. We choose to compute negativity [35] as the measure of entanglement. For a bipartite pure state |ϕ AB in a Ground state vs GHZ state for g1 = g2 and different qudit sizes d = 2, 3, 4 (top to bottom at g1 = 0). Fidelity between the Hamiltonian ground state (obtained through numerical diagonalization) and the GHZ state Eq. (13) as a function of the coupling strength for (a) Ω1 = 0.15ω, Ω2 = 0.1ω, where |ϕ AB ϕ| TB is the partial transpose of |ϕ AB ϕ| and · 1 is the trace norm. To extract bipartite pair-wise entanglement in a tripartite system, we use the reduced density matrix of |ϕ ABC on subsystem A ⊗ B by tracing over subsystem C: ρ AB = tr C |ϕ ABC ϕ|. Figure 4 displays the density plot of the negativity as a function of the coupling strengths g 1 , g 2 in the qubit (Fig. 4a), qutrit ( Fig. 4b), ququart ( Fig. 4c) cases, for Ω 1 = Ω 2 = 0.1ω. Note that the negativity is clearly symmetric under the exchange g 1 ↔ g 2 in the qubit case (since Ω 1 = Ω 2 ), and it becomes gradually more asymmetric by increasing the number of levels. The bipartite entanglement between the two artificial atoms has a nontrivial response to the coupling strengths between subsystems. In particular the negativity is maximum for intermediate values of the couplings, and it is strongly suppressed at large couplings. The position of the maximum depends on the number of levels in the qudit and reads: For a better visualization, in Fig. 4d we consider cuts of Figs. 4a-c at a fixed value of g 2 = 0.2ω. The negativity first rises to a maximum with increased coupling strength g 1 , before decaying to zero exponentially (as we will discuss below). This phenomenon is reported in Ref. 44 where an approximate expression is derived to explain the curve at weaker coupling. Here, we obtain an analytical expression for the decaying curve. In addition, our approach demonstrates that the entanglement suppression is a consequence of the structure of the entanglement encoded in the ground state [see Eq. (13)]: the tripartite GHZ state at large coupling limit(g 1 , g 2 → ∞), hinders the bipartite entanglement obtained after tracing out one of the subsystems, that asymptotically vanishes. This property of the GHZ states results in a counter-intuitive implication: the strong coupling can destroy entanglement between the two quantum systems connected by the resonator. Now we show that the approximate expression of the ground state, i.e., the GHZ state of Eq. (13), leads to an accurate prediction for the entanglement at large coupling; we can easily calculate the negativity for those states. Indeed, the corresponding reduced density matrix ρ ′ with resonator degree-of-freedom traced out is a (2d × 2d) matrix and has only four non-zero matrix elements: with K = exp{−2[g 1 +(d−1)g 2 ] 2 /ω 2 }. Therefore the analytical expression for the negativity of the ground state in Eq. (13) is These approximate expressions are displayed (dashed) in Fig. 4. Note that the approximation describes very accurately the exponential decay of the negativity at large coupling. We close the section by noting that the negativity measure is not a sufficient test of entanglement for systems with dimensions beyond 2×3. Under such circumstances, a state with zero negativity could possibly be a PPT or "bound entangled" state, which is argued to be metrologically useful [45][46][47].

IV. DYNAMICS
In this section, we discuss the dynamical evolution of the coupled system. We consider two complementary cases: the quench dynamics starting from the non interacting state and the adiabatic preparation of the GHZ state.

A. Quench dynamics
We start by discussing the dynamics of the system under non-adiabatic switching of the interaction. We consider the system initially prepared in the ground state of the uncoupled HamiltonianĤ 0 , i.e. |ψ 0 = |g 0 0 . For simplicity, we consider an instantaneous quench of the coupling constant to the final values g 1 = g 2 = 0.3ω. In this case, the time-evolved state reads Due to the non-adiabatic control, the state of the system is different from the ground state of the interacting Hamiltonian after the quench, and evolves in time.
Figures 5a-c display the time evolution of the fidelity between the initial state |ψ 0 and the time evolved state |ψ(t) = e −iĤt |ψ 0 in the qubit, qutrit, and ququart cases, respectively, for Ω 1 = 0.12ω and Ω 2 = 0.1ω. We note that the dynamics of |ψ 0 involves multiple frequencies, and displays a clear dependence on the number of levels of the qudit. An approximate description of the evolution can be obtained by neglecting the free terms of the atoms in the Hamiltonian. Namely, we set Ω 1,2 = 0, and expand |ψ 0 in the eigenstates basis ofH [see Eq. (10)], i.e., |ψ 0 = σmN σ m N σm |g 0 0 |σ m N σm . For the qubit case, we obtain ± α 2N ± /N !, with α ± = (g 1 ± g 2 )/ω. For the qutrit and the ququart case the summation in Eq. (18) involves d terms. Such approximated dynamics is displayed [48] in Figs. 5a-c using dashed curves. We note that the approximate expressions capture the initial decrease in fidelity and the amplitude of its oscillations, as well as the main frequency components of the evolution. We carried out the frequency analysis of the dynamical response through the fast Fourier transform of the curves displayed in Figs. 5a-c [49]. The results are reported in Figs. 5d-f. In the qubit case (d = 2, Fig. 5d), the analytical expression of Eq. , and the mean photon numberâ †â (panel c) for the qubit, qutrit and ququart cases. Consider first the qubit population σ z (t) : in all the three cases, the evolution is non-harmonic, and the number of relevant frequencies increases by increasing the number of levels in the qudit. In particular, for precise modeling of the dynamic in this regime different Fock number states must be included in the calculation and neglecting the off-diagonal terms in the basis ofH would provide a poor description of the physics (see the above discussion).
Consider for instance, the time evolution ofσ z in the qubit case (blue curve in Fig. 6d). We can work within the approximation exploited above for the strong coupling regime. Namely, we approximateĤ withH of Eq. (10), and we write |ψ 0 in the basis ofH. By performing the calculation, it is possible to derive a simple expression for σ z when g 1 = g 2 (as in the plot of Fig. 6d), namely This approximate expression is displayed in Fig. 6a (dashed). In agreement with the results displayed for the fidelity in Fig. 5a, the approximation leading to Eq. (19) gives good estimates for both the amplitude and the main frequency of the oscillations. This is confirmed by the Fourier analysis in Fig. 6d. Similar considerations apply to the time evolution of the expectation value of Ĵ z d . Since the qubit system is symmetric under the exchange 1 ↔ 2, the approximated evolutionĤ ∼H can be readily obtained from Eq. (19); recalling thatĴ z d=2 =σ z /2, we obtain Ĵ z (t) ∼ −1/2 cos(4g 2 1 t/ω − N ωt) [the minus sign comes from our sign convention in Eq. (1)]. As a consequence, the Fourier transform of Ĵ z d=2 (t) , shown in Fig. 6e, perfectly reproduces the one displayed in Fig. 6d, up to a factor 2.
Notably, the dynamics of the mean photon number â †â is much more regular, as displayed in Fig. 6c; for visualization purposes, the various curves are offset by the quantity d − 2. The plot clearly displays two relevant features: i) there is a frequency mode which is independent of the number of levels in the qudit; ii) the amplitude of the oscillations increases with d. These features can be discussed retaining the approximationĤ ∼H. In particular, by applying the Baker-Hausdorff expansion to the operator e iHtâ †â e −iHt , we derive ψ(t)|â †â |ψ(t) = 4[g 2 1 + (d − 1)g 2 2 ] sin 2 (ωt/2). (20) The validity of this approximation is investigated in Fig. 6c (dashed curves). While we find a good agreement for the qubit and the qutrit cases, the deviations are more relevant for the ququart. This result is confirmed by the Fourier transform of the curves, displayed in Fig. 6f; note that the curves have a vertical offset of (d − 2)/2 for visualization purposes. As discussed above, all the curves share a harmonic component with frequency ω, captured by the adiabatic approximation Eq. (20). This mode represents the primary frequency component in the qubit (d = 2, blue) and qutrit cases (d = 3, green). The situation is more complex for the ququart (d = 4, purple), with several sidebands and an enhanced low-frequency oscillation ∼ 0.1ω.
We conclude this section by stressing that the presence of additional levels is beneficial to inducing non-adiabatic photon generation in the ultrastrong coupling regime.

B. Adiabatic state preparation
In this section we demonstrate how the state in Eq. (13) can be prepared with high fidelity by adiabatic evolution:Ĥ The initial state and the final state correspond to the ground state of the HamiltoniansĤ in =Ĥ(t = 0) and H, respectively, and µ(t) is a function that goes from 0 to 1 when t goes from 0 to the final evolution time t f ; µ(t) is chosen to be a linear function in our simulations. Here we propose two simple schemes to prepare the hybrid GHZ state: I. switch on the couplings at fixed frequencies; II. change the frequencies at fixed couplings. In both schemes, the final HamiltonianĤ is set to be the Hamiltonian of the qubit-resonator-qudit system in Eq. (1). I. In the first approach, the system is initialized without coupling terms g 1 (t = 0) = g 2 (t = 0) = 0, i.e.,Ĥ in = During the adiabatic evolution, the coupling terms are gradually switched on to the final value g f = 0.5ω, reading g 1 (t) = g 2 (t) = µ(t)g f , see the inset in Fig. 7a. The main panel of Fig. 7a displays the evolution of the fidelity between the instantaneous state of the system under the time-dependent HamiltonianĤ(t) and the expected GHZ state at the final time, obtained by setting g 1 = g 2 = 0.5ω in Eq. (13). The different curves corresponds to the qubit, qutrit and ququart cases. In all the cases, the fidelity is minimum at t = 0 and grows monotonically with time, reaching the unit value (within numerical accuracy) for t = t f = 500/ω. Note that for a given time the fidelity is maximum in the qubit case, and typically decreases by increasing the number of levels in the qudit.
II. Here, we keep fixed the coupling terms g 1 , g 2 , while tuning the characteristic frequencies of the artificial atoms. Specifically, the system is initialized , with the initial transition frequencies satisfying Ω ′ 1 ≫ Ω 1 and Ω ′ 2 ≫ Ω 2 . In the adiabatic preparation, the artificial atoms frequencies Ω 1 (t), Ω 2 (t), are linearly reduced to the final values Ω 1 = Ω 2 = 0.1ω (see the inset of Fig. 7b). The corresponding evolution of the fidelity with the final GHZ state is shown in Fig. 7b. As in the previous case, the fidelity grows monotonically to 1 as time approaches t f in all the cases. Notably, the presence of additional levels in the qudit is displayed to be beneficial to the process: for a given t > t f /2, the fidelity is larger in the qutrit and ququart case with respect to the qubit case.

V. CONCLUSION
We formulated and studied a quantum Rabi type model describing the interaction between a two-level and a multi-level system mediated by a single mode bosonic field (in quantum technology such bosonic field is realized by a resonator). In the weak and in the strong coupling limits, we devised two different analytical schemes. In the weak-coupling limit, the effective Hamiltonian is obtained through a suitable Schrieffer-Wolff transformation. Assuming the qudit degenerating in a two-level system, the spectrum of the effective Hamiltonian has been obtained exactly. In the strong coupling limit, the ground state of the system is provided by a tripartiteentangled state of the GHZ type. A known feature of the GHZ states is that they do not allow bipartite entanglement between the three partners. As a result, qubit and qudit cannot be highly entangled, and the correlation between them drops exponentially with increasing coupling in the ultrastrong coupling regime. Such analysis is corroborated by the study of the negativity providing the sufficient condition for the qubit-qudit entanglement. We analyzed the system dynamics both under quenching and adiabatic switching of the couplings between the atoms and the resonator. The non-adiabatic nature of the quantum quench leads to the generations of photons in the resonator, with a magnified effect by increasing the number of levels in the qudit. By this procedure, the GHZ states can be adiabatically prepared with high fidelity. Both the analysis of the spectrum and the dynamics indicate that the interaction is effectively mag-nified by the number of the levels of the qudit. The study of the dynamics gives preliminary information on the qubit-qudit coherent state transfer. Our work provides relevant information for applications in quantum technology, particularly for hybrid quantum networks design [21,22,24,50]. The proposed scheme for the GHZ preparation may be implemented in cQED platforms, where both ultrastrong and deep-strong coupling regimes have been studied theoretically [51][52][53] and implemented experimentally [11,12,[54][55][56]. Very recently, the ultrastrong coupling regime has been reached in hybrid semiconducting-superconducting technology [57], too.
The generalization of the Rabi model to the qudit system can be applied to discuss superconducting circuits based on transmon qubits. Indeed, the multi-level nature of the energy spectrum cannot be ignored in these elements due to the weak anharmonicity. These circuits have been extensively investigated, even in combination with semiconducting qubits [58]. On the theoretical side, our scheme can be investigated in a more general setting. A certainly interesting direction to go would be studying the system dynamics in presence of dissipation [59].

ACKNOWLEDGMENTS
We thank Yvonne Gao, Dariel Mok, and Andrea Iorio for fruitful discussions. LCK is supported by the Ministry of Education and the National Research Foundation of Singapore. WJ Fan would like to acknowledge the support from NRF-CRP19-2017-01.

Appendix A: Low coupling spectrum
Here we report the low-coupling approximation of the effective Hamiltonian Eq. (5), obtained by neglecting the (Ĵ ± d ) 2 terms, and performing a RWA-like approximation. For the qutrit d = 3, For the ququart d = 4,

Appendix B: Adiabatic approximation
To diagonalize Eq. (1), we generalize the adiabatic approximation, which was first adopted in [42] to find solution to the single spin quantum Rabi model. In the adiabatic limit Ω 1 , Ω 2 ≪ ω, the Hamiltonian is approximated as: where the free energy terms of the atoms have been dropped from Eq. (1). We consider eigenstates ofH in the form |σ m N σ,m = |σ ⊗ |m ⊗ |N σ,m where σ =↑, ↓ are the eigenstates of σ x withσ x |↑, ↓ = ±|↑, ↓ , |m are the eigenstates of the qudit spin operator (Ĵ + d +Ĵ − d ). Next we find the eigenvalues ofH:

Qutrit case
The eigenstates |m of the qutrit operator (Ĵ + 3 +Ĵ − 3 ) with eigenvalues E m = 0, ±2, can be already obtained through diagonalization: The eigenstates |N σ,m satisfy the set of six eigenvalue equations: Repeating the steps of the previous section, we obtain the eigenstates: with eigenvalues

Ququart case
The eigenstates |m of the ququart operator (Ĵ + 4 +Ĵ − 4 ) are obtained as: Repeating the steps of the previous section, we obtain the eigenstates:

Qudit case
The above explicit results for the eigenvalues and eigenstates can be also obtain in the general qudit case by applying two unitary transformations to the adiabatic limit Hamiltonian in Eq. (B1). The transformations are of the Lang-Firsov type [60]: The transformed Hamiltonian can be written as The eigenvalues are found immediately from those of of σ x ,Ĵ + d +Ĵ − d , andâ †â , namely σ = ±1, m = −d, −d + 2, . . . , d−2, d, and N = 0, 1, . . ., respectively. The eigenstates can be correspondingly given as |σ, m, N , where after the unitary transformations the oscillator Fock state is independent of the qubit and qudit state. In the original basis, the eigenstates are found by applying to these eigenstates the above unitary transformations which, after acting on eigenstates ofσ x andĴ + d +Ĵ − d , reduce to displacement operators acting on the oscillator states, with the displacement depending on σ and m.
Appendix C: Perturbation correction to the energy and state We are interested in resolving the degeneracy in the subspace (denoted by D) spanned by {| ↑, +, 0 ↑,+ , | ↓, −, 0 ↓,− } with degenerate energy E 0 ↑,+ = E 0 ↓,− . The first order perturbation equation reads [61]: n + E (1) n |ψ (0) n , (C1) and first order energy correction: As we show below, all the first order corrections are zero. Since the energy degeneracy in the subspace of interest D is not lifted by first order corrections, we need to find another good basis that diagonalize a new matrix M : which turns out to be a 2 × 2 symmetric matrix, with M 00 = M 11 and M 01 = M 10 . Diagonalizing the matrix M , we obtain the new basis: The energy degeneracy in the subspace with lowest energy is lifted at the second order in perturbation theory. Therefore, the degeneracy-lifted energies are given by The first order correction to the ground state are obtained through the formula which produces the higher order terms in Eq. (12).
Since the Hamiltonian can be exactly diagonalized with t = 0, we can evaluate the leading contribution by treating H ′ =Ĥ d=3 −Ĥ d=3 (t = 0) as perturbation of H d=3 (t = 0). Namely, we rewrite the Hamiltonian in the previous basis, reading (keeping terms up to u 2 ) Since t is not explicitly present in the diagonal of the rotated matrix, the first order correction is zero. The second order correction in t can be computed with degenerate second order perturbation theory, giving a 2 × 2