Mean photon number dependent variational method to the Rabi model

TL;DR

This paper introduces a mean-photon-number dependent variational method for the Rabi model, effectively spanning all coupling regimes and improving upon existing approximations for both ground and excited states.

Contribution

It develops a novel variational approach that adapts to photon number, unifying and extending previous methods for the Rabi model across coupling regimes.

Findings

Accurately evaluates ground and excited states in the Rabi model.

Recovers previous approximations in their valid regimes.

Provides improved results over existing methods for excited states.

Abstract

We present a mean-photon-number dependent variational method, which works well in whole coupling regime if the photon energy is dominant over the spin-flipping, to evaluate the properties of the Rabi model for both the ground state and the excited states. For the ground state, it is shown that the previous approximate methods, the generalized rotating-wave approximation (only working well in the strong coupling limit) and the generalized variational method (only working well in the weak coupling limit), can be recovered in the corresponding coupling limits. The key point of our method is to tailor the merits of these two existing methods by introducing a mean-photon-number dependent variational parameter. For the excited states,our method yields considerable improvements over the generalized rotating-wave approximation. The variational method proposed could be readily applied to the more complex models, for which an analytic formula is difficult to be formulated.

Equations

$$H=\omega a^{†}a+\frac{\Omega }{2}{\sigma }_x+g({\sigma }_{-}+{\sigma }_{+})(a+a^{†}),$$

$$H_{\text{RWA}}=\omega a^{†}a+\frac{\Omega }{2}{\sigma }_x+g({\sigma }_{-}{a}^{†}+{\sigma }_{+}a).$$

$$E_{\mathrm{RWA}}^{(\pm ,N)}=(N-\frac{1}{2})\omega \pm \sqrt{\frac{(\omega -\Omega {)}^2}{4}+N{g}^2},$$

$$\tilde{H}=\omega a^{†}a-\frac{g^2}{\omega }+\frac{\Omega }{2}{\sigma }_{x}F(\lambda )+\frac{i\Omega }{2}{\sigma }_{y}G(\lambda ).$$

$${\tilde{H}}_{\text{AA}}=\omega a^{†}a-\frac{g^2}{\omega }+\frac{\Omega f_0(\lambda ,a^{†}a)}{2}{\sigma }_x,$$

$$\begin{array}{rl}& E_{\text{AA}}^{\pm ,N}=N\omega -\frac{g^2}{\omega }\pm \frac{\Omega f_0(\lambda ,N)}{2},\\ & |{\tilde{\Psi }}_{\text{AA}}^{\pm ,N}⟩=|{\pm }_x,N⟩,\end{array}$$

$$\begin{array}{rl}|{\Psi }_{\text{AA}}^{\pm ,N}⟩& =U^{†}|{\pm }_x,N⟩\\ & =e^{\lambda (a^{†}-a)}|{+}_z,N⟩\pm e^{-\lambda (a^{†}-a)}|{-}_z,N⟩.\end{array}$$

$${\tilde{H}}_{\text{GRWA}}={\tilde{H}}_{\text{AA}}+\frac{\Omega }{2}[{\sigma }_{-}{a}^{†}{f}_1(\lambda ,a^{†}a)+\text{h.c.}].$$

$$\tilde{H}_{\text{GRWA}}^{\text{BLOCK}}=\left(\begin{array}[]{cc}E_{\text{AA}}^{+,N-1}&h^{\prime}_{N-1_{+},N_{-}}\\ h^{\prime}_{N_{-},N-1_{+}}&E^{-,N}_{\text{AA}}\\ \end{array}\right),$$

$$h_{N-1_{+},N_{-}}^{\prime }=h_{N_{-},N-1_{+}}^{\prime }=\frac{1}{2}\Omega \sqrt{N}{f}_1(\lambda ,N).$$

$$|{\tilde{\Psi }}_{\text{GRWA}}^{\pm ,N}⟩$$

$$|{\Psi }_{\text{GRWA}}^{\pm ,N}⟩=U^{†}|{\tilde{\Psi }}_{\text{GRWA}}^{\pm ,N}⟩,$$

$$E_{\text{GVM}}^{(0)}={\lambda }^2\omega +2g\lambda -\frac{\Omega }{2}{f}_0(\lambda ,0),$$

$$\lambda =-\frac{g}{\omega }\frac{1}{1+\frac{\Omega }{\omega }{e}^{-2g^2/(\omega +\Omega {)}^2}}.$$

$$H=\left(\begin{array}[]{ccccc}E_{0}&h_{0_{-},0_{+}}&0&h_{0_{-},1_{+}}&\cdots\\ h_{0_{+},0_{-}}&E_{0}&h_{0_{+},1_{-}}&0&\cdots\\ 0&h_{1_{-},0_{+}}&E_{1}&h_{1_{-},1_{+}}&\cdots\\ h_{1_{+},0_{-}}&0&h_{1_{+},1_{-}}&E_{1}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right),$$

$$H_{\text{AA}}=\left(\begin{array}[]{ccccc}E_{0}&h_{0_{-},0_{+}}&0&0&\cdots\\ h_{0_{+},0_{-}}&E_{0}&0&0&\cdots\\ 0&0&E_{1}&h_{1_{-},1_{+}}&\cdots\\ 0&0&h_{1_{+},1_{-}}&E_{1}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right).$$

