{'fraction_non_alphanumeric': 0.0426413530719102, 'fraction_numerical': 0.009793950352205523, 'mean_word_length': 4.208554051402786, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective Hamiltonians in the Quantum Rabi Problem

We revisit the theoretical description of the ultrastrong light-matter interaction in terms of exactly solvable effective Hamiltonians. A perturbative approach based on polaronic and spin-dependent squeezing transformations provides an effective Hamiltonian for the quantum Rabi model up to the second order in the expansion parameter. The model consistently includes both rotating and counter-rotating terms, going therefore beyond the rotating wave approximation. Analytical and numerical results show that the proposed Hamiltonian performs better than the Bloch-Siegert model when calculating operator averages (e.g.\, the mean photon number and number of excitations). This improvement is due to a refined calculation of the dressed states within the present model. Regarding the frequency shift induced by the qubit-photon interaction, we find a different sign from the Bloch-Siegert value. This influences the eigenstates structure in a non-trivial way and ensures the correct calculation of the number of excitations associated to a given dressed state. As a consistency check, we show that the exactly solvable independent boson model is reproduced as a special limit case of the perturbative Hamiltonian.


I. INTRODUCTION
The physics of strongly coupled qubit-photon hybrid systems has become a steadily developing field [1,2] which challenges the existing theoretical methods. In particular, as the field-matter coupling strength increases, the Jaynes-Cummings model [3], which neglects the so-called off-resonant transitions within the rotating wave approximation (RWA), becomes less accurate.
More general models taking fully into acount the interaction between a two-level system (e.g. a superconducting qubit or a quantum dot) and a boson mode (cavity photon or transmission line resonator) are defined by the class of Hamiltonians where the quasi-spin notation is used for the qubit and a, a † are bosonic operators of the single-mode field. The corresponding energies are ε and ω respectively, and g 1 and g 2 are the coupling strengths for the resonant (rotating) and off-resonant (counter-rotating) processes. We take = 1 throughout. The Jaynes-Cummings (JC) Hamiltonian is then recovered from Eq.(1) by taking g 2 = 0, while the quantum Rabi model (QRM) [4,5] corresponds to g 1 = g 2 = g: Generally, the problems defined by Eq.(1) are not exactly soluble, one notable exception being the quantum Rabi Hamiltonian whose spectrum was calculated [6][7][8] as zeroes of recursively defined new special functions or in terms of confluent Heun functions [9]. The mathematics is quite involved both for the eigenvalues and for eigenvectors. In view of this situation, in the literature one prefers to devise effective Hamiltonians [2] to obtain approximate eigenvalues and eigenvectors, in suitably restricted domains of the parameter space defined by {ε, ω, g}. Also, simple, analytic effective descriptions are strongly preferred when including additional interactions to describe optical, transport [10,11] or relaxation processes [12][13][14].
A useful test for these approximations consists in inspecting limit cases in which exact solutions are known for the full Hamiltonian. One such case is the JC model mentioned above, in which the RWA brings the Hamiltonian to a 2 × 2 block-diagonal form. A second one is the independent boson model (IBM) [15] which is extensively used for quantum optical and transport calculations involving carrier interaction with photons, phonons [16] or vibrons [17]. In the present case the IBM limit is obtained by taking ε = 0 in the QRM, Eq. (2), and the corresponding Hamiltonian is diagonalized by a displacement of the bosonic oscillator centers brought in by the Lang-Firsov, or polaronic unitary transform.
In deriving approximate effective Hamiltonians for the quantum Rabi problem one usually employs a limited expansion of various unitary transforms in terms of a suitable small parameter. An example is provided by the dispersive Hamiltonian [18,19] obtained by a unitary transform aiming at diagonalizing the qubit-boson interaction, expanded to second order in the parameter λ = g/|ε − ω|, in the weak coupling limit λ ≪ 1. The resulting effective Hamiltonian reproduces the JC spectrum to second order in λ, i.e. within the domain of parameters in which it was derived. Moreover, it recovers the exact result in the IBM limit.
Another popular effective Hamiltonian is the Bloch-Siegert (BS) one [20,21]. It treats the QRM problem up to second order in the parameter λ = g/(ε + ω) and brings it to a renormalized form of the JC Hamiltonian. Two succesive unitary transforms are used in the process and higher order off-resonant terms are discarded as rapidly oscillating, in the spirit of the RWA. Supposedly, the validity domain of this approach is roughly g/ω 0.3 and corresponds to the so-called perturbative ultra strong coupling regime (pUSC) [2]. The perturbative criteria for the BS Hamiltonian was recently used [21] to establish a spectral classification of the field-matter coupling regimes in the QRM. On the experimental side the Bloch-Siegert shift was recently observed in strongly coupled qubit-oscillator systems [22]. An alternative improvement of the JC model relies on the Van Vleck perturbation theory to the second order in the qubit-field coupling [23].
Surprisingly, a simple analysis shows that the IBM limit is not recovered by the BS Hamiltonian, despite the fact that taking ε = 0 does not preclude the expansion parameter, now equal to g/ω, to remain small. Still, the BS model is tacitly accepted as the main candidate to describe the ultra-strong coupling regime.
In view of this situation, in what follows we show how to obtain an effective Hamiltonian for the QRM, in the same perturbative conditions as for the BS model, i.e. up to the second order in λ, having the same simple, renormalized JC structure and which additionally meets the test of recovering the exact IBM result for ε = 0. Moreover, it will become obvious that the procedure is purely perturbative, i.e. it does not have to rely on discarding off-resonant terms.
As the qubit energy ε approaches resonance the spectra of the two effective models become more similar. Nevertheless, we show that important differences remain in the dressed states and are inherited by the expectation values provided by these states. Especially the excitation number turns out to be sensitive to these differences.
The rest of the paper is organized as follows. The perturbative effective Hamiltonian is derived in Section II. A detailed comparison of this perturbative Hamiltonian with the Bloch-Siegert model is presented and discussed in Section III. Their spectra, average photon number, population inversion and number of excitations are obtained and compared both analytically and numerically in subsections III A, III B, and III C. We also compare our result with the so called generalized rotating wave Hamiltonian (GRWA) obtained by Irish [25] in the deep strong coupling case (subsection III D). This regime is defined by large values of the coupling, g/ω ≈ 1 and beyond, but the numerical evidence [26][27][28] shows that this approach gives good results at lower couplings too. Therefore comparison against it is relevant. Besides, the GRWA Hamiltonian meets the IBM limit test. Section IV is left for conclusions.

II. AN IBM-COMPATIBLE EFFECTIVE HAMILTONIAN
We consider unitary transforms of the form e S , with S † = −S and their expansion in iterated commutators with their generators S. Then any operator A transforms as Now let us rewrite the Rabi Hamiltonian of Eq.(2) as in which we separate the free qubit contribution and the independent boson model part. The linear part in H IBM can be absorbed in the quadratic term by a shift of the bosonic operators, produced by the generator with the value of λ to be chosen later. Using it is clear that the canonical transformation amounts to a simple shift As a result, the transformed Hamiltonian becomes Choosing λ = g/ω removes the interaction term and diagonalizes the IBM Hamiltonian.
Therefore the IBM problem alone is exactly solvable. The solution consists of the oscillator keeping its frequency ω, but shifted by g/ω, the shift direction being decided by the eigenstate of σ x The presence of the free qubit part in Eq.(4) spoils this simple scenario. Both the frequency and the shift are modified, and we show below that in these circumstances a different choice for the value of λ is more appropriate.
Obviously, σ x is left unchanged by the IBM transform. In contrast, the change in σ z is rather complex. One can easily check that: This means that successive commutators generate alternatively the σ z and iσ y operators, and each step brings in an additional factor −2 λ(a † − a). The σ z operator collects the even terms, which add up to a hyperbolic cosine function, while the iσ y terms generate a hyperbolic sine function. The net result is Up to now no approximation is involved. A simpler Hamiltonian, valid to the second order in λ, is obtained by expanding cosh(x) ≈ 1 + x 2 /2 and sinh(x) ≈ x. Taking into account that σ x = σ † + σ and iσ y = σ † − σ, this expansion leads to an effective Hamiltonian H ′ of the form Now it is clear that an appropriate choice for λ is which cancels out the off-resonant term linear in the bosonic operators. This is the same small parameter as the one used in the Bloch-Siegert formalism [20]. The effective Hamiltonian has now the expression Finally, one has to handle the last, quadratic term. To this end one considers a second unitary transform, with the generator This σ z -dependent squeezing transform [30] leads to in which again succesive commutators generate alternatively the same two operators. One gets a Bogolyubov-type transform Different Hamiltonian pieces change accordingly, for instance the photon number operator and From Eqs. (18), (19) above the quantities in H ′ S ′ containing (a † − a) 2 are identified as: The terms contained in X z cancel out by choosing η so that The condition makes η of second order in λ and, neglecting higher contributions, one has In this case the term X 0 in Eq. (21) is of fourth order in λ and is discarded. The terms in the Hamiltonian Eq. (14) not considered up to now are the free qubit energy, which commutes with S ′ , and the JC interaction 2λ ε(a † σ + aσ † ), which is already a first order quantity, and therefore its transformation by S ′ ∼ λ 2 brings in changes of at least third order. Therefore both are kept unchanged.
Finally, expanding systematically the hyperbolic functions to second order in λ (i.e. first order in η) and spelling out the expressions of λ and η in terms of the model parameters, leaves us with H ′ S ′ ≈ H P where the perturbative IBM-compatible effective Hamiltonian is given by This expression of the effective Hamiltonian H P is the main result of the paper.
It is now obvious that for ε = 0, H P recovers exactly the IBM result, Eq. (9) and, moreover, this holds for all values of the coupling constant g.

III. COMPARING EFFECTIVE HAMILTONIANS
The Bloch-Siegert and the perturbative effective Hamiltonian given by Eq. (24) have both the structure of a renormalized JC model The difference comes from the parameters Ω, γ and δ. They are given by [21] for H = H BS , and by for H = H P . The JC structure of Eq.(25) reduced the problem to 2 × 2 blocks, indexed by the photon number n. In the corresponding basis {|n, ↑ , |n + 1, ↓ } the n-th block has the form where σ 0 is the 2 × 2 unit matrix and we introduced with ∆ n = ε − ω + 2 Ω (n + 1) acting as an effective, n-dependent detuning. The angle θ n is defined by the relations cos θ n = ∆ n /(2R n ) and sin θ n = γ √ n + 1/R n . Note that Ω/2 + δ turns out to take the same value g 2 /(ω + ε) in both effective Hamiltonians considered, therefore the middle points of their spectra coincide.
It is known that the spectrum ofĥ n is E n,± = α n ± R n with the eigenvectors (i.e. dressed states) given by the linear combination where (u + , v + ) = (cos θ n /2, sin θ n /2) for the higher state |ϕ n,+ and (u − , v − ) = (sin θ n /2, − cos θ n /2) for the lower state |ϕ n,− . To complete the picture we also calculated give the energy of the ground state |0, ↓ as E 0 = −ε/2 − (Ω/2 + δ). By examining Eqs. (26) and (27), it is clear that the differences between the two models are more pronounced for small values of ε. As already mentioned, at ε = 0, H P recovers the exact IBM result, while H BS does not.
On the other hand at resonance, ε = ω, the parameters of the two models coincide, except for the sign of Ω. The latter has no bearing on the spectra, which therefore are the same for the two models, but it modifies the structure of the eigenstates. Indeed, at resonance the sign of Ω carries over to a sign change in cos θ n which in turn leads to exchanging cos θ n /2 and sin θ n /2. In other words, even though at resonance the spectra of H BS and H P coincide, because of the sign difference in ∆ n the attribution of the eigenvalues is flipped: for example, the weight of the two basis vectors |n, ↑ and |n + 1, ↓ in the upper state |ϕ n,+ of the BS Hamiltonian is the same as for the lower state |ϕ n,− of the perturbative Hamiltonian. Further consequences of this fact are analyzed in subsection III C.
To resume the situation, the two models predict different spectra and states, as suggested by these two examples. Intermediate values for ε ∈ [0, ω] are expected to lead to similar conclusions. Let us also remark that Ω is interpreted as a bosonic frequency shift which depends on the qubit state, and detecting the former gives information about the latter [2]. Therefore the sign difference means that this way of reading the qubit state works in opposite directions in the predictions of the two models.
Still, these arguments are merely heuristic. Eigenstates are not directly comparable, H P and H BS having been obtained by different unitary transforms, i.e. the two effective Hamiltonians work in different bases. The relevant comparison is on basis-independent quantities, to which we turn below.

A. Energy spectra
At resonance, ε = ω, all parameters in Eq.(28) are exactly the same for the two cases, except for sign of ∆ n , and thus the spectra are identical. On the other hand, in the limit ε = 0, H P becomes exact. Therefore, it is expected that for the whole subresonant range, ε ω, our approach is an improvement over the BS efective Hamiltonian.  To illustrate this fact we compare in Fig. 1 the spectra of the two effective models (BS and perturbative), and those produced by numerical diagonalization of the full Rabi Hamiltonian. For the latter a cutoff at N ph = 30 photons is found to be sufficient. We use ω = 1 units throughout, since the problem is homogeneous of degree 1 in the parameters, i.e. it is defined essentially by only two dimensionless parameters ε/ω and g/ω.
We selected the pair of spectral branches E 4,± for three values of the qubit energy ε. However, in the IBM limit (i.e for ε = 0, see Fig. 1(a)) two additional branches are added, to illustrate that the coinciding spectra of the Rabi and perturbative Hamiltonians are doubly degenerate, while the Bloch-Siegert spectrum displays a degeneracy lifting as the qubit-photon coupling strenght increases. Note that for the sake of clarity the spectral branches of H P are represented by red circles.
To understand these differences one has to recall that in the limit ε = 0, H P becomes exact, while the BS spectra are correct only to λ 2 accuracy, as it is easy to check on the expression for E n,± . With increasing ε the spectral difference between the models diminishes, see Fig. 1(b),(c), as predicted by the analytic expression of H P and H BS until, for ε = 1 they coincide (not shown). In the process the perturbative values remain closer to the exact Rabi result than the BS ones.

B. Operator averages. Photon population
When calculating the expectation value of an operator on a given eigenstate of the effective Hamiltonian, one has to take into acccount the unitary transforms that lead to it. The operators should be 'rotated' by the same procedure as the Hamiltonian.
In the case of H P , we use again the two generators S = λ(a † −a)σ x and S ′ = 1 2 η(a † 2 −a 2 )σ z , with the parameters given in Eqs. (13) and (23). The observable we are interested in is the photon numbern = a † a. From Eq.(7) we obtain by the first transform The second step involves S ′ and is apparently more complicated but it gets simpler by limiting the expansion to second order in λ and taking into account that η ∼ λ 2 . For instance, it is clear that the second term in Eq.(32) remains unchanged. Indeed, changes are at least of first power of η and with the prefactor λ they are third order quantities. For the first term a † a one uses Eq.(18) expanded to second order in λ, with the net result that after both transforms one haŝ This is the operator to be averaged on the eigenstates of the form given in Eq.(31). Things turn out to be further simplified by noting that operators like a 2 , aσ and their hermitic conjugates do not contribute to the expectation value, since they do not connect vectors belonging to the same two-dimensional subspace.
To conclude, one has n P,± = ϕ n,± | [a † a − λ(a † σ + aσ † ) + λ 2 ] |ϕ n,± where the index P recalls that the average is calculated w.r.t. the dressed states of H P . Obviously, n P,± depends on the specific two-dimensional subspace n as well but we drop this additional index for the simplicity of writing. Replacing u ± and v ± leads to the values A negative sign of Ω means a lower cos θ n value, and this encourages a higher photon population in the upper state. For the ground state one obtains n P,0 = λ 2 . The prediction of the Bloch-Siegert Hamiltonian for the same quantity is evaluated by the same procedure. Some differences arise though, beginning with the first generator which now readsS = λ(a † σ † − aσ) with the result that now The second transform has the same formS ′ = 1 2 η(a † 2 − a 2 )σ z with η ∼ λ 2 . As before, the terms containing λ remain unchanged up to the second order, and the first term a † a is modified by a combination of a 2 and a † 2 , which do not contribute to the expectation value. The new feature is that no contribution comes from the second term in Eq.(36) above. Therefore one is lead to n BS,± = ϕ n,± | [a † a + λ 2 σσ † − λ 2 a † aσ z ] |ϕ n,± = u 2 ± (n − λ 2 n) + v 2 ± (n + 1 + λ 2 (n + 2)) .
As a result one has n BS,± = n + 1 2 For the ground state one finds n BS,0 = λ 2 . The numerical results presented in Fig. 2 complement our analysis. We calculated the mean photon number N corresponding to the eigenstates of the fourth two-dimensional Jaynes-Cummings subspace for the exact and the effective Hamiltonians. Again, we find that in the IBM limit (see Fig. 2(a)) the perturbative Hamiltonian reproduces the exact expectation value of the Rabi model (again, we used red circles to highlight the results of the perturbative model). Note that for ε = 0 one has S ′ = 0 too. Thus one is left with only the first transform Eq.(32), which is exact. Thus not only the Hamitonian, but the 'rotated' form of then operator is exact. In contrast, the results obtained from the Bloch-Siegert model show noticeable differences starting around the qubit-photon coupling strenght g ∼ 0.15. In particular, the mean photon number corresponding to |ϕ 4,+ decreases as g increases. These differences are even more pronounced for higher two-dimensional subspaces. As in the case of the spectrum, the photon number expectation values given by the two models become closer as ε increases, with the perturbative ones remaining nearer to the exact Rabi values (see Figs. 2(b),(c)).

C. Operator averages. Population inversion and number of excitations
More clearcut differences exist between the two models than suggested by the above results, especially at resonance. We show this by performing the same comparison of expectation values, this time for the population inversion operator σ z and for the number of excitationsN x =n + 1/2(1 + σ z ).
For H P the change of σ z under the first unitary transform is obtained as in the derivation of Eqs.(10), (11).
The expression remains unchanged by the second transform, which commutes with σ z . The terms containing λ are left unchanged, as discussed in the previous cases. Also, as before, the operators a 2 , a † 2 and the off-resonant ones a † σ † , aσ do not contribute to the expectation value. We are left with which by straightforward calculations leads to Specifically for the two states one obtains σ z P,± = ± cos θ n [1 − 4λ 2 (n + 1)] ±2λ sin θ n √ n + 1 + 2λ 2 .
A negative value for cos θ n hints to a higher population inversion on the lower state. For the ground state one finds σ z P,0 = −1 + 2λ 2 . Finally, we address the same problem for H BS . The transform of σ z underS = λ(a † σ † − aσ) is complicated, but for the second order in λ the quantities [S, σ z ] and [S, [S, σ z ]] are sufficient. It is easy to check that This leads to the result Specifically, for the two states one has σ z BS,± = ± cos θ n (1 − 2λ 2 n) + 2λ 2 (1 ∓ cos θ n ) = ± cos θ n [1 − 2λ 2 (n + 1)] + 2λ 2 .
This is not a mere coincidence, and relies essentially on the fact that the generators of the respective unitary transforms differ byS := S −S = λ(a † σ − a σ † ), which commutes withN x [24].
Again, the interesting situation is provided by the resonant case. There one can see analytically that both models predict exactly the same results for a given spectral subspace. The curves for N ± x are superimposed since, except for the sign, the cosine values coincide. But because of this sign difference, inherited from the opposite signs of Ω P and Ω BS , the assignment of N x expectation values to the energy branches is switched, N ± x,BS = N ∓ x,P . More precisely, in the BS model N − x , corresponding to the lower branch, takes the higher value, while N + x , stemming from the upper branch, is smaller. In the perturbative case, it is the other way around. This is an important difference, expressing the distinct structure of the eigenstates, and which cannot gradually disappear as λ gets smaller. To decide which attribution is the correct one we rely on comparison with the result of the full Rabi Hamiltonian.
The outcome is shown in Fig. 3, in which the values of the two models are plotted, along with the numerical calculation using the dressed states provided by diagonalization of the Rabi Hamiltonian. The left panel corresponds to the eigenstates |ϕ 4,− and the right one to |ϕ 4,+ . One can check by inspection of Figs. 3 (a) and (b) that N ± x,P = N ∓ x,BS , and that comparison with the Rabi solutions decides in favor of the perturbative results.
Finally we compare in Fig. 4 the number of excitations for the same values of the qubit energy considered in previous discussions. In the IBM-limit N ± x,P follow qualitatively the exact curves but also deviate from them as g increases. (Let us stress that for ε = 0 the exact result is analytically available in terms of Laguerre polynomials [15], but here we stick to Eq.(47), in which the operators are rotated only to second order). On the other hand, the excitation numbers of the BS model are quite far from the N ± x,R for all values of ε. In particular, at ε = 0.66 the branches of N ± x,BS exhibit an incorrect crossing point around g = 0.235. Moreover, the BS attribution of the excitation number to a given spectral branch is still flipped and the inequality N + x,BS < N − x,BS still holds. On the contrary, the branches N ± x,P follow closely the exact ones in the sub-resonant regime too. The success of the model in the DSC regime relies on the fact that all terms arising from the transformed qubit energy become small at large g, due to the exponential prefactor. Suprisingly, its predictions work quite well at lower couplings too. This means that at least in certain circumstances the off-resonant terms are negligible indeed. On the contrary, in situations like ε ≈ 2ω this is not the case any more. Terms like a † 2 σ now describe resonant transitions in the system and thus are no longer negligible. In contrast, these processes are included in the present perturbative approach by its second unitary transform S ′ .
The situation is illustrated in Fig. 5 where we compare a part of the low-energy spectra of the Rabi, GRWA and perturbative Hamiltonians in the resonant (ε = ω) and supra-resonant (ε > ω) regimes. The eigenvalues of the GRWA Hamiltonian [25] are given analytically by Eq.(52). Clearly, the eigenvalues of both GRWA and perturbative Hamiltonians reasonably follow the Rabi spectrum in the resonant case along the whole range of the pUSC regime (see Fig. 5(a)). However, Fig. 5(b) reveals that in the supra-resonant regime the GRWA spectrum is no longer accurate in the pUSC regime while the errors of the perturbative model are quite small.

IV. CONCLUSIONS
The effective Hamiltonian given by Eq.(24) was derived from the quantum Rabi model systematically to second order in the small parameter λ = g/(ε + ω). No other approximations (like discarding terms as off-resonant) are involved. In this sense the proposed model is purely perturbative. The expansion parameter λ is the same as in the Bloch-Siegert theory. In contrast to the latter, the present formalism recovers exactly the IBM Hamiltonian in the limit ε → 0. The proposed model allows analytical calculations of various operator averages which are used for a detailed comparison with the predictions of the BS approach.
Even in the cases when the spectra of two models are close, we find that the main difference from the BS results is seen in the structure of the dressed eigenstates. The source of this dissimilarity is traced back to the different sign of the frequency shift induced by the qubit-field coupling. This is furher reflected in the expectation values of relevant observables, like the excitation number, and their attribution to the appropriate eigenstates.
The energy spectrum and the dressed state structure associated to the full quantum Rabi model are well approximated by the IBM-compatible effective Hamiltonian in the resonant and off-resonant perturbative USC regime. Our analysis shows that recovering the independent-boson model in the limit ε → 0 is a valuable test for effective Hamiltonian approaches to the quantum Rabi model.