$$E_{\text{AA}}^{\pm ,N}$$

$$|{\Psi }_{\text{AA}}^{\pm ,N}⟩$$

$$H=\left(\begin{array}[]{cccccc}E_{\text{AA}}^{-,0}&0&0&h^{\prime}_{0_{-},1_{+}}&h^{\prime}_{0_{-},2_{-}}&\cdots\\ 0&E_{\text{AA}}^{+,0}&h^{\prime}_{0_{+},1_{-}}&0&0&\cdots\\ 0&h^{\prime}_{1_{-},0_{+}}&E^{-,1}_{\text{AA}}&0&0&\cdots\\ h^{\prime}_{1_{+},0_{-}}&0&0&E_{\text{AA}}^{+,1}&h^{\prime}_{1_{+},2_{-}}&\cdots\\ h^{\prime}_{2_{-},0_{-}}&0&0&h^{\prime}_{2_{-},1_{+}}&E_{\text{AA}}^{-,2}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right),$$

$$H_{\text{GRWA}}=\left(\begin{array}[]{cccccc}E_{\text{AA}}^{-,0}&0&0&0&0&\cdots\\ 0&E_{\text{AA}}^{+,0}&h^{\prime}_{0_{+},1_{-}}&0&0&\cdots\\ 0&h^{\prime}_{1_{-},0_{+}}&E^{-,1}_{\text{AA}}&0&0&\cdots\\ 0&0&0&E_{\text{AA}}^{+,1}&h^{\prime}_{1_{+},2_{-}}&\cdots\\ 0&0&0&h^{\prime}_{2_{-},1_{+}}&E_{\text{AA}}^{-,2}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right).$$

$$|G_A⟩=\frac{1}{\sqrt{2}}(|{+}_z,\lambda ⟩-|{-}_z,-\lambda ⟩),$$

$$E_0=⟨G_A|H|G_A⟩={\lambda }^2\omega +2g\lambda -\frac{\Omega }{2}{e}^{-2{\lambda }^2},$$

$$⟨a^{†}a{⟩}_0=⟨G_A|a^{†}a|G_A⟩={\lambda }^2.$$

$$|{\tilde{\Psi }}_0⟩=(1+K{)}^{-1/2}(|{\tilde{G}}_A⟩+{\sum _{\{\pm ,N\}}}^{\prime }{c}_{\pm ,N}|{\pm }_x,N⟩),$$

$$P(\lambda )=(1+K{)}^{-1/2}.$$

$${\lambda }_{\text{eff}}=\frac{K_B{\lambda }_A+K_A{\lambda }_B}{K_A+K_B},$$

$$|{\Psi }_A^{\pm ,N}⟩=|{\Psi }_{\mathrm{GRWA}}^{\pm ,N}(\lambda )⟩.$$

$$|{\tilde{\Psi }}_{\pm ,N}⟩=\frac{1}{\sqrt{1+K}}(|{\tilde{\Psi }}_A^{\pm ,N}⟩+{\sum _{\{\pm ,M\}}}^{\prime }{c}_{\pm ,M}|{\tilde{\Psi }}_A^{\pm ,M}⟩),$$

$$E_A^{\pm ,N}$$

$$⟨a^{†}a{⟩}_A^{\pm ,N}$$

$$F^{\pm ,N}=R_{N,\pm }^2[(\frac{\sqrt{N-1}{S}_{N-1,+}}{E_{+,N-1}^A-E_{\pm ,N}^A}{)}^2+(\frac{\sqrt{N-1}{S}_{N-1,-}}{E_{-,N-1}^A-E_{\pm ,N}^A}{)}^2]$$

$$+S_{N,\pm }^2[(\frac{\sqrt{N}{R}_{N+1,+}}{E_{+,N+1}^A-E_{\pm ,N}^A}{)}^2+(\frac{\sqrt{N}{R}_{N+1,-}}{E_{+,N+1}^A-E_{\pm ,N}^A}{)}^2].$$

Full text

I Introduction

The Rabi model describes a two-level system interacting with a single-mode bosonic field rabi . It plays a fundamental role in quantum optics cavityQED , quantum information raimond , and condensed matter physics holstein . Although it has been intensively explored, only in recent years has the integrability of this model been formulated Braak2011 . However, this analytic achievement is not the end of the study on this model, oppositely, it has triggered more theoretical and experimental studies Solano2011 ; Romero2012 ; Restrepo2014 . Explicitly, the Rabi model has been experimentally simulated in optical waveguide Crespi2012 , superconducting circuit system circuit1 ; circuit2 ; StrongCouplingExps2 ; ultra , and solid-state semiconductorGunter2004 ; Deveaud2007 ; Cristofolini2012 ; Carusotto2013 ; wen , which provides a perfect test bed to explore the physics of light-matter interaction in the deep strong coupling regime. Another significance of the analytic achievement is that it supplies some insight to understand the involved physics, e.g., the vacuum induced Berry phase Larson2012 , and the quantum phase transition in the related multi-model Rabi model, i.e., the spin-boson model Leggett_RMP ; weiss_book , for which an exact solution is quite difficult to obtain, and a well-established approximate method is desirable.

For decades of study on the Rabi model, besides the numerical treatment casanova , there exist many approximate analytic methods hausinger ; firstorder ; qhchen . The most famous approximation is the rotating-wave approximation (RWA) jc . Working in the near-resonance and weak-coupling regime, the RWA neglects the counter-rotating terms in the interaction and results in the Jaynes-Cummings (J-C) model jc . It has served as a basic starting point in understanding many quantum phenomena involved in light-matter interaction jc_review , because most of the practical quantum optical experiments work in the weak coupling regime raimond ; mabuchi . However, in circuit quantum electrodynamics system, the neglected counter-rotating term becomes important due to the strong circuit1 ; circuit2 or the ultra-strong coupling ultra between the bosonic field and the two-level system. To treat the strong coupling, Irish et al. adiabatic proposed an adiabatic approximation (AA) in the limit that the frequency of the field is much larger than the one of the two-level system. After working in the displaced oscillator basis, it takes the frequency of the two-level system as perturbation and results in a truncated Hamiltonian with the interaction effects collected in a renormalization factor to the frequency of the two-level system. In 2007 grwa , the AA was improved by considering the RWA-type interaction in the reformulated Hamiltonian in the displaced oscillator basis. This scheme was named as generalized RWA (GRWA). Although the GRWA works well in a quite broad parameter regime, especially in the strong coupling regime, it does not work well in the weak coupling regime, especially for the positive detuning case. In addition, the mean photon number predicted by the GRWA is independent of the frequency of the two-level system, which is actually not true. As an improvement, a generalized variational method (GVM) GVMground ; GVM has been introduced, where the displacement of the displaced oscillator basis is determined by minimizing the ground state energy. Indeed, the GVM evidently improves the GRWA in weak coupling regime with positive detuning, and yields a frequency dependent ground state mean photon number. However, for strong coupling and intermediate coupling regimes, the GVM is no longer applicable. Moreover the GVM is limited to the ground state.

Obviously, the merit of the GRWA and the AA comes from the introduction of the displaced oscillator basis, which captures the essential physics in the large coupling regime. However, its disadvantage lies in fixing the displacement, which leads to a frequency independent mean photon number of the obtained ground state. On the contrary, the GVM frees the displacement, but it does not introduce the displaced oscillator basis and has been excessively simplified in the analytic treatment. In the present work, we combine the merits of the GRWA and the GVM to obtain a novel analytic method. We start from the GRWA formula but further introduce a mean photon number dependent variational method to determine the displacement. As a result, our approximation method is applicable in both weak and strong coupling regimes. In the weak coupling regime, it recovers the result of the GVM and in the strong coupling regime it recovers the GRWA. In the intermediate coupling, it provides a natural crossover from the GVM to the GRWA. This variational method is not only valid for the ground state, but also for the excited states. To show the merit of the our method, we focus on the energy spectrum and mean photon number of the Rabi model and compare the result with that obtained by the GVM and the GRWA, taking the exact numerical result as a benchmark.

The paper is organized as follows. In Sec. II we introduce the Rabi model and give a review to the previous approximate methods for self containing and also for convenience of later discussions. In Sec. III we present our method and make some detailed comparisons with the results obtained by the previous methods. Finally, Sec. IV is devoted to conclusions and discussions.

II The model and some previous methods

The Hamiltonian of the Rabi model reads

[TABLE]

where $a$ and $a^{{\dagger}}$ are the annihilation and creation operators of the quantized single-mode bosonic field with frequency $\omega$, $\sigma_{x}$ is the Pauli matrix for the two-level system with level splitting $\Omega$, and $\sigma_{\pm}=(\sigma_{z}\mp i\sigma_{y})/2$ are the transition operators between the two levels, and $g$ is the coupling strength. Here, for convenience of comparison we follow the notations in Ref.grwa to use spin-flipping $\sigma_{x}$ for the level-splitting term instead of $\sigma_{z}$ commonly used in quantum optics scully . However, these two notations can be transformed into each other by a rotation on the two-level system. According to the tuning relationship between the two-level system and the field, the model takes three cases: resonance ($\omega=\Omega$), positive detuning ($\omega<\Omega$) and negative detuning ($\omega>\Omega$). Throughout the paper we take $\Omega$ as unit of energy.

Essentially, the existing approximate methods can be formulated in two ways: One is to truncate Eq. (1) into J-C-like exactly solvable form, and the other is to expand Eq. (1) on a proper basis and then truncate the obtained matrix into the block-diagonal form. In the following, we reformulate these approximations in the two ways in order to compare their performance.

II.1 Truncated Hamiltonian

RWA: Neglecting the counter-rotating terms $\sigma_{-}a+\sigma_{+}a^{\dagger}$ in Eq. (1) yields the RWA Hamiltonian

[TABLE]

This is the J-C Hamiltonian jc , which is exactly solvable. Its eigen solution reads

[TABLE]

with the ground eigen-energy $E_{\text{RWA}}^{(0)}=-\frac{\Omega}{2}$, which is just the J-C energy ladder Fink2008 . 2. 2.

AA: Performing a unitary transformation $U=e^{\lambda\sigma_{z}(a-a^{{\dagger}})}$ with $\lambda=-\frac{g}{\omega}$ to Eq. (1), one obtains $\tilde{H}=UHU^{{\dagger}}$ with adiabatic

[TABLE]

Here $F(\lambda)=\sum_{k=0}^{\infty}[a^{{\dagger}2k}f_{2k}(\lambda,a^{{\dagger}}a)+\text{h.c.}]$, $G(\lambda)=\sum_{k=0}^{\infty}[a^{{\dagger}2k+1}f_{2k+1}(\lambda,a^{{\dagger}}a)-\text{h.c.}]$, and

$f_{m}(\lambda,x)=\frac{(-2\lambda)^{m}e^{-2\lambda^{2}}(x+m)!}{x!}L_{x}^{m}(4\lambda^{2})$ with $L_{x}^{m}$ being the associated Laguerre polynomial (see Appendix A). In the small $\Omega$ case[$\Omega\ll(\omega,g)$], keeping only the zero-th order term of $a$ and $a^{\dagger}$ in $F(\lambda)$ is a good approximation, which leads to

[TABLE]

whose eigensolution can be evaluated readily as

[TABLE]

with $|\pm_{x}\rangle$ being the eigenstates of $\sigma_{x}$ and $|N\rangle$ being the Fock state. After the inverse transformation, through representing the $|\pm_{x}\rangle$ by the original $|\pm_{z}\rangle$ basis, one gets the eigen-state under the AA:

[TABLE] 3. 3.

GRWA: Going beyond the AA, one further considers the zeroth order term in G($\lambda$) involving one-excitation terms. Only considering the “energy-conserving” one-excitation terms, one arrives at the GRWA Hamiltonian grwa

[TABLE]

On the basis of $|\pm_{x},N\rangle$, Eq. (8) is block-diagonalized with $2\times 2$ subblocks

[TABLE]

which gives a pair of eigen-vectors $\{R_{N,\pm},S_{N,\pm}\}$. The off-diagonal entries are defined by

[TABLE]

Thus, the eigenstates of Eq. (8) read as

[TABLE]

The states to the original Hamiltonian (1) are obtained by the inverse transformation:

[TABLE]

while the ground state $|\tilde{\Psi}_{\text{GRWA}}^{(0)}\rangle=|-_{x},0\rangle$ is the same as that of AA. 4. 4.

GVM: Different from the above two methods, the parameter $\lambda$ here is not fixed but is optimized by minimizing the ground-state energy GVMground

[TABLE]

which results in the equation to determine $\lambda$ as $g+\omega\lambda+\lambda e^{-2\lambda^{2}}=0$. Since it cannot be solved analytically, Zhang et al. GVMground took the following approximate solution

[TABLE]

Below we address the conditions under which the above methods work well. The RWA is valid in the very weak coupling regime ($g\ll\Omega,\omega$) and under the near-resonance ($\omega\sim\Omega$) conditions. Beyond the usual strong coupling regime, namely, in the strong coupling limit, the RWA is no longer valid but the AA shows its advantage. For either large $\omega$ or large $g$, the term of displaced oscillator is dominant in (4) and the $\Omega$ terms can be treated as perturbation. Thus the validity of the AA lies in strong coupling limit ($g\gg\omega$) or negative detuning ($\omega>\Omega$) regime. Because the GRWA further keeps all one-excitation “energy-conserving” terms unincorporated in the AA, its applicable range for the excited states is extended to the regime that covers those of both RWA and AA, which is nearly the whole parameter regime. The reason can be due to “the fundamental similarity between the standard RWA and AA model: both involved calculating the energy splitting due to an interaction between two otherwise degenerate basis states”, as clearly stated in Ref. [grwa, ]. However, the validity regime of the GRWA could be further broadened if the following aspects can be properly treated. First, the ground-state energy of the GRWA is the same as the AA, no improvement has been obtained. Second, its energy spectrum requires a more accurate calculation for small ratio of $\omega/\Omega$ in the weak coupling regime, especially for the ground state. Third, it predicts an incorrect $\Omega$-independent mean photon number due to the fixed $\lambda$. The GVM improves the accuracy of the ground-state energy and its mean photon number behavior captures the $\Omega$-dependent property in the weak coupling regime, especially for the positive detuning case. However, since an oversimplified analytic treatment has been applied, the results of the GVM becomes even worse than the GRWA in the strong coupling limit regime.

II.2 Basis Formulation

Truncating the Hamiltonian in AA and GRWA can be understood in an alternative way by discarding the remote off-diagonal elements of the Hamiltonian matrix on certain basis grwa ; adiabatic . Here we reformulate the AA and the GRWA based on this idea.

Choosing the basis $|N_{\pm}\rangle=e^{-\lambda\sigma_{z}(a-a^{{\dagger}})}|\pm_{z},N\rangle$ with $\lambda=-g/\omega$, Eq. (1) can be rewritten as

[TABLE]

with $E_{N}=\omega N$ and $h_{N_{\alpha},M_{\beta}}=\langle N_{\alpha}|H|M_{\beta}\rangle$. Discarding the remote off-diagonal elements leads to a $2\times 2$ block-diagonal matrix

[TABLE]

By diagonalizing Eq. (16), one finds the eigen solution

[TABLE]

which matches well with Eq. (7) obtained under the AA.

Irish et al. further used the eigenstates in Eq. (18) as basis to expand the Hamiltonian (1), which reads

[TABLE]

with $h^{\prime}_{N_{\alpha},M_{\beta}}=\langle\Psi^{(\alpha,N)}_{\text{AA}}|H|\Psi^{(\beta,M)}_{\text{AA}}\rangle$. Then dropping the remote off-diagonal matrix elements gives rise to

[TABLE]

Based on this form, the energy spectra and the eigenstates can be readily solved, which are consistent with those obtained under the GRWA, i.e., Eq. (11).

III Mean photon-number dependent variational method

III.1 Method description and improvements for the ground state

From the above analysis on the previous approximations, we can see that truncating the Hamiltonian matrix into block-diagonalized form in a completed orthogonal basis is equivalent to truncating the Hamiltonian operator expansions. The better a basis is chosen, the more information an approximation obtains. Moreover, a proper basis can even be considered as an approximate state. First, let us focus on the ground state. For the existing approximate methods, the ground state takes the following form

[TABLE]

where $|\pm\lambda\rangle=e^{\pm\lambda(a^{\dagger}-a)}|0\rangle$ are the coherent states. For the AA/GRWA, $\lambda=-g/\omega$, while for the GVM, $\lambda$ is approximately given by Eq. (14). This motivates us to take Eq. (21) as our trial state but completely free the parameter $\lambda$. The reasons for taking the form of Eq. (21) as our trial state are as follows. On one hand, it is seen that Eq. (21) can reproduce the previous known results from different approximations like AA/GRWA and GVM. In particular, the Irish’s scheme is valid in nearly the whole coupling regime. On the other hand, it is motivated from the competition nature between the displaced oscillator and spin-flipping in the model. It is noticed that if the spin-flipping term $\frac{\Omega}{2}\sigma_{x}$ is neglected, the model reduces to a displaced oscillator Hamiltonian with two degenerated ground states $|+_{z},-g/\omega\rangle$ and $|-_{z},g/\omega\rangle$. Considering an infinitesimal spin-flipping, the degeneracy is lifted and their linear combination, namely, $\frac{1}{\sqrt{2}}(|+_{z},-g/\omega\rangle-|-_{z},g/\omega\rangle)$ is just the trial state Eq.(21) with $\lambda=-g/\omega$. When increasing the spin-flipping, the competition between the displaced oscillator and the spin-flipping motivates us to free the $\lambda$ as a variational parameter. In addition, $|G_{A}\rangle$ is also eigenstate of parity operator $\Pi=-\sigma_{x}(-1)^{a^{{\dagger}}a}$, which commutes with the model Hamiltonian. Thus this approximate ground state has a definite parity. With the assumed ground state Eq. (21), the energy and the mean photon number of the ground state is easy to obtain:

[TABLE]

Obviously, the parameter $\lambda$ is optimal if the projection $P(\lambda)=\langle G_{A}(\lambda)|\Psi_{0}\rangle$ is exactly equal to one, where $|\Psi_{0}\rangle$ is the exact ground state. Unfortunately, a simple expression of $|\Psi_{0}\rangle$ is unknown (though a series expression of $|\Psi_{0}\rangle$ can be given but it is quite useless due to the infinite series form). We here adopt an approximate but accurate enough form for $|\Psi_{0}\rangle$. Note that a unitary transformation $U=e^{\lambda\sigma_{z}(a-a^{{\dagger}})}$ can recast $|G_{A}\rangle$ to $|\tilde{G}_{A}\rangle=U|G_{A}\rangle=|-_{x},0\rangle$, which can be regarded as the zero-th order approximation for $|\tilde{\Psi}_{0}\rangle=U|\Psi_{0}\rangle$. Taking the complete orthogonal basis of $|\tilde{G}_{A}\rangle$ into consideration, we can construct the perturbative corrections of $|\tilde{\Psi}_{0}\rangle$. For perturbative calculation, we choose the basis of AA $|\pm_{z},N\rangle$ to be the complete orthogonal basis of $|\tilde{G}_{A}\rangle$. Note that it should work better to choose a more accurate basis, e.g., the GRWA basis, to expand the perturbation. But here, the choice of the AA or the GRWA basis makes little difference in the ground state calculation. According to perturbation theory, we can expand $|\tilde{\Psi}_{0}\rangle$ to the first order of the perturbation as

[TABLE]

where $c_{\pm,N}=\frac{\langle\pm_{x},N|\Delta\tilde{H}|\tilde{G}_{A}\rangle}{\langle\pm_{x},N|\tilde{H}_{0}|\pm_{x},N\rangle-E_{0}}$, $K={\sum\limits_{\{\pm,N\}}}^{\prime}c_{\pm,N}^{2}$ and $\tilde{H}_{0}=\tilde{H}_{\rm{AA}}$ with $\Delta\tilde{H}=\tilde{H}-\tilde{H}_{0}$. Note that the primed summation excludes the ground state itself with the label $\{-,0\}$, and $c_{\pm,N}$ will vanish if state $|\pm_{x},N\rangle$ has a different parity from $|\tilde{G}_{a}\rangle$. Then we can calculate

[TABLE]

For given $\omega,~{}\Omega$ and $g$, if $K$ is minimized by choosing optimal $\lambda$, then the obtained $|G_{A}\rangle$ would be optimal ground state. This process can be done numerically.

Before presenting the numerical results, it is useful to discuss some limit cases analytically. First, our result can recover the ones under the AA/GRWA in the strong coupling limit. In this limit $g$ is much larger than $\Omega$, the two-level splitting $\frac{\Omega}{2}\sigma_{x}$ in Eq. (1) can be safely neglected. Then one can verify that Eq. (25) takes the form $P_{\text{SC}}(\lambda)=[1+(\omega\lambda+g)^{2}/(E_{\text{AA}}^{+,1}-E_{0})^{2}]^{-{1\over 2}}$, which has an optimal value only when $\lambda=-g/\omega$. It corresponds exactly to the result under the AA/GRWA. Second, our result can reduce to the one under the GVM in the weak coupling limit. In this case, one can calculate $P_{\text{WC}}(\lambda)\simeq[1+(\omega\lambda+\Omega\lambda+g)^{2}/(E_{\text{AA}}^{+,1}-E_{0})^{2}]^{-{1\over 2}}$, which has the optimal value when $\lambda=-{g\over\omega+\Omega}$. This recovers the result under the GVM in Ref. GVMground .

For other cases, the analytic evaluation of the optimization on $P(\lambda)$ or $K$ is difficult. We should resort to numerical evaluation of the expression of $K$, which has much less numerical work than that of exact numeric. Figure 1 shows $K$ as the function of $\lambda$ in different $\omega$ and $g$. The leading four components in the summation of $K$ are also plotted. We can see the following characters: (i) When $\omega/\Omega$ is sufficiently large [see Fig. 1(a)], all the components except $c_{+,1}^{2}$ are negligible. Thus, $c_{+,1}^{2}$ is a good substitution of $K$ for the minimization. (ii) When $\omega$ is comparable to $\Omega$, $c_{\pm,N}^{2}$ for $N>1$ becomes sizable [see Fig. 1(b)]. However, $K$ still has only one minimum, which means that $c_{+,1}^{2}$ still can act as a substitution of $K$ for the minimization. (iii) When $\omega/\Omega$ is small [see Fig. 1(c)], $c_{\pm,N}^{2}$ for $N>1$ become important and $K$ shows two minimums. Therefore, none of $c_{\pm,N}^{2}$ can be taken as a substitution of $K$ for the minimization. The two minimums of $K$ should be considered equally. (iv) For $\omega/\Omega$ sufficiently small [see Fig. 1(d)], the series of $c_{\pm,N}^{2}$ lose convergency and a multi-minimum structure of $K$ appears. This complicated structure indicates that the coherent state form of the trial wavefunction Eq. (21) cannot capture the physics dominated by the spin-flipping and our scheme is no longer valid in this regime.

Figure 2 shows the optimal $\lambda$ and the corresponding $P(\lambda)$ for the negative-detuning ($\omega=2.0$), the resonance ($\omega=1.0$), and the positive-detuning ($\omega=0.5$) cases. For the two-minimum situation in the positive-detuning [see Fig. 2(a) and Fig. 2(d)] case, the optimal $\lambda$ can be evaluated effectively as

[TABLE]

where $\lambda_{A}$ and $\lambda_{B}$ are the two minimum positions of $\lambda$, and $K_{A}$ and $K_{B}$ are their corresponding $K$. We find that $\lambda_{\text{eff}}$ is dependent of the coupling strength, which is quite different from the fixed $\lambda$ result under the AA/GRWA grwa . Furthermore, $\lambda/g$ approaches to $-0.67$ in the weak coupling limit, which is consistent with the analytic result of $\lambda/g\rightarrow\frac{-1}{\omega+\Omega}$ obtained under the GVM. And $\lambda/g$ approaches to $-2.0$ in the strong coupling limit, which is consistent with the analytic result of $-\frac{1}{\omega}$ obtained under the AA/GRWA. In the whole parameter range, $P(\lambda)$ shows little deviation from $1$, which indicates that our obtained $|G_{A}\rangle$ is almost the same as the exact ground state. For the one-minimum situation in the resonance [see Fig. 2(b) and (e)] and in the negative-detuning [see Fig. 2(c) and (f)] cases, where the optimal $\lambda$ is determined by minimizing $c_{+,1}^{2}$, it is interesting to find that the change scope of $P$ converges closer and closer to one. It means that our scheme performs better and better with the increase of the $\omega$.

Figure 3 shows the ground-state energy $E_{0}$ and the mean photon number $\langle a^{\dagger}a\rangle_{0}$ as a function of the coupling strength $g$ for different detuning cases evaluated by different methods. We stress that although $E_{0}$ obtained by various methods in the negative detuning case is almost the same [see Fig. 3(a)], the mean photon number obtained by different methods behaves quite differently, as shown in Fig. 3(d). The GVM works better than the AA/GRWA in the weak coupling regime, while it gets worse in the intermediate coupling regime. However, our result obtained by optimizing $c_{+,1}^{2}$ matches well with the exact one even in the whole coupling regime. The improvement of our scheme to $E_{0}$ becomes more obvious with the decrease of $\omega$. In resonance case, we can see from Fig. 3(b) that the result from the AA/GRWA has a clear deviation from the exact value in the weak coupling regime and the one by the GVM shows a dramatic deviation in the strong coupling regime, while our result is consistent with the exact one almost in the whole coupling regime. For mean photon number in Fig. 3(e), our result is obviously more accurate than those obtained by the other methods. With a further decrease of $\omega$, the AA/GRWA and the GVM become worse and worse, but our results remain its good performance in evaluating $E_{0}$ and $\langle a^{\dagger}a\rangle_{0}$, as shown in Figs. 3(c) and (f).

Thus, our result reproduces the result of the GVM in weak coupling regime and the one of the AA/GRWA in the strong coupling limit regime, respectively. Our method tailors the advantages of the existing GVM and AA/GRWA methods. Nevertheless, with further decreasing $\omega$, the $K(\lambda)$ shows a multiple-minimum structure [see Fig. 1(d)], and the performance of our scheme also gets inaccurate. This indicates that the coherent state basis is no longer a good starting point in this case where the spin-flipping becomes dominant.

III.2 Applications to excited states

Our variational method can be also applied to the excited states. The GRWA-form excited state is adopted to be the trial state, but with unfixed parameter $\lambda$

[TABLE]

This trial state possesses a definite parity. As in the ground state, the perturbation scheme is again employed to determine the optimal value of $\lambda$. For convenience, we still discuss in the transformed representation. The zero-th order Hamiltonian is $\tilde{H}_{0}=\tilde{H}_{\rm GRWA}$, and the perturbation is $\Delta\tilde{H}=\tilde{H}-\tilde{H}_{\rm GRWA}$. $\lambda$ is determined by maximizing the projection $P(\lambda)=\langle\tilde{\Psi}_{A}^{\pm,N}(\lambda)|\tilde{\Psi}_{\pm,N}\rangle$, where $|\tilde{\Psi}_{\pm,N}\rangle$ is the exact excited state corresponding to $|\tilde{\Psi}_{A}^{\pm,N}\rangle$. $|\tilde{\Psi}_{\pm,N}\rangle$ can be evaluated perturbatively as

[TABLE]

where $c_{\pm,M}=\frac{\langle\tilde{\Psi}_{A}^{\pm,M}|\Delta\tilde{H}|\tilde{\Psi}_{A}^{\pm,N}\rangle}{E^{A}_{\pm,M}-E^{A}_{\pm,N}}$ and $K={\sum\limits_{\{\pm,M\}}}^{\prime}c_{\pm,M}^{2}$ with $E^{A}_{\pm,N}=\langle\tilde{\Psi}_{A}^{\pm,N}|\tilde{H}_{0}|\tilde{\Psi}_{A}^{\pm,N}\rangle$. Here, similarly to (24), the primed summation excludes the trial state itself with the label $\{\pm,N\}$. Then $\lambda$ can be calculated by optimizing $P(\lambda)$. The corresponding energy and mean photon number are

[TABLE]

In the large $g$ limit, $P(\lambda)$ approaches $[1+F^{\pm,N}(\omega\lambda+g)^{2}]^{1\over 2}$, where

[TABLE]

Then the optimal $\lambda$ takes $-g/\omega$, which recovers the result of the GRWA. In the small $g$ limit, $P(\lambda)$ approaches $[1+F^{\pm,N}(\omega\lambda+\Omega\lambda+g)^{2}]^{1\over 2}$. Then the optimal $\lambda$ reads $\lambda=\frac{-g}{\omega+\Omega}$. In both of the two limits, the chosen $\lambda$ is independent of excitation label $\{\pm,N\}$. It means the series of the approximate states hold the orthogonality. For other coupling cases, where $\lambda$ can be extracted by simple numerics, one can expect that $\lambda$ depends on excitation label $\{\pm,N\}$. Exactly speaking, differences in $\lambda$ would come to break the orthogonality, this arises from the simplicity of the trial state (27) we have adopted. Despite this small price, it is worth using such a simple trial state to gain quite many improvements in the physical properties, such as the energy spectrum and photon number.

To compare different methods for the excited states we illustrate by the example around resonance, i.e. $\omega=\Omega$, which is the most typical case. Since in the energy spectrum the level crossing occurs amongst the excited states with certain parities [see Fig. 4(a) $\&$ (d)], it is inconvenient to order the excited states in terms of energy. We order the excited states according to the label sequence of the GRWA basis. In the following, the first and second excited states are taken as examples. Since the GRWA modifies the AA in the excited cases, its improvement to the AA is remarkable and performs well in a quite broad regime. Thus, if one is viewing from a large scale of the coupling strength, the energy spectrum calculated by the GRWA nearly recovers the exact results, just as what our scheme performs [see Fig. 4(a) $\&$ (d)]. However, in the more detailed scales, the outcome of our method is more consistent with the exact one than the GRWA, especially in weak coupling regime [Fig. 4(b) $\&$ (e)]. For the mean photon number, although for the first excited state both the GRWA and our result are fairly accurate and thus show little difference in comparison with the exact one [Fig. 4(c)], for the second excited state the dramatic improvements over the GRWA from our variational method can be seen [Fig.4(f)]. For the second excited state, the GRWA does not capture the concave feature of the photon number in small $g$, while our result coincides with the exact one well in almost the whole coupling regime except some small discrepancy in a narrow window of intermediate coupling regime.

Besides the tuning case, we also check the validity of our method by considering a set of experiment-related parameters in Ref.ultra , which reads $\Omega=(4.20\pm 0.02)GHz$, $\omega/2\pi=(8.13\pm 0.01)GHz$ and $g/2\pi=(0.82\pm 0.03)GHz$. This is a large detuning case and the detuning is also much larger than the coupling strength since $\omega=12.16\Omega$ and $g^{*}=1.227\Omega$. In this case, we calculate the first and second exited state energies. Referring to the numerically exact results, Fig.5 presents a comparison between the results by our method and those by the GRWA, in which a significant improvement is seen. It should be mentioned that, since the AA is modified by the GRWA and the GVM is limited to the ground state, they are not included in the above comparisons.

IV Conclusions and discussions

We have introduced a mean photon number dependent variational method to evaluate the properties of the Rabi model. Our scheme combines the advantages of the existing AA/GRWA and GVM approximations. For the ground state, the trial state is the superposition of two coherent states with opposite displacements, and the key parameter $\lambda$ is determined by maximizing the projection of the assumed state and the exact one, which has been approximated by a perturbation theory. In the weak coupling regime our result is in agreement with that of the GVM which is accurate in this regime but deviates from the exact one in the strong coupling regime. On the other hand, in the strong coupling regime, our result is consistent with that obtained by the AA/GRWA which works well in this regime but deviates from the exact one in the weak coupling regime. In the intermediate regime, our method not only provides a natural crossover from the AA/GRWA to the GVM but also yields an obvious improvement over all of them. It is shown that the improvements for the mean photon number are even more substantial than the energy. Thus our method is valid in whole coupling regime with not sufficiently small frequency of the bosonic field. In the small limit of the frequency of the bosonic field, both our method and the existing AA/GRWA and GVM work no longer well, which indicates that the position-displaced oscillator basis is no longer a good trial state and one should explore new starting point in this regime.

Although most variational methods limit to the ground state, our variational scheme can be also applied to the excited states. For the excited states, the deviation of the GRWA in the weak coupling regime is still considerable. In contrast, the validity of our scheme for the whole coupling regime still remains. The quantitative deviation of the GRWA energy in the weak coupling regime and the qualitative missing of concave feature for the mean photon number in the GRWA are well rectified in our scheme.

In short, our variational scheme efficiently improves several previous widely-used approximations such as the AA, the GRWA and the GVM, with better qualitative and quantitative descriptions on the physics of the model. Despite that the integrability and exactly analytical expressions of energy spectra have been obtained for the Rabi model in Ref.Braak2011 , series expansion form of its wavefunction is still inconvenient to calculate the physical variables in the model. On the contrary, our method directly starts from the wavefunction assumption and emphasizes its physics meaning. For example, the ground state form of our wavefunction is not only directly related the mean photon number but also useful for discussion of nonclassical states preparation nori . In particular, our method to evaluate properties of the ground state and the low excited states might be applicable to the multi-mode Rabi model, i.e., the so-called spin-boson model, where a novel quantum phase transition characterized by the low level energies is intensively studied recently Vojta2012 ; Tong2011 .

Acknowledgements

We greatly appreciate Gang Chen and Lixian Yu for useful discussions. The work is partly supported by the programs for NSFC, PCSIRT (Grant No. IRT1251), the national program for basic research and the Fundamental Research Funds for the Central Universities of China.

Appendix A Unitary transformation on the Hamiltonian

A unitary operator $U=e^{\lambda\sigma_{z}(a-a^{{\dagger}})}$ transforms the model Hamiltonian $H=\omega_{0}a^{{\dagger}}a+\frac{1}{2}\Omega\sigma_{x}+\lambda\sigma_{z}(a^{{\dagger}}+a)$ into $\tilde{H}$ in the new representation. With the formula

[TABLE]

where $C_{n}=B$ if $n=0$ and $C_{n+1}=[A,C_{n}]$ otherwise, one can obtain

[TABLE]

The terms of $\cosh[2\lambda(a-a^{{\dagger}})]$ and $\sinh[2\lambda(a-a^{{\dagger}})]$ can be expanded in powers of $a$ and $a^{{\dagger}}$ according to formula

[TABLE]

Below we will use the associated Laguerre function defined by

[TABLE]

and and the Laguerre function

[TABLE]

For the factor $(a^{{\dagger}})^{m}a^{n}$, one has

[TABLE]

where

[TABLE]

Here $\hat{N}=a^{{\dagger}}a$ is particle number operator. Set $\nu=-2\lambda$.

[TABLE]

For $m-n=2k\geq 0$,

[TABLE]

For $m-n=-2k<0$,

[TABLE]

By the definition of the function

[TABLE]

one can expand Eq.(43) as

[TABLE]

By the same way, one has

[TABLE]

Substitute Eq.(43) and Eq.(44) into the transformed Hamiltonian Eq.(33), the expansion Eq.(4) is obtained